Impr ov ed accuracy of plane-w av e electr omagnetic modelling by application of the total and scattered ﬁeld decomposition and perfectly matched layers

and magnetic ﬁelds by using Dirichlet BC can be expected to partly cancel out in the magnetotelluric transfer functions, for example the impedance tensor. In this work, we quantify this cancellation effect. The inaccuracy is less than typical measurement errors in the vast majority of apparent resistivity and phase data, even, when the primary ﬁelds are strongly inaccurate.


I N T RO D U C T I O N
The plane-wave electromagnetic induction method is a mainstream geophysical technique to image the structure of the Earth.Working in the far-field zone, its sources are either magnetospheric and ionospheric currents (magnetotellurics), the global distribution of lightning discharges (audio-magnetotellurics) or distant radio transmitters (radio-magnetotellurics).Depending on the signal frequencies, plane-wave electromagnetic induction methods have been widely used in engineering and environmental geophysics (e.g.Tezkan 1999 ), groundwater prospecting (e.g.Kalscheuer et al. 2015 ), mineral exploration, oil and gas exploration and in imaging the deep structure of the Earth (e.g.Chen et al. 1996 ;Nelson et al. 1996 ;Becken et al. 2011 ;Dong et al. 2016 ;Zhang et al. 2016 ).Since the distributions of conductivity, dielectric permittivity and magnetic permeability in the Earth can be inferred from electromagnetic induction data, these are considered proxies for geological interpretation.
Accuratel y simulating plane-w ave electromagnetic induction responses is challenging, as the sources are not located inside the modelling domain.Therefore, the precise representation of the electromagnetic fields generated by the remote sources is critical.Conceptually, the source field can be represented by splitting up the total field in two equi v alent w ays.In the first, the total field is the sum of a plane-wave incident field originating from the source and a scattered field generated by interaction of the incident field with the Earth (Ward & Hohmann 1988 ).In the second, the background of a 2-D model is defined as the mean 1-D distribution of material parameters, in which 2-D structures of interest are embedded.Here, the sum of the incident field and the field part scattered by the backg round for m the normal field , which is assumed 1-D and anal yticall y computable, and the field part scattered at 2-D anomalies is called anomalous field .Thus, the total field is the sum of the normal and anomalous fields.Therefore, two approaches are available to simulate magnetotelluric responses.The total-field approach (e.g.Franke et al. 2007 ;Nam et al. 2007 ;Li & Pek 2008 ;Ren et al. 2013 ;Usui 2015 ;Jahandari & Farquharson 2017 ) formulates the Helmholtz equation in terms of the total field and impresses an approximate solution to the total field (typically the normal field) weighted by mesh and material parameters as an equi v alent source at the boundary.This procedure means, that inhomogeneous Dirichlet boundary conditions (BC) for the total field are used.The anomalous-field approach (e.g.Wannamaker et al. 1987 ;Key & Weiss 2006 ;Kordy et al. 2016 ) takes the anomalous field as a variable and impresses an equi v alent source in form of the normal electric field times the anomalous conductivity (or admittivity) in the anomalous model part.This approach involves the approximation of the anomalous field at the boundaries by homogeneous Dirichlet BC.Generally, the use of approximate solutions as BC may lead to unreliable and inaccurate solutions.Fur ther more, for realistic models with profound lateral parameter changes, where the model background is not strictly 1-D, the impressed electromagnetic field cannot be anal yticall y obtained.
When BC do not exactly account for the modelled solution, socalled boundary effects occur.Since we consider the electromagnetic field in frequency domain, we do not formally analyse the dif fusi ve or w ave-like field propagation.Nevertheless, these boundary effects can physically be understood as the scattering of the unaccounted part of the field back into the modelling domain so we refer to them as boundary reflections.In the total-field approach, the used inhomogeneous Dirichlet BC are supposed to excite the plane-wave incident field, account for the scattered field, and simulate an infinite domain beyond its boundary.This is achieved by prescribing the anal yticall y computable normal field at the boundary.Thus, the prescribed normal field is accounted for, which means that it is not reflected, whereas the anomalous field is not accounted for and therefore reflected.Boundary reflections are supposedly minimized by enlarging the domain several skin depths so that the anomalous field decays suf ficientl y tow ards the boundaries (e.g.Jones & Price 1970 ;Jones & Pascoe 1971 ).A conserv ati ve minimum extent toward the side boundaries is the horizontal adjustment distance (e.g.Ranganayaki & Madden 1980 ;Fainberg & Singer 1987 ;Weaver & Dawson 1992 ).The model depth is chosen such that the electromagnetic field is zero at the bottom.Realistic magnetotelluric targets, ho wever , exhibit profound lateral changes in material distribution, for example subduction zones, which include dipping layers, topography and ocean-continent contacts.Typically, to allow for different stratifications at the left and right sides of the model, the corresponding analytic 1-D solutions are prescribed at the respective side, while the field values are interpolated along the top, for example linearly (Pek & Verner 1997 ;Li 2002 ), using cubic splines (Franke et al. 2007 ) or a cosine taper (Wannamaker et al. 1987 ).These approaches assume that the total field varies in a specific way along the top, which may differ from the actual 2-D solution, leading to reflections from the top boundary and thus unreliable and inaccurate solutions for the electric and magnetic fields at the magnetotelluric stations.In the anomalous-field approach, equi v alent problems arise.
Using Dirichlet BC for laterally var ying models, boundar y reflections cannot be avoided by enlarging the modelling domain, because the anomalous field in air does not decay toward the boundary.Thus, the modelled electromagnetic fields may be dominated by interference of outgoing and reflected anomalous fields and depend more on the placement of the boundaries than on the model itself, making them highly unreliable.Those parts of the erroneous boundary reflections that are common to the electric and magnetic fields may cancel out in magnetotelluric transfer functions, such as the impedance tensor.The limited usefulness of comparing modelling results for transfer functions has long been realized.Thus, we focus our analysis on comparing electromagnetic field components.
Efficient absorbing boundary methods such as perfectly matched layers (PML) are an alternative to Dirichlet BC.Introduced by Berenger ( 1994 ), PML add an artificial layer to the boundary which damps the fields, such that no boundary reflections occur.Since PML damp the solution to zero, they can directly replace homogeneous Dirichlet BC in the anomalous-field approach (Rivera-Rios et al. 2019 ;Lei et al. 2022 ), but not inhomogeneous Dirichlet BC in the total-field approach.Ho wever , applying PML to the anomalousfield approach is limited to cases with a 1-D background.
In this work, we implement a numerical method called total and scattered field decomposition (TSFD), which enables the use of PML in plane-wave electromagnetic induction modelling.The TSFD divides the modelling domain into two regions (Jin & Riley 2008 , chapter 5): The total-field region (TFR), where the problem is formulated in terms of total-field values, contains the area of interest and, depending on the case, parts of the outer model areas and boundaries.The scatter ed-field r e gion (SFR), w here the prob lem is formulated in terms of scattered-field values, includes the boundaries, where PML may be applied to absorb the scattered field.At the TSFD interface between these regions, the field incident to the TFR is impressed.
Though, to our knowledge, no studies have implemented the TSFD in geophysical modelling, a few antenna applications assess the method in time domain, both using finite differences (Merewether et al. 1980 ;Umashankar & Taflove 1982 ;Taflove et al. 1995 ;Martin & Pettersson 2000 ;Watts 2003 ) and finite elements (FE; Ledfelt 2001 ;Yao & Wang 2003 ;Riley et al. 2006 ), and in frequency domain (Kirsch & Monk 2002 ;Rumpf 2012 ).In most cases, the scatterer is surrounded by free space, which is in contrast to geophysical models, where the Earth occupies half of the modelling domain and boundaries.Such a setup is addressed by a FE time-domain scheme by Lou et al. ( 2005 ), though only for models with a homogeneous half-space background.Here, the TSFD interface surrounding the area of interest runs partly through air and partly through Earth, intersecting material changes.
We develop two different TSFD setups for specific geophysical settings.The horizontal TSFD is constructed for models with laterall y v ar ying backg rounds.SFR and TFR are separated along a horizontal TSFD interface in air, where the plane-wave incident field is excited (Fig. 1 a).Dirichlet BC prescribe zero at the bottom and the 1-D scattered and normal field solutions of the left and right model sections at the edges to the respective lateral SFR and TFR boundary nodes, while PML at the top allow the 2-D solution to adjust freely without adverse boundary effects (Table 1 ).
To completely remove boundary reflections at high frequencies, we develop the surrounding TSFD , which extends the scheme by Lou et al. ( 2005 ) to allow for horizontally layered backgrounds rather than only homogeneous half-spaces (Fig. 1 b).Here, the normal field is impressed at the TSFD interface and the solution variable in the SFR is the anomalous field ( cf .Table 1 ).
A feature of the TSFD, that remains largely undiscussed in previous studies, is a considerable reduction in computational problem size compared to the standard approach.PML allow to truncate the domain close to the area of interest and thereby save mesh nodes, but  simultaneously add nodes in the boundary layers.In terms of computational cost, PML are meaningful when more nodes are saved than added.
In the quasi-static case, displacement currents are neglected and the skin depth is small compared to the w avelength.Thus, onl y a coarse discretization is necessary for sufficient accuracy (e.g.Pek & Verner 1997 ;Li 2002 ;Key & Weiss 2006 ;Franke et al. 2007 ;Li & Pek 2008 ).In unstructured FE meshes, cell sizes flexibly gro w to w ards the boundaries, so onl y fe w nodes are located in the outer model area, and PML cannot save more unknowns than they add.In the upper radio-magnetotelluric frequency range of 10 kHz-1 MHz, ho wever , displacement currents may need to be accounted for (Kalscheuer et al. 2007 ).Thus, the skin depth is larger and the wavelength is smaller than in the quasi-static case at the same frequency.Hence, as compared to an erroneous quasi-static modelling approach, a larger domain size and a denser discretization are necessary to correctly capture the electromagnetic field.Moreover, the pollution error limits the numerical accuracy of 2-D FE solutions for high-frequency scattering problems (Bayliss et al. 1985 ;Babuska & Sauter 1997 ;Gillman et al. 2015 ).Thus, to retain numerical accuracy, the required discretization per wavelength is even finer at higher frequencies.Shrinking such a relati vel y large and densely discretized domain saves more nodes than PML add.Thereby, the TSFD can substantially reduce the computational problem size at high frequencies.Especially the surrounding TSFD is expected advantageous, as it shrinks the domain from all sides, whereas the horizontal TSFD only reduces the thickness of the air layer.
Our objectives are (1) to accurately solve for electromagnetic fields of laterally varying models by avoiding systematic errors related to boundary reflections, (2) to completely remove boundary reflections in cases where the surrounding TSFD is applicable, (3) to reduce the computational problem size at high frequencies and (4) to e v aluate the modelling errors made by using Dirichlet BC.We introduce the numerical incorporation of the two TSFD setups into a 2-D plane-wave FE frequency-domain modelling code and demonstrate their respective benefits based on a fat dyke model, a quarter-space model, a realistic subduction model and a highfrequency model of two buried anomalies.

