# Multidimensional surrogate modelling for Airborne TDEM data

Wouter Deleersnyder1,2 David Dudal1,3 Thomas Hermans2
1KU Leuven Campus Kortrijk - KULAK  Department of Physics  Etienne Sabbelaan
53  8500 Kortrijk  Belgium.
E-mail<EMAIL_ADDRESS>
2Ghent University  Department of Geology  Krijgslaan 281 - S8  9000 Gent
Belgium
3Ghent University  Department of Physics and Astronomy  Ghent  Krijgslaan 281
- S9  9000 Gent  Belgium

(Received 2023 XXXX XX; in original form 2023 XXXX XX)

###### Abstract

The computational resources required to solve the full 3D inversion of time-
domain electromagnetic data are immense. To overcome the time-consuming 3D
simulations, we construct a surrogate model, more precisely, a data-driven
statistical model. It is trained on 3D simulation data and predicts the
approximate output much faster. Given the computational cost related to the
simulations, there are limitations in the number of training samples that can
be generated. In addition, certain applications require a wide range of
parameters to be sampled, such as the electrical conductivity parameters in a
saltwater intrusion case. This chapter is therefore limited to a two-layer
model. We construct a surrogate model that predicts the discrepancy between a
1D two-layered subsurface model and a deviation of the 1D assumption. The
latter response is quickly computed with a semi-analytical 1D forward model.
The results are encouraging even with few training samples, but obtaining a
high accuracy is difficult with relatively simple data fit models. We propose
to view the performance in terms of learning gain, representing the gain from
the surrogate model whilst still acknowledging a residual discrepancy.

###### keywords:

Forward modelling, Machine Learning, Surrogate modelling, Electromagnetics,
Airborne

##  1 Introduction

Airborne Electromagnetic (AEM) mapping is an increasingly used method to image
near-surface geological features over large areas via the bulk electrical
resistivity. Common applications are, for example, in mineral exploration
(Macnae and Milkereit,, 2007), contamination (Pfaffhuber et al.,, 2017),
hydrogeological mapping (Mikucki et al.,, 2015), and saltwater intrusion
(Siemon et al.,, 2019; Deleersnyder et al., 2023b, ). Part of the AEM systems
collects time-domain electromagnetic (TDEM) data, but those systems require
significantly more computer power to run 3D modelling simulations, compared to
frequency-domain electromagnetic (FDEM) data. 3D FDEM forward modelling is
becoming applicable for smaller areas (Ansari et al.,, 2017; Cockett et al.,,
2015; Haber et al.,, 2004; Werthmüller et al.,, 2019), while 3D TDEM forward
modelling remains impractical (Börner et al.,, 2015; Haber et al.,, 2007). For
example, Engebretsen et al., (2022) could carry out a 2D TDEM inversion with
418 soundings, which took 30 hours, 97 GB of memory for only 16 iterations.
Auken et al., (2017) argues that the computational demand for a full 3D AEM
survey is still too large so that only few companies are capable of realizing
such inversions.

By consequence, quasi-2D/3D inversion methods, using a much faster and
approximate (1D) forward model, methods are standard for large-scale mappings.
The 1D forward model assumes cylindrical symmetry and is thus limited to
horizontally layered conductivity models. Using that forward model in a
quasi-2D/3D method may yield erroneous interpretations and the results may
better not be interpreted quantitatively. Deleersnyder et al., (2022) provide
an appraisal tool which verifies whether an erroneous interpretation may occur
as a result of the approximate 1D forward model, and indicates areas that
should be reinterpreted with a full 3D forward model.

To overcome the time-consuming simulations, a common data-driven approach is
to construct surrogate models. It is a statistical model that can approximate
the simulation output. A single evaluation of the trained statistical model is
much faster than a simulation. The statistical relation is built from pairs of
full simulation data and its inputs. Surrogate modelling is thus a case of
supervised machine learning. Surrogate modelling and machine learning may
replace 3D forward modelling on a mesh during a 3D inversion, sometimes coined
Fast Forward Modelling (FFM).

The existing literature of using Machine Learning (ML) and 1D forward
modelling is already significant. Bording et al., (2021) use a conventional
Neural Network (NN) for ground-based TDEM data and claims a factor 13 speed up
compared to the accurate 1D forward models with 98% accuracy within the 3%
relative error limit. Other examples of 1D surrogate models are found in
surface wave dispersion (Hou et al.,, 2019), seismic and marine controlled-
source EM data (Puzyrev and Swidinsky,, 2021) and seismic full-waveform
inversion (Wu et al.,, 2018). For AEM modelling, Asif et al., (2022) train a
NN for 1D forward modelling, including the flight altitude as additional
parameter, which adds complexity to the learning problem. Additionally, two
separate NNs were trained, one to predict the forward response and one for the
partial derivatives. A speed-up with a factor 50 is claimed. Wu et al., (2023)
have trained an alternative NN with a long short-term memory (LSTM)
architecture with similar ‘functionalities’ and claim a factor 2700 speed-up.
The step towards 3D modelling has already been made for gravity data (Zhang et
al.,, 2021) or frequency domain magnetotellurics (Conway et al.,, 2019; Peng
et al.,, 2023), which suffer less from the computational burden to generate
training data.

Other work shows that NNs can be used for a direct inversion (Noh et al.,,
2020; Feng et al.,, 2020; Bai et al.,, 2020; Li et al.,, 2020; Puzyrev and
Swidinsky,, 2021; Laloy et al.,, 2021). This avoids the potentially time-
consuming optimization step of the inverse problem, but one loses control over
the specific tunable aspects of the inversion process. The trained NN needs to
handle the non-uniqueness of the geophysical problem, as well as the noise. By
training a separate forward/surrogate model, one retains control over the
inversion: The type of regularization and the regularization parameter(s) can
be easily modified. This allows to use the trained surrogate forward model in
flexible regularization schemes where, for example, the sharpness can be tuned
(Deleersnyder et al., 2023b, ; Klose et al.,, 2022).

While much work has been done in exploring the surrogate modelling and Machine
Learning for 1D forward problems in TDEM problems, the step towards
multidimensional forward modelling, to our knowledge, has not yet been made.
We examine what is involved to train a 3D surrogate forward model. This work
is rather atypical as it will not provide conclusive answers to questions as:
Are the 1D methods straightforwardly applicable to multidimensional problems?
Will we be able ever to create enough training data? Yet, we show encouraging
results in a limited parameter space (two layers) that can be useful for
future multidimensional surrogate modelling endeavours. Likewise, it is a call
for ‘insightful’ learning, where physics and data are paving the way to a
successful, potentially ever generalizable surrogate model. For example,
rather than predicting the 3D data directly, we predict the relative error
between the 3D and 1D data (Deleersnyder et al., 2023a, ). The computation of
the 1D approximate forward model is fast and explains most of the variability
in the geophysical data. The discrepancy, originating from the
multidimensional nature of the subsurface, is predicted with machine learning
tools. We focus on a minimal, well-defined context (a central-loop airborne
system flying at 40 m altitude, covering an area with high conductivity
contrasts, such as when imaging a fresh saltwater interface) and examine what
can be learned from it. Also, we have opted for more classical ML methods,
which are more insightful, rather than using (advanced) neural networks. Also,
demonstrating the performance of a surrogate model for a simplified case
without advanced ML leaves room to use more advanced ML in more extensive
contexts (with more parameters).

We first discuss the data generation in Section 2.1 in terms of low and high
fidelity data, the 1D forward model and 3D simulations respectively. In
Section 2.1.2, we discuss the data transformations used in this work. Inspired
by the characteristics of the transformed data (discussed in Section 3.1) and
preliminary tests, a Gaussian Process Regression (GPR) is picked as surrogate
model (Section 2.2). In Section 3.2, the performance of the trained surrogate
models is described, suggesting successful routes for future work.

##  2 Methods

###  2.1 Data generation

The objective of this work is to construct an ML-based correction model on an
incomplete, but fast evaluable low-fidelity forward model. For TDEM data we
use a 1D semi-analytical model (Hunziker et al.,, 2015), which assumes
horizontal layers, to generate Low Fidelity (LF) data. A Python, open-source
implementation is used (Werthmüller,, 2017) . The correction model is trained
using High Fidelity (HF) data, which must be generated by a forward model that
can accommodate multidimensional structures. The only viable method is through
numerical simulations, described in Section 2.1.1. Data transformations,
before training, are described in Section 2.1.2.

The surrogate correction model predicts the target with the following listed
features:

  * •

The depth $d$ (in meter) of the interface between the two layers,
$d\in[10,50]$,

  * •

the angle $\theta$ (in ∘) of the dipping layer (0° being purely horizontal),
$\theta\in[-75,75]$ and,

  * •

the Electrical Conductivity of the first ($EC1$) and second ($EC2$) layer (in
Sm-1), $EC1,EC2\in[0,0.75]$..

To test the generalizability of this approach, the prior (or the training
data) is not specifically adapted to a specific application. The data points
are sampled uniformly within the feature space within the specified
boundaries.

####  2.1.1 Finite Volume Method simulations to generate High Fidelity data

The finite volume method (FVM) is a numerical discretization technique for
representing and evaluating partial differential equations by conversion to
algebraic equations. The FVM evaluates exact expressions for the average value
of the quantities over some (finite) volume and uses this discretization to
construct an approximate solution. In this work, we resort to the SimPEG
package (Cockett et al.,, 2015; Heagy et al.,, 2017), which is an open-source
project in Python which implements the Maxwell’s equations and the most common
geophysical electromagnetic sources, and the parallel, direct PARADISO-solver
(Schenk and Gärtner,, 2004).

The High Fidelity mesh in this work is a 3D regular rectangular mesh with 288
000 cells. The core ranges from $-50$ to $50$ m in the lateral orientations
and from $-140$ to $60$ m in the vertical orientation with cell sizes of
$5$-by-$5$ m. Outside the core region, cells are progressively increased with
a rate of 20% up to a cell size of $200$ m. The total mesh extent is 2,35 km
in the lateral orientation and 2,80 km in the vertical orientation. The
accuracy of the $HF$ data is illustrated on a two horizontally layered case
with 0.1 and 0.4 Sm-1 and an interface at 50 m depth. The $HF$ and $LF$ data
are shown in Figure 1. It is unavoidable to have a small discrepancy. In this
case the relative discrepancy averages 0.19% per time gate and a maximum of
0.8% at late time gates; see Figure 1.

Figure 1: Illustrative comparison of Low Fidelity ($LF$) and High Fidelity
($HF$) data with horizontal layers for 0.1 and 0.4 Sm-1 and an interface at 50
m depth illustrate the accuracy of the $HF$ mesh. The relative discrepancy
between the $HF$ and $LF$ data is also shown; it indicates larger errors at
later times.

####  2.1.2 Dataset and transformations

The training and test dataset consists of 5090 and 566 samples, respectively.
The target variable $\frac{\partial\mathbf{B}}{\partial t}$ is sampled
logarithmically from $10^{-5}$ to $10^{-2}$ seconds with 100 time steps. The
target output is thus a vector of 100 data points, calling for a
dimensionality reduction, which is described in Section 2.2.1. Before the
dimensionality reduction, the magnetic data is transformed to a relative error
$RE$ between the $LF$ and $HF$ data.

$RE=\frac{HF-LF}{LF}$ (1)

By predicting the error, solely the residual multidimensional variability has
to be learned. Transforming the target to a relative error, rather than an
absolute error, ensures that the target outputs are within the same range, as
the original magnetic data (and consequently the absolute noise) spans
multiple orders of magnitude.

The sensitivity of the forward response does not vary linearly with changes in
the conductivity. Bording et al., (2021) suggests the following normalization
for a better optimization:

$\tilde{EC}=a+\frac{(b-a)\left(\log_{10}EC-\log_{10}(\min
EC\right)}{\log_{10}\max(EC)-\log_{10}\min(EC)},$ (2)

where $\tilde{EC}$ is the normalized electrical conductivity. The minimum and
maximum correspond to the minimum and maximum electrical conductivity in the
training dataset. The values are transformed within the interval $[a,b]$,
which is $[-1,1]$ in our case.

Before training and after the dimensionality reduction of the target, all
features and target outputs are standardized. This is a common requirement for
many machine learning estimators.

###  2.2 Surrogate modelling

####  2.2.1 Dimensionality Reduction

The target consists of 100 variables. However, the smoothness of the
underlying physics implies that neighbouring variables are correlated, i.e.,
the $HF$ data varies smoothly throughout time, which calls for a
dimensionality reduction technique. Dimensionality reduction techniques
transform high-dimensional data into a lower dimensional space, retaining the
most descriptive features, depending on the technique considered.
Dimensionality reduction is often used in machine learning, as it combats the
computational cost and is able to control overfitting (Lopez-Alvis et al.,,
2021).

In this work, the target variables are consecutively transformed to the
functional and principal component space. This is called Functional Component
Analysis (FCA) (Silverman,, 1996; Satija and Caers,, 2015). First, the target
is represented in a functional data space via Functional Data Analysis (FDA)
(Ramsay et al.,, 2006), which represents the original data as a linear
combination of differentiable mathematical functions:

$y(t)\simeq\sum_{k}^{K}c_{k}\phi_{k}(t),$ (3)

where $\phi_{k}(t)$ is a set of basis functions. The true, underlying basis of
our target is unknown. Therefore, we view each functional datum (each sample),
using a non-parametric representation, meaning that each value of the 100
target variables are stored in a finite grid of points. The underlying physics
is smooth and thus the functional data are smooth and one may use
interpolation to evaluate between the grid. With FDA, we represent that data
in a parametric fashion, using the linear combination of Eq. (3). Common
systems rely on the Fourier basis, monomial basis, finite element basis etc.
We work with the B-spline basis. It can represent smooth functions in a sparse
fashion and it is computationally efficient.

We can now reduce the dimensionality, by representing the non-parametric
(original) data in a parametric fashion, i.e. coefficients of the
corresponding spline basis, where the number of coefficients $K<100$ can still
relatively accurately represent the original data.

On those functional data coefficients, Principal Component Analysis (PCA) is
performed, providing a representation in terms of the principal vector basis.
It consists of a linear mapping which maximizes the variance of the data in
that principal vector basis or principal components. Those principal
components are linearly uncorrelated and are ordered by decreasing variance.
By excluding the principal components which explain little variance, the
dimensionality can be further reduced. We denote the number of retained
components with $L$.

Note that FCA is different than the method coined Functional Principal
Component Analysis (FPCA), which retains modes in functional variation of
longitudinal data (Hall et al.,, 2006).

####  2.2.2 Gaussian Process Regression

In Gaussian Process Regression (GPR), the concept of a Gaussian Process (GP)
is used for making predictions. Williams and Rasmussen, (2006) defines a GP as
a collection of random variables, any finite number of which have a joint
Gaussian distribution and is fully specified by second order statistics: a
mean and a covariance function. In a GPR context, the latter is referred to as
the _kernel_. The kernel encodes the type of structure that is found in
samples drawn from the GP prior by defining the “similarity” of two datapoints
combined with the assumption that similar datapoints should have similar
target values. When making a prediction, we evaluate the samples at a given
point, but we must first incorporate the knowledge from the training data
(Scheidt et al.,, 2018; Duvenaud,, 2014; Williams and Rasmussen,, 2006).

Our proposed model is a GPR and relies on the Radial Basis Function (RBF) in
combination with a White Noise kernel. The RBF is a squared-exponential kernel

$k_{\text{RBF}}(x_{i},x_{j})=\exp\left(-\frac{d(x_{i},x_{j})^{2}}{2\ell^{2}}\right),$
(4)

where $\ell$ is the length-scale of the kernel and $d(\cdot,\cdot)$ the
Euclidean distance. It is a default kernel with appealing properties (it is
stationary, i.e., it does not directly depend on the value of $x$, but the
distance between two points, and it has infinitely many derivatives)
(Duvenaud,, 2014). The RBF generates smooth samples in the GP prior and is
generally a good choice if no discontinuities are expected, as in our case.
Increasing the length scale leads to more slowly varying GPs.

The White Noise kernel is included,

$k_{\text{WN}}(x_{i},x_{j})=\left\\{\begin{array}[]{ll}\sigma^{2}&\mbox{if
}x_{i}=x_{j}\\\ 0&\mbox{else}\end{array}\right.$ (5)

meaning that the GPR will explicitly learn the variance $\sigma^{2}$ of the
noise, rather than to solely rely on the Tikhonov regularization (smoothing)
which is assumed on the training point’s covariance matrix111As the covariance
matrix may be ill-conditioned. (controlled by parameter $\alpha$ within the
GPR, it works similarly as the White Noise Kernel, but the variance is not a
free parameter). The main origin of the noise in the dataset is due to the
discretization of the subsurface model for the $HF$ fidelity data. For
example, when the interface at depth $d$ does not coincide with the edges of a
cell, the depth is rounded to coincide with the edge of that cell, causing an
apparent random noise. We are not interested in the noise and want to prevent
overfitting, such that the fitted length scales of the RBF kernel follow the
slowly varying, general trend, rather than the short-scale variations due to
the noise. The latter would limit the interpolating power of the GPR.

We train a separate GPR for each functional component, without using a
multiple-output GPR. The latter could exploit output correlations, however, we
train in the functional component space, where PCA was used to create
uncorrelated functional components. In Section 4, we discuss the simplicity of
using single-output GPRs over one multiple-output GPR which has the advantage
of considering correlation between responses. This means that we will optimize
for each functional component separately. The length scales and the noise-
level will potentially be different for each component.

##  3 Results

###  3.1 Insights from transformed data

We can easily apply the interpretation of the diffusing currents to the
obtained relative error (RE) curves in Figure 2. The peaks in the RE can be
explained by the bulk of the current distribution reaching the interface,
where the subsurface deviates from the horizontally layered models. If we
perform a sweep over the angle, we find that larger angles yield larger maxima
in the RE. A sweep over the depth shows, not unexpectedly, that a deeper
interface yields a maximum at later times, because the diffusing currents
travel longer to reach the discrepancy at the interface. The EC of the first
layer has a large impact on the ‘arrival time’ of these diffusing currents. We
know that electrical currents diffuse faster in lower conductive media, which
is confirmed in the RE curves. The EC in the second layer does not impact the
time of the peaks in RE, but the figure nicely illustrates that the contrasts
between two electrical conductivities are more significant to explain the
value of the RE peaks, rather than the actual electrical conductivities
itself.

Figure 2: Relative error (RE) curves in a parameter study. In each plot, only
one parameter is varied at a time. The fixed parameters were
$\theta=15^{\circ}$, $d=50$ m, EC1 = 100 mSm-1, EC2 = 400 mSm-1.

We need to check to what extent the target can represent the original relative
error after dimensionality reduction. With $L<100$, the number of functional
components will no longer be able to represent all relative errors $RE$
accurately. For this analysis, the dimensionality reduction was trained on the
training dataset and its performance, expressed in Mean Absolute Error (MAE),
tested on the test dataset. The MAE is computed between the original data
($RE$) and the back transform of the reduced data $\tilde{RE}$,

$MAE={\frac{1}{n}\sum_{i=1}^{n}\left|RE_{i}-\tilde{RE}_{i}\right|},$ (6)

where $n$ is the number of samples in the training dataset. The minimum MAE in
this analysis will be the infimum for the surrogate model performance.

With the FDA, the dimensionality is reduced to $K=10,20$ and 30 dimensions,
meaning that we use $K$ B-spline basis functions to represent our original RE.
Then, $L$ functional components are retained after PCA on the coefficients of
the FDA. The MAE is shown in Figure 3 for each $L\leq K$. Also, the MAE for
all $K$ in $[1,30]$, while retaining all functional components, is shown
(black dotted line). This MAE is identical to the MAE as if you would not
perform PCA on the FDA coefficients.

From Figure 3, we observe that an increase in the number of B-spline functions
$K$ can better represent the original relative error RE, as expected. However,
this effect becomes insignificant for $K>20$. Note that the minimum number of
B-splines is 4. This does not hold for PCA. Secondly, we learn that many
functional components are required to represent the original data in a reduced
dimensional space accurately. If $K$ is low, the results suggests to retain
almost all $L$ functional components. For larger $K$, one can gain quite a bit
of dimensionality reduction by using PCA. For example, for $K=30$, one can
easily further halve the dimensions without much loss of accuracy.

We do not decide on the optimal number of B-spline basis functions $K$ and the
number of functional components $L$ using this analysis. It will serve as a
guide to reduce the optimization for these parameters during the training of
the GPR. The MAE of the _predicted_ relative error in our original space will
be the only criterion to determine $K$ and $L$.

Figure 3: Mean Absolute Error (MAE) introduced by working on a reduced
dimensional space (rather than inaccurate prediction of the estimator), for
different retained functional components $L$ of the $K$ B-spline functions.
MAE for $L>30$ is approximately constant and therefore not shown.

###  3.2 Performance

A. B.

Figure 4: The accuracy (in terms of Mean Absolute Error (MAE) ) of the trained
surrogate model on the test dataset, varying the dimensionality reduction
dimensions. A. The accuracy for varying the number of functional basis
functions $K$ and retaining all functional components ($K=L$) and B. The MAE
and the (PCA) explained variance ratio for varying the number of retained
functional components $L$ for $K=16$.

A. B.
C. D.

Figure 5: The predicted and the true Relative Error (RE) for four samples from the test dataset (see Table 1). Table 1: The features from four samples from the test dataset shown in Figure 5. |  $d$ (m) |  $\theta$ (∘) |  $EC1$ (S/m) |  $EC2$ (S/m)
---|---|---|---|---
Case A | 50 | 50.5 | 0.47 | 0.80
Case B | 29 | -16 | 0.69 | 0.89
Case C | 11.6 | -2.33 | 0.91 | 0.77
Case D | 37.5 | 8.3 | 0.016 | 0.72

The performance of the trained surrogate model is measured with the MAE in
terms of (original) data ($RE$) and the back transform of the reduced data
$\tilde{RE}$. The surrogate model is trained for the number of basis functions
$K$ in the functional data analysis and the number of functional components
$L$ to retain. The parameter $\alpha$ controls the Tikhonov regularization
(smoothing) and during optimization, we experienced that a higher value for
$\alpha$ was required for functional components that explain less variance
(apparently, the White Noise Kernel could not handle the apparent noise
properly). A value of $10^{-5}$ was used for the first four functional
components, $10^{-3}$ for the next two and $10^{-2}$ for $l>6$.

Based on our insights from the dimensionality reduction in Figure 3, we have
optimized the surrogate model for ten different $K$s. The MAE of the predicted
data for varying $K$ is shown in Figure 4A, where the maximum number of
functional components $L=K$ is retained. Starting from $K=8$, the prediction
accuracy starts to stabilize with $K=16$ as minimum. In Figure 4B, we show the
MAE in terms of number of functional components $L$ for $K=16$. Starting from
$L=8$, not much is gained from adding more functional components. Note that
this type of optimization leads to a different choice than relying on the
(cumulative) explained variance ratio, which is often done for PCA. For the
remainder of this text, we proceed with $K=16$ and $L=8$.

To illustrate the reduction in the computational burden, the $HF$ data
production requires one hour and 20 minutes of computation time on one 36-core
High Performance Computing node (2 x 18-core Intel Xeon Gold 6140 (Skylake @
2.3 GHz)) and requires 150 GB of RAM. Evaluating the correction model takes
two milliseconds on a laptop. The initialization of a pre-trained GPR requires
30 seconds for a training dataset with 5000 samples (this initialization time
grows with the number of training samples).

In Figure 5, we show some samples illustrating some typical behaviour from the
test dataset. Figure 5A and B are quite successful, the MAE is 0.00455 and
0.0028 respectively, which is in the order of best MAE given the
dimensionality reduction (see Figure 4). The sample from Figure 5C is an
example where the angle $\theta$ is small ($-2.33^{\circ}$). The trend in the
predicted relative error $\tilde{RE}$ is good, however the absolute values are
too small. The surrogate model thus underestimates the true relative error.
The sample in Figure 5D is a curious one, as it does not really resemble the
typical shape that is observed in the other samples and in the sweeps in
Figure 2. Yet, it cannot be ascribed to the dimensionality reduction: the MAE
of the relative error of the back transformed reduced representation is only
0.00040. The reason can be ascribed to the underlying physics, where the
dissipated currents move much more rapidly in more resistive media. This calls
for a more thorough sampling of this part of the parameter space.

A. B.

Figure 6: A. The learning gain using the surrogate model as correction to the
$LF$ model. The black line represents the ordered total MAE of the test
samples. The green area shows the total reduced MAE between the 1D ($LF$) and
3D ($HF$) data and the orange area shows the residual MAE of the RE. The cases
where the surrogate model makes the prediction of the $HF$ data worse are
shown in red. B. The first and second PC coefficients on the feature space.
Samples in the test dataset with the 10% largest MAE are generally found
further from the centre and indicate that the large MAE is a result of being
at on the periphery of the trained parameter space.  Figure 7: The learning
curve suggests that more data could increase the overall accuracy.

We can also view the performance per test sample separately in Figure 6A. The
test samples are ordered from small discrepancies between $LF$ and $HF$ data
to large discrepancies. The overall green part is what one gains by using the
correction/surrogate model to correct the 1D forward $LF$ model. The bottom
orange parts indicate the residual error, which remains after applying the
surrogate model correction. In some cases, the surrogate model corrects the
$LF$ data such that the predicted $HF$ data is worse than when only the $LF$
data was used (red bars), this is in 6% of the cases. Those cases correspond
to low angles, small EC contrasts between layers or a low EC of the first
layer and combinations thereof. No particular cause could be identified. In
Figure 6B we apply PCA to the four input parameters or features. By plotting
the first two coefficients, we show the sampled parameter space in reduced
dimensions. Samples far from the center (the origin) are extreme values in the
sampled data set. The 10% predictions with the largest MAE are usually in the
periphery, suggesting that the limiting extrapolating capabilities of the GPR
in combination with the RBF kernel is the cause of the large MAE. Additional
sampling in those parameter regions could improve the MAE.

Figure 7 shows the learning curve of the training. It shows the MAE of the GPR
accuracy for various training dataset sizes. The MAE was each time computed on
the same test dataset of the above correction model. It is of course
impossible to make predictions for the course of the curve for more training
samples, but it is reasonable to assume that the accuracy can be further
improved by providing more training data.

##  4 Discussion

To our knowledge, we lack a measure to determine the overall performance of
such a correction model. The effect of a modelling error, e.g., on an
inversion model, is still relatively unknown. The approach to look at the
learning gain, i.e. to look at what we gain with such correction model
compared to only relying on a $LF$ model, requires further validation.
However, such a performance measure will be needed. Based on the results
presented, we can conclude that an outstanding accuracy is not achievable with
these straightforward ML techniques, even in this simplified case with two
layers. Other surrogate model approaches, such as projection based method
(Peherstorfer et al.,, 2018), could provide a way out, but require further
research. A second key would be to adjust the prior or feature space to the
specific application, which will require significant $HF$ data for each
project and may not be the most elegant and sustainable approach.

GPR provides relatively simple data fit models, meaning that it interpolates
between observations. It is versatile, due to the option for different
kernels. The RBF (default) kernel works in most smooth cases. It is difficult
to obtain general models, however, and extrapolation is only possible $\ell$
(see Eq. (4)) units away from the trained prior (Duvenaud,, 2014). If one
could identify or construct a more closely related underlying kernel, a more
general GPR can be constructed with some extrapolating power. Note that GPR’s
are not sparse and thus use the whole dataset to perform one prediction. If
one would have the resources to produce millions of training samples,
alternatives to GPR, such as neural networks, should be considered.

We restate our intention to construct a surrogate model without resorting to
advanced ML techniques. In Section 2.2, we have explicitly opted to train a
single-output GPR on each functional component, rather than using a multiple
output GPR. The latter turns out not to be so trivial (Wang and Chen,, 2015;
Liu et al.,, 2018; Constantinescu and Anitescu,, 2013). A truly multi-output
GPR should not only describe the correlation between data points, but also the
correlation between responses (Wang and Chen,, 2015). Hence, it could leverage
information from one output to provide a more accurate prediction for another
output, rather than modeling them individually (Liu et al.,, 2018). However,
multiple-output GPR is still a field of active research and existing codes are
sparse. Constantinescu and Anitescu, (2013) point out the difficulty of
finding good covariance functions, which turn out to be much more challenging
than for single-output GPs.

##  5 Conclusion

We have learned from this simplified multidimensional surrogate modelling case
that it is possible to construct a surrogate model with basic ML concepts
relatively easily. The performance is promising with relatively few training
samples. To further increase the accuracy, significantly more training data is
ultimately needed.

Next steps in future work should consider (1) active learning, to more
knowingly sample the parameter space which maximizes the accuracy gain per
extra training sample, especially because the 3D forward simulations are so
expensive, (2) generating data on a ‘Medium fidelity’ mesh to increase the
number of training samples (while the test dataset should still be generated
on an accurate high-fidelity mesh), (3) examine more sparse dimensionality
reduction techniques (fitting functions to the RE curves did not turn out to
perform better, test were not presented in this work) and (4) examine
projection-based reduction methods instead of data fitting ML models.
Eventually, as a community, we should define meaningful targets to evaluate
the performance of different surrogate models.

## acknowledgments

The research leading to these results has received funding from FWO (Fund for
Scientific Research, Flanders, grant 1113020N and 1113022N) and the King
Baudouin Foundation Ernest du Bois prize 2022. The resources and services used
in this work were provided by the VSC (Flemish Supercomputer Center), funded
by the Research Foundation - Flanders (FWO) and the Flemish Government.

## DATA AVAILABILITY

Data archiving is underway, and it will be stored on Zenodo. Data has been
uploaded as a temporary measure to the submission portal.

## References

  * Ansari et al., (2017) Ansari, S., Farquharson, C., and MacLachlan, S. (2017).  A gauged finite-element potential formulation for accurate inductive and galvanic modelling of 3-D electromagnetic problems.  Geophysical Journal International, 210(1):105–129.
  * Asif et al., (2022) Asif, M. R., Foged, N., Maurya, P. K., Grombacher, D. J., Christiansen, A. V., Auken, E., and Larsen, J. J. (2022).  Integrating neural networks in least-squares inversion of airborne time-domain electromagnetic data.  Geophysics, 87(4):E177–E187.
  * Auken et al., (2017) Auken, E., Boesen, T., and Christiansen, A. V. (2017).  A review of airborne electromagnetic methods with focus on geotechnical and hydrological applications from 2007 to 2017.  Advances in Geophysics, 58:47–93.
  * Bai et al., (2020) Bai, P., Vignoli, G., Viezzoli, A., Nevalainen, J., and Vacca, G. (2020).  (quasi-) real-time inversion of airborne time-domain electromagnetic data via artificial neural network.  Remote Sensing, 12(20):3440.
  * Bording et al., (2021) Bording, T. S., Asif, M. R., Barfod, A. S., Larsen, J. J., Zhang, B., Grombacher, D. J., Christiansen, A. V., Engebretsen, K. W., Pedersen, J. B., Maurya, P. K., et al. (2021).  Machine learning based fast forward modelling of ground-based time-domain electromagnetic data.  Journal of Applied Geophysics, 187:104290.
  * Börner et al., (2015) Börner, R.-U., Ernst, O. G., and Güttel, S. (2015).  Three-dimensional transient electromagnetic modelling using rational Krylov methods.  Geophysical Journal International, 202(3):2025–2043.
  * Cockett et al., (2015) Cockett, R., Kang, S., Heagy, L. J., Pidlisecky, A., and Oldenburg, D. W. (2015).  SimPEG: An open source framework for simulation and gradient based parameter estimation in geophysical applications.  Computers & Geosciences, 85:142–154.
  * Constantinescu and Anitescu, (2013) Constantinescu, E. M. and Anitescu, M. (2013).  Physics-based covariance models for Gaussian processes with multiple outputs.  International Journal for Uncertainty Quantification, 3(1).
  * Conway et al., (2019) Conway, D., Alexander, B., King, M., Heinson, G., and Kee, Y. (2019).  Inverting magnetotelluric responses in a three-dimensional earth using fast forward approximations based on artificial neural networks.  Computers & geosciences, 127:44–52.
  * Deleersnyder et al., (2022) Deleersnyder, W., Dudal, D., and Hermans, T. (2022).  Novel Airborne EM Image Appraisal Tool for Imperfect Forward Modeling.  Remote Sensing, 14(22):5757.
  * (11) Deleersnyder, W., Dudal, D., and Hermans, T. (2023a).  Machine learning assisted fast forward 3D modelling for time-domain electromagnetic induction data–lessons from a simplified case.  In EGU General Assembly 2023, Location: Vienna, Austria.
  * (12) Deleersnyder, W., Maveau, B., Hermans, T., and Dudal, D. (2023b).  Flexible quasi-2D inversion of time-domain AEM data, using a wavelet-based complexity measure.  Geophysical Journal International, 233(3):1847–1862.
  * Duvenaud, (2014) Duvenaud, D. (2014).  Automatic model construction with Gaussian processes.  PhD thesis, University of Cambridge.
  * Engebretsen et al., (2022) Engebretsen, K. W., Zhang, B., Fiandaca, G., Madsen, L. M., Auken, E., and Christiansen, A. V. (2022).  Accelerated 2.5-D inversion of airborne transient electromagnetic data using reduced 3-D meshing.  Geophysical Journal International, 230(1):643–653.
  * Feng et al., (2020) Feng, B., Zhang, J.-f., Li, D., and Bai, Y. (2020).  Resistivity-depth imaging with the airborne transient electromagnetic method based on an artificial neural network.  Journal of Environmental and Engineering Geophysics, 25(3):355–368.
  * Haber et al., (2004) Haber, E., Ascher, U. M., and Oldenburg, D. W. (2004).  Inversion of 3D electromagnetic data in frequency and time domain using an inexact all-at-once approach.  Geophysics, 69(5):1216–1228.
  * Haber et al., (2007) Haber, E., Oldenburg, D. W., and Shekhtman, R. (2007).  Inversion of time domain three-dimensional electromagnetic data.  Geophysical Journal International, 171(2):550–564.
  * Hall et al., (2006) Hall, P., Müller, H.-G., and Wang, J.-L. (2006).  Properties of principal component methods for functional and longitudinal data analysis.  The annals of statistics, 34(3):1493–1517.
  * Heagy et al., (2017) Heagy, L. J., Cockett, R., Kang, S., Rosenkjaer, G. K., and Oldenburg, D. W. (2017).  A framework for simulation and inversion in electromagnetics.  Computers & Geosciences, 107:1–19.
  * Hou et al., (2019) Hou, S., Angio, S., Clowes, A., Mikhalev, I., Hoeber, H., and Hagedorn, S. (2019).  Learn to invert: Surface wave inversion with deep neural network.  In 81st EAGE Conference and Exhibition 2019 Workshop Programme, pages 1–5. European Association of Geoscientists & Engineers.
  * Hunziker et al., (2015) Hunziker, J., Thorbecke, J., and Slob, E. (2015).  The electromagnetic response in a layered vertical transverse isotropic medium: A new look at an old problem.  Geophysics, 80(1):F1–F18.
  * Klose et al., (2022) Klose, T., Guillemoteau, J., Vignoli, G., and Tronicke, J. (2022).  Laterally constrained inversion (LCI) of multi-configuration EMI data with tunable sharpness.  Journal of Applied Geophysics, 196:104519.
  * Laloy et al., (2021) Laloy, E., Linde, N., and Jacques, D. (2021).  Approaching geoscientific inverse problems with vector-to-image domain transfer networks.  Advances in Water Resources, 152:103917.
  * Li et al., (2020) Li, R., Zhang, H., Zhuang, Q., Li, R., and Chen, Y. (2020).  BP neural network and improved differential evolution for transient electromagnetic inversion.  Computers & Geosciences, 137:104434.
  * Liu et al., (2018) Liu, H., Cai, J., and Ong, Y.-S. (2018).  Remarks on multi-output Gaussian process regression.  Knowledge-Based Systems, 144:102–121.
  * Lopez-Alvis et al., (2021) Lopez-Alvis, J., Laloy, E., Nguyen, F., and Hermans, T. (2021).  Deep generative models in inversion: The impact of the generator’s nonlinearity and development of a new approach based on a variational autoencoder.  Computers & Geosciences, 152:104762.
  * Macnae and Milkereit, (2007) Macnae, J. and Milkereit, B. (2007).  Developments in broadband airborne electromagnetics in the past decade.  In Proceedings of Exploration, volume 7, pages 387–398.
  * Mikucki et al., (2015) Mikucki, J. A., Auken, E., Tulaczyk, S., Virginia, R., Schamper, C., Sørensen, K., Doran, P., Dugan, H., and Foley, N. (2015).  Deep groundwater and potential subsurface habitats beneath an Antarctic dry valley.  Nature communications, 6(1):1–9.
  * Noh et al., (2020) Noh, K., Yoon, D., and Byun, J. (2020).  Imaging subsurface resistivity structure from airborne electromagnetic induction data using deep neural network.  Exploration Geophysics, 51(2):214–220.
  * Peherstorfer et al., (2018) Peherstorfer, B., Willcox, K., and Gunzburger, M. (2018).  Survey of multifidelity methods in uncertainty propagation, inference, and optimization.  Siam Review, 60(3):550–591.
  * Peng et al., (2023) Peng, Z., Yang, B., Liu, L., and Xu, Y. (2023).  Rapid surrogate modeling of magnetotelluric in the frequency domain using physics-driven deep neural networks.  Computers & Geosciences, page 105360.
  * Pfaffhuber et al., (2017) Pfaffhuber, A. A., Lysdahl, A. O., Sørmo, E., Skurdal, G. H., Thomassen, T., Anschütz, H., and Scheibz, J. (2017).  Delineating hazardous material without touching - AEM mapping of Norwegian alum shale.  First Break, 35(8).
  * Puzyrev and Swidinsky, (2021) Puzyrev, V. and Swidinsky, A. (2021).  Inversion of 1D frequency-and time-domain electromagnetic data with convolutional neural networks.  Computers & Geosciences, 149:104681.
  * Ramsay et al., (2006) Ramsay, J. O. S., Kotz, C., Read, N., Balakrishnan, B., Vidakovic, and Johnson, N. (2006).  Functional Data Analysis, volume 1.  John Wiley & Sons, Ltd.
  * Satija and Caers, (2015) Satija, A. and Caers, J. (2015).  Direct forecasting of subsurface flow response from non-linear dynamic data by linear least-squares in canonical functional principal component space.  Advances in Water Resources, 77:69–81.
  * Scheidt et al., (2018) Scheidt, C., Li, L., and Caers, J. (2018).  Quantifying uncertainty in subsurface systems, volume 236.  John Wiley & Sons.
  * Schenk and Gärtner, (2004) Schenk, O. and Gärtner, K. (2004).  Solving unsymmetric sparse systems of linear equations with PARDISO.  Future Generation Computer Systems, 20(3):475–487.
  * Siemon et al., (2019) Siemon, B., van Baaren, E., Dabekaussen, W., Delsman, J., Dubelaar, W., Karaoulis, M., and Steuer, A. (2019).  Automatic identification of fresh – saline groundwater interfaces from airborne electromagnetic data in Zeeland, the Netherlands.  Near Surface Geophysics, 17(1):3–25.
  * Silverman, (1996) Silverman, B. W. (1996).  Smoothed functional principal components analysis by choice of norm.  The Annals of Statistics, 24(1):1–24.
  * Wang and Chen, (2015) Wang, B. and Chen, T. (2015).  Gaussian process regression with multiple response variables.  Chemometrics and Intelligent Laboratory Systems, 142:159–165.
  * Werthmüller et al., (2019) Werthmüller, D., Mulder, W. A., and Slob, E. C. (2019).  emg3d: A multigrid solver for 3D electromagnetic diffusion.  Journal of Open Source Software, 4(39):1463.
  * Werthmüller, (2017) Werthmüller, D. (2017).  An open-source full 3D electromagnetic modeler for 1D VTI media in Python: empymod.  GEOPHYSICS, 82(6):WB9–WB19.
  * Williams and Rasmussen, (2006) Williams, C. K. and Rasmussen, C. E. (2006).  Gaussian processes for machine learning, volume 2.  MIT press Cambridge, MA.
  * Wu et al., (2023) Wu, S., Huang, Q., and Zhao, L. (2023).  A deep learning-based network for the simulation of airborne electromagnetic responses.  Geophysical Journal International, 233(1):253–263.
  * Wu et al., (2018) Wu, Y., Lin, Y., and Zhou, Z. (2018).  InversionNet: Accurate and efficient seismic waveform inversion with convolutional neural networks.  In 2018 SEG International Exposition and Annual Meeting. OnePetro.
  * Zhang et al., (2021) Zhang, L., Zhang, G., Liu, Y., and Fan, Z. (2021).  Deep learning for 3-D inversion of gravity data.  IEEE Transactions on Geoscience and Remote Sensing, 60:1–18.


We first review the problem of vanishing reflection at grazing incidence, and
it’s relation to so–called ‘half-bound states’. Consider a plane wave incident
from the left of a dielectric profile $\epsilon(x)$, the profile becoming
homogeneous as $|x|\to\infty$. For zero reflection, the wave must be of the
form $\exp(ik_{x}x)$ on the far left and far right of the profile. Outside the
slab the Helmholtz equation becomes, in the limit of grazing incidence,
Laplace’s equation $d^{2}\phi/dx^{2}=0$ with solutions $a+bx$, where $a$ and
$b$ are constants. But if there is to be zero reflection in this limit,
outside the slab we must have $\exp(ik_{x}x)\to 1+ik_{x}x$. As $k_{x}\to 0$
the wave equals unity everywhere outside the slab. Therefore if there is to be
zero reflection in the limit of grazing incidence the Helmholtz equation must
have two independent solutions, one which is equal to a constant on both sides
of the slab, and the other proportional to the $x$ coordinate, with the same
proportionality constant on both sides of the slab.

Provided the permittivity profile is symmetric about its centre (so that both
$\phi(x)$ and $\phi(-x)$ are solutions), and that one of the solutions to the
TE field is constant outside the slab, we thus ensure zero reflection at
grazing incidence (although non–symmetric potentials can also exhibit less
than complete reflection at grazing nogami1996). In quantum mechanics a state
which is constant outside of a potential where it is ‘bound’ is not
normalizable, and is known as a ‘half-bound state’ (see e.g. page 280 of
newton2002 for a classification of the bound states of potentials). These
states are bound states with an energy that is at the edge of the continuum.
In the reflection from a dielectric profile we can understand these states as
wave–guide modes whose decay constant is close to zero. From this point of
view we are providing a method for designing waveguides with one mode with
vanishing decay constant outside the guide.

The simplest example of a system supporting a half-bound state is a
homogeneous slab waveguide of width $L$, with a constant permittivity
$\epsilon(x)=\epsilon$ that spans a region $x\in[-L/2,L/2]$. Re–introducing
the coordinate $x$, (eq:TEEQ) becomes

$\left[\frac{d^{2}}{dx^{2}}+k_{0}^{2}\left(\epsilon(x)-\epsilon_{b}\right)+k_{x}^{2}\right]\phi(x)=0,$
(1)

The solutions within the two regions are given by $\phi(x)=\exp(\pm{\rm
i}k_{x}x)$ ($|x|>L/2$) and $\phi(x)=\exp(\pm{\rm
i}\sqrt{k_{0}^{2}(\epsilon-\epsilon_{b})+k_{x}^{2}}x)$ ($|x|<L/2$). If the
phase accumulated inside the slab is an integer number of $\pi$, the field on
the right of the slab is $\pm$ the field on the left, and the boundary
conditions can be fulfilled without any reflected wave. This condition is
given by

$\exp\left[\pm{\rm
i}\sqrt{k_{0}^{2}(\epsilon-\epsilon_{b})+k_{x}^{2}}L\right]=\pm 1=\exp({\rm
i}n\pi),$ (2)

a condition which can be fulfilled with a slab of length
$L=n\pi/\sqrt{k_{0}^{2}(\epsilon-\epsilon_{b})+k_{x}^{2}}$, which is the
standard condition for a transmission resonance. If we choose the angle of
zero reflection to be grazing incidence, then the slab waveguide must have
width

$L=\frac{n\pi}{k_{0}\sqrt{\epsilon-\epsilon_{b}}}$ (3)

In this case the field outside the slab reduces to a constant at grazing
incidence, which is the half-bound state just mentioned. fig:Homog shows the
numerically computed reflectivity as a function of angle for a homogenous slab
of permittivity $\epsilon=3$ and background permittivity $\epsilon_{b}=1$
whose length satisfies eq:SlabLength. For the numerical integration we apply
the odeint function from the SciPy library [2020SciPy-NMeth] to integrate Eq.
(1). As the figure shows, the TE reflectivity is zero as we approach
$|\gamma|=|\sin(\theta)|=1$.

[width=0.5valign=m]Figures/HomogRefEpsField.pdf

Figure : (top) Reflectivity as a function of angle $\gamma=\sin(\theta)$ for a
homogeneous slab (permittivity profile inset) in free space
($\epsilon_{b}=1$), whose length obeys eq:SlabLength for $\epsilon=3$,
$k_{0}=1$ and $n=2$. (bottom) the field of a TE wave indicent onto the slab
(boundaries indicated as vertical dashed lines), in the cases where the
initial wave is at normal incidence (blue) and close to grazing (orange). At
grazing incidence, the field decays to a non-zero constant outside the region
where the dielectric exists, indicating that the slab supports a HBS.

##  Factorisation Method

In this section, we use the factorisation method usually applied in
supersymmetric quantum mechanics [Cooper1995] to design a family of HBS
supporting dielectric profiles, which do not reflect TE waves at grazing
incidence. We note that our approach connects with recent work in quantum
mechanics by Ahmed et al. where zero energy reflectionless potentials have
been investigated in several cases [ahmed2017b, dijk2017, ahmed2020], and
where half-bound state supporting potentials have been subject to
supersymmetric transformations to generate e.g. potentials supporting no bound
states at all [ahmed2017]. Here we start by writing the Helmholtz equation (1)
as a product of first order operators

$\left(-\frac{d}{dx}+k_{0}\alpha(x)\right)\left(\frac{d}{dx}+k_{0}\alpha(x)\right)\phi(x)=\hat{a}^{\dagger}\hat{a}\phi(x)=0,$
(4)

where $\hat{a}^{\dagger},\hat{a}$ are analogous to raising and lowering
operators in quantum mechanics (although their commutator is not proportional
to a constant). This factorization can only be carried out for certain forms
of the permittivity $\epsilon(x)$. Assuming grazing incidence $k_{x}=0$,
expanding out the product in (4), and comparing to eq:TEHHDimensionless gives
the relationship between the function $\alpha(x)$, and the material parameter
$\epsilon(x)$

$\epsilon(x)-\epsilon_{b}=\frac{1}{k_{0}}\frac{d\alpha(x)}{dx}-\alpha^{2}(x).$
(5)

Having factorized the Helmholtz equation in the form (4) we can immediately
find the solution where $\hat{a}\phi(x)=0$

$\phi_{1}(x)=\phi_{1}(-L/2)\exp\left[-k_{0}\int_{-L/2}^{x}\alpha(x^{\prime})dx^{\prime}\right].$
(6)

As described earlier, the second solution does not need to be considered. By
inspection of eq:Alpha, one can see that if $\alpha(x)\rightarrow 0$ as
$|x|\rightarrow\infty$, then $\epsilon(x)=\epsilon_{b}$ as required. For
vanishing reflection of grazing incidence waves, the amplitude of the wave
outside of the slab must be the same on either side. This requires
$\int_{-L/2}^{L/2}\alpha(x)dx=0$, a condition which we enforce by choosing
$\alpha(x)$ antisymmetric around $x=0$. Unlike the case of a true bound state
where, for a given $\alpha(x)$, one of the kernels of $\hat{a}$ and
$\hat{a}^{\dagger}$ will diverge at infinity, here the functions satisfying
$\hat{a}\phi=0$ and $\hat{a}^{\dagger}\phi=0$ _both_ represent half bound
states. This means that the operators can be ordered either way in (4) and
both corresponding profiles will support HBS and have vanishing reflection at
grazing incidence. To be precise, both the profile given in Eq. (5) and that
where $\epsilon(x)-\epsilon_{b}=-k_{0}^{-1}d\alpha/dx-\alpha^{2}(x)$ have zero
reflection at grazing incidence. With these conditions, an entire family of
these profiles can be generated. Two such examples are given in
fig:SUSY1,fig:SUSY2. To generate these HBS supporting profiles, we can choose
any odd function $\alpha(x)$, with a value and derivative that goes to zero at
$|x|=L/2$. The permittivity profiles derived from eq:Alpha then give us
designs for dielectric materials that don’t reflect TE waves at grazing
incidence. Interestingly, with $\alpha(x)$ being an odd function, the
permittivity profiles necessarily have regions where
$\epsilon(x)<\epsilon_{b}$. This is immediately evident from (5): after
increasing in magnitude away from $x=-L/2$, $\alpha(x)$ must return to zero at
$x=+L/2$. Therefore at some point the function $\alpha(x)$ must turn around,
having finite value but zero gradient. At this point
$\epsilon(x)-\epsilon_{b}=-\alpha^{2}<0$. Close to grazing the wavevector in
the direction of propagation, $k_{x}$ becomes complex. In e.g. the WKB method,
at the turning point where $\epsilon=\epsilon_{b}$, the wave is matched to a
combination of Airy functions. In the case of a single turning point this
gives rise to complete reflection [heading2013]. The profiles found here have
several turning points and run counter to this intuition, allowing complete
transmission, despite the regions of exponential decay. Nevertheless the
intuition is not completely wrong. The larger the region where $\epsilon(x)$
dips below $\epsilon_{b}$, the larger the reflection is _close_ to grazing
incidence, and the narrower the range of angles close to $|\sin(\theta)|$
where the reflection vanishes. This effect can be seen in fig:SUSY1, where the
comparatively large region of $\epsilon(x)<\epsilon_{b}$, leads to a mirror-
like behaviour above $\gamma\sim 0.8$, and then a very narrow angular range of
low reflectivity around grazing incidence. Conversely in fig:SUSY2 the region
where $\epsilon(x)<\epsilon_{b}$ is relatively small, resulting in a larger
region of low reflectivity close to grazing incidence. It is an interesting
question whether those profiles with a very sharp drop in reflectivity close
to grazing could be realised with sufficient precision (using e.g. dielectric
multilayers or metamaterials) to make this narrow angular region of low
reflectivity observable.

[width=0.5valign=m]Figures/SimonRefEpsField.pdf

Figure : Reflectivity as a function of angle $\gamma=\sin(\theta)$ (top) for a
permittivity ($\epsilon_{b}=1$) generated using (5) with both permittivity and
$\alpha(x)$ inset. In this example there are large regions where
$\epsilon(x)<1$ which would be expected to cause extremely strong reflection
as the angle of incidence increases. The onset of this strong reflection is
evident beyond $\gamma\sim 0.8$ after which the reflectivity approaches unity.
Rather surprisingly as the angle becomes close to grazing, this strong
reflection rapidly decreases to zero. Note that a dip in reflectivity as we
approach grazing incidence is evident for both polarizations. The TE field
(bottom) is shown in the case of normal (blue) and grazing incidence (orange),
demonstrating that—in accordance with our design—at grazing incidence there is
a mode that is constant outside the slab.

[width=0.5valign=m]Figures/DeanRefEpsField.pdf

Figure : As in fig:SUSY1, for a different choice of $\alpha(x)$ profile. In
this case, $\alpha(x)$ is chosen such that the region where $\epsilon(x)<1$ is
reduced. Accordingly the reflectivity reduces to zero over a greater angular
range than in fig:SUSY1, as we approach grazing incidence. Again notice that
the TM polarization also exhibits a dip in reflectivity close to $|\gamma|=1$.