# MAGIC: Modular Auto-encoder for Generalisable Model Inversion with Bias
Corrections

Yihang She
University of Cambridge
<EMAIL_ADDRESS>
&Clement Atzberger
Mantle Labs
<EMAIL_ADDRESS>
&Andrew Blake
University of Cambridge
<EMAIL_ADDRESS>
&Adriano Gualandi
University of Cambridge
<EMAIL_ADDRESS>
&Srinivasan Keshav
University of Cambridge
<EMAIL_ADDRESS>

###### Abstract

Scientists often model physical processes to understand the natural world and
uncover the causation behind observations. Due to unavoidable simplification,
discrepancies often arise between model predictions and actual observations,
in the form of systematic biases, whose impact varies with model completeness.
Classical model inversion methods such as Bayesian inference or regressive
neural networks tend either to overlook biases or make assumptions about their
nature during data preprocessing, potentially leading to implausible results.
Inspired by recent work in inverse graphics, we replace the decoder stage of a
standard autoencoder with a physical model followed by a bias-correction
layer. This generalisable approach simultaneously inverts the model and
corrects its biases in an end-to-end manner without making strong assumptions
about the nature of the biases. We demonstrate the effectiveness of our
approach using two physical models from disparate domains: a complex radiative
transfer model from remote sensing; and a volcanic deformation model from
geodesy. Our method matches or surpasses results from classical approaches
without requiring biases to be explicitly filtered out, suggesting an
effective pathway for understanding the causation of various physical
processes. The code is available on https://github.com/yihshe/MAGIC.git.

##  1 Introduction

The quantitative and mathematical modeling of natural phenomena is a
foundational task in science. Inverting a physical model, to deduce underlying
causal variables from observations, is a powerful technique that has been used
widely (Adler and Holder, 2021; Lahoz and Schneider, 2014; Zhdanov, 2015).
Prior attempts at model inversion have struggled with the issue of model
incompleteness, where outcomes predicted by a model differ from actual
outcomes (Groetsch and Groetsch, 1993; Zhdanov, 2015). These approaches either
overlook model biases or make assumptions about them during data
preprocessing, leading to significant error in inferred variable values.

Two major recent advances encourage us to revisit this long-standing problem.
First, modern sensor technologies have made available extensive data-sets of
measured observations (Claverie et al., 2018). Second, although standard
statistical approaches in machine learning may be unable to yield
interpretable representations of that data (Chen et al., 2016), in recent
years auto-encoders have been used to uncover disentangled representations in
low dimensions (Locatello et al., 2019). However these approaches lack an
effective inference mechanism to obtain variables with clearly defined
physical meanings (Kumar et al., 2018).

It is notable that researchers in graphics have explored incorporating a
differentiable renderer into an auto-encoder to recover interpretable
representations from visual data (Loper and Black, 2014; Kato and Harada,
2019). Similar approaches have also been used in some specific physics
problems (Lu et al., 2020; Zérah et al., 2022) but it is unclear whether these
approaches can be generalized. This is because auto-encoder gradients need to
backpropagate through the physical model but existing physical models are
often written in legacy code, and not designed to be differentiable in the way
that deep learning demands.

A further issue is that existing autocoder-based approaches rarely compensate
for systematic biases in physical models. Ideally, a physical model’s outputs
exhibit no systematic bias when compared with measured observations. However,
real-world models may exhibit particularly prominent systematic biases because
they only approximately represent a complex, and often noisy, physical system
(Abramowitz et al., 2007; Zhdanov, 2015; Widlowski et al., 2013). This
requires bias correction, with the complexity of the bias correction depending
on the incompleteness of the physical model. Unfortunately, the addition of a
bias correction layer can modify the causal relationship between the physical
model’s causal variables and the corrected output. Thus, there is a fine
balance between the complexity of bias correction and interpretability of the
causal variables extracted by inversion of the bias corrected model.

To address these issues, we first demonstrate how to augment a deterministic
physical model with a bias correction layer that learns a transformation from
model output to observed outcomes. This greatly increases model prediction
accuracy. Second, we use an autoencoder-like network to invert this augmented
model to extract causal variables from observations. We do this in two
disparate physical science domains: a radiative transfer model (RTM)
(Atzberger, 2000) that computes spectral signatures from biophysical variables
and is not inherently differentiable; and a volcanic deformation model (Mogi,
1958) that computes surface deformations due to magma chamber activity. We
demonstrate the generalisability of our approach by making the complex RTM
fully differentiable and stable during training, using the power of Large
Language Models (LLMs), an approach that could be applied more generally to
inverting legacy implementations of complex physical models. Finally, we use
two physical models with different levels of completeness to investigate the
complexity of the corresponding bias correction.

By detailed comparison of learning outcomes against baselines derived from
classical methods, we demonstrate that auto-encoders can be used as a
generalizable approach for simultaneous model inversion and bias correction.
This suggests a pathway for a more accurate understanding of the causation of
diverse physical processes, by means of end-to-end learning.

##  2 Related Work

###  2.1 Disentangled representation learning

A key aspect of our approach is on interpretability and disentanglement. This
relates to the concept of disentangled representation learning, which seeks to
uncover lower-dimensional, semantically meaningful factors of variation from
high-dimensional data (Locatello et al., 2019). In the paradigm of
unsupervised learning of disentangled representations, Variational
Autoencoders (VAEs) (Kingma and Welling, 2013) and their variants, such as
$\beta$-VAE (Higgins et al., 2017), are based on probabilistic modeling. These
approaches allow disentangled variables to emerge from learning (Kumar et al.,
2018). Alternatively, Generative Adversarial Networks (GANs) (Goodfellow et
al., 2014) have also demonstrated success in learning interpretable and
disentangled representations for image synthesis, as seen in InfoGAN (Chen et
al., 2016) and StyleGAN (Karras et al., 2019). However, as generative models,
the disentangled variables typically emerge only during learning and lack an
effective inference mechanism (Kumar et al., 2018). Moreover, as Locatello et
al. (2019) suggested, inductive biases are essential for unsupervised
approaches to learning disentangled representations. For physical models, the
situation is simplified because the forward model is deterministic and
benefits from having a predetermined set of interpretable variables, namely
the physical input parameters to the model.

###  2.2 Inverse problems in physical sciences

Across various domains in natural sciences, forward numerical models are
developed based on physical principles and inverting them has led to
significant applications e.g. electrical impedance tomography in healthcare
(Adler and Holder, 2021), data assimilation in climate science (Lahoz and
Schneider, 2014), and understanding the Earth’s structure in geophysics
(Zhdanov, 2015). We test our approach on two domains (see section 3.2). First,
Radiative Transfer Models (RTMs) from remote sensing (Suits, 1971; Atzberger,
2000) belong to a family of canopy reflectance models that simulate spectral
bidirectional forest reflectance, given inputs of important forest
characteristics (Widlowski et al., 2013; Li and Strahler, 1985; Suits, 1971;
Atzberger, 2000). Second, volcanic deformation models from geodesy simulate
surface deformation due to magma chamber activity (Segall, 2013; Mogi, 1958)
and come from a broader family of geophysical models where inversion is a key
problem (Zhdanov, 2015). Classical inversion approaches include Bayesian
inference, numerical optimization, and look-up tables (Groetsch and Groetsch,
1993; Goel, 1988; Combal et al., 2003). More recently, neural network
regressors have also been used (Gong, 1999; Schlerf and Atzberger, 2006;
Mousavi et al., 2024).

###  2.3 Autoencoder as an approach to inversion

We do not claim originality simply for using auto-encoders to invert physical
models; that has been extensively investigated in computer science, especially
in inverse graphics (Kato et al., 2020). A seminal contribution is OpenDR
(Loper and Black, 2014), which proposes differentiable rendering to invert a
graphics renderer by learning with backpropagation. Since then, auto-encoders
have been a popular choice to incorporate a differentiable renderer into self-
supervised learning pipelines (Kato et al., 2020), and have been applied to 3D
reconstruction, and human pose estimation from single-view images (Kato and
Harada, 2019; Pavlakos et al., 2018). This approach has been applied to a few
scientific domains (Lu et al., 2020; Zérah et al., 2022), and recently,
specifically for a canopy-reflectance model (Zérah et al., 2023). However none
of these studies addresses systematic bias, and that is what is especially
distinctive about our approach.

##  3 Methods

###  3.1 Learning the inverse of physical models

To infer physical variables $Z$ from measured observations $X$, we use an
auto-encoder architecture with bias correction to train an encoder that maps
measured observations to physical variables (Fig. 1). The physical process
that generates the observations $X$ is written:

$X=\mathbf{F}(Z)+\mathbf{B}(Q)+\eta$ (1)

where $\mathbf{F}$ is the modelled process acting on physical variables $Z$;
$\mathbf{B}(Q)$ is the systematic bias between $\mathbf{F}(Z)$ and $X$
resulting from unmodeled processes $\mathbf{B}$, and depending on a collection
of other variables $Q$; and $\eta$ represents random noise.

A classical auto-encoder $\mathbf{M}_{\text{A}}$ consists of an encoder
$\mathbf{E}_{\text{A}}$ and a decoder $\mathbf{D}_{\text{A}}$ (eq. 2).

$X_{\text{A}}=\mathbf{D}_{\text{A}}(\mathbf{E}_{\text{A}}(X))$ (2)

Given a sufficiently sophisticated network architecture,
$\mathbf{M}_{\text{A}}$ can minimize Mean Squared Error
$\mathcal{L}_{\text{MSE}}(\theta_{\mathbf{E}_{\text{A}}})$ between $X$ and
$X_{\text{A}}$, but the latent codes $Z_{\text{A}}$ computed as
$\mathbf{E}_{\text{A}}(X)$ need not be semantically interpretable in terms of
well defined physical properties. Replacing $\mathbf{D}_{\text{A}}$ in eq. 2
with the physical model $\mathbf{F}$ however ensures that an encoder
$\mathbf{E}_{\text{B}}$ embeds $X$ as interpretable physical variables
$Z_{\text{B}}=\mathbf{E}_{\text{B}}(X)$, which are then passed to $\mathbf{F}$
to reconstruct observations (autoencoder $\mathbf{M}_{\text{B}}$):

$X_{\text{B}}=\mathbf{F}(\mathbf{E}_{\text{B}}(X))$ (3)

Ideally, when $\mathbf{F}$ is perfectly built to account fully for the
physical processes, there is no need for $\mathbf{B}(Q)$ in eq. 1. In that
case, $\mathcal{L}_{\text{MSE}}(\theta_{\mathbf{E}_{\text{B}}})$ should be on
a par with $\mathcal{L}_{\text{MSE}}(\theta_{\mathbf{E}_{\text{A}}})$. In
reality the rigidity and oversimplification of $\mathbf{F}$ gives appreciable
reconstruction mismatch between $X$ and $X_{\text{B}}$ so we add additional
non-linear bias-correction layers $\mathbf{C}$ after $\mathbf{F}$ (autoencoder
$\mathbf{M}_{\text{C}}$, also see Fig. 1):

$X_{\text{C}}=\mathbf{C}(\mathbf{F}(\mathbf{E}_{\text{C}}(X))).$ (4)

The correction layer $\mathbf{C}$ increases the flexibility of $\mathbf{F}$
and decreases
$\mathcal{L}_{\text{MSE}}(\theta_{\mathbf{E}_{\text{C}}},\theta_{\mathbf{C}})$
to improve the learning outcome of $Z_{\text{C}}$. The risk, though, is that
input variables $Z_{\text{C}}$ are no longer anchored as strongly, especially
if with appreciable flexibility in $\mathbf{C}$, allowing $Z_{\text{C}}$ to
entangle and drift away from physical interpretations. A solution is to limit
the complexity of $\mathbf{C}$ to two linear transformations with a single
ReLU activation $\sigma$ between them:

$X_{\text{C}}=W_{2}\sigma(W_{1}\mathbf{F}(\mathbf{E}_{\text{C}}(X))+b_{1})+b_{2}.$
(5)

Here, $W_{1}$ expands the input to four times its original dimension, and
$W_{2}$ reduces it back to its original size. This design mirrors the feed-
forward head in Transformers (Vaswani et al., 2017), which has allowed just
enough flexibility in non-linear transformations — exactly what we need (see
section 4).

Figure 1: Learning the inverse, end-to-end, in an autoencoder that includes
correction layers $\mathbf{C}$.

###  3.2 Physical models for inversion

INFORM: The INFORM radiative transfer model from remote sensing (denoted by
$\mathbf{F}_{\text{RTM}}$), developed by Atzberger (2000), consists of a set
of sub-models designed to simulate the reflectance of forest canopies given
parameters representing their biophysical and geometric properties. Despite
decades of research, the complexity of forest structures leads to
discrepancies between simulated and measured spectra. These biases stem from
over-simplification of the physics and significantly affect accuracy of
inferred variables with classical approaches (Widlowski et al., 2013). We aim
to estimate certain set of the biophysical variables, whose plausible ranges
are known (Tab. 1), while fixing other variables at typical values for the
study area (see section A.3). In practice, $\mathbf{F}_{\text{RTM}}$ is
implemented with NumPy and comprises approximately 1700 lines of code, and is
not explicitly differentiable.

Mogi: The Mogi model (denoted by $\mathbf{F}_{\text{Mogi}}$) is a model used
in geodesy to describe the displacement field at the surface of the Earth
resulting from a spherical pressure source, typically a magma chamber, at
depth (Mogi, 1958). The inverse problem for $\mathbf{F}_{\text{Mogi}}$
involves estimating the location and volume change of magma chamber activity
from observed surface displacements (Tab. 2 and section A.5). Compared with
$\mathbf{F}_{\text{RTM}}$, $\mathbf{F}_{\text{Mogi}}$ is built on simple
assumptions, with the complexity of the Earth’s crust being crudely
approximated by an elastic half-space. Despite its simplicity, the Mogi model
is able to capture the first-order patterns of surface deformation associated
with volcanic sources. The challenge is that, the volcanic deformation can be
small if compared to other sources of deformation (e.g. hydrological loading),
posing challenges to the bias correction (see section A.6.2).

Table 1: Biophysical variables of $\mathbf{F}_{\text{RTM}}$ to learn. Bound
ranges of $Z_{\mathrm{h}}$ and $Z_{\mathrm{cd}}$ are marked with *, as they
will be inferred from $Z_{\mathrm{fc}}$ and fixed stem density using
allometric equations (Jucker et al., 2017) to mitigate the ill-posed nature of
the inverse problem (see section A.3).

Sub-model | Leaf Model | Canopy Model | Forest Model
---|---|---|---
Variable |  Structure Parameter |  Chlorophyll A+B |  Water Content |  Dry Matter |  Leaf Area Index |  Undergrowth LAI |  Fractional Coverage |  Tree Height |  Crown Diameter
Acronym | N | cab | cw | cm | LAI | LAIu | fc | h | cd
Min | 1 | 10 | 0.001 | 0.005 | 0.01 | 0.01 | 0.1 | * | *
Max | 3 | 80 | 0.02 | 0.05 | 5 | 1 | 1 | * | *

Table 2: Geophysical variables of $\mathbf{F}_{\text{Mogi}}$ to learn. The
source center $(x_{m},y_{m})$ is assumed to lie within the GNSS stations’
boundaries (section 3.4), with depth and volume ranges matching those used in
(Walwer et al., 2016).

Variable | Epicenter-X ($\mathrm{km}$) | Epicenter-Y ($\mathrm{km}$) | Depth ($\mathrm{km}$) | Volume Change ($10^{6}~{}\mathrm{m^{3}}$)
---|---|---|---|---
Acronym | $x_{m}$ | $y_{m}$ | $d$ | $\Delta V$
Min | -9.33 | -5.80 | 2 | -10
Max | 14.35 | 7.62 | 20 | 10

###  3.3 End-to-end learning with physical models

Differentiability is required for auto-encoder workflow and the implementation
of $\mathbf{F}_{\text{Mogi}}$ in PyTorch is straightforward as it only
involves a few lines of equations (section A.5). The RTM
$\mathbf{F}_{\text{RTM}}$ is complex and all its NumPy operations need
conversion to PyTorch, identifying functions that are not differentiable, and
constructing backward function for each to obtain approximated derivatives.

Recent developments in LLMs would suggest using GPT-4 (OpenAI, 2023) so that
common operations get automatically converted to PyTorch, and others are
flagged for custom treatment. In this way, we converted 1,742 lines of code in
less than a week (see section A.1). In practice, canopy-radiation interaction
needs exponential, logarithm, and square root functions which can lead to
numerically unstable derivatives. This is handled by bypassing instability
points during backpropagation (see section A.2).

###  3.4 Datasets

The dataset for $\mathbf{F}_{\text{RTM}}$ consists of 17692 individual spectra
$X_{\text{S2}}=\\{x_{\text{S2},i},i=1,\ldots,11\\}$, extracted from 1283
individual sites across Austria observed by Sentinel-2 (Drusch et al., 2012).
Spectra $X_{\text{S2}}$ form a temporal sequence covering April to October
2018 — both coniferous and deciduous forests — consisting of 12 species.
Spectra comprise 11 bands (see section A.4).

For $\mathbf{F}_{\text{Mogi}}$, the dataset comes from 12 Global Navigation
Satellite System (GNSS) stations around the Akutan Volcano, one of the most
active volcanoes in the Aleutian Islands (Alaska, USA). Every day, each GNSS
station records its position at the surface in east, north and vertical
directions (denoted
$X_{\text{GNSS}}=\\{x_{\text{e,i}},x_{\text{n,i}},x_{\text{v,i}},i=1,\ldots,12\\}$).
In total, we have observations from 6525 consecutive days from 2006 until
early 2024 (http://geodesy.unr.edu, (Blewitt et al., 2018)) (see section A.6).

For both problems, collecting ground truth data is largely impractical. The
interpretation of results has to rely on plausibility checks on variable
distributions and temporal evolution patterns. For $\mathbf{F}_{\text{RTM}}$,
we use prior knowledge of tree species and temporal data to validate
biophysical estimates. For $\mathbf{F}_{\text{Mogi}}$, geophysical estimates
of a magma source and reconstructed transient GNSS signals are validated
against existing literature (Ji and Herring, 2011; Walwer et al., 2016).

###  3.5 Baseline

In addition to comparing $\mathbf{M}_{\text{A}}$, $\mathbf{M}_{\text{B}}$, and
$\mathbf{M}_{\text{C}}$, we have also compared with classical approach
baselines. Further implementation details can be found in section A.8.

For $\mathbf{F}_{\text{RTM}}$, we trained a neural network regressor
$\mathbf{M}_{\text{D, RTM}}$ using sampled variables and spectra from
$\mathbf{F}_{\text{RTM}}$ — a classical approach for inversion (Gong, 1999) as
a baseline. Model $\mathbf{M}_{\text{D, RTM}}$ shares the same architecture as
the encoder of $\mathbf{M}_{\text{C, RTM}}$. It is first trained to regress
the sampled variables from the synthetic spectra and then applied to the real
spectra $X_{\text{S2}}$ to predict physical variables ($Z_{\text{D, RTM}}$).
To obtain the spectral reconstruction $X_{\text{D, RTM}}$, $Z_{\text{D, RTM}}$
is passed to $\mathbf{F}_{\text{RTM}}$.

For $\mathbf{F}_{\text{Mogi}}$, we compared our inversion results with
classical results (Walwer et al., 2016) which preprocessed the GNSS time
series over 2006-2014 as follows. First they used Singular Spectrum Analysis
to detect likely transient signals by inspecting the processed time series
across stations and directions. Then, they performed a grid search of
parameters to invert $\mathbf{F}_{\text{Mogi}}$ by fitting its estimates to
the extracted transient signals. They fixed the location of the Mogi source
and estimated the depth and volume change for periods likely due to magma
movement. Notably they detected significant volcanic inflation of Akutan in
early 2008, also reported by Ji and Herring (2011). Thus for validation, we
have saved the GNSS time series from 2006 to 2009 as a test set to cover the
2008 event.

##  4 Results

###  4.1 Effects of bias correction

The systematic bias of $\mathbf{F}_{\text{RTM}}$ is evident comparing the
input spectra $X_{\text{S2}}$ with the reconstruction $X_{\text{S2,B}}$ from
$\mathbf{M}_{\text{B, RTM}}$ trained without bias correction (Fig. 2(a)).
Further illustrations of bias corrections are detailed in section A.9.

Similarly $\mathbf{F}_{\text{Mogi}}$ shows systematic bias comparing
$X_{\text{GNSS}}$ with the reconstruction from $\mathbf{M}_{\text{B, Mogi}}$
(Fig. 2(b)). Bias correction using $\mathbf{M}_{\text{C, Mogi}}$ achieves some
improvement in reconstruction of GNSS signals (Tab. 4) though less
dramatically than $\mathbf{M}_{\text{C, RTM}}$, as indicated by the $R^{2}$
values (Figs. 2(a) and 2(b)). The RTM is a more complex model, presenting
greater challenges in inversion, but imposes a strong prior and captures
physical processes effectively. In contrast, $\mathbf{F}_{\text{Mogi}}$ makes
simpler assumptions — considering only a point-like pressure source for
surface deformation — forming a less informative bottleneck in the decoder and
so poses a greater challenge for bias correction. However our correction layer
is still able to capture most of the systematic bias and improves the learned
inverse.

Table 3: RTM: The reconstruction loss of models under different steps of the inversion. $\mathcal{L}_{\mathrm{MSE}}$ under $\mathbf{M}_{\mathrm{C,RTM}}$ is comparable to the classical auto-encoder $\mathbf{M}_{\mathrm{A,RTM}}$, and is reduced by more than an order of magnitude compared to the regressor baseline $\mathbf{M}_{\mathrm{D,RTM}}$. Model | Ablation | $\mathcal{L}_{\mathrm{MSE,train}}$ | $\mathcal{L}_{\mathrm{MSE,val}}$ | $\mathcal{L}_{\mathrm{MSE,test}}$
---|---|---|---|---
$\mathbf{M}_{\text{A, RTM}}$ | w/o $\mathbf{F}_{\text{RTM}}$, w/o $\mathbf{C}_{\text{RTM}}$ | 0.0193 | 0.0219 | 0.0191
$\mathbf{M}_{\text{B, RTM}}$ | w/ $\mathbf{F}_{\text{RTM}}$, w/o $\mathbf{C}_{\text{RTM}}$ | 0.0875 | 0.0833 | 0.0856
$\mathbf{M}_{\text{C, RTM}}$ | w/ $\mathbf{F}_{\text{RTM}}$, w/ $\mathbf{C}_{\text{RTM}}$ | 0.0210 | 0.0235 | 0.0217
$\mathbf{M}_{\text{D, RTM}}$ | N/A | - | - | 0.6676
Table 4: Mogi: The reconstruction loss under different steps. $\mathcal{L}_{\mathrm{MSE}}$ for $\mathbf{M}_{\mathrm{C,Mogi}}$ outperforms $\mathbf{M}_{\mathrm{B,Mogi}}$ and matches $\mathbf{M}_{\mathrm{D,Mogi}}$, but the bias correction is less pronounced compared to the RTM model — a more complete model with less complex biases to address (see Fig. 2). Model | Ablation | $\mathcal{L}_{\mathrm{MSE,train}}$ | $\mathcal{L}_{\mathrm{MSE,val}}$ | $\mathcal{L}_{\mathrm{MSE,test}}$
---|---|---|---|---
$\mathbf{M}_{\text{A, Mogi}}$ | w/o $\mathbf{F}_{\text{Mogi}}$, w/o $\mathbf{C}_{\text{Mogi}}$ | 0.4389 | 0.9119 | 0.6324
$\mathbf{M}_{\text{B, Mogi}}$ | w/ $\mathbf{F}_{\text{Mogi}}$, w/o $\mathbf{C}_{\text{Mogi}}$ | 0.6998 | 0.9721 | 0.8095
$\mathbf{M}_{\text{C, Mogi}}$ | w/ $\mathbf{F}_{\text{Mogi}}$, w/ $\mathbf{C}_{\text{Mogi}}$ | 0.4803 | 0.9338 | 0.7110

(a) RTM: Reconstruction accuracy illustrated for 2 of the spectral bands.

(b) Mogi: Reconstruction accuracy of vertical displacements illustrated for 2
of the GNSS stations.

Figure 2: Reconstruction accuracy without bias correction module displays
clear bias. Blue: $\mathbf{M}_{\mathrm{C}}$. Red: $\mathbf{M}_{\mathrm{B}}$.
After bias correction, reconstruction for Mogi is not as linear as for RTM,
likely due to its simplicity which acts as an information bottleneck in the
decoder, highlighting the trade-off between model simplicity and bias
correction complexity.

Bias correction also enhances the learning of variables in latent space for
$\mathbf{F}_{\text{RTM}}$ and for $\mathbf{F}_{\text{Mogi}}$ — compare
retrieved variables with and without a bias correction (Fig. 3) — for
$\mathbf{F}_{\text{RTM}}$ the distributions of several variables are distorted
to extremes if bias is not corrected (Fig. 3(a)).

Similarly for $\mathbf{F}_{\text{Mogi}}$ bias pushes the center of the Mogi
source to extremes (Fig. 3(b)) if not corrected, and although $Z_{\text{d}}$
and $Z_{\text{$\Delta$V}}$ follow unimodal or bimodal distributions, they have
implausible value ranges without bias correction. Both Ji and Herring (2011)
and Walwer et al. (2016) have suggested that the depth $Z_{\text{d}}$ less
than $10~{}\text{km}$ and the volume change $Z_{\text{$\Delta$V}}$ range
rarely exceed $-5\cdot 10^{6}~{}\text{m}^{3}$ to $5\cdot
10^{6}~{}\text{m}^{3}$ for the evaluated period at the Akutan volcano (Walwer
et al., 2016). Indeed, if $Z_{\text{d}}$ exceeds $10~{}\text{km}$, then the
volume change required to reproduce the observed displacement would be
unreasonably large. It is possible that, without bias correction, the model
attempts to replicate a common mode signal putting the source at greater depth
and thus generating a large wavelength signal at the surface.

(a) Distributions of variables: $Z_{\text{B, RTM}}$ vs. $Z_{\text{C, RTM}}$

(b) Distributions of variables: $Z_{\text{B, Mogi}}$ vs. $Z_{\text{C, Mogi}}$

Figure 3: Bias correction tightened distributions of variables. Without bias
correction, variable distributions (red) are implausible, breaking out from
preset bounded ranges. With bias correction (blue) learned distributions are
more aligned with the physical sense.

(a) Co-distributions: $Z_{\mathrm{D,RTM}}$

(b) Co-distributions: $Z_{\mathrm{C,RTM}}$

Figure 4: Pairwise co-distributions of variables illustrated for
$Z_{\text{LAI}}$. Red: coniferous forest. Blue: deciduous forest. Our model
can learn distinct physical patterns.

###  4.2 Evaluation of physical variables

We clustered RTM variables by forest type. Coniferous and deciduous forests
have distinct distributions and the co-distribution of variables
$Z_{\text{LAI}}$ appears positively correlated with $Z_{\text{cw}}$, and
negatively with $Z_{\text{LAIu}}$ under $\mathbf{M}_{\text{C, RTM}}$ (Fig.
4(b)). Without bias correction no correlation patterns are visible (Fig.
4(a)). This makes sense: higher leaf area at the canopy level
($Z_{\text{LAI}}$) indicates plant growth and hence higher water content
($Z_{\text{cw}}$) and chlorophyll content ($Z_{\text{cab}}$) but reduces light
penetration, reducing plant growth in the understory, hence lower
$Z_{\text{LAIu}}$. Further illustrations of physical variables are detailed in
section A.10.

Figure 5: Temporal variations of inferred physical parameters
$Z_{\mathrm{C,RTM}}$ illustrated for $Z_{\mathrm{LAI}}$, $Z_{\mathrm{LAIu}}$,
and $Z_{\mathrm{fc}}$. $\mathbf{M}_{\mathrm{C,RTM}}$ effectively captures
distinct, temporally smooth, and plausible variations for different forest
types, while those from $\mathbf{M}_{\mathrm{D,RTM}}$ appear to be less smooth
(see section A.10.1).

Temporal variation exhibits more consistent and clearer patterns for some
variables (Fig. 5). Understory plant growth ($Z_{\text{LAIu}}$) declines from
April due to reduced light penetration from increased canopy growth
($Z_{\text{LAI}}$). Coniferous forests exhibit different patterns than
deciduous ones, with higher fractional coverage ($Z_{\text{fc}}$) and less
variation over time, as would be expected. Temporal variations learned by the
$\mathbf{M}_{\text{D, RTM}}$ appear less consistent (see section A.10.1).

The 7 inferred variables are clustered by species and Jeffereys-Matusita (JM)
distance, varying between 0 and 2, is calculated to assess separability of the
12 species. As shown in Fig. 6, our model achieves clear separation between
different forest types (deciduous vs. coniferous), and improved separability
within forest type, suggesting a degree of disentanglement of inferred
variables. Note that the coniferous species ‘Larix Decidua’ (European Larch)
appears close to the deciduous group, but actually that makes sense as it is
one of the few conifers that sheds its needles in the fall.

(a) JM distance based on $Z_{\text{D, RTM}}$

(b) JM distance based on $Z_{\text{C, RTM}}$

Figure 6: Pairwise Jeffreys-Matusita (JM) distance between species based on
the learned variables. Red: coniferous species. Blue: deciduous species.
Species distributions from different forest groups become more disentangled
after bias correction.

For the Mogi model, classical signal processing of the time series and Mogi
inversions for 2006–2009 show a transient deformation event due to the
inflation of a magma source in early 2008 (Ji and Herring, 2011; Walwer et
al., 2016). It lasts about half a year and is characterized by the uplift and
radial extension of the surface around the Mogi source. This is consistent
with the temporal variation of $Z_{\text{$\Delta$V}}$ learned by our model.
The magnitude of $Z_{\text{$\Delta$V}}$ and the estimated values of
$Z_{\text{d}}$ (ranging from 6 to 9 km) and the epicenter (located roughly
beneath Akutan volcano), which remain relatively constant over time, align
with previous estimates Walwer et al. (2016) which estimated Mogi parameters
using grid search to fit the model to the processed GNSS signals.

Examining two GNSS stations, AV08 and AV12, located on opposite sides of the
volcano, we compare their GNSS displacements as reconstructed by
$\mathbf{F}_{\text{Mogi}}$ under the $\mathbf{M}_{\text{C, Mogi}}$ model. The
model successfully reconstructs the transient signals of the 2008 inflation
and the subsequent slow deflation, characterised by specific directional
movements, with consistency across stations and directions that align with Ji
and Herring (2011) and Walwer et al. (2016). Note that classical methods
(Walwer et al., 2016) required the GNSS signals to be filtered for assumed
biases, and transient signals identified by comparing processed signals across
GNSS stations, followed by grid search to invert the Mogi model. In contrast,
our auto-encoder approach allows $\mathbf{M}_{\text{C, RTM}}$ to capture
transient signals, estimate Mogi parameters, and correct biases in a single
stage without relying on additional data filtering or strong bias assumptions.
Furthermore, this inference is achieved through a single forward pass of
$\mathbf{M}_{\text{C, RTM}}$, which is trained on datasets with no overlap
with the test period.

Figure 7: Time series of inferred physical parameters. The red curve is fitted
by a Kalman filter to indicate the temporal variations more clearly. The
volume change $Z_{\Delta V}$ successfully captured a major volcanic inflation
event in 2008, reaching a peak of approximately $5\cdot
10^{6}~{}\mathrm{m^{3}}$ at a depth of just less than 10 km, consistent with
previous findings (Walwer et al., 2016). Figure 8: Time series of GNSS
displacements illustrated for two GNSS stations, located at the opposite side
of the Mogi source. Orange: $X_{\mathrm{GNSS}}$. Blue: Mogi reconstructions
using $Z_{\mathrm{C,Mogi}}$. Red: Mogi reconstructions fitted by a Kalman
filter. Our model successfully reconstructs the 2008 volcanic inflation
transient signals, characterized by radial expansion and vertical uplift, with
displacement magnitudes of 10 mm, aligned with prior findings (Ji and Herring,
2011; Walwer et al., 2016).

##  5 Conclusion and Limitations

We have developed an end-to-end approach to learn the inverse of deterministic
physical models while correcting for biases. We demonstrated the
generalisability of our approach by applying it to physical models from very
different domains, each with distinct complexities, differentiability issues,
and data dynamics. Extensive analysis shows that our approach has improved
reconstruction loss and strongly suggests the plausibility of the recovered
physical variables.

Our work suffers from some limitations. First, our approach does not extend to
physical models with stochastic outcomes, where the same inputs can lead to
different outputs, a much more complex inversion problem. Second, we have only
studied two representative physical models. However, we anticipate that our
approach can be applied to models in other physical domains. Finally, we use
the same bias correction independent of the physical model. While we found
this to be effective, given the trade-off between model completeness and
complexity necessary for bias correction, future work could use neural
architecture search (Ren et al., 2021) to identify the optimal bias correction
layer for a specific inverse problem.

Acknowledgement: This work was supported by the UKRI Centre for Doctoral
Training in Application of Artificial Intelligence to the study of
Environmental Risks (reference EP/S022961/1) and Cambridge Centre for Carbon
Credits. We would also like to thank Markus Immitzer from Mantle Labs for
sharing the Sentinel-2 data with us.

## References

  * Adler and Holder [2021] Andy Adler and David Holder.  _Electrical impedance tomography: methods, history and applications_.  CRC Press, 2021.
  * Lahoz and Schneider [2014] William A Lahoz and Philipp Schneider.  Data assimilation: making sense of earth observation.  _Frontiers in Environmental Science_ , 2:16, 2014.
  * Zhdanov [2015] Michael S Zhdanov.  _Inverse theory and applications in geophysics_ , volume 36.  Elsevier, 2015.
  * Groetsch and Groetsch [1993] Charles W Groetsch and CW Groetsch.  _Inverse problems in the mathematical sciences_ , volume 52.  Springer, 1993.
  * Claverie et al. [2018] Martin Claverie, Junchang Ju, Jeffrey G Masek, Jennifer L Dungan, Eric F Vermote, Jean-Claude Roger, Sergii V Skakun, and Christopher Justice.  The harmonized landsat and sentinel-2 surface reflectance data set.  _Remote sensing of environment_ , 219:145–161, 2018.
  * Chen et al. [2016] Xi Chen, Yan Duan, Rein Houthooft, John Schulman, Ilya Sutskever, and Pieter Abbeel.  Infogan: Interpretable representation learning by information maximizing generative adversarial nets.  _Advances in neural information processing systems_ , 29, 2016.
  * Locatello et al. [2019] Francesco Locatello, Stefan Bauer, Mario Lucic, Gunnar Rätsch, Sylvain Gelly, Bernhard Schölkopf, and Olivier Bachem.  Challenging common assumptions in the unsupervised learning of disentangled representations, 2019.
  * Kumar et al. [2018] Abhishek Kumar, Prasanna Sattigeri, and Avinash Balakrishnan.  Variational inference of disentangled latent concepts from unlabeled observations, 2018.
  * Loper and Black [2014] Matthew M Loper and Michael J Black.  Opendr: An approximate differentiable renderer.  In _Computer Vision–ECCV 2014: 13th European Conference, Zurich, Switzerland, September 6-12, 2014, Proceedings, Part VII 13_ , pages 154–169. Springer, 2014.
  * Kato and Harada [2019] Hiroharu Kato and Tatsuya Harada.  Learning view priors for single-view 3d reconstruction.  In _Proceedings of the IEEE/CVF conference on computer vision and pattern recognition_ , pages 9778–9787, 2019.
  * Lu et al. [2020] Peter Y Lu, Samuel Kim, and Marin Soljačić.  Extracting interpretable physical parameters from spatiotemporal systems using unsupervised learning.  _Physical Review X_ , 10(3):031056, 2020.
  * Zérah et al. [2022] Yoël Zérah, Silvia Valero, and Jordi Inglada.  Physics-guided interpretable probabilistic representation learning for high resolution image time series.  _IEEE Transactions on Geoscience and Remote Sensing_ , 2022.
  * Abramowitz et al. [2007] Gab Abramowitz, Andy Pitman, Hoshin Gupta, Eva Kowalczyk, and Yingping Wang.  Systematic bias in land surface models.  _Journal of Hydrometeorology_ , 8(5):989–1001, 2007.
  * Widlowski et al. [2013] J-L Widlowski, B Pinty, M Lopatka, C Atzberger, D Buzica, Michaël Chelle, M Disney, J-P Gastellu-Etchegorry, M Gerboles, N Gobron, et al.  The fourth radiation transfer model intercomparison (rami-iv): Proficiency testing of canopy reflectance models with iso-13528.  _Journal of Geophysical Research: Atmospheres_ , 118(13):6869–6890, 2013.
  * Atzberger [2000] Clement Atzberger.  Development of an invertible forest reflectance model: The infor-model.  In _A decade of trans-european remote sensing cooperation. Proceedings of the 20th EARSeL Symposium Dresden, Germany_ , volume 14, pages 39–44, 2000.
  * Mogi [1958] Kiyoo Mogi.  Relations between the eruptions of various volcanoes and the deformations of the ground surfaces around them.  _Earthq Res Inst_ , 36:99–134, 1958.
  * Kingma and Welling [2013] Diederik P Kingma and Max Welling.  Auto-encoding variational bayes.  _arXiv preprint arXiv:1312.6114_ , 2013.
  * Higgins et al. [2017] Irina Higgins, Loic Matthey, Arka Pal, Christopher Burgess, Xavier Glorot, Matthew Botvinick, Shakir Mohamed, and Alexander Lerchner.  beta-vae: Learning basic visual concepts with a constrained variational framework.  In _International conference on learning representations_ , 2017.
  * Goodfellow et al. [2014] Ian Goodfellow, Jean Pouget-Abadie, Mehdi Mirza, Bing Xu, David Warde-Farley, Sherjil Ozair, Aaron Courville, and Yoshua Bengio.  Generative adversarial nets.  _Advances in neural information processing systems_ , 27, 2014.
  * Karras et al. [2019] Tero Karras, Samuli Laine, and Timo Aila.  A style-based generator architecture for generative adversarial networks.  In _Proceedings of the IEEE/CVF conference on computer vision and pattern recognition_ , pages 4401–4410, 2019.
  * Suits [1971] Gwynn H Suits.  The calculation of the directional reflectance of a vegetative canopy.  _Remote Sensing of Environment_ , 2:117–125, 1971.
  * Li and Strahler [1985] Xiaowen Li and Alan H Strahler.  Geometric-optical modeling of a conifer forest canopy.  _IEEE Transactions on Geoscience and Remote Sensing_ , (5):705–721, 1985.
  * Segall [2013] Paul Segall.  Volcano deformation and eruption forecasting.  _Geological Society, London, Special Publications_ , 380(1):85–106, 2013.
  * Goel [1988] Narendra S Goel.  Models of vegetation canopy reflectance and their use in estimation of biophysical parameters from reflectance data.  _Remote sensing reviews_ , 4(1):1–212, 1988.
  * Combal et al. [2003] B Combal, Frédéric Baret, M Weiss, Alain Trubuil, D Mace, A Pragnere, R Myneni, Y Knyazikhin, and L Wang.  Retrieval of canopy biophysical variables from bidirectional reflectance: Using prior information to solve the ill-posed inverse problem.  _Remote sensing of environment_ , 84(1):1–15, 2003.
  * Gong [1999] P Gong.  Inverting a canopy reflectance model using a neural network.  _International Journal of Remote Sensing_ , 20(1):111–122, 1999.
  * Schlerf and Atzberger [2006] Martin Schlerf and Clement Atzberger.  Inversion of a forest reflectance model to estimate structural canopy variables from hyperspectral remote sensing data.  _Remote sensing of environment_ , 100(3):281–294, 2006.
  * Mousavi et al. [2024] S Mostafa Mousavi, Gregory C Beroza, Tapan Mukerji, and Majid Rasht-Behesht.  Applications of deep neural networks in exploration seismology: A technical survey.  _Geophysics_ , 89(1):WA95–WA115, 2024.
  * Kato et al. [2020] Hiroharu Kato, Deniz Beker, Mihai Morariu, Takahiro Ando, Toru Matsuoka, Wadim Kehl, and Adrien Gaidon.  Differentiable rendering: A survey.  _arXiv preprint arXiv:2006.12057_ , 2020.
  * Pavlakos et al. [2018] Georgios Pavlakos, Luyang Zhu, Xiaowei Zhou, and Kostas Daniilidis.  Learning to estimate 3d human pose and shape from a single color image.  In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , pages 459–468, 2018.
  * Zérah et al. [2023] Yoël Zérah, Silvia Valero, and Jordi Inglada.  Physics-constrained deep learning for biophysical parameter retrieval from sentinel-2 images: inversion of the prosail model.  _Available at SSRN 4671923_ , 2023.
  * Vaswani et al. [2017] Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Łukasz Kaiser, and Illia Polosukhin.  Attention is all you need.  _Advances in neural information processing systems_ , 30, 2017.
  * Jucker et al. [2017] Tommaso Jucker, John Caspersen, Jérôme Chave, Cécile Antin, Nicolas Barbier, Frans Bongers, Michele Dalponte, Karin Y van Ewijk, David I Forrester, Matthias Haeni, et al.  Allometric equations for integrating remote sensing imagery into forest monitoring programmes.  _Global change biology_ , 23(1):177–190, 2017.
  * Walwer et al. [2016] Damian Walwer, Eric Calais, and Michael Ghil.  Data-adaptive detection of transient deformation in geodetic networks.  _Journal of Geophysical Research: Solid Earth_ , 121(3):2129–2152, 2016.
  * OpenAI [2023] OpenAI.  Gpt-4 technical report, 2023.
  * Drusch et al. [2012] Matthias Drusch, Umberto Del Bello, Sébastien Carlier, Olivier Colin, Veronica Fernandez, Ferran Gascon, Bianca Hoersch, Claudia Isola, Paolo Laberinti, Philippe Martimort, et al.  Sentinel-2: Esa’s optical high-resolution mission for gmes operational services.  _Remote sensing of Environment_ , 120:25–36, 2012.
  * Blewitt et al. [2018] Geoffrey Blewitt, William Hammond, et al.  Harnessing the gps data explosion for interdisciplinary science.  _Eos_ , 99, 2018.
  * Ji and Herring [2011] Kang Hyeun Ji and Thomas A Herring.  Transient signal detection using gps measurements: Transient inflation at akutan volcano, alaska, during early 2008.  _Geophysical Research Letters_ , 38(6), 2011.
  * Ren et al. [2021] Pengzhen Ren, Yun Xiao, Xiaojun Chang, Po-Yao Huang, Zhihui Li, Xiaojiang Chen, and Xin Wang.  A comprehensive survey of neural architecture search: Challenges and solutions.  _ACM Computing Surveys (CSUR)_ , 54(4):1–34, 2021.
  * Gualandi et al. [2016] Adriano Gualandi, Enrico Serpelloni, and Maria Elina Belardinelli.  Blind source separation problem in gps time series.  _Journal of Geodesy_ , 90(4):323–341, 2016.

###### Contents

  1. 1 Introduction
  2. 2 Related Work
    1. 2.1 Disentangled representation learning
    2. 2.2 Inverse problems in physical sciences
    3. 2.3 Autoencoder as an approach to inversion
  3. 3 Methods
    1. 3.1 Learning the inverse of physical models
    2. 3.2 Physical models for inversion
    3. 3.3 End-to-end learning with physical models
    4. 3.4 Datasets
    5. 3.5 Baseline
  4. 4 Results
    1. 4.1 Effects of bias correction
    2. 4.2 Evaluation of physical variables
  5. 5 Conclusion and Limitations
  6. A Appendix
    1. A.1 Making a fully differentiable physical model for deep learning
      1. A.1.1 Conversion of $\mathbf{F}_{\mathrm{RTM}}$ to PyTorch assisted by GPT-4
      2. A.1.2 Validating the PyTorch implementation of $\mathbf{F}_{\mathrm{RTM}}$
    2. A.2 Bypassing instability points in the training loop
      1. A.2.1 Numerical instability of $\mathbf{F}_{\mathrm{RTM}}$
      2. A.2.2 Forward pass
      3. A.2.3 Backpropagation
      4. A.2.4 Stabilized training process
    3. A.3 Biophysical variables of $\mathbf{F}_{\mathrm{RTM}}$
    4. A.4 Sentinel-2 dataset: spectral bands, temporal, and species information
    5. A.5 Mathematical formulations of $\mathbf{F}_{\mathrm{Mogi}}$
    6. A.6 GNSS dataset
      1. A.6.1 GNSS stations around the Akutan volcano
      2. A.6.2 Assumed components of GNSS time series
    7. A.7 Smoothness term added to MSE loss for Mogi inversion
    8. A.8 Implementation details
    9. A.9 Extended results for bias correction
      1. A.9.1 Bias correction for $\mathbf{F}_{\mathrm{RTM}}$
      2. A.9.2 Bias correction for $\mathbf{F}_{\mathrm{Mogi}}$
    10. A.10 Extended results for the evaluation of physical variables
      1. A.10.1 Retrieved biophysical variables of $\mathbf{F}_{\mathrm{RTM}}$
      2. A.10.2 Reconstructed GNSS signals based on $Z_{\mathrm{C,Mogi}}$

##  Appendix A Appendix

###  A.1 Making a fully differentiable physical model for deep learning

####  A.1.1 Conversion of $\mathbf{F}_{\mathrm{RTM}}$ to PyTorch assisted by
GPT-4

The original modules of $\mathbf{F}_{\mathrm{RTM}}$ are extracted from a
vegetation app called EnMap. It is implemented using Numpy arrays and
operations. To track the computational graph and enable gradient
backpropagation, it was reimplemented using PyTorch operations. However,
$\mathbf{F}_{\mathrm{RTM}}$ is not inherently differentiable, and as a
complicated physical model, converting it in PyTorch manually was challenging.
In light of the recent development of large language models (LLMs), we decided
to utilize GPT-4 [OpenAI, 2023] to assist in the conversion from NumPy to
PyTorch. While GPT-4 reduced a significant amount of repetitive work, we ran
into the following limitations:

GPT-4 is useful for converting common operations from NumPy into PyTorch. We
used it to successfully convert 1,742 lines of code across various scripts
from the original implementation to PyTorch (e.g. Listing LABEL:lst:np and
Listing LABEL:lst:torch). However, some operations e.g. the exponential
integral function from SciPy, have no equivalent in PyTorch (Listing
LABEL:lst:torch). To calculate their derivative, we implemented our own
backward functions to define derivatives.

Additionally, GPT-4 only accepts text input of limited length, which means it
can only convert a script snippet piece by piece. As the conversation
lengthens, it made mistakes. For instance, it converted numpy.radians to
torch.radians, which does not exist in PyTorch. Thus, it was still important
to carefully check the conversion, run unit tests, and work on compiling all
the scripts together. Interestingly, GPT allows for prompt engineering, which
means we can provide feedback to refine its conversions, such as requesting it
to highlight uncertain parts for review in subsequent conversions.

⬇

1n = PD_refractive

2k = (np.outer(Cab, PD_k_Cab) + np.outer(Car, PD_k_Car) + np.outer(Anth,
PD_k_Anth) +

3 np.outer(Cbrown, PD_k_Brown) + np.outer(Cw, PD_k_Cw) + np.outer(Cm,
PD_k_Cm)) / N[:, np.newaxis]

4

5ind_k0_row, ind_k0_col = np.where(k == 0)

6

7if len(ind_k0_row) > 0:

8 k[ind_k0_row, ind_k0_col] = np.finfo(float).eps

9trans = (1 - k) * np.exp(-k) + (k ** 2) * exp1(k)

Listing 1: Example code of $\mathbf{F}_{\mathrm{RTM}}$ initially implemented
using NumPy. It features an exponential integral function $\mathrm{exp1}$ in
line 9. Additionally, it’s important to highlight that in line 3, the
structure parameter $N$ functions as a denominator in the calculation,
therefore, it cannot be zero.

⬇

1n = PD_refractive

2k = (torch.outer(Cab, PD_k_Cab) + torch.outer(Car, PD_k_Car) +
torch.outer(Anth, PD_k_Anth) +

3 torch.outer(Cbrown, PD_k_Brown) + torch.outer(Cw, PD_k_Cw) + torch.outer(Cm,
PD_k_Cm)) / N.unsqueeze(-1)

4

5ind_k0_row, ind_k0_col = torch.where(k == 0)

6

7if len(ind_k0_row) > 0:

8 k[ind_k0_row, ind_k0_col] = torch.finfo(float).eps

9trans = (1 - k) * torch.exp(-k) + (k ** 2) *
torch.from_numpy(exp1(k.numpy()))

Listing 2: Conversion of $\mathbf{F}_{\mathrm{RTM}}$ to PyTorch by GPT. It is
followed by feedback on $\mathrm{exp1}$ to raise user’s attention. Here is an
excerpt: "Please note that the $\mathrm{exp1}$ function is not available in
PyTorch as of my last training cut-off in September 2021. We can import it
from the scipy.special library, but we have to convert the tensor to numpy
before using it and convert the result back to a tensor. It’s not ideal in
terms of performance, but I’m afraid there’s no other option for now. Please
check if there is an updated version of PyTorch that includes this function."

This experiment was done in May 2023, when GPT was still in its early phase of
development. It has had several features updated since then. For instance, it
can now accept longer text inputs. Nonetheless, we believe the limitations
discussed above still largely hold true and are worth considering when
converting legacy implementations of other complex physical models.

####  A.1.2 Validating the PyTorch implementation of
$\mathbf{F}_{\mathrm{RTM}}$

To validate the implementation of $\mathbf{F}_{\mathrm{RTM}}$ in PyTorch
against its original version in NumPy, we randomly sampled 10,000 sets of
variables, passed them through each implementation, and compared their outputs
(Tab. 5). The mismatch rate of the outputs is 0.457% when the absolute
tolerance is set to 1e-5. The maximum absolute difference over all 130000
output reflectances from the two implementations is 3.050e-5 wherein the
physical unit of reflectance here is 1. This indicates that
$\mathbf{F}_{\mathrm{RTM}}$ implemented in PyTorch is substantially equivalent
to its original version and, therefore can be used in subsequent tasks. The
PyTorch implementation for the RTM has been released on GitHub
(https://github.com/yihshe/MAGIC.git) as a public Python library to facilitate
its integration into other deep-learning workflows.

Table 5: Unit test of the outputs of our PyTorch implementation of
$\mathbf{F}_{\mathrm{RTM}}$ against the original NumPy version. The total
number of elements to compare is 130000 as we simulate the same set of 10000
samples for both implementations, each with 13 spectral bands of Sentinel-2.

Total Elements | Mismatched Elements | Absolute Tolerance | Mismatch Ratio | Max Absolute Difference
---|---|---|---|---
130000 | 594 | 1e-5 | 0.457% | 3.050e-5

###  A.2 Bypassing instability points in the training loop

####  A.2.1 Numerical instability of $\mathbf{F}_{\mathrm{RTM}}$

With $\mathbf{F}_{\mathrm{RTM}}$ now implemented in PyTorch, we can
backpropagate gradients through this physical model. However, the model is
sensitive to input variable ranges and is non-differentiable at certain
points. Furthermore, complex operations like exponential, logarithm, and
square root, often lead to numerical instability during derivative
computation, even though these operations are theoretically differentiable. To
solve these issues, we have developed the following methods to stabilize
training when using a complex physical model like $\mathbf{F}_{\mathrm{RTM}}$.

####  A.2.2 Forward pass

By definition, input variables $Z_{\mathrm{RTM}}$ are non-negative values with
physical meanings (Tab. 1). In practice, $\mathbf{F}_{\mathrm{RTM}}$ can fail
in the forward pass when $Z_{\mathrm{RTM}}$ from other value ranges are
provided (e.g. structure parameter $Z_{\mathrm{N}}$ in Listing LABEL:lst:np).
The same concept applies to $\mathbf{F}_{\text{Mogi}}$. Thus, for both models,
we choose to initially infer normalized variables $\Lambda\in[0,1]$ and
subsequently map these to variables $Z\in[Z_{\text{min}},Z_{\text{max}}]$ in
their original scale, according to their value ranges (Tabs. 1 and 2).

$Z=(Z_{\text{max}}-Z_{\text{min}})\cdot\Lambda+Z_{\text{min}}$ (6)

This approach not only facilitates a successful forward pass of physical
models but also narrows down the search space for optimization.

####  A.2.3 Backpropagation

Despite encoding the forward pass successfully, backpropagation remained
unstable for $\mathbf{F}_{\mathrm{RTM}}$. For successful learning, it is
necessary to backpropagate gradients through the physical model. The
aforementioned caveats of $\mathbf{F}_{\mathrm{RTM}}$ resulted in
backpropagation failures, leading to NaN gradient and then NaN loss. During
experiments to stabilize the backpropagation, we noticed that NaN gradients,
when they initially appear, are always in association with specific sets of
$\mathbf{F}_{\mathrm{RTM}}$ variables. Ideally, one would apply constraints to
specific $\mathbf{F}_{\mathrm{RTM}}$ operations causing these NaN gradients.
However, working out such constraints seemed challenging unless we understood
the semantics of $\mathbf{F}_{\mathrm{RTM}}$ and its differentiability
thoroughly, a complex task.

We, therefore, introduced a simple workaround: since it is the NaN gradients
that cause subsequent learning failure, we replace them with a small random
value, sampled from a uniform [0,1] distribution and scaled by 1e-7, whenever
these gradients first appear. The specific workflow of the algorithm is
described in algorithm 1

Algorithm 1 Gradient stabilization

Calculate gradients

Initialize list $grads$, containing gradients of all model parameters where
gradients exist and contain any NaN values

if $grads$ is not empty then

Set a small constant epsilon equal to 1e-7

for each v in $grads$ do

Generate random values of the same shape as v, scaled by epsilon

Create a mask where v is NaN or equals 0

Replace values in v where the mask is True with corresponding random values

end for

end if

Update gradients

####  A.2.4 Stabilized training process

Our methods effectively stabilized the learning of
$\mathbf{M}_{\mathrm{B,RTM}}$ and $\mathbf{M}_{\mathrm{C,RTM}}$ in the
autoencoder framework, despite instability of $\mathbf{F}_{\mathrm{RTM}}$
during backpropagation. As can be seen in Fig. 9, both training and validation
losses of $\mathbf{M}_{\mathrm{C,RTM}}$ converge over the epochs. Therefore,
by preventing the propagation of NaN gradients to earlier layers and allowing
the normal progression of the forward pass, the stabilizer has effectively
enabled the training process to bypass instability points in the optimization
space and allowed a continuous search for optimal points.

Figure 9: Stabilized training of $\mathbf{M}_{\mathrm{C,RTM}}$. Our algorithm
to update the gradients has overcome the instability of
$\mathbf{F}_{\mathrm{RTM}}$ during backpropagation and allowed the convergence
of training loss.

###  A.3 Biophysical variables of $\mathbf{F}_{\mathrm{RTM}}$

The input of $\mathbf{F}_{\mathrm{RTM}}$ consists of biophysical variables of
three hierarchical levels. Note that the original $\mathbf{F}_{\mathrm{RTM}}$
does not have fractional coverage $Z_{\mathrm{fc}}$ as an input variable. We
include $Z_{\mathrm{fc}}$ as one of the seven variables to be learned directly
by our model, which will be used to infer crown diameter $Z_{\mathrm{cd}}$ and
height $Z_{\mathrm{h}}$ based on derived equations. As the fractional coverage
is jointly defined by the stem density and crown diameter within each unit
hector (or 10,000 $m^{2}$), $Z_{\mathrm{cd}}$ can be derived given
$Z_{\mathrm{fc}}$ and $Z_{\mathrm{cd}}$ using eq. 7. Furthermore, to derive
$Z_{\mathrm{h}}$, we fit an allometric equation (Fig. 10, eq. 8) using the
samples of temperate forests from the global allometric database [Jucker et
al., 2017]. $R^{2}$ of the derived equation between $Z_{\mathrm{h}}$ and
$Z_{\mathrm{cd}}$ is 0.383 (eq. 8).

Table 6: Overview of the biophysical variables of $\mathbf{F}_{\mathrm{RTM}}$.
These variables can be attributed to three hierarchical levels. 7 variables
are learned directly. * $Z_{\mathrm{cd}}$ and $Z_{\mathrm{h}}$ will be
inferred from $fc$ using eq. 7 and eq. 8, respectively.

Group | Variable | Acronym | To Learn | Default Value | Sample Range
---|---|---|---|---|---
Min | Max
Background | Soil brightness factor | psoil | ✗ | 0.8 | - | -
Leaf Model | Structure Parameter | N | ✓ | - | 1 | 3
Chlorophyll A+B | cab | ✓ | - | 10 | 80
Water Content | cw | ✓ | - | 0.001 | 0.02
Dry Matter | cm | ✓ | - | 0.005 | 0.05
Carotenoids | car | ✗ | 10 | - | -
Brown Pigments | cbrown | ✗ | 0.25 | - | -
Anthocyanins | anth | ✗ | 2 | - | -
Proteins | cp | ✗ | 0.0015 | - | -
Cabon-based Constituents | cbc | ✗ | 0.01 | - | -
Canopy Model | Leaf Area Index | LAI | ✓ | - | 0.01 | 5
Leaf Angle Distribution | typeLIDF | ✗ | Beta Distribution | - | -
Hot Spot Size | hspot | ✗ | 0.01 | - | -
Observation Zenith Angle | tto | ✗ | 0 | - | -
Sun Zenith Angle | tts | ✗ | 30 | - | -
Relative Azimuth Angle | psi | ✗ | 0 | - | -
Forest Model | Undergrowth LAI | LAIu | ✓ | - | 0.01 | 1
Stem Density | sd | ✗ | 500 | - | -
Fractional Coverage | fc | ✓ | - | 0.1 | 1
Tree Height | h | ✓ | * | * | *
| Crown Diameter | cd | ✓ | * | * | *

$Z_{\mathrm{cd}}=2\cdot\sqrt{\frac{Z_{\mathrm{fc}}\cdot 10000}{\pi\cdot
Z_{\mathrm{sd}}}}$ (7)

Figure 10: Allometric relationship between canopy height $Z_{\mathrm{h}}$ and
crown diameter $Z_{\mathrm{cd}}$. We fit an allometric equation between
$Z_{\mathrm{h}}$ and $Z_{\mathrm{cd}}$ (eq. 8). $R^{2}$ of the fitted equation
is 0.383.

$Z_{\mathrm{h}}=\exp\left(2.117+0.507\cdot\ln(Z_{\mathrm{cd}})\right)$ (8)

###  A.4 Sentinel-2 dataset: spectral bands, temporal, and species information

The detailed information of the spectral bands of $X_{\mathrm{S2}}$ can be
viewed in Tab. 7. Statistics of the $X_{\mathrm{S2}}$ can be viewed in Tab. 8.
These spectra were sampled from individual sites covering a time sequence of
14 timestamps. Note that our current approach to inverting
$\mathbf{F}_{\mathrm{RTM}}$ does not integrate temporal information into the
learning process, although temporal variations are evaluated. Such temporal
information could serve as useful prior knowledge to boost the model’s
performance, especially to ensure the consistency of temporal variations,
which we will consider in future works.

Table 7: Sentinel-2 bands to use. VNIR stands for Visible and Near Infrared.
SWIR stands for Short Wave Infrared.

Band | B2 | B3 | B4 | B5 | B6 | B7 | B8 | B8a | B9 | B11 | B12
---|---|---|---|---|---|---|---|---|---|---|---
Resolution | 10 m | 10 m | 10 m | 20 m | 20 m | 20 m | 10 m | 20 m | 60 m | 20 m | 20 m
Central Wavelength | 490 nm | 560 nm | 665 nm | 705 nm | 740 nm | 783 nm | 842 nm | 865 nm | 940 nm | 1610 nm | 2190 nm
Description | Blue | Green | Red | VNIR | VNIR | VNIR | VNIR | VNIR | SWIR | SWIR | SWIR

Table 8: Information of the Sentinel-2 dataset $X_{\mathrm{S2}}$

Total Number of Spectra | Number of Individual Sites | Number of Dates | Number of Species
---|---|---|---
17962 | 1283 | 14 | 12

All samples from $X_{\mathrm{S2}}$ cover both coniferous and deciduous forests
consisting of 5 and 7 species, respectively:

  * •

Coniferous forest: Pseudotsuga Menziesii, Picea Abies, Pinus Nigra, Larix
Decidua, Pinus Sylvestris.

  * •

Deciduous forest: Prunus Spp, Fagus Sylvatica, Carpinus Betulus, Quercus Spp,
Acer Pseudoplatanus, Fraxinus Excelsior, Alnus Glutinosa.

###  A.5 Mathematical formulations of $\mathbf{F}_{\mathrm{Mogi}}$

$\mathbf{F}_{\mathrm{Mogi}}$ describes the displacement field on the surface
resulting from a spherical pressure source, typically a magma chamber, at
depth. The inverse problem of $\mathbf{F}_{\mathrm{Mogi}}$ deals with
estimating the location and volume change of magma chamber activities given
the observed displacements at ground stations. In this project,
$\mathbf{F}_{\mathrm{Mogi}}$ to use is further simplified by assuming only a
point pressure source. Specifically, for $N$ ground stations located at
$\\{x_{i},y_{i},z_{i},i=1,\cdots,N\\}$, the displacement field
$X_{\mathrm{GNSS}}=\\{x_{e,i},x_{n,i},x_{v,i},i=1,\cdots,N\\}$ at the surface
along the east, north, and vertical directions due to the subsurface
volumetric change of a Mogi source is estimated by:

$\displaystyle\hat{x}_{e,i}$
$\displaystyle=\frac{\alpha(x_{i}-Z_{\mathrm{x_{m}}})}{R_{i}^{3}},$ (9)
$\displaystyle\hat{x}_{n,i}$
$\displaystyle=\frac{\alpha(y_{i}-Z_{\mathrm{y_{m}}})}{R_{i}^{3}},$ (10)
$\displaystyle\hat{x}_{v,i}$
$\displaystyle=\frac{\alpha(z_{i}-Z_{\mathrm{d}})}{R_{i}^{3}}$ (11)

where $Z_{\mathrm{x_{m}}}$ and $Z_{\mathrm{y_{m}}}$ are the horizontal
coordinates of the Mogi source, and $Z_{\mathrm{d}}$ is its depth. $\alpha$ is
the scaling factor that includes the volume change $Z_{\mathrm{\Delta V}}$ and
the material properties (Poisson’s ratio of the medium), given by:

$\alpha=\frac{(1-\nu)Z_{\mathrm{\Delta V}}}{\pi}$ (12)

and $R_{i}$ is the radial distance from the Mogi source point to the station
at ($x_{i}$, $y_{i}$, $0$) (the equations are valid for a half-space, so
elevation is neglected):

$R_{i}=\sqrt{(x_{i}-x_{m})^{2}+(y_{i}-y_{m})^{2}+Z_{\mathrm{d}}^{2}}$ (13)

For the inversion problem, we assume the following parameters are known:

  * •

$\nu$, Poisson’s ratio, usually set at 0.25, a typical value for the Earth’s
crust.

  * •

$\\{x_{i},y_{i}\,i=1,\cdots,N\\}$, locations of the GNSS stations.

Thus, our goal is to estimate the location
$(Z_{\mathrm{x_{m}}},Z_{\mathrm{y_{m}}},Z_{\mathrm{d}})$ and volume change
$Z_{\mathrm{\Delta V}}$ of the Mogi source, given the measured displacements
recorded by the GNSS stations
$X_{\mathrm{GNSS}}=\\{u_{e,i},u_{n,i},u_{v,i},i=1,\cdots,N\\}$.

###  A.6 GNSS dataset

####  A.6.1 GNSS stations around the Akutan volcano

Surface displacement is monitored via continuous GNSS (Global Navigation
Satellite System) stations that produce daily position time series
(http://geodesy.unr.edu, Blewitt et al. [2018]. For the Akutan volcano, we
have 12 GNSS stations within the region (Fig. 11). From these stations, we
have compiled daily time-series GNSS displacement data spanning from 2006 to
2024, which we will use to test our inversion approach. It is common for GNSS
stations to occasionally miss data due to various issues such as power outages
or instrument replacements. Approximately 20% of our raw time series data is
missing, which, although significant, is still considered good in terms of
GNSS observations. We have employed a variational Bayesian Independent
Component Analysis [Gualandi et al., 2016] to fill up these data gaps using
existing observations to make sure we have a complete daily time series.

Our model takes the displacement field observed at each single day as a
training sample. However, we can check the plausibility of the inversion by
evaluating the temporal variations of the retrieved variables and
reconstructed signals comparing them against existing literature [Ji and
Herring, 2011, Walwer et al., 2016].

Figure 11: Locations of the 12 GNSS stations around the Akutan Volcano.

####  A.6.2 Assumed components of GNSS time series

In addition to deformations due to a Mogi source, the surface displacement
observed by a GNSS station can be influenced by other physical processes and
noise sources. In simplified form, a displacement $u(t)$ of a station along
certain direction at time $t$ can be modeled as [Walwer et al., 2016]:

$u(t)=S_{\mathrm{linear}}(t)+S_{\mathrm{annual}}(t)+S_{\mathrm{semiannual}}(t)+S_{\mathrm{eruption}}(t)+\epsilon_{\mathrm{pink}}+\epsilon_{\mathrm{white}}$
(14)

Where:

$\displaystyle S_{\mathrm{linear}}(t)=q+m\cdot t$ (15) $\displaystyle
S_{\mathrm{annual}}(t)=A_{1}\cdot\sin\left(\frac{2\pi
t}{T_{1}}\right)+B_{1}\cdot\cos\left(\frac{2\pi t}{T_{1}}\right)$ (16)
$\displaystyle S_{\mathrm{semiannual}}(t)=A_{2}\cdot\sin\left(\frac{2\pi
t}{T_{2}}\right)+B_{2}\cdot\cos\left(\frac{2\pi t}{T_{2}}\right)$ (17)
$\displaystyle S_{\mathrm{eruption}}(t)=C\cdot Z_{\mathrm{\Delta V,t}}$ (18)

$S_{\mathrm{linear}}$ represents the linear trend, often due to tectonic plate
movement; $S_{\mathrm{annual}}$ and $S_{\mathrm{semiannual}}$ represent the
annual and semiannual seasonality of the signal, potentially due to
hydrological processes; and $S_{\mathrm{eruption}}$ accounts for the
deformation due to Mogi source since the eruption event happening at $t_{s}$
and depends on the volume change $Z_{\mathrm{\Delta V,t}}$ at time step $t$.
$\epsilon_{\mathrm{pink}}$ and $\epsilon_{\mathrm{white}}$ are pink and white
noises, respectively, commonly observed in natural phenomena.

The volcanic deformation $S_{\mathrm{eruption}}(t)$ can be small (usually at
$\mathrm{mm}$ level) if compared to other sources of deformation in eq. 14
($\mathrm{mm}$ to $\mathrm{cm}$ level), making the bias correction both an
essential ingredient to get accurate results and a challenging task due to the
incompleteness of $\mathbf{F}_{\mathrm{Mogi}}$.

###  A.7 Smoothness term added to MSE loss for Mogi inversion

Compared to $X_{\mathrm{S2}}$, $X_{\mathrm{GNSS}}$ exhibits stronger temporal
dynamics as a daily time series. While our model takes the displacement field
of each day as input, enhancing the smoothness of the Mogi reconstructions
along the temporal dimension will align the results more closely with the
expectation that the deformation from the volcanic source
($S_{\mathrm{eruption}}$) should demonstrate temporal continuity, rather than
being influenced by the patterns of systematic biases from other components in
eq. 14. To enhance the temporal continuity of the Mogi reconstructions, we
structure the data such that each batch contains mini-sequences sampled from
$T$ consecutive days. We have also introduced an additional smoothness term to
$\mathbf{M}_{\mathrm{C,Mogi}}$, which is expressed as follows:

$\mathcal{L}_{\text{total}}(\theta_{\mathbf{E}_{\mathrm{C}}},\theta_{\mathbf{C}})=\mathcal{L}_{\text{MSE}}(\theta_{\mathbf{E}_{\mathrm{C}}},\theta_{\mathbf{C}})+\lambda\sum_{j=1}^{B}\sum_{t=1}^{T-1}\|x_{\mathrm{\mathbf{F},C,j}}^{(t+1)}-x_{\mathrm{\mathbf{F},C,j}}^{(t)}\|^{2}$
(19)

where the sequence size $T$ is the number of time steps in each sequence, and
the batch size $B$ is the number of sequences in each batch. In our case,
$T=5$ and $B=12$, leading to a total of 60 samples per batch, matching the
training setup used for the classical auto-encoder
$\mathbf{M}_{\mathrm{A,Mogi}}$. The parameter $\lambda$ is set at 0.01 to
ensure minimal dramatic changes within the 5-day window, reducing noise and
ensuring continuity in the Mogi reconstruction (Fig. 12). The same smoothness
term is applied during the training process of $\mathbf{M}_{\mathrm{B,Mogi}}$.

Figure 12: The smoothness term has enhanced the temporal continuity of Mogi
reconstruction illustrated for AV08. Orange: $X_{\mathrm{GNSS}}$. Blue: Mogi
reconstructions using $Z_{\mathrm{C,Mogi}}$.

###  A.8 Implementation details

Both $X_{\text{S2}}$ and $X_{\text{GNSS}}$ are further split into the train,
validation, and test sets. For $X_{\text{S2}}$, samples from the same
individual site are placed in the same set to avoid data leakage. For
$X_{\text{GNSS}}$, we have saved the GNSS time series from 2006 to 2009 as a
test set to cover the 2008 event for validation. The batch size is set to 64
for loading $X_{\mathrm{S2}}$ and 60 for loading $X_{\mathrm{GNSS}}$. Training
loss is MSE and for $\mathbf{M}_{\text{B, Mogi}}$ and $\mathbf{M}_{\text{C,
Mogi}}$ with an added regularizer to smooth the reconstructed time series from
$\mathbf{M}_{\text{C, Mogi}}$ (eq. 19). We used the Adam optimization, setting
the initial learning rate at 0.001 and weight decay at 0.0001. The maximum
number of epochs is 100 and the learning rate is reduced by a factor of 10
after 50 epochs; training stops if $\mathcal{L}_{\text{MSE}}$ on the
validation set does not improve after an additional 10 epochs. All experiments
were conducted on a system equipped with an AMD EPYC 7702 64-Core Processor
and 1TB of RAM. Training times varied depending on the model complexity.
Specifically, models that require backpropagation through
$\mathbf{F}_{\mathrm{RTM}}$ (such as $\mathbf{M}_{\mathrm{B,RTM}}$ and
$\mathbf{M}_{\mathrm{C,RTM}}$) took approximately 3 hours to train. In
contrast, training times for all other models did not exceed half an hour.

###  A.9 Extended results for bias correction

####  A.9.1 Bias correction for $\mathbf{F}_{\mathrm{RTM}}$

Figure 13: Reconstruction accuracy of $X_{\mathrm{B,RTM}}$ displays clear bias
without bias correction module. Figure 14: Reconstruction accuracy of
$X_{\mathrm{D,RTM}}$ displays clear bias without bias correction module,
showing a similar pattern of bias as Fig. 14. Figure 15: Reconstruction
accuracy of $X_{\mathrm{C,RTM}}$ significantly improved after bias correction.
Figure 16: Reconstructed spectral signatures based on $Z_{\mathrm{C,RTM}}$ and
clustered by species, which display distinct patterns between forest types.
Notably, the bias correction for ‘Larix decidua’ displays a pattern similar to
those of deciduous forest, which is consistent with our findings in Fig. 6(b).

####  A.9.2 Bias correction for $\mathbf{F}_{\mathrm{Mogi}}$

Figure 17: Reconstruction accuracy of $X_{\mathrm{B,Mogi}}$ illustrated for
vertical displacements, which displays clear bias and is more complex than
those shown by $\mathbf{F}_{\mathrm{RTM}}$. Figure 18: Reconstruction accuracy
of $X_{\mathrm{C,Mogi}}$ illustrated for vertical displacements. The bias
correction has improved the reconstruction but its as not as linear as for
RTM, highlighting the trade-off between model simplicity and bias correction
complexity. Figure 19: Reconstructed GNSS signals and corrected biases
illustrated for north displacements for the first half of 12 GNSS stations.
Orange: measured observations $X_{\mathrm{S2}}$. Blue: Mogi output
$X_{\mathbf{F}_{\mathrm{Mogi,C}}}$ (first column); Corrected biases
$B(Q)_{\mathrm{C,Mogi}}=X_{\mathrm{C,Mogi}}-X_{\mathbf{F}_{\mathrm{Mogi,C}}}$
(second column); Corrected output $X_{\mathrm{C,Mogi}}$ (third column).
$X_{\mathbf{F}_{\mathrm{Mogi,C}}}$ clearly captures the transient signals in
2008, characterized by the radial extension along horizontal directions. In
contrast, $B(Q)_{\mathrm{C,Mogi}}$ mainly captures the seasonal components and
noises as suggested in eq. 14, and its temporal variations across stations do
not contain prominent transient signals, which is aligned with our
expectations about the seasonal components and random noises. On the left part
of the black line are validations on the test set, and on the right part are
validations on the training set. The predictions are consistent throughout the
training and test sets. Figure 20: Reconstructed GNSS signals and corrected
biases illustrated for north displacements for the second half of 12 GNSS
stations. Orange: measured observations $X_{\mathrm{S2}}$. Blue: Mogi output
$X_{\mathbf{F}_{\mathrm{Mogi,C}}}$ (first column); Corrected biases
$B(Q)_{\mathrm{C,Mogi}}=X_{\mathrm{C,Mogi}}-X_{\mathbf{F}_{\mathrm{Mogi,C}}}$
(second column); Corrected output $X_{\mathrm{C,Mogi}}$ (third column).
$X_{\mathbf{F}_{\mathrm{Mogi,C}}}$ clearly captures the transient signals in
2008, characterized by the radial extension along horizontal directions. In
contrast, $B(Q)_{\mathrm{C,Mogi}}$ mainly captures the seasonal components and
noises as suggested in eq. 14, and its temporal variations across stations do
not contain prominent transient signals, which is aligned with our
expectations about the seasonal components and random noises. On the left part
of the black line are validations on the test set, and on the right part are
validations on the training set. The predictions are consistent throughout the
training and test sets.

###  A.10 Extended results for the evaluation of physical variables

####  A.10.1 Retrieved biophysical variables of $\mathbf{F}_{\mathrm{RTM}}$

Figure 21: Distributions of variables are tightened up by bias correction
illustrated for $Z_{\mathrm{C,RTM}}$ vs. $Z_{\mathrm{B,RTM}}$. Without bias
correction the variable distributions (red) are implausible, tending to break
out beyond the preset bounded ranges. With bias correction (blue), the model
is able to learn more plausible distributions. Figure 22: Application of the
regressor baseline $\mathbf{M}_{\mathrm{D,RTM}}$ to $X_{\mathrm{S2}}$ leads to
implausible parameter distributions of $Z_{\mathrm{D,RTM}}$ that tend to break
out of the present parameter ranges. Figure 23: $Z_{\mathrm{C,RTM}}$ learned
by our model are plausible and distinguish between forest types. Figure 24:
Pairwise co-distributions of $Z_{\mathrm{D,RTM}}$ learned by the regressor
baseline show no clear physical patterns. Red: coniferous forest. Blue:
deciduous forest. Figure 25: Pairwise co-distributions of $Z_{\mathrm{C,RTM}}$
learned by our model show distinct and plausible physical patterns. Red:
coniferous forest. Blue: deciduous forest. Figure 26: Temporal variations of
inferred physical parameters $Z_{\mathrm{D,RTM}}$ are less consistent over
time. Figure 27: Temporal variations of inferred physical parameters
$Z_{\mathrm{C,RTM}}$. $\mathbf{M}_{\mathrm{C,RTM}}$ effectively captures
distinct, temporally smooth and plausible variations for different forest
types. Table 9: Statistics of $Z_{\mathrm{D,RTM}}$ learned by
$\mathbf{M}_{\mathrm{D,RTM}}$. We have also used $Z_{\mathrm{D,RTM}}$ to
calculate the JM distance between species distributions, which are less
disentangled compared to $Z_{\mathrm{C,RTM}}$.

Species | N | cab | cw | cm | LAI | LAIu | fc
---|---|---|---|---|---|---|---
Pseudotsuga menziesii | 1.43 $\pm$ 0.54 | 38.31 $\pm$ 14.87 | 1.79 $\pm$ 0.25 | 1.91 $\pm$ 1.02 | 2.15 $\pm$ 1.24 | 0.72 $\pm$ 0.34 | 0.61 $\pm$ 0.19
Picea abies | 1.36 $\pm$ 0.44 | 39.49 $\pm$ 13.24 | 1.87 $\pm$ 0.17 | 2.15 $\pm$ 0.98 | 1.72 $\pm$ 1.25 | 0.84 $\pm$ 0.25 | 0.58 $\pm$ 0.13
Pinus nigra | 1.56 $\pm$ 0.55 | 39.50 $\pm$ 15.81 | 1.86 $\pm$ 0.20 | 2.47 $\pm$ 1.13 | 2.32 $\pm$ 1.42 | 0.54 $\pm$ 0.37 | 0.62 $\pm$ 0.19
Larix decidua | 1.50 $\pm$ 0.50 | 38.26 $\pm$ 13.85 | 1.61 $\pm$ 0.42 | 1.43 $\pm$ 0.69 | 1.90 $\pm$ 1.04 | 0.56 $\pm$ 0.40 | 0.67 $\pm$ 0.21
Pinus sylvestris | 1.56 $\pm$ 0.58 | 42.34 $\pm$ 15.68 | 1.90 $\pm$ 0.11 | 1.99 $\pm$ 1.01 | 2.40 $\pm$ 1.25 | 0.67 $\pm$ 0.37 | 0.59 $\pm$ 0.18
Prunus spp | 1.90 $\pm$ 0.84 | 43.22 $\pm$ 22.17 | 1.35 $\pm$ 0.57 | 1.60 $\pm$ 1.04 | 1.79 $\pm$ 0.82 | 0.46 $\pm$ 0.34 | 0.71 $\pm$ 0.25
Fagus sylvatica | 1.98 $\pm$ 0.76 | 46.75 $\pm$ 20.22 | 1.71 $\pm$ 0.42 | 1.21 $\pm$ 0.79 | 2.17 $\pm$ 0.92 | 0.33 $\pm$ 0.37 | 0.82 $\pm$ 0.23
Carpinus betulus | 1.96 $\pm$ 0.71 | 47.55 $\pm$ 19.40 | 1.67 $\pm$ 0.36 | 1.36 $\pm$ 0.84 | 1.91 $\pm$ 0.91 | 0.55 $\pm$ 0.40 | 0.75 $\pm$ 0.21
Quercus spp | 1.91 $\pm$ 0.69 | 48.59 $\pm$ 18.28 | 1.73 $\pm$ 0.41 | 1.46 $\pm$ 0.83 | 2.09 $\pm$ 0.94 | 0.55 $\pm$ 0.40 | 0.75 $\pm$ 0.21
Acer pseudoplatanus | 2.21 $\pm$ 0.75 | 51.34 $\pm$ 20.31 | 1.72 $\pm$ 0.45 | 1.51 $\pm$ 0.79 | 2.01 $\pm$ 0.99 | 0.27 $\pm$ 0.33 | 0.86 $\pm$ 0.18
Fraxinus excelsior | 2.16 $\pm$ 0.71 | 51.67 $\pm$ 18.05 | 1.76 $\pm$ 0.37 | 1.93 $\pm$ 1.17 | 2.05 $\pm$ 1.07 | 0.44 $\pm$ 0.39 | 0.78 $\pm$ 0.21
Alnus glutinosa | 2.22 $\pm$ 0.72 | 55.01 $\pm$ 18.09 | 1.78 $\pm$ 0.34 | 2.03 $\pm$ 0.97 | 2.04 $\pm$ 1.03 | 0.46 $\pm$ 0.41 | 0.76 $\pm$ 0.20

Table 10: Statistics of $Z_{\mathrm{C,RTM}}$ learned by
$\mathbf{M}_{\mathrm{C,RTM}}$. We have also used $Z_{\mathrm{C,RTM}}$ to
calculate the JM distance between species distributions and our model can
learn disentangled representations.

Species | N | cab | cw | cm | LAI | LAIu | fc
---|---|---|---|---|---|---|---
Pseudotsuga menziesii | 1.65 $\pm$ 0.16 | 39.96 $\pm$ 5.33 | 1.00 $\pm$ 0.24 | 2.36 $\pm$ 0.23 | 2.66 $\pm$ 0.43 | 0.49 $\pm$ 0.07 | 0.66 $\pm$ 0.04
Picea abies | 1.53 $\pm$ 0.16 | 43.33 $\pm$ 4.89 | 0.99 $\pm$ 0.24 | 2.37 $\pm$ 0.21 | 2.70 $\pm$ 0.43 | 0.50 $\pm$ 0.07 | 0.68 $\pm$ 0.04
Pinus nigra | 1.75 $\pm$ 0.16 | 42.71 $\pm$ 5.63 | 0.96 $\pm$ 0.19 | 1.97 $\pm$ 0.28 | 2.73 $\pm$ 0.37 | 0.48 $\pm$ 0.06 | 0.64 $\pm$ 0.04
Larix decidua | 1.86 $\pm$ 0.13 | 36.52 $\pm$ 7.50 | 0.88 $\pm$ 0.25 | 2.28 $\pm$ 0.33 | 2.34 $\pm$ 0.48 | 0.49 $\pm$ 0.08 | 0.59 $\pm$ 0.05
Pinus sylvestris | 1.71 $\pm$ 0.14 | 39.87 $\pm$ 5.18 | 0.99 $\pm$ 0.18 | 2.24 $\pm$ 0.24 | 2.68 $\pm$ 0.34 | 0.48 $\pm$ 0.05 | 0.66 $\pm$ 0.03
Prunus spp | 1.93 $\pm$ 0.12 | 34.71 $\pm$ 8.65 | 0.80 $\pm$ 0.30 | 2.26 $\pm$ 0.37 | 2.13 $\pm$ 0.54 | 0.50 $\pm$ 0.10 | 0.53 $\pm$ 0.05
Fagus sylvatica | 1.98 $\pm$ 0.17 | 29.76 $\pm$ 10.71 | 0.94 $\pm$ 0.33 | 2.24 $\pm$ 0.53 | 2.32 $\pm$ 0.69 | 0.43 $\pm$ 0.11 | 0.54 $\pm$ 0.06
Carpinus betulus | 1.90 $\pm$ 0.11 | 33.38 $\pm$ 9.65 | 0.99 $\pm$ 0.31 | 2.37 $\pm$ 0.48 | 2.42 $\pm$ 0.61 | 0.45 $\pm$ 0.10 | 0.56 $\pm$ 0.05
Quercus spp | 1.87 $\pm$ 0.11 | 33.46 $\pm$ 8.55 | 1.06 $\pm$ 0.33 | 2.35 $\pm$ 0.48 | 2.56 $\pm$ 0.63 | 0.44 $\pm$ 0.10 | 0.59 $\pm$ 0.06
Acer pseudoplatanus | 1.99 $\pm$ 0.14 | 30.82 $\pm$ 10.09 | 0.99 $\pm$ 0.34 | 2.12 $\pm$ 0.49 | 2.45 $\pm$ 0.68 | 0.42 $\pm$ 0.12 | 0.54 $\pm$ 0.05
Fraxinus excelsior | 1.90 $\pm$ 0.12 | 34.73 $\pm$ 8.27 | 1.02 $\pm$ 0.32 | 2.19 $\pm$ 0.39 | 2.57 $\pm$ 0.56 | 0.45 $\pm$ 0.10 | 0.58 $\pm$ 0.05
Alnus glutinosa | 1.89 $\pm$ 0.13 | 35.02 $\pm$ 8.37 | 0.99 $\pm$ 0.27 | 2.22 $\pm$ 0.30 | 2.52 $\pm$ 0.48 | 0.46 $\pm$ 0.08 | 0.59 $\pm$ 0.05

####  A.10.2 Reconstructed GNSS signals based on $Z_{\mathrm{C,Mogi}}$

Figure 28: Time series of GNSS displacements illustrated for the first half of
12 GNSS stations, located around the Mogi source. Orange: $X_{\mathrm{GNSS}}$.
Blue: Mogi reconstructions using $Z_{\mathrm{C,Mogi}}$. Red: Mogi
reconstructions fitted by a Kalman filter. Our model successfully reconstructs
the 2008 volcanic inflation transient signals, characterized by radial
expansion and vertical uplift, with displacement magnitudes of 10 mm, aligned
with prior findings [Ji and Herring, 2011, Walwer et al., 2016]. Figure 29:
Time series of GNSS displacements illustrated for the second half of 12 GNSS
stations, located around the Mogi source. Orange: $X_{\mathrm{GNSS}}$. Blue:
Mogi reconstructions using $Z_{\mathrm{C,Mogi}}$. Red: Mogi reconstructions
fitted by a Kalman filter. Our model successfully reconstructs the 2008
volcanic inflation transient signals, characterized by radial expansion and
vertical uplift, with displacement magnitudes of 10 mm, aligned with prior
findings [Ji and Herring, 2011, Walwer et al., 2016].


# Span of the Jones polynomials of certain v-adequate virtual links

Minori Okamura and Keiichi Sakai  Faculty of Mathematics, Shinshu University,
3-1-1 Asahi, Matsumoto, Nagano 390-8621, Japan<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract.

It is known that the Kauffman-Murasugi-Thistlethwaite type inequality becomes
an equality for any (possibly virtual) adequate link diagram. We refine this
condition. As an application we obtain a criterion for virtual link diagram
with exactly one virtual crossing to represent a properly virtual link.

##  1\. Introduction

The Kauffman-Murasugi-Thistlethwaite inequality [4, 9, 10], KMT inequality for
short, is known as an effective tool to estimate, and in some cases determine,
the minimal crossing number of (classical) links in terms of the span of the
Jones polynomial (or equivalently of the Kauffman bracket polynomial). The KMT
inequality is a strict inequality for some links, and the defect is closely
related to the Euler characteristics of the _Turaev surface_ [11] (also known
as the _atom_ [8]). Thus the KMT inequality has a refined form involving the
Euler characteristics of the Turaev surfaces (Theorem 4.1; see also [1, 2]).
Moreover it is known that the refined KMT inequality becomes an equality for
_adequate_ link diagrams. The proofs of these results seem easier than that of
the original KMT inequality, and in fact the refined inequality provides a
simple proof of the Tait conjecture [11].

In this paper we refine this sufficient condition for (possibly virtual) link
diagram under which condition the refined KMT inequality becomes an equality.
We moreover introduce the notion of _v-adequate link diagrams_ to be diagrams
obtained by virtualizing exactly one real crossing of some adequate diagrams,
and as an application we prove that the refined KMT inequality becomes an
equality for v-adequate link diagrams with certain condition. This means that
we can determine the span of the Jones polynomials of an adequate diagram and
a v-adequate one obtained from the former, and this allows us to show that one
of the span of the Jones polynomials of these links cannot be divided by four.
We therefore obtain a recipe for producing properly virtual links. The
“certain condition” is valid for v-adequate diagrams derived from classical
adequate diagrams, and hence we indeed have infinitely many examples to which
our criterion is applied. These examples can be seen as generalizations of
Kishino’s result [6].

This paper is organized as follows. In §2 we review the Kauffman bracket and
the Turaev surface of (possibly virtual) link diagrams. We introduce the
notion of (pseudo-) adequacy in §3, and we prove our main theorems in §4.

##  2\. Preliminaries

We follow the conventions in [5] for virtual link diagrams. Two virtual link
diagrams $D$ and $D^{\prime}$ are said to be equivalent if $D$ can be
transformed into $D^{\prime}$ by a finite sequence of Reidemeister moves and
virtual Reidemeister moves [5, Figure 2].

###  2.1. The Kauffman bracket polynomial

The transformations of a virtual link diagram in a neighborhood of a real
crossing as in Figure 2.1 are called respectively _A-splice_ and _B-splice_.
We draw a dotted line at each spliced real crossing as in Figure 2.1, which we
call a _connecting arc_.

B-spliceA-splice Figure 2.1. A-splice and B-splice, and connecting arcs

Let $D$ be a virtual link diagram. A _state_ $s$ of $D$ is a map from the set
of the real crossings of $D$ to the set $\\{\mathrm{A},\mathrm{B}\\}$. We
denote by $\mathcal{S}$ the set of the states of $D$. For $s\in\mathcal{S}$,
let $D(s)$ be the virtual link diagram obtained from $D$ by performing
$s(p)$-splice at each real crossing $p$ of $D$. Define the _weight_ of
$s\in\mathcal{S}$, denoted by $\left\langle D/s\right\rangle$, by

$\left\langle D/s\right\rangle\coloneqq
A^{\alpha(s)-\beta(s)}(-A^{2}-A^{-2})^{\sharp D(s)-1}\in\mathbb{Z}[A,A^{-1}],$

where $\alpha(s)\coloneqq\sharp s^{-1}(A)$, $\beta(s)\coloneqq\sharp
s^{-1}(B)$ and $\sharp D(s)$ is the number of components of $D(s)$. The
_Kauffman bracket polynomial_ $\left\langle D\right\rangle$ of $D$ is defined
by

$\left\langle D\right\rangle\coloneqq\sum_{s\in\mathcal{S}}\left\langle
D/s\right\rangle.$

For a Laurent polynomial $f\in\mathbb{Z}[A,A^{-1}]$ we denote the maximal
(resp. the minimal) degree of $f$ by $\overline{\deg}\,f$ (resp.
$\underline{\deg}\,f$), and define $\mathrm{span}(f)$ as

$\mathrm{span}(f)\coloneqq\overline{\deg}\,f-\underline{\deg}\,f.$

It is well known that $\mathrm{span}\left\langle D\right\rangle$ is invariant
under the generalized Reidemeister moves. For a virtual link $L$ we define
$\mathrm{span}(L)$ as $\mathrm{span}\left\langle D\right\rangle$ for any
diagram $D$ representing $L$. The following property is well known.

######  Lemma 2.1.

If $L$ is a classical link, then $\mathrm{span}(L)$ can be divided by $4$.

###  2.2. The Turaev surface

Recall the definition of the Turaev surface of a virtual link diagram $D$ with
$c(D)>0$. We replace all the real crossings of $D$ with disks and join these
disks with bands, each of which corresponds to an arc of $D$, so that the
resulting surface $F^{\prime}(D)$ is orientable and is almost embedded into
$\mathbb{R}^{2}$, except for the neighborhoods of virtual crossings (Figure
2.2). The surface $F^{\prime}(D)$ with boundary has been introduced in [3] and
is called the _fat frame_ of $D$.

real crossing$F^{\prime}(D)$virtual crossing$F^{\prime}(D)$ Figure 2.2. Fat
frame $F^{\prime}(D)$

Around the real crossings of $D$, we color $\partial F^{\prime}(D)$ in
checkerboard manner so that the components corresponding to curves obtained by
A-splice (resp. B-splice) is colored by red (resp. blue), as in Figure 2.3.

Figure 2.3. Checkerboard coloring of $F^{\prime}(D)$ around the real crossings
of $D$

We twist each band of $F^{\prime}(D)$ if necessary so that the checkerboard
colorings of $\partial F^{\prime}(D)$ near the real crossings extend
consistently to whole boundary, as in Figure 2.4.

Figure 2.4. A link diagram with unorientable Turaev surface

The resulting surface with boundary is called the _twisted fat frame_ of $D$,
and denoted by $F(D)$. Now let $T_{D}$ be the closed surface obtained by
attaching disks along $\partial F(D)$, and call $T_{D}$ the Turaev surface of
$D$.

######  Remark 2.2.

By definition the Turaev surface of a trivial link diagram is a union of
$S^{2}$.

######  Remark 2.3.

$T_{D}$ is orientable for any classical link diagram [11]. This is not
necessarily the case for virtual diagrams; see Figure 2.4 for example. It is
not hard to see that $F(D)=F^{\prime}(D)$, and hence $T_{D}$ is orientable, if
$D$ is alternating.

######  Definition 2.4.

For a virtual link diagram $D$, let $s_{A}$ (resp. $s_{B}$) be the state of
$D$ that maps every real crossing of $D$ to A (resp. B).

######  Remark 2.5.

By construction, the red boundary of $\partial F(D)$ coresspond to $D(s_{A})$,
and the blue boundary of $\partial F(D)$ coresspond to $D(s_{B})$. Thus we
have a one-to-one corespondence between the components of $D(s_{A})\cup
D(s_{B})$ and those of $\partial F(D)$.

######  Lemma 2.6.

For a virtual link diagram $D$, let $\chi(D)$ be the Euler characteristic of
$T_{D}$. Then we have $\sharp D(s_{A})+\sharp D(s_{B})=\chi(D)+c(D)$.

###### Proof.

$T_{D}$ can be decomposed into a cell complex; 0-cells are in one-one
correspondence to the real crossings of $D$, 1-cells correspond to the arcs of
$D$, and 2-cells corresspond to the disks attached to $F(D)$ along the
boundary. The number of arcs is equal to $2c(D)$, since $D$ can be seen as a
4-valent graph whose vertices are the real crossings of $D$ and whose edges
are the arcs of $D$. Since the number of 2-cells is $\sharp D(s_{A})+\sharp
D(s_{B})$ as mentioned in Remark 2.5, we have

$\displaystyle\chi(D)$ $\displaystyle=c(D)-2c(D)+\sharp D(s_{A})+\sharp
D(s_{B})$ $\displaystyle=-c(D)+\sharp D(s_{A})+\sharp D(s_{B}).\qed$

##  3\. Pseudo-adequate daigrams and v-pseudo-adequate diagrams

######  Definition 3.1.

A virtual link diagram $D$ is said to be _pseudo-adequate_ if, for any
$s\in\mathcal{S}$,

  1. (a)

$\sharp D(s_{A})\geq\sharp D(s_{A}(1))$ and

  2. (b)

$\sharp D(s_{B})\geq\sharp D(s_{B}(1))$,

where $s_{A}(1)$ (resp. $s_{B}(1)$) is a state that sends one real crossing to
B (resp. A) and all the other crossings to A (resp. B); for notations see §4.

######  Remark 3.2.

A virtual link diagram $D$ is said to be _A-adequate_ (resp.  _B-adequate_) if
the condition (a) (resp. (b)) in Definition 3.1 holds after replacing “$\geq$”
with “$>$”, and _adequate_ if $D$ is A-adequate and B-adequate [7]. $D$ is
adequate if and only if the four components of $\partial F(D)$ near each real
crossing $p$ are all distinct.

Clearly an adequate diagram is pseudo-adequate. Any pseudo-adequate classical
diagram is adequate, but this is not the case for virtual diagrams; the
diagram $H_{n}$ in Figure 3.1 is pseudo-adequate but not adequate, because
$\sharp H_{n}(s_{A})=\sharp H_{n}(s_{A}(1))=1$ and $\sharp H_{n}(s_{B})=\sharp
H_{n}(s_{B}(1))=1$.

$\cdots$$H_{n}$$\cdots$$H_{n}(s_{A})$$\cdots$$H_{n}(s_{A}(1))$$\cdots$$H^{\prime}_{n}$
Figure 3.1. $H_{n}$ is pseudo-adequate but not adequate, and $H_{n}$ is
v-pseudo-adequate obtained from $H^{\prime}_{n}$

A _virtualization of a real crossing_ is the replacement of the real crossing
with a virtual crossing.

######  Definition 3.3.

A _v-adequate_ (resp.  _v-pseudo-adequate_) diagram is a diagram $D$ obtained
from an adequate (resp. a pseudo-adequate) diagram $D^{\prime}$ by
virtualizing one real crossing of $D^{\prime}$.

######  Example 3.4.

The diagram $H_{n}$ in Figure 3.1 is v-pseudo-adequate; indeed $H_{n}$ can be
obtained from $H^{\prime}_{n}$ by virtualizing one of its real crossings, and
$H^{\prime}_{n}$ is pseudo-adequate because $\sharp
H^{\prime}_{n}(s_{A})=\sharp H^{\prime}_{n}(s_{B})=2$ while $\sharp
H^{\prime}_{n}(s_{A}(1))$ and $\sharp H^{\prime}_{n}(s_{B}(1))$ are $1$ or
$2$.

######  Definition 3.5.

For $s\in\mathcal{S}$, a connecting arc $\gamma$ in $D(s)$ is said to be
_admissible_ if

  1. (1)

$\gamma$ connect two distinct components of $D(s)$, or

  2. (2)

both the endpoint of $\gamma$ are on a single component of $D(s)$ and any
orientation of the component looks as in Figure 3.2.

(1)(2) Figure 3.2. Admissible connecting arcs

##  4\. The refined KMT inequality

######  Theorem 4.1 ([1], see also [2, 8]).

Let $D$ be a virtual link diagram representing a virtual link $L$, and
$\chi(D)$ be the Euler characteristic of $T_{D}$. Then we have

$\mathrm{span}(L)\leq 4c(D)+2(\chi(D)-2).$

It is known that the equality holds in Theorem 4.1 for adequate diagrams. This
condition can be refined as follows:

######  Theorem 4.2.

The equality holds in Theorem 4.1 for pseudo-adequate diagrams.

The above Theorem 4.2 deduces our main theorem.

######  Theorem 4.3.

Let $D$ be a v-pseudo-adequate virtual link diagram obtained from a pseudo-
adequate diagram $D^{\prime}$ by replacing a real crossing $p$ with a virtual
crossing $vp$. Let $C_{A}$ and $C_{B}$ be the components of $D(s_{A})$ and
$D(s_{B})$ including $vp$ respectively. If any connecting arcs of $D(s_{A})$
(resp. $D(s_{B})$) whose endpoints are both on $C_{A}$ (resp. $C_{B}$) are
admissible, then $\mathrm{span}\left\langle
D\right\rangle=4c(D)+2(\chi(D)-2)$.

######  Corollary 4.4.

Let $D$ be a v-adequate diagram obtained from an adequate classical diagram
$D^{\prime}$. Then $D$ does not represent any classical link.

###### Proof.

Any connecting arcs of $D(s_{A})$ and $D(s_{B})$ whose endpoints are both on
$C_{A}$ and $C_{B}$ (Theorem 4.3) are admissible, because both $C_{A}$ and
$C_{B}$ contain exactly one virtual crossing; see Figure 3.2. Theorem 4.3
therefore implies $\mathrm{span}\left\langle
D\right\rangle=4c(D)+2(\chi(D)-2)$.

On the other hand, because $D^{\prime}$ is adequate (and hence pseudo-
adequate), Theorem 4.2 implies $\mathrm{span}\left\langle
D^{\prime}\right\rangle=4c(D^{\prime})+2(\chi(D^{\prime})-2)$. Since moreover
$c(D^{\prime})=c(D)+1$, $\sharp D^{\prime}(s_{A})=\sharp D(s_{A})+1$ and
$\sharp D^{\prime}(s_{B})=\sharp D(s_{B})+1$, we obtain
$\chi(D^{\prime})=\chi(D)+1$ by Lemma 2.6. Thus

$\displaystyle\mathrm{span}\left\langle D^{\prime}\right\rangle$
$\displaystyle=4(c(D)-1)+2(\chi(D)-1-2)=4c(D)+2(\chi(D)-2)-6$
$\displaystyle=\mathrm{span}\left\langle D\right\rangle-6.$

Since $\mathrm{span}\left\langle D^{\prime}\right\rangle$ is divisible by $4$
by Lemma 2.1, $\mathrm{span}\left\langle D\right\rangle$ cannot be divided by
$4$. Thus $D$ cannot represent any classical link. ∎

######  Remark 4.5.

In his master’s thesis, Kishino [6] proved that the $\mathrm{span}\left\langle
D\right\rangle=4c(D)-2$ if a diagram $D$ is obtained from a proper alternating
classical diagram, and hence $D$ does not represent any classical link.
Corollary 4.4 generalizes his result; any proper alternating diagram is
adequate. Theorem 4.3 has another application to proper alternating virtual
diagrams; see [3].

######  Example 4.6.

The diagram $H_{n}$ in Figure 3.1 is pseudo-adequate and v-pseudo-adequate.
Moreover the endpoints of any connecting arcs of $H_{n}(s_{A})$ (resp.
$K_{n}(s_{B})$) are admissible. Thus both Theorems 4.2 and 4.3 can be applied
to deduce $\mathrm{span}\left\langle
H_{n}\right\rangle=4c(H_{n})+2(\chi(H_{n})-2)$. Since $c(H_{n})=n$ and we can
see that $\sharp H_{n}(s_{A})=\sharp H_{n}(s_{B})=1$, Lemma 2.6 tells us that
$\chi(H_{n})=2-n$ and hence $\mathrm{span}\left\langle H_{n}\right\rangle=2n$.
In particular, if $n$ is odd, then $\mathrm{span}\left\langle
H_{n}\right\rangle$ cannot be divided by $4$ and therefore $H_{n}$ cannot
represent any classical links.

The rest of this paper is devoted to the proofs of Theorems 4.2 and 4.3.

For arbitrarily chosen real crossings $p_{1},\dots,p_{j}$ ($1\leq j\leq
c(D)$), let $s_{A}(j)$ (resp. $s_{B}(j)$) be a state of $D$ that maps
$p_{1},\dots,p_{j}$ to B (resp. A) and the other real crossings to A (resp.
B). Note that any $s\in\mathcal{S}$ other than $s_{A}$ and $s_{B}$ can be
expressed as $s_{A}(j)$ or $s_{B}(j)$. By definition

(4.1) $\displaystyle\overline{\deg}\,\left\langle D/s_{A}(j)\right\rangle$
$\displaystyle=c(D)-2j+2\sharp D(s_{A}(j))-2,$ (4.2)
$\displaystyle\underline{\deg}\,\left\langle D/s_{B}(j)\right\rangle$
$\displaystyle=-c(D)+2j-2\sharp D(s_{B}(j))+2.$

######  Lemma 4.7.

If $j\geq 1$, then we have

$\overline{\deg}\,\left\langle
D/s_{A}(j)\right\rangle\leq\overline{\deg}\,\left\langle
D/s_{A}\right\rangle,\quad\underline{\deg}\,\left\langle
D/s_{B}(j)\right\rangle\geq\underline{\deg}\,\left\langle
D/s_{B}\right\rangle.$

Thus $\mathrm{span}\left\langle
D\right\rangle\leq\overline{\deg}\,\left\langle
D/s_{A}\right\rangle-\underline{\deg}\,\left\langle D/s_{B}\right\rangle$.

###### Proof.

Figure 4.1 shows

(4.3) $\displaystyle\sharp D(s_{A}(j-1))-1\leq\sharp D(s_{A}(j))\leq\sharp
D(s_{A}(j-1))+1.$

A-spliceB-
splice$D(s_{A}(j-1))$$D(s_{A}(j-1))$$D(s_{A}(j-1))$$D(s_{A}(j))$$D(s_{A}(j))$$D(s_{A}(j))$$p$$p$$p$A-spliceB-
spliceA-spliceB-splice Figure 4.1. Proof of (4.3)

By (4.3) we inductively deduce

(4.4) $\sharp D(s_{A}(j))\leq\sharp D(s_{A}(j-1))+1\leq\sharp
D(s_{A}(j-2))+2\leq\dots\leq\sharp D(s_{A})+j.$

Similarly, we have

(4.5) $\sharp D(s_{B}(j))\leq\sharp D(s_{B})+j.$

By definition

(4.6) $\displaystyle\overline{\deg}\,\left\langle D/s_{A}\right\rangle$
$\displaystyle=c(D)+2\sharp D(s_{A})-2,$ (4.7)
$\displaystyle\underline{\deg}\,\left\langle D/s_{B}\right\rangle$
$\displaystyle=-c(D)-2\sharp D(s_{B})+2.$

By (4.4) and (4.5),

$\displaystyle\overline{\deg}\,\left\langle D/s_{A}(j)\right\rangle$
$\displaystyle=c(D)-2j+2\sharp D(s_{A}(j))-2$ $\displaystyle\leq c(D)+2\sharp
D(s_{A})-2,$ $\displaystyle\underline{\deg}\,\left\langle
D/s_{B}(j)\right\rangle$ $\displaystyle=-c(D)+2j-2\sharp D(s_{B}(j))+2$
$\displaystyle\geq-c(D)-2\sharp D(s_{B})+2.$

These inequalities and (4.6), (4.7) imply

$\overline{\deg}\,\left\langle D\right\rangle\leq\overline{\deg}\,\left\langle
D/s_{A}\right\rangle\quad\text{and}\quad\underline{\deg}\,\left\langle
D\right\rangle\geq\underline{\deg}\,\left\langle D/s_{B}\right\rangle.$

Thus we have $\mathrm{span}\left\langle
D\right\rangle\leq\overline{\deg}\,\left\langle
D/s_{A}\right\rangle-\underline{\deg}\,\left\langle D/s_{B}\right\rangle$. ∎

###### Proof of Theorem 4.1.

Lemma 2.6 and Lemma 4.7 imply

(4.8) $\displaystyle\mathrm{span}\left\langle D\right\rangle$
$\displaystyle\leq\overline{\deg}\,\left\langle
D/s_{A}\right\rangle-\underline{\deg}\,\left\langle D/s_{B}\right\rangle$
$\displaystyle=2c(D)+2(\sharp D(s_{A})+\sharp D(s_{B}))-4$
$\displaystyle=2c(D)+2(\chi(D)+c(D))-4$ (4.9)
$\displaystyle=4c(D)+2(\chi(D)-2).\qed$

###### Proof of Theorem 4.2.

Pseudo-adequacy of $D$ implies

(4.10) $\sharp D(s_{A}(1))\leq\sharp D(s_{A}),\quad\sharp
D(s_{B}(1))\leq\sharp D(s_{B}).$

Thus (4.4) and (4.5) can be sharpend in this case as

(4.11) $\displaystyle\sharp D(s_{A}(j))\leq\sharp D(s_{A})+j-1,\quad\sharp
D(s_{B}(j))\leq\sharp D(s_{B})+j-1.$

(4.6), (4.7) and (4.11) imply

$\displaystyle\overline{\deg}\,\left\langle D/s_{A}(j)\right\rangle$
$\displaystyle=c(D)-2j+2\sharp D(s_{A}(j))-2$ (4.12) $\displaystyle\leq
c(D)+2\sharp D(s_{A})-4,$ $\displaystyle\underline{\deg}\,\left\langle
D/s_{B}(j)\right\rangle$ $\displaystyle=-c(D)+2j-2\sharp D(s_{B}(j))+2$ (4.13)
$\displaystyle\geq-c(D)-2\sharp D(s_{B})+4.$

These estimations deduce $\overline{\deg}\,\left\langle
D/s_{A}(j)\right\rangle<\overline{\deg}\,\left\langle D/s_{A}\right\rangle$
and $\underline{\deg}\,\left\langle
D/s_{B}(j)\right\rangle>\underline{\deg}\,\left\langle D/s_{B}\right\rangle$
for any choices of $p_{1},\dots,p_{j}$ ($j\geq 1$), and hence we have
$\overline{\deg}\,\left\langle D\right\rangle=\overline{\deg}\,\left\langle
D/s_{A}\right\rangle$ and $\underline{\deg}\,\left\langle
D\right\rangle=\underline{\deg}\,\left\langle D/s_{B}\right\rangle$. Thus the
inequality in (4.8) becomes an equality. ∎

###### Proof of Theorem 4.3.

In general, it is not hard to see that a connecting arc in $D(s)$
($s\in\mathcal{S}$) corresponing to a real crossing $p$ is admissible if and
only if $\sharp D(s(1))\leq\sharp D(s)$ where $s(1)$ is the state of $D$
obtained from $s$ by changing splice at $p$. Thus if all the connecting arcs
of $D(s_{A})$ and those of $D(s_{B})$ are admissible, then $D$ is pseudo-
adequate.

If $D$ is such a diagram as in Theorem 4.3, all the connecting arcs of
$D(s_{A})$ and $D(s_{B})$ are admissible; because $D^{\prime}$ is
(pseudo-)adequate, the connecting arcs one of whose endpoints is not on
$C_{A}$ nor $C_{B}$ are admissible. By assumption the other arcs are also
admissible. ∎

## Acknowledgements

The authors express their sincere appreciation to Professor Naoko Kamada and
Professor Shin Satoh for their fruitful comments and encouragements. KS is
partially supported by JSPS KAKENHI Grant Number 16K05144.

## References

  * [1] Y. Bae, H. R. Morton, _The spread and extreme terms of Jones polynomials_ , J. Knot Theory Ramifications 12 (3) (2003), 359-373
  * [2] O. T. Dasbach, D. Futer, E. Kalfagianni, X.-S. Lin, N. Stoltzfus, _The Jones polynomial and graphs on surfaces_ , J. Combin. Theory Ser. B 98 (2) (2008), 384-399
  * [3] N. Kamada, _Span of the Jones polynomial of an alternating virtual link_ , Algebr. Geom. Topol. 4 (2004), 1083-1101
  * [4] L. Kauffman, _State models and the Jones polynomial_ , Topology 26 (1987), 395-407
  * [5] L. Kauffman, _Virtual knot theory_ , European J. Combin. 20 (1999), 663-690
  * [6] T. Kishino, _On classification of virtual links whose crossing numbers are equal to or less than 6_ (in Japanese), Master Thesis, Osaka City University (2000)
  * [7] W. B. R. Lickorish and M. B. Thistlethwaite, _Some links with non-trivial polynomials and their crossing numbers_ , Comment. Math. Helv. 63 (1988), 527-539
  * [8] V. O. Manturov and D. P. Ilyutko, _Virtual Knots. The state of the Art_. Series on Knots and Everything, 51, World Scientific Publishing Co. Pte. Ltd., Hackensack, NJ (2013)
  * [9] K. Murasugi, _Jones polynomials and classical conjectures in knot theory_ , Topology 26 (1987), 187-194
  * [10] M. B. Thistlethwaite, _A spanning tree expansion of the Jones polynomial_ , Topology 26 (1987), no. 3, 297-309
  * [11] V. G. Turaev, _A simple proof of the Murasugi and Kauffman theorems on alternating links_ , Enseign. Math. (2) 33 (1987), no. 3-4, 203-225