2-D Helmholtz equation
Our code defines the y -z plane as the modelling plane in a righthanded coordinate system with the geo-electric strike in x -direction.It accounts for variations of electric conductivity, magnetic permeability and dielectric permittivity and for both conduction and displacement currents.The electromagnetic induction problem is described by the Helmholtz equation, which for the total field in 2-D and in frequency domain (time-dependence e −i ωt ) is We distinguish two modes transverse electric (TE): transverse magnetic (TM): (2) Here, E x [V m −1 ] and H x [A m −1 ] denote electric and magnetic fields, respecti vel y, in strike direction.The impedivity is ˆ z = −iωμ, and the admittivity is ˆ y = σ − iωε, with angular frequency ω [rad s -1 ], magnetic permeability μ [H m −1 ], dielectric permittivity ε [F m −1 ] and electric conductivity σ [S m −1 ].The modelling domain is divided into a mesh of unstructured triangles with constant material parameters within each element.The Helmholtz equation is discretized using node-based FE following Jin ( 2014 , chapter 4).

Dirichlet BCs
In all cases considered, here, we assume the plane-wave source field to be vertically incident.The Dirichlet BC in our code prescribe the total-field solution for a given background model.An incident field of (1 + 0 i ) Vm -1 (TE mode) or (1 + 0 i ) Am -1 (TM mode) is assumed at the top nodes (i.e.high up in the air).Using Wait's recursion formula (Wait 1953 ), the total field for the background model is computed at the lateral edges.Though the incident field at the top has a fixed value of (1 + 0 i ) Vm -1 or (1 + 0 i ) Am -1 , the total field at the top depends on the respective 1-D model that is solved.The model depth is chosen such, that the total field can be approximated to be zero at the bottom of the domain.
If the background of the 2-D model is a horizontally layered halfspace, the 1-D sections at the left and right boundary are identical and the same total-field solution is prescribed at the left and right boundaries.Thus, the total-field solution is constant along the top boundary and equal to that of the top left and top right nodes.
If the 2-D model changes from the left to the right, we compute two analytical solutions separately for the model sections at the left and right boundaries, and apply these solutions at the respective boundary of the 2-D model.As the prescribed total field at the top left and top right corners of the 2-D domain has dif ferent v alues, we interpolate linearly between these two values at the top boundary.The prescribed boundary values along the left and right boundary match the model.At the top boundary, ho wever , the interpolation only approximates the true 2-D solution and is, thus, imperfect.In the lack of better alternatives, this inhomogeneous BC has long been considered the most adequate solution.In this study, we quantify the error related to the associated boundary reflections.

Total and scattered field decomposition
The discretized form of eq. ( 1) is written as a system of equations Ku = r , derived using the weak formulation of the FE problem.The system matrix K contains discretization and model parameters, the vector r the BC, and the vector u the solution at the nodes.Depending on the approach, the solution u in eq. ( 1) is total-field approach, u a anomalous-field approach, mixed: u s in SFR u t in TFR horizontal TSFD, mixed: u a in SFR u t in TFR surrounding TSFD. (3) Here, superscripts denote total ( t ), scattered ( s ), incident ( i ), anomalous ( a ) and normal fields ( n ) and u t = u s + u i = u a + u n .We define u i as the known incident plane-wave source field, u s as the scattering response of the whole Earth, u n as the 1-D normal solution to a horizontally layered background and u a as the anomalous field generated by 2-D anomalies within that background ( cf .Table 1 ).Since the horizontal TSFD is formulated in terms of total and scattered fields (in TFR and SFR, respecti vel y), the solution has a jump of exactly the incident field across the TSFD interface.For the surrounding TSFD, the problem is formulated in terms of total and anomalous fields (in TFR and SFR, respecti vel y), so the solution has a jump of exactly the normal field across the TSFD interface.Elements are completely located in either of the two regions (Fig. 1 c).Nodes at the TSFD interface are defined as total field nodes.Thus, SFR elements touching the TSFD interface have corner nodes in both regions and must incorporate the jump between the two fields.This is achieved by adding the appropriate field jump to the nodes at the TSFD interface and by subtracting it from the SFR nodes adjacent to the TSFD interface.This excites the incident (horizontal TSFD) or normal field (surrounding TSFD) in the TFR, and simultaneously prevents it from entering the SFR (Jin & Riley 2008 ).Moving the known incident or normal field contributions of the involved nodes from the solution vector to the right-hand side vector r , we obtain where u r k = u i k for horizontal TSFD and u r k = u n k for surrounding TSFD.

Perfectly matched layers
Our PML is based on the expression by Sj ögreen & Petersson ( 2005 ).A sequence of layers is added at the boundaries, which successi vel y damps the fields towards zero.In horizontal (top and bottom boundaries) and vertical (left and right boundaries) PML, the damping is governed by the absorption coefficients respecti vel y.Scaling factor α depends on the total thickness of the PML and exponent β governing the damping should be between 2 and 3. Subscripts y and z denote directions.Coordinates y node and z node are those of nodes in the PML and y PML and z PML denote the location of the interface between PML and SFR.At the outermost nodes of the PML, homogeneous Dirichlet BC are applied.Since PML have a finite thickness, eq. ( 1) has to be solved at the nodes inside the PML, using modified coefficients Outside the PML, η y and η z are zero, and the modified equation reduces to the original Helmholtz equation.We use the elemental averages of the terms in brackets in the expressions for g y and g z , but nodal values of η y and η z in the expressions for h .(f-h) relative errors of the solutions computed using Dirichlet BC and horizontal TSFD to the semi-analytic solution for the original model.Note, that the station distribution is not entirely symmetric around the origin.Weaver et al. ( 1985Weaver et al. ( , 1986 ) (Section 3.1 ) for which a semi-analytic solution exists.Secondly, we analyse the impact of the Dirichlet BC on the solution for a quarter-space model (Section 3.2 ) with increasing conductivity contrasts between the two quarter spaces.This model is an extreme case, where 1-D Dirichlet BC fail to correctly prescribe boundary v alues.Thirdl y, we use the horizontal TSFD, to solve a realistic subduction model, comprising strong lateral conducti vity v ariations, complex topo graph y and bath ymetry (Section 3.3 ).To quantify the reduction of the computational problem size with the two TSFD setups at high frequencies, we study the model of two rectangular anomalies embedded in a layered half-space (Section 3.4 ).The models are created using FacetModeller (Leli èvre et al. 2018 ) and the mesh generator triangle (Shewchuk 1996 ).The semi-analytic solutions for the fat dyke model and the quarter-space model are obtained with the code by Weaver et al. ( 1985Weaver et al. ( , 1986 ) ).We found that the correctness of the semi-analytic solution in TEmode strongly depends on the chosen modelling parameters and the receiver spacing.We ensured that the semi-analytic solution used for reference in this work is of significantly higher accuracy than the accurac y achiev ed b y an y finite-element solution, so that the semi-analytic solution can be relied upon as a reference solution.We therefore refer to the deviations between the finite-element solution relative to the reference solution as 'errors'.Specific cases, where the accuracy of the semi-analytic solution fails, are clearly indicated.For more details, see Appendix A .

Fat d yk e model
We verify all boundary methods implemented in our code with the fat dyke model by Weaver et al. ( 1985Weaver et al. ( , 1986 ) ) for which a semianalytic solution can be computed.The model consists of three blocks of different conductivities overlying a perfect conductor.As a first example, we use a modified version of the model, where the left and right block have the same conductivities (Fig. 2 a).This model has a 1-D background and with the boundary suf ficientl y far away from the anomaly (the dyke), the Dirichlet BC solution can be expected to be suf ficientl y accurate.As a second example, we use the original model from Weaver et al. ( 1985 , 1986 ), which has slightl y dif ferent conducti vities in the left and right blocks (Fig. 2 e).In this case, one would expect the horizontal TSFD solution to be superior to the Dirichlet solution.The surrounding TSFD is not suitable for this setup.
We solve the two models with the different boundary methods.The errors of the computed solutions relative to the semi-analytic solutions are shown for the E x -component and the real part of the The imaginary part of B x is zero, thus no meaningful relative error can be calculated, but its absolute deviation from the semi-analytic solution is of the order of 10 −4 .Note, that all fields are normalized with the B -field solution at a large horizontal offset to the left as in Weaver et al. ( 1985Weaver et al. ( , 1986 ) ).Since the uncertainty of good-quality field data is about 2 per cent, the modelled solutions should be at least as reliable as that.The results show that all boundary methods yield suf ficientl y accurate results for both models.For the modified model, the horizontal and surrounding TSFD solutions are slightly more accurate than the Dirichlet solution, presumably, because some small boundary effects still occur using Dirichlet BC, w hich are eradicated, w hen PML are applied instead.A slight asymmetry with respect to the profile midpoint at y = 0 m is observed in the numerical E x -solutions, which, accordingly, results in an asymmetric distribution of the relative errors ( cf .Figs 2 b and c) of no more than 0.1 and 0.2 per cent, for the real and imaginary parts, respecti vel y.This may be an effect of the unstructured mesh, and though manual mesh refinement did not result in a higher level of symmetry, goal-oriented refinement possibly would.The results for the original model are slightly less accurate than for the modified model, but still well within the limit of acceptable accuracy, proving that if the sections at the left and the right of the model do not differ strongly, both the Dirichlet option and the horizontal TSFD option of our code perform well.
Figure 3. Quarter-space model: (a) Sketch of the model, (b-g) relative errors ε R [per cent] of the solutions computed using Dirichlet BC (b, d, f) and horizontal TSFD (c, e, g, cf .Fig. 1 a) to the semi-analytic solution for increasing conductivity contrasts between the quarter spaces (of panel a).Note that the right-hand plot in panel (f) is a zoom-in into the left-hand plot for a better visibility of the individual curves.For σ 1 0.0005 S m −1 , the accuracies of both the semi-analytic and finite-element solutions for E x seem to deteriorate, such that their difference decreases, leading to an apparent (but misleading) decrease of the relative errors in (b) and (d).

Quarter-space model
Here, we analyse how the solution behaves, if the sections at the left and right of the model differ more strongly.For this, we use the quarter-space model (Fig. 3 a, details in Appendix B2 ) and stepwise increases of the conductivity contrast between the two quarter spaces.We solve the model with Dirichlet BC and horizontal TSFD at a frequency of 0.1 Hz and successi vel y lower the conductivity in the left quarter space from 0.1 to 0.000001 S m −1 , while keeping it fixed at 0.1 S m −1 in the right quarter space.For this example, the semi-analytic code for the TE mode does not provide numerical solutions for conductivities below 0.0001 S m −1 ( cf .-e) varies.This can be explained by the fact, that Re( B x ) ought to be constant along the entire profile, while E x varies with y .The Dirichlet solution for B x clearly loses accuracy with increasing conductivity contrasts ( cf .Fig. 3 f), whereas the horizontal TSFD solution for B x (Fig. 3 g) does not show this dependence on the conductivity contrast and remains more accurate at higher contrasts.For E x , the semi-analytic solution is not available for the two largest contrasts.For the remaining cases, a less clear, but systematic dependence of the Dirichlet solution on the strength of the conductivity contrast is observed, but only up to a contrast of two decades.For higher contrasts (from σ 1 = 0.0005 S m −1 ), the relative error of our solution drops, which is probably indicative of an increasingly inaccurate semi-analytic solution, rather than an increasingly accurate numerical solution.Nevertheless, we find the horizontal TSFD solution to be generally more accurate than the Dirichlet solution.For the Dirichlet solution to be as accurate as possible, a large distance to the outer boundaries is required.For the horizontal TSFD, the distance to the PML at the top can be much smaller, which results in a considerably thinner air layer, even including the PML.This does, ho wever , not significantly affect the number of unknowns, as the discretization in the air can be coarse, but about the same number of nodes are saved by shrinking the air layer as are added by the PML.This results in an improved solution accuracy of the horizontal TSFD compared to Dirichlet BC, while keeping the number of nodes approximately the same.
In the next step, we examine how the solution behaves for a large conductivity contrast at different frequencies.We set the conductivities to 0.0001 S m −1 in the left and to 0.1 S m −1 in the right quarter space and solve at frequencies of 10 −5 −10 4 Hz.In order to make the results at the stations at the different frequencies comparable, the interstation spacing and length of the magnetotelluric profile are set depending on frequency and skin depths in the two quarter spaces.Thus, the model scales with frequency.Unfortunatel y, the semi-anal ytic code (E-polarization) is limited as to where the stations can be placed and how they can be spaced and the chosen distribution cannot be realized with it.Also, the semianalytic solution for E x for a quarter-space model cannot be obtained for low frequencies, when conductivities are also low (details in Appendix A ).At higher frequencies, when displacement currents are non-negligible, the semi-analytic solution will be incorrect as well, so for the upper and lower end of the investigated frequency range, no reference solution is available.Even without a reference solution, the reliability of the horizontal TSFD solution can be tested by moving the PML boundaries close to the area of interest.If the PML would not be full y absorbing, boundary ef fects would occur and these would profoundl y af fect the solution at the stations.This effect would be dependent on the distance of the PML to the profile.We therefore solve the model several times with the PML at the top placed at different heights (2, 3 and 10 δ above the surface) and e v aluate the maximum deviation between the three solutions as Here, u j denotes the solution at station j normalized by the solution obtained at a large horizontal offset from the anomalies and receivers that is held fixed in different simulations (normalization sim. to Weaver et al. 1985Weaver et al. , 1986 ) ). Superscripts refer to cases with the PML at different heights (1: 10 δ, 2: 3 δ, 3: 2 δ).Additionally, we calculate the station-average of the maximum deviation per frequency.The mesh is fine enough, such that the solution at the stations has converged and (major) mesh effects are excluded.We find, that the top boundary can be placed close to the area of interest without significantly ( d ≤ 2 per cent) affecting the solution at the stations (Figs 4 a and b).This proves that the PML is wellabsorbing, meaning that the horizontal TSFD solution is largely unaf fected b y boundary reflections.As the only difference between the three models is the domain height, the 2 per cent deviation could result from small boundary reflections at the Dirichlet BC at the sides.These would differ, as the side boundaries have different v ertical e xtents in the three cases.Another reason could be an imperfect PML.In our tests, we found that slightl y dif ferent scaling coefficients α are optimal for different electromagnetic field components.In practice, ho wever , one w ould use one value for α for all computations and this α is then a compromise that optimizes all components.Nonetheless, the variation of the solution is suf ficientl y small at all frequencies, indicating that the horizontal TSFD yields satisfactory results throughout the entire magnetotelluric frequency range.We find that the Dirichlet solution does not deviate more than 2 per cent from the horizontal TSFD solution (Figs 4 c and d) with PML placed 10 δ above the surface throughout the frequency range.This is in line with the findings of the previous experiment, where the horizontal TSFD and Dirichlet BC were found to yield similar results for the conductivity contrast investigated here.The previous experiment showed, that the accuracy of the Dirichlet BC deteriorated at conductivity contrasts of four decades or higher ( cf .F ig. 3 f), w hereas the horizontal TSFD solution was more accurate ( cf .Fig. 3 g).This experiment suggests that the TSFD solution has an advantage over the Dirichlet solution throughout the frequency range, yielding more reliable results at high conductivity contrasts at all frequencies.Thus, the horizontal TSFD can be used to obtain reference solutions for models, where no other reference solution is available.

Subduction model
The subduction model for land and seafloor magnetotelluric data is approximated from the inversion result by Worzewski et al. ( 2011 ) (Fig. 5 , details in Appendix B3 ).The aim of this experiment is to e v aluate the developed TSFD for a complex model.The frequency range is 10 −5 −10 4 Hz.At frequencies ≥1 Hz (TE) and ≥0.1 Hz (TM), seafloor stations are not considered, because the field decays to near zero at the seafloor.We solve the model with horizontal TSFD (Fig. 1 a) and place the PML at 2, 3 and 10 δ above the area of interest.The stationaverage of the maximum deviation between the horizontal TSFD solutions of different heights (Fig. 6 a) is small, meaning that the solution at the stations is not significantly affected by the placement of the PML.Thus, as for the quarter-space model, we can use the horizontal TSFD solution as a reference to e v aluate the Dirichlet solution (Fig. 6 b).For the Dirichlet solution, we observe a stronger difference from the horizontal TSFD solution at higher frequencies.A likely reason for this are stronger boundary effects from the top at higher frequencies.
In the apparent resistivities ρ a and phases , which are computed from the 2-D impedance transfer function E i H j as the effect of the BC on the solution is largely cancelled out because the division of electric by magnetic fields partly removes modelling errors that propagate from one type of field to the other by the linearity of Maxwell's equations.For the purpose of demonstrating this effect, we apply three different options for Dirichlet BC: (i) An analytical 1-D solution to the left and right sections at the respectiv e v er tical boundar y and linear interpolation along the top.
(ii) An analytical 1-D solution to the right model section at all boundaries.
(iii) An analytical 1-D solution to the left model section at all boundaries.
The first option is the same as used throughout this paper, which was already shown to yield relati vel y accurate results, whereas the two additional options are wrong at one side and the top of the computational domain and cannot be expected to function equally well as the interpolation option.Indeed, the error of the electric field values obtained with the two additional Dirichlet BC relative to the horizontal TSFD solution is considerably higher than that of the interpolation variant (Fig. 7 a).This error is solely related to the prescribed boundary values and amounts up to 100 per cent at high frequencies.Nevertheless, for the transfer functions, the errors of the different Dirichlet cases relative to the horizontal TSFD case are significantly smaller than for the primary quantities ( E x , H x ).Assuming a relative error of 2 per cent for the impedance transfer functions gives a limit of 4 per cent relative error for the apparent resistivities and of 1.14 • absolute error for the phases.These limits are exceeded only around the highest peak of the topography (Figs 7 b  and c), which proves that the significant modelling error that we (purposefully) added to Dirichlet cases 2 and 3 largely cancels out in the impedances.Ho wever , it also sho ws that the effect of the BC on these quantities is not completely negligible, and this observation is valid for all Dirichlet BC cases.This demonstrates that the horizontal TSFD yields more reliable results than the standard Dirichlet total-field approach at exposed locations of the topography.
Although in many cases the standard Dirichlet BC based on linear interpolation yields suf ficientl y accurate solutions, the tests for both quarter-space and subduction models indicate that this approach systematically affects the solutions for models with a large conductivity contrast between the left and the right section and for complicated models at higher frequencies, respecti vel y, e ven when the outer boundaries are placed far away from the area of interest.In general, the horizontal TSFD yields more reliable results.Hence, in the absence of a reference solution for a particular model, the TSFD can be used to provide one such reference.

Anomaly model
The fourth model consists of two anomalies in a layered half-space, solved at a frequency of 300 kHz (Fig. 8 a, details in Appendix B4 ).
The air layer at such high frequency must be densely discretized, making the computational prob lem unreasonab ly large.The skin effect in air is negligible, such that the anomalous field is reflected by Dirichlet BC, even when the distance to the outer boundaries is large.Using standard Dirichlet BC, the solution at the stations therefore strongly depends on the distance to the boundaries, and does not converge, even for large domains.
The horizontal TSFD (Fig. 1 a) solution is independent of the distance to the PML so the top may be placed close to the area of interest.The solution converges with increasing distance to the side boundaries, making the air layer thin, but wide.With the surrounding TSFD (Fig. 1 b), all boundaries may be placed close to the area of interest, without adversely affecting the solution.To ensure this, we use one surrounding TSFD solution for a large modelling domain as a reference solution ( cf .Appendix B4 ).We then gradually decrease the domain extent of the surrounding TSFD models and calculate the station-averaged error of those increasingly smaller cases relative to the large surrounding TSFD case.The solution is hardly af fected at all, e ven for the smallest surrounding TSFD case which contains 2820 unknowns (Fig. 8 b).Fur ther more, we calculate the station-av eraged relativ e error of various horizontal TSFD cases for models of decreasing width, again using the large surrounding TSFD case as a reference.The smallest horizontal TSFD case, which yields a suf ficientl y accurate solution (i.e.all components below 2 per cent) contains 11 880 unknowns (Fig. 8 b).Only the imaginar y par t of the B x -component does ne ver full y converge, but deviates no more than 5 per cent from the reference, if the domain is suf ficientl y large.On the contrary, the surrounding TSFD yields the most accurate solution and requires the fewest unknowns, making it the computationally most ef ficient.The adv antage of both TSFD setups over the Dirichlet approach, ho wever , is increased solution reliability and accuracy.

D I S C U S S I O N
The crucial part of the TSFD is the PML at the boundary, so the PML parameters must be carefully chosen to guarantee full absorption.An imperfect PML directly affects the solution.When Dirichlet BC cannot correctly prescribe the solution at the boundary, replacing them by well-working PML yields more accurate and reliable results.Since the PML consists of several layers, it adds unknowns to the computational problem.In terms of computational cost, it only makes sense to apply PML, if the discrete problem can be shrunk more than it is enlarged by the PML.For complicated models with strong 2-D scattering as well as high frequencies, we have shown that using TSFD and PML is far superior to Dirichlet BC in terms of solution accuracy.Thus, there are general inherent advantages to our novel approach.At high frequencies, the surrounding TSFD is most cost-efficient, accurate and reliable, but its application is limited to models with 1-D backgrounds.Since the TSFD interface intersects material changes, the field incident to the TFR is not only the plane-wave source field, but also the response of the layers crossed by the TSFD interface.As the incident field at the TSFD interface is prescribed, it must be anal yticall y computable so only 1-D background models can be allowed.The horizontal TSFD is the most versatile setup for any 2-D model at all frequencies.It prevents reflections at the top boundary.If the side boundaries are placed suf ficientl y far outw ards, artificial scattered fields, that are generated by reflections at the lateral boundaries, are largely absorbed by the PML at the top boundary.
For specific, simpler 2-D models, semi-analytic solutions are available as references to evaluate the accuracy of the developed code.Unfortunatel y, semi-anal ytic solutions cannot be obtained for challenging models, which may incorporate non-uniform station distribution, strong conductivity contrasts or topography, and over a broad frequency range.We used such models to verify our code, but for the complicated test models, we used our own horizontal TSFD solution as reference to test the Dirichlet solution.Being obtained with the same code, both solutions might contain partly identical numerical errors.Ho wever , reference solutions obtained with an independent code would also be affected by boundary reflections, as all other available codes use Dirichlet BC.Therefore, the computed TSFD solutions are the best reference presently available.In this study we choose a limit of 2 per cent measurement uncertainty in the electromagnetic fields.A higher accuracy of the TSFD solutions is achie v able if the PML are adapted to be more absorbing.This would mean adding more layers and thereby more unknowns to the computational problem, but yield a higher accurac y.A natural further dev elopment of this study is the extension into 3-D, where the TSFD is expected to be even more beneficial.

C O N C L U S I O N S
Our objective is to remove the systematic impact of boundary reflections from the total-field approach to magnetotelluric forward modelling for rele v ant geolo gic models b y application of two novel TSFD setups.The surrounding TSFD is suitable for models with horizontally layered backgrounds.At high frequencies, it is shown to significantly reduce the computational problem size and simultaneously yield more accurate and reliable solutions for the electromagnetic field than both the standard Dirichlet approach and the horizontal TSFD.This setup is advantageous in radio-magnetotelluric settings for accurately imaging shallow targets in resistive environments.The horizontal TSFD yields higher solution accuracy and reliability for target areas with backgrounds that exceed horizontally layered half-spaces in complexity.This is the case in many relevant geological settings, such as passive continental margins, subduction zones, mountain chains, topographic plateaus and fault zones.Our setup allows for the first time to solve such models without limiting solution accuracy by the applied BC.This improved solution accuracy is ef fecti ve for targets with both low and high conductivity contrasts to the surrounding geological units.For high contrasts, the anomalous field is stronger leading to stronger boundary reflections, when Dirichlet BC are used.Hence, our combination of TSFD and PML approaches to magnetotelluric forward modelling leads to superior solution accuracy for the electromagnetic field under general premises.Owing to the linearity of Maxwell's equations, modelling errors in the electric and magnetic field propagate to modelling errors in the magnetic and electric field, respecti vel y.Thus, the errors introduced to the electric and magnetic fields partly cancel another in the transfer functions used as input data in inversion.Whereas it has long been assumed that this cancellation is significant, it is first with the introduction of our combined TSFD and PML scheme that it can be accurately quantified.Considering the high level of inaccuracy of electric and magnetic field components generated by Dirichlet BC it is surprising (and reassuring to the practitioner) to see that the inaccuracy in the transfer functions is smaller than typical measurement errors in the majority of data considered in our subduction zone e xample.Nev ertheless, with more sophisticated field instrumentation and time-series processing becoming av ailable, highl y accurate forw ard modelling codes, such as our TSFD and PML scheme, are needed to harness the full potential in detecting small-scale or low-contrast anomalies, which have responses at the brink of measurement accuracy.

A C K N O W L E D G M E N T S
The computations were enabled by resources in projects SNIC 2018/8-167, SNIC 2020/15-265 and SNIC 2021/22-883 provided by the Swedish National Infrastructure for Computing (SNIC) at UPP-MAX, partially funded by the Swedish Research Council through g rant ag reement no.2018-05973.This work w as co-financed b y The Swedish Foundation for International Cooperation in Research and Higher Education (STINT) under contract number CH2017-7233 and the National Natural Science Foundation of China (41922027).We thank Peter Leli èvre for implementing the option to set boundary markers in FacetModeller.

D ATA AVA I L A B I L I T Y
The software used for the simulations in the study is available at Zenodo via 10.5281/zenodo.6700659under the GNU LGPL v3.0.

R E F E R E N C E S A P P E N D I X A : A C C U R A C Y O F T H E S E M I -A N A LY T I C S O L U T I O N
As the code to compute reference solutions for the fat dyke and the quarter-space models is not full y anal ytic, but semi-anal ytic (computed by an approximative scheme based on analytic formulations), we analyse its accuracy to ensure its reliability as a reference.The semi-analytic solutions of the TM and TE modes are obtained by the B-polarization formulation (Weaver et al. 1985 ) and the E-polarization formulation (Weaver et al. 1986 ), respecti vel y.The solutions can be obtained for a three-segment model of thickness t (in Weaver et al. 1985Weaver et al. , 1986 , t is called d ) overlying a perfect conductor, with the half width a of the middle block.If a = 0 km, the model consists of only two blocks and if t → ∞ , it becomes a quarter-space model.The conductivities of the blocks are denoted σ 1, 2, 3 from left to right.
For testing the accuracy, we have to compare the results to a known (and highly accurate) solution.For a homogeneous halfspace, the solution at the surface can be calculated completely anal yticall y.We obtain the semi-analytic solution with the threesegment option ( a = 10 km) and with the two-segment option for a homogeneous half-space by setting σ 1 = σ 2 = σ 3 = 0.1 S m −1 and t = 200 km, which corresponds to 40 skin depths at a frequency of 0.1 Hz and can be regarded as t → ∞ .The solution of the Epolarization is computed at nodes between y = −v 1 and y = v 1 with a spacing of d 1 .In the inner part of the profile between y = −v 0 and y = v 0 a finer spacing of d 0 can be chosen.We set the integration limits to v 0 = 10 km and v 1 = 110 km and test different spacings d 0 and d 1 (Fig. A1 ).
The homogeneous half-space solution is 1-D and thus constant along the profile.Ho wever , the semi-analytic solution in Epolarization varies in the vicinity of the contacts between the blocks, even though all blocks are set to the same conductivity, so that, physically, no contacts exist.We generally find, that the solution accuracy depends on the node spacing and that the error is larger for the quarter-space model than for the three-segment model.If the spacing is too large, the solution jumps spuriously.With finer spacing, the inaccuracy is characterized by less severe jumps and a smoother variation.With even finer spacing, the inaccuracies around the contacts become small enough to be neglected.
We repeat the test for the quarter-space model with a reduced conductivity of 0.05 S m −1 in the left block.To keep the approximative t → ∞ we set t = 280 km, all other parameters remain the same.When the conductivities of the blocks are not identical, we have no reference to compare the semi-analytic results to.Nevertheless, the curves seem to be dependent on the spacing and for some cases there are clearly jumps in the solution (Fig. A2 ), which contradict the continuity of the electromagnetic field components.
To explain this, we have to take a closer look at how the semianalytic solution is obtained.At each observation point, the general solution is calculated based on the Green's functions for each block.This general solution is the sum of the particular solution, associated with the block on which the current observation point is located, and the auxiliary solutions associated with the two other blocks ( cf .eq. 14 in Weaver et al. 1986 ).The Green's functions consist of an integration along y , which is limited by −∞ and a (block 1), −a and a (block 2) a and ∞ (block 3).The integral limits that are −∞ and ∞ in the analytical expression become finite ( −v 0 and v 0 ), when e v aluated numericall y.Each of these inte grals involv es a Fourier cosine series expansion that sums over the vertical wave number from zero to ∞ .This sum has to be cut at a finite index.The results can be correct, if the successi vel y summed terms decrease in size, such that the series expansion con verges.Ho wever , turning infinite bounds to finite ones can result in inaccuracies.As also noted in Weaver et al. ( 1986 ), the (complicated) analytic expressions for the auxiliary solutions exhibit logarithmic singularities at observation points located exactly at the contacts to the neighbouring blocks (e.g. for an observation point at y = −a on block 2, the auxiliary solution from block 1 becomes singular, vice versa for an observation point at y = −a on block 1, the auxiliary solution from block 2 becomes singular and likewise at the contact between blocks 2 and 3).Singularity extraction techniques (addition of an analytic term to the singular integrand and compensation for this addition to the integrand by subtracting a corresponding integral) are used to remove this singularity, and, as stated in Weaver et al. ( 1986 ) and as our tests show, the solution at the contacts is not singular.Nevertheless, in the vicinity of the 'removed' singularity, the solution clearly deteriorates towards the contact ( cf .F igs A1 , A2 ), w hich suggests, that the solution is less stab le near the contacts, if the integration increments are not suf ficientl y small.Another source for inaccuracies could be round-off errors or similar numerical effects, or the applied numerical integration method (Gauss-Laguerre and Simpson quadratures).Since the observed inaccuracy depends on the spacing and placement of the integration (or quadrature) points, it is likely an effect due to the numerical integration.
As a result of the integration method, the fields can be obtained at locations which are multiples of the spacing d 0 (in the inner part of the profile) or d 1 (in the outer part of the profile), but not at arbitrary locations.
An important note regards the calculation of the two-segment model.Internall y, the semi-anal ytic code changes the half width of the middle block from zero to a = 1 km and assigns the middle and right block the same conductivity.Then, the origin is shifted 1 km to the right.So, essentially, a three-segment model is solved internally, but the right segment consist of two segments of identical conductivity and the user-intended contact between the blocks is located at y = 0 km.Here, attention has to be paid to the desired station locations, because the entire profile, including the input integration limits of v 0 and v 1 are shifted by + 1 km.
To conclude, the integration limits v 0 and v 1 and increments d 0 and d 1 have to be chosen carefully.For computing the reference solutions for this work, we used v 1 = 170 km, v 0 = 10 km, d 1 = 2 km and d 0 = 0.5 km, as these gave the most accurate results in our test.
Another notable limitation regards the maximum depth that can be set in the E-polarization code, which is rele v ant when solving a quarter-space model, where t → ∞ .In that case, t has to be set as a multiple (e.g. 5 or 10) of the skin depth.Assuming (for example) a frequency of 0.1 Hz, the maximum t that can be chosen is 900 km (otherwise, the code yields not-a-number as results).This corresponds to just above six skin depth in a medium with a conductivity 0.0001 S m −1 .For lower conductivities, the perfect conductor at depth t would affect the results, and the considered model would not be a quarter-space model.In the B-polarization, this problem does not occur.

A P P E N D I X B : D E TA I L E D M O D E L D E S C R I P T I O N S
Here, we provide full details of the tested models and the performed computations.

B1 Fat d yk e model
The fat dyke model was used once in a modified form (Table B1 ,

B3 Subduction model
The subduction zone model (Fig. 5 ) is approximated from the inversion result obtained by Worzewski et al. ( 2011 ).The conductivity of the air is σ = 10 −8 Sm −1 .The conductivity of the Earth ranges between σ = 2 × 10 −4 and 1 S m −1 , and the sea water has a conductivity of σ = 3.3 S m −1 .The relative magnetic permeability and relati ve dielectric permitti vity are both one and constant throughout the model.The area of interest of the subduction model extends from y = 0 to 400 km and from z = −10 to 350 km.The resolution of the FE mesh is suf ficientl y fine for the solution to be independent of the mesh.Especially around the stations, the mesh is strongly refined.The permitted interior angles of the triangles are between 30 • and 120 • .
When Dirichlet BC are applied at all outer boundaries, these are placed at a distance of 10 skin depths in the most resistive Earth structure away from the area of interest.The most resistive Earth structure has a conductivity of σ = 2 × 10 −4 S m −1 .The used frequencies are 10 −5 −10 4 Hz.In Table B5 , the skin depth δ and the padding distance of 10 δ are listed for each frequency f .
When the model responses are solved for with the horizontal TSFD ( cf .Fig. 1 a), the distances between area of interest and the left, right and bottom outer boundaries are 10 δ.PML are applied at the top, and the model is solved three times, with the top boundary at different heights.The interface between domain and PML is placed 2, 3 and 10 δ above the top of the area of interest.The TSFD interface is in each case placed at exactly half this distance (1, 1.5 and 5 δ).The PML in our computations al wa ys contains 10 layers and has a total thickness of t PML = 2 δ, which is listed in Table B5 for each individual frequency.The damping of the PML is governed by scaling factor α, which is listed in Table B5 as well, and exponent β, which is al wa ys set to β = 2 in our computations.

B4 Anomaly model
The anomaly model (shown in Fig. 8 a) comprises two rectangular anomalies of σ = 0.1 S m −1 embedded in a layered half-space with a 10 m thick top layer of σ = 10 −3 S m −1 , and an underlying half-space of σ = 10 −4 S m −1 .The air has a conductivity of σ = 10 −8 S m −1 .The relative magnetic permeability and relative dielectric permittivity are both one and constant throughout the model.The first anomaly has an extent of −10 ≤ y ≤ 10 m and 0 ≤ z ≤ 10 m and the second one of −30 ≤ y ≤ −20 m and 5 ≤ z ≤ 25 m.A profile of 21 stations is spread out between y = −100 and 100 m at the Earth's surface at z = 0 m.The model responses were solved for at a frequency of 300 kHz, where the skin depth in the underlying half-space is δ ≈ 100 m.The resolution of the FE mesh is suf ficientl y fine for the solution to be independent of the mesh.Especially around the stations and in the air layer, the mesh is dense.The permitted interior angles of the triangles are between 30 • and 120 • .
For the traditional approach using Dirichlet BC at all outer boundaries, the solution does not converge with increasing padding distance towards the boundaries, and due to the fine discretization in the air, the computational problem becomes unreasonably large.
For the horizontal TSFD ( cf .

Figure 1 .
Figure1.Sketches of (a) horizontal TSFD for laterall y v arying models, with a horizontal interface between total (TFR) and scattered field region (SFR), and PML at the top, (b) surrounding TSFD for models with anomalies in a horizontally layered background with PML at all boundaries and (c) FE mesh at the TSFD interface with elements sharing nodes in the TFR and SFR.Different colours in (a) and (b) represent different conductivity structures.

Figure 2 .
Figure 2. Fat dyke model by Weaver et al. ( 1985 , 1986 ): (a, e) Sketches of the modified and original versions of the model (not true to scale).(b-d) relative er rors ε R [per cent] of the solutions computed using Dirichlet BC, horizontal TSFD ( cf .Fig. 1 a) and sur rounding TSFD ( cf .Fig. 1 b) to the semi-analytic solution for the modified model and(f-h) relative errors of the solutions computed using Dirichlet BC and horizontal TSFD to the semi-analytic solution for the original model.Note, that the station distribution is not entirely symmetric around the origin.
Appendix A ).The relative errors of our Dirichlet BC and horizontal TSFD solutions to the semi-analytic solution are shown in Figs 3 (b)-(g) and allow a number of observations.The relative error of the B x -component (Figs 3 f and g) remains almost constant along the profile, whereas the relative error of E x (Figs 3 b

Figure 4 .
Figure 4. Quarter-space model (Fig. 3 a, with σ 1 = 0.0001 S m −1 ): (a) Maximum deviation d (eq.7 ) between three different horizontal TSFD cases (PML at 2, 3 and 10 δ above area of interest, cf .Fig. 1 a) at a frequency of 10 Hz, (b) station-averaged deviation between the three different horizontal TSFD cases, (c) relative error of the Dirichlet solution to horizontal TSFD solution with PML at a height of 10 δ at 10 Hz and (d) station-averaged relative errors of the Dirichlet solution relative to the horizontal TSFD solution with PML at a height of 10 δ.

Figure 5 .
Figure 5. Subduction model: The 398 stations are not plotted, b ut distrib uted at 1 km spacing along the seafloor and the Earth's surface.

Figure A1 .Figure A2 .
Figure A1.Semi-analytic TE-mode solution for a homogeneous half-space, obtained with the three-segment model (left-hand columns) and two-segment model (right-hand columns) with equal conductivities in all segments and different spacings d 0 and d 1 between the integration points (in the rows, see text).The analytic solution is shown in red for comparison.d 0 = 0.5 km d 1 = 2 km d 0 = 0.5 km d 1 = 10 km d 0 = 10 km d 1 = 10 km Figs 2 (a)-(d)) and in the original form by Weaver et al. ( 1985 ) (Table B2 , Figs 2 (e)-(h)).
Fig. 1 a), the bottom boundary is placed at z = 500 m.At the top, a 10-layer PML with a total thickness of 100 m is applied with the PML-domain interface at different heights between z = −20 and −500 m.The PML coefficients are α = 3250 and β = 2 and the TSFD interface is placed at z = −10 m.The distances to the left and right boundaries are varied.The smallest case in terms of computational problem size with a suf ficientl y accurate solution contains 11 880 unknowns.This case has a horizontal domain extent of −2650 ≤ y ≤ 2650 m and a PMLdomain interface height of z = −100 m.When this height is smaller, the horizontal domain extent must become even larger for sufficient accuracy, leading to larger problem sizes.When solving the model with the surrounding TSFD ( cf .Fig. 1 b), PML are applied at all outer boundaries.All PML have a total thickness of 100 m and contain 10 layers with α = 3250 and β = 2.For the smallest case of sufficient accuracy, the interfaces between PML and domain are located at z = −100 m (top), z = 50 m (bottom), y = −110 m (left) and y = 110 m (right).The TSFD interfaces are placed at z = −10 m (top), z = 40 m (bottom), y = −105 m (left) and y = 105 m (right).This case has a computational problem size of 2820 unknowns.See Fig. 8 for a comparison of the solution accuracy and the problem size of the smallest horizontal TSFD and surrounding TSFD cases.Sufficient accuracy is reached, when the mean deviation of the solution at the stations relative to a reference remains below 2 per cent.As a reference we use a solution that is converged towards by all computations as the domain extent is increased.It is obtained using the surrounding TSFD setup with the interfaces between PML and domain located at z = −100 m (top), z = 50 m (bottom), y = −3000 m (left) and y = 3000 m (right).The TSFD interfaces are located at z = −10 m (top), z = 40 m (bottom), y = −105 m (left) and y = 105 m (right).Note, that for the horizontal TSFD setup the imaginar y par t of the B x -component does ne ver full y converge, but deviates no more than 5 per cent from the reference, if the domain is suf ficientl y large.

Table 1 .
Overview: solution variables, BC, types of fields prescribed at outer boundary and field impressed at TSFD interface for different modelling approaches and TSFD schemes ( cf .Fig.1).