Lateral variations of electrical conductivity in the lower mantle constrained by Swarm and CryoSat-2 missions

The electrical conductivity is an important geophysical parameter connected to the thermal, chemical, and mineralogical state of the Earth’s mantle. In this paper, we apply the previously developed methodology of forward and inverse EM induction modeling to the latest version of satellite-derived spherical harmonic coefficients of external and internal magnetic field, and obtain the first 3-D mantle conductivity models with contributions from Swarm and CryoSat-2 satellite data. We recover degree 3 conductivity structures which partially overlap with the shape of the large low-shear velocity provinces in the lower mantle.


Introduction
The electrical conductivity is an important geophysical parameter connected to the thermal, chemical, and mineralogical state of the Earth's mantle. A traditional technique to study the distribution of electrical conductivity in deep regions of the Earth is the electromagnetic induction (EMI) method. On the global planetary scale, the natural sources of electromagnetic energy with the largest potential to probe the deep Earth structure are present in the magnetosphere, the ionosphere, and in the Earth' oceans. An excellent overview of recent advances in the area of global EMI forward and inverse modeling is given by Kuvshinov (2015).
In the last two decades, the low-orbit satellite missions dedicated to measurements of the Earth's magnetic fields, such as CHAMP, Ørsted, SAC-C, and Swarm, have provided invaluable insight into the structure and dynamics of various physical processes from the dynamo in the Earth's core to the magnetosphere (Stolle et al. 2018). Recently, these datasets have been extended by careful processing and calibration of platform magnetometer records from other missions, such as CryoSat-2 .
In the area of global EM induction, spherically symmetric (1-D) conductivity models based on satellite data were obtained (Kuvshinov and Olsen 2006;Velímský 2010;Püthe et al. 2015;Civet et al. 2015;Grayver et al. 2017, among others). The global three-dimensional (3-D) inverse problem represents a significant challenge in terms of input data processing, efficient implementation of the forward solver, and balanced regularization. So far, only ground observatory data were inverted in terms of 3-D mantle conductivity models (Kelbert et al. 2009;Tarits and Mandéa 2010;Semenov and Kuvshinov 2012;Sun et al. 2015). The main challenge in the 3-D inversion of satellite data is the separation of the primary, inducing components, and the secondary, induced fields observed by a small number of fast moving platforms with sufficient spatial and temporal resolution. One possible path that we follow here, is to isolate the time series of external and internal spherical harmonic (SH) coefficients by the process of comprehensive inversion (Sabaka et al. 2013). These time series are provided as a Swarm mission Level 2 product by the European Space Agency. Since the rather critical assessment of the early versions of this product by Martinec et al. (2018), the situation has improved significantly (Sabaka et al. 2018(Sabaka et al. , 2020. In this paper, we apply the methodology developed by Velímský (2013) and Maksimov and Velímský (2017) to the latest version of satellite-derived SH coefficients, and obtain the first 3-D mantle conductivity models with contributions from Swarm and CryoSat-2 satellite data.

Methods
Detailed descriptions and tests of the forward solver, the inverse solver, and the sensitivity calculations in the context of global 3-D EM induction problem in the time domain are presented in Velímský and Martinec (2005), Velímský (2013), and Maksimov and Velímský (2017), respectively. Our approach to the 3-D mantle conductivity inversion builds on the process of comprehensive inversion (CI). Here we provide only a short summary of the method, and refer the reader to the detailed description in the works of the original authors (Sabaka et al. 2013(Sabaka et al. , 2018(Sabaka et al. , 2020, and in particular in the paper by Kuvshinov et al. (2020) in this issue.

The forward problem
In the spherical Earth G of radius a with 3-D distribution of electrical conductivity σ (r) > 0 , and under magnetoquasistatic approximation, the magnetic field B(r; t) is governed by the EM induction equation, Here r = (r, �) = (r, ϑ, ϕ) stands for the position vector described by the radial coordinate r, colatitude ϑ , longitude ϕ , t is time, and µ 0 is the magnetic permeability of vacuum. In the surrounding insulating atmosphere, the magnetic field is fully described by a scalar magnetic potential U (r; t) , satisfying the homogeneous Laplace equation, On the surface r = a , continuity of the magnetic field is required, The solution of the Laplace Eq. (2) can be written in terms of SH series, (1) (2) ∇ 2 U (r; t) = 0, for r ≥ a.
(3) B(r = a, �; t) = − ∇U (r, �; t)| r=a . (4) where G (e) jm (t) and G (i) jm (t) are the SH coefficients describing, respectively, the external (primary, inducing) and internal (secondary, induced) magnetic field. Here we use fully normalized, real spherical harmonics Y jm (�) . The spherical harmonic-finite element approach to global time-domain EM induction introduced by Velímský and Martinec (2005) uses the external boundary condition. The time series of external field coefficients G (e) jm (t i ) is prescribed in the discretized linear system of equations corresponding to formulation (1, 2, 3). The SH expansion is truncated at a finite degree j (e) max , and constant sampling in time is assumed, For a given conductivity distribution σ (r) , the method integrates Eq. (1) using the Crank-Nicolson scheme. Besides the complete solution B(r; t) , it provides also the series of internal field coefficients G (i) jm (t i ) , using the same time sampling, and truncated at j (i) max ≥ j (e) max .

Observatory, Swarm and CryoSat-2 data
The CI is based on fitting of a complex parameterization of individual magnetic field constituents into observatory and satellite data, including the along-track and crosstrack gradients. In the first phase of the CI process, the models of the core field and its secular variations, the lithospheric field, the ionospheric, magnetospheric, and oceanic tidal fields, are constructed from quiet-time signals. The EMI process is accounted for in the parameterization of the magnetospheric fields by means of RC index (Olsen et al. 2014), using a 1-D conductivity model, and in the parameterization of ionospheric fields by a transfer matrix Q(ω) , based on a 1-D mantle conductivity with a 2-D surface layer representing the conductivity contrast between the oceans and continents (Sabaka et al. 2013). The rotation angles between the common reference frame of each spacecraft, and the reference frame of individual magnetometers are also co-estimated. The data selection is restricted to night-side, quiet-time signals, defined by the limiting values of Kp index, and rate of change of the Dst index (Sabaka et al. 2013(Sabaka et al. , 2018(Sabaka et al. , 2020. For our purpose, the second phase of the CI processing is crucial. Here, the observatory and satellite data from all times, including the disturbed periods, are assembled. The observatory dataset comprises 5 or 6 years of hourly mean vector measurements of magnetic field from 161 ground observatories (Macmillan and Olsen 2013). The Swarm dataset consists of the vector magnetic field measurements from Swarm Alpha and Bravo satellites at 1-min sampling, covering the same period. Finally, the (5) Page 3 of 12 Velímský and Knopp Earth, Planets and Space (2021) 73:4 CryoSat-2 dataset includes 5 years of calibrated magnetic field vectors measured by the platform magnetometers at 1-min sampling ). The core, lithospheric, and in most cases also the ionospheric signals obtained by the CI are subtracted, and the resulting residua are further down-weighted, or completely removed in the polar areas. The SH coefficients of external and internal field up to a specified degree, j (e) max , j (i) max , and order are then obtained by SH analysis of data in bins of selected duration τ , typically several hours. Distinct values of τ are used for the dipole terms (1, 0), and the remaining spherical harmonics. Table 1 gives an overview of various datasets from the CI, and respective settings. The particular choices of j (e) max , j (i) max , and τ are subject to a trade-off, governed by the spatiotemporal resolution of the datasets. Increasing the width of the time bin would allow to increase formally also the SH truncation degree, but the field variations on shorter time scales would then yield spurious spatial variations. Moreover, the smaller spatial scales of magnetospheric and induced fields are both reduced at the altitude of low-orbit satellites by downward, and upward continuation, respectively ). Although the 3-D mantle conductivity inversion was carried out for all of the datasets, we limit the graphic presentation of the results to only two cases in this paper. The dataset 0513 is based on 5 years of data, including the data from the CryoSat-2 platform magnetometers. The dataset 0604 extends the length of the time series to 6 years, however, without the inclusion of CryoSat-2 data.
Three additional remarks deserve the attention of the reader. First, the external and internal field SH coefficients provided by the CI use the traditional Schmidt's semi-normalization, and must be converted to fully normalized values using formulas (5, 6) from Velímský (2013). Second, the higher sampling rate of the dipolar coefficients is not exploited in our modeling, as the forward and inverse solvers require the same constant time step for all spherical harmonics. Therefore, the dipolar coefficients are downsampled by averaging to the same sampling rate as the remaining coefficients. The renormalization and time-averaging can be combined into a single linear operator Ŵ, where g represents a vector arrangement of the CI-based, Schmidt semi-normalized SH coefficients in the original time sampling, and G is a similar arrangement of the fully normalized, time-averaged harmonics. Figure 1 shows the time series of the external and internal field coefficients for the 0513 dataset, after such processing.
Thirdly, the CI also provides the inverse covariance matrices D MMA , relating the information on the (6) G = Ŵ g, uncertainties of the SH coefficients. Without going into the detailed structure of Ŵ , we can use it also to obtain the inverse covariance matrix of the fully normalized time-averaged coefficients as In the synthetic tests performed in Velímský (2013) only the diagonal terms of the inverse covariance matrix were used. Here we take the next step by including the terms relating the SH coefficients of different degree and order (j, m), (j ′ , m ′ ) that belong to the same time bin of length τ . The elements of D then satisfy where \ denotes the integer division operator. Since i and i ′ are used as the outer indices, the inverse covariance matrix has block-diagonal structure. Figure 2 shows an example of one such block for the 0513 dataset, corresponding to time level 2018.0. Note that the diagonal terms dominate. Although the matrix is positive definite, small negative weights do appear off the diagonal.

The inverse problem
The purpose of the inverse problem is to reconstruct the conductivity distribution σ (r) from the time series of external and internal field coefficients G (e,obs) obtained by processing of observatory and satellite data; year the obs supersript now differentiates the observationbased time series from the predictions of the forward solver. Because of inherent non-uniqueness of the inverse problem, the class of admissible conductivity models is restricted by regularization. In our case, the inversion is defined over a finite-dimensional manifold M . Since the electrical conductivity is positive-valued, a natural choice is to parameterize its (decimal) logarithm where σ 0 = 1 S/m is used to preserve correct dimensionality, and the negative sign stems from the internal representation of the model by means of resistivity, rather than conductivity. The radial discretization is defined on K layers, where ξ k (r) = 1 in the kth layer r k ≤ r ≤ r k+1 , and zero otherwise. The lateral resolution is governed by the truncation degree j (9) is used for all r < a − δ . In the uppermost layer of thickness δ = 13 km , the electrical conductivity is not included in the model vector m , but fixed to an a priori distribution, taking into account the conductivities of igneous rocks, seawater, continental, coastal, and ocean sediments, along the lines of Everett et al. (2003).
Following Velímský (2013), we solve the regularized inverse problem by finding the minimum of a penalty function over M . Here, χ 2 (m) is the data misfit, measuring the distance between observed data and prediction of the forward problem. Let G (i) jm (m; t i ) be the solution of the forward problem (1, 2, 3), excited by the series of external field coefficients G (e) jm (t i ) = G (e,obs) jm (t i ) , and using the conductivity model σ (m; r) obtained from a particular realization of m ∈ M using Eq. (9). The data misfit is then defined as  Table 1 Overview

of various datasets for which the 3-D conductivity inversion is attempted
The datasets vary by the truncation degree and order of the external field, time sampling (the first value corresponds to the external and internal dipole terms 1/0, the second value to all other terms), the length of series, the inclusion of observatory, Swarm, and CryoSat-2 data (O, S, and C, respectively), and data selection. DW55GM: data polewards of dipole latitude 55 • down-weighted by sin(ϑ GM )/100 ; DW50, DW55: data polewards of dipole latitude 50 • or 55 • down-weighted by sin(ϑ GC )/100 ; R50, R55: data polewards of dipole latitude 50 • or 55 • removed, NS: only night-side data, when the Sun is at least 10 • below horizon are used; NI: the comprehensive model of the ionosphere is not subtracted from data. ϑ GM and ϑ GC stand for geomagnetic (dipole) and geocentric colatitudes, respectively where Note that in Velímský (2013), only the simplified formula with diagonal covariance matrix was used. The regularization function R 2 (m) in the present work constrains the L 2 norm of the first spatial derivative of the log-conductivity over the volume G δ , which excludes the uppermost Kth layer, where h k = r k − r k−1 . Thus, the fixed near-surface conductivity map, and its interface with the underlying mantle do not contribute to the regularization value. Other options were also explored (e.g., constraining the second derivatives of log-conductivity, or using different weighting of radial and lateral parts of the smoothing operators), but the results do not not convey any substantially different information, and are not presented here.
The penalty function F (m; ) defined in Eq. (10) is minimized for a series of 15 regularization coefficients = (5, 2, 1, 0.5, 0.2, 0.1, . . . , 0.0001) , providing a series of models The first run for = 5 is started from m = 0 , corresponding to a homogeneous sphere with unit conductivity, and each subsequent run is initiated at the minimum obtained for the nearest larger . The minimization algorithm is based on the limited-memory quasi-Newton method with Broyden-Fletcher-Goldfarb-Shanno formula (Press et al. 1992, Section 10.7) used for gradual updates of the approximate Hessian. The gradient of the data misfit in the model space ∇ m χ 2 (m) is evaluated using the adjoint approach (Velímský 2013;Maksimov and Velímský 2017). Finally, the data misfits χ 2 (m( )) are plotted versus the regularizations R 2 (m( )) , in the form of L-curve (Hansen 1992). The optimal regularization parameter ˜ and corresponding model m(˜ ) are then chosen visually near the maximum inflection point, thus balancing the data misfit and the regularization function. Figure 3 presents the L-curves obtained for the two inversion runs described in detail below.

Sensitivity analysis
We use the methodology described by Maksimov and Velímský (2017) to calculate the data Hessian matrix in the final, optimally regularized model m(˜ ), for each dataset. The efficient algorithm based on the combination of the forward, adjoint, forward scatter, and adjoint scatter problems allows the assembly of the entire matrix at the cost equivalent to 2M + 2 forward runs. The Hessian matrix provides a useful insight into the sensitivity of data to individual model parameters. Figure 7 shows the data Hessians for the models presented in this paper, and detailed discussion is presented in the next section. Note, that in the minimum of penalty function F (m) , only the full Hessian matrix H F , which includes also the second derivatives of regularization term, is positive definite. The data Hessian H χ can have small negative values on the diagonal, but we find it more instructive for discussion of data sensitivity.

Inversion settings
The input settings used for the CUP20-OSC5 and CUP20-OS6 runs are summarized in Table 2, together with the respective diagnostics results. The radial parameterization of the conductivity model is selected closely to the previous 3-D inversions by Semenov and Kuvshinov (2012, further referenced as ETH12) and Sun et al. (2015, further referenced as OSU15). The layer interfaces are placed at the depths of 40, 250, 410, 520, 670, 900, 1200, 1600, 2400, and 2891.2 km. The lateral resolution of the forward solver is set to j max = 8 , and the layers from the inverse model are further sub-discretized into 192 radial nodes in total, comprising also the near-surface layer, and the core with constant conductivity of 2 × 10 5 S/m.

Results
We start this section by comparison of the 1-D, spherically symmetric part of the recovered 3-D models. It is governed only by the ρ k 00 components of the model vector m . Figure 4 shows the 1-D profiles for all 3-D inversion runs of various CI datasets. The CUP20-OSC5 and CUP20-OS6 models are marked, respectively, by thick red and blue Page 6 of 12 Velímský and Knopp Earth, Planets and Space (2021) 73:4 lines. In addition, we show by thick black line the 1-D conductivity profile MIN1DM_0501, which was obtained by joint regularized inversion of satellite-derived C-responses and M 2 tidal signals (Grayver et al. 2017). Although this approach uses different datasets, model parameterization and regularization, and the forward modeling is based on the integral-equation method in the frequency domain, there is remarkable agreement between the 1-D profiles below 250 km. Our models fail to resolve the low-conductivity layer below the lithosphere-asthenosphere boundary, but this is expected as the periods shorter than 6 hrs are filtered out from our input datasets. Use of alternative settings of the data processing in the CI has the largest impact on the 1-D conductivity profiles in the upper mantle. Below 900 km, the spherically symmetric parts of all models are in good agreement. Figures 5, 6 compare the lateral cross-sections of CUP20-OSC5, CUP20-OS6, OSU15-j3, and ETH12-j3 models in the upper and lower mantle, respectively. In order to compare only the large-scale features of all models, we have low-pass filtered the OSU15 and ETH12 models by least-squares fitting of the SH series up to degree j ρ max = 3 from Eq. (9) into the grid values of log-conductivity in each layer, hence the -j3 appendage. The lateral variations in the ETH12-j3 model start only below the depth of 410 km.
We can observe some systematic features in the CUP20-OSC5, CUP20-OS6, OSU15-j3, and ETH12-j3 models throughout the transition zone of the upper   3 Trade-off plots between data misfit and regularization (the L-curves) for CUP20-OS6 (left) and CUP20-OCS5 (right) inversions, respectively. The red numbers denote the respective values of the regularization parameter

Table 2 The input settings, and the diagnostics outputs (regularization parameter, data misfit, and regularization value) of the 3-D runs presented in this paper
Parameter   Table 1 are shown by thin lines of various colors without further distinction Page 7 of 12 Velímský and Knopp Earth, Planets and Space (2021) 73:4 mantle. All four models prefer increased conductivity below Eurasia, the southern Pacific and the Indian Ocean, and lower conductivity below the western Pacific. The largest disagreement between the models is below South America. Here both CUP20 models, and to a lesser extent also the ETH12-j3 model, infer a conductivity increase from the Pacific coast, across the continent, and further into the South Atlantic, while the OSU15-j3 model shows rather weaker lateral dependence, and in the opposite direction. The ETH12-j3 model is significantly more conductive in the lowest part of the transition zone, and shows less lateral variations. In Fig. 6, we summarize the results of the four conductivity models for the lower mantle. Note that the color scale changes with respect to the one used in Fig. 5. However, it remains consistent across the models. In the uppermost part of the lower mantle we observe good agreement between both CUP20 models and the ETH12-j3 model in the position of the lateral conductivity heterogeneities, although in the latter case more suppressed by regularization. Increased conductivity is observed below the South America, and Eurasia, and there is also a less distinct conductive patch below the southern Pacific. This feature also occurs in the OSU15-j3 model. On the other hand, a large negative conductivity anomaly is present below both Americas only in the OSU15-j3 model.
Proceeding below the depths of 1000 km, the OSU15-j3, and the ETH12-j3 models show almost no large-scale lateral conductivity variations. However, in the case of CUP20 models, it is interesting to compare the lower mantle conductivity features with the distribution of the large low shear-wave velocity provinces (LLSVP) observed by seismic tomography (Lekic et al. 2012;Doubrovine et al. 2016;Hosseini et al. 2018). A conductive feature similar in its shape and size to the African LLSVP, Page 8 of 12 Velímský and Knopp Earth, Planets and Space (2021) 73:4 although partially shifted to the east, is present in both of our models, spanning from the Eurasia to the southern Indian Ocean. A second conductive feature overlaps well with the Pacific LLSVP, although it seems to extend further below the South America, which is not observed in the tomography models. A resistive area spanning the entire depth of the lower mantle in the CUP20 models is located below the Australia. This feature is consistent with the inversion of monthly mean values from geomagnetic observatories by Tarits and Mandéa (2010, not shown here). A quantitative interpretation of the conductivity models in terms of lower mantle temperature and composition is beyond the scope of this paper. In general, where the conductivity increase coincides with lower seismic velocities, the temperature is probably the main control parameter (Verhoeven et al. 2009). However, other effects, such as the iron spin transitions and MORB presence can influence the lower mantle conductivity in both ways (Lin et al. 2007;Ohta et al. 2010).
Finally, the sensitivity of the datasets used in the inversions to individual model parameters of CUP20-OS6 and CUP20-OSC5 is assessed by studying the respective Hessians in Fig. 7. In both cases, the largest sensitivities are observed below the transition zone, and slowly decreasing down to 2400 km. A second, less distinct area of large second derivatives, is present just above the transition zone, especially for the CUP20-OS6 model. However, given the 6 hrs sampling rate of the external and internal SH coefficients, we do not believe that the data has significant resolution above the transition zone. The data sensitivity to individual model parameters is also influenced by the thickness of individual layers in a non-linear way. We have tested it in several inversion runs with alternative parameterization, using layers of constant thickness of 100 km (not shown here). In those cases, the relative Cross-sections of electrical conductivity models CUP20-OS6, CUP20-OSC5, OSU15-j3, and ETH12-j3 in the lower mantle Page 9 of 12 Velímský and Knopp Earth, Planets and Space (2021) 73:4 amplitudes of Hessian elements above the transition zone are reduced. Throughout the models, there are significant correlations and anticorrelations between the SH coefficients of the same degree and order across neighboring layers, which quickly fall off for more distant layers. Within each layer, the diagonal sensitivities dominate, and they decrease with increasing SH degree.

Conclusions and outlook
Together with the paper by Kuvshinov et al. (2020), our work represents the first attempt to obtain the 3-D image of mantle conductivity using satellite data. The small number of satellites is the main limiting factor in these studies. As we deal with the transient signals originating in the magnetosphere, and their induced counterparts, the separation of temporal and spatial variations from observations on a moving platform becomes a challenging problem with some necessary trade-offs. So far, the process of CI has been able to recover the time series of external and induced field coefficients up to degree and order 3, with time resolution starting at 6 hrs. Hence, the 3-D conductivity inversion must be focused on the lower mantle, and only large-scale, averaged features can be reconstructed from the data. The 3-D conductivity models discussed in the preceding section have recovered robust features which have only small sensitivity to the changes in model parameterizations and datafiltering settings in the CI. The inclusion of CryoSat-2 platform magnetometer data had only small influence on the inversion results, it did allow for a slightly larger reduction of the data misfit in the CUP20-OSC5 inversion. However, it is not sufficient to significantly increase the spatial resolution of the inversion. It remains an open question, whether inclusion of additional platform magnetometer data with Swarm-derived calibration will enable us to progress to SH degree 5 or beyond in the near future.
The previous global 3-D inversions used the concept of local observatory responses. Hence, more details can be recovered in densely covered areas (e.g., below Europe) at the risk of obtaining spurious variations in poorly covered areas, unless corrected by regularization. Such approach is not suitable for satellite data, where global descriptions of the inducing and induced signals, or corresponding global transfer functions are needed. A combination of both approaches is one possible way for future studies.
Moreover, the small-scale induced fields in the lower mantle, carrying the information on the small-scale conductivity variations, are subject to geometric attenuation and EM shielding throughout the Earth, just like the considerably stronger small-scale fields of the geodynamo.
While the 3-D EMI sounding of the upper mantle can benefit from additional energy sources in the ionosphere (Kuvshinov and Koch 2015;Guzavina et al. 2019) and oceans (Grayver et al. 2017;Zhang et al. 2019), the choice of the EM excitation of the lower mantle is more limited. The inversion of the geomagnetic jerks (Nagao et al. 2003) is burdened by the principal ambiguity in the recovery of source field and conductivity, and similar argument can be used for the exploitation of the lengthof-day variations through the electromagnetic coremantle coupling (Holme 2000).
The EMI sounding by magnetospheric sources thus remains the only viable approach to provide a geophysical constraint on the lower mantle conductivity distribution, and thus contribute to the knowledge of the thermochemical and mineralogical state of the mantle. The main challenge for the future lies in the sparsity of observations, and data processing methods to overcome it. The 3-D forward EMI problem requires a boundary condition, external or other, to be specified everywhere on the surface for sufficiently long time series, or large periods. The source model thus needs to be interpolated in space and time/frequency from the available datapoints to fill the gaps, while suppressing incoherent contributions of non-inductive origin. Several promising data-driven approaches are being developed. One direction is the use of an a priori conductivity model in the construction of the source model. The source then can be constrained by using only the horizontal component of magnetic field which is less sensitive to lateral conductivity variations (Kuvshinov 2008;Grayver et al. 2020). The principal component analysis, as applied for ionospheric sources by Sun et al. (2015), could be also exploited here. Another path calls for a more detailed pre-processing and cleaning of satellite data in the polar areas by accounting for the polar electrojets and field-aligned currents by means of optimized model parameterization (Martinec et al. 2018). These signals, while not able to induce significant currents inside the Earth due to their geometry, are a strong source of noise at high latitudes. Finally, the physics-based models of the magnetosphere Fig. 7 Hessians (second derivatives of data misfit) obtained for the final models CUP20-OS6 (top) and CUP20-OCS5 (bottom), respectively. The parameters on the horizontal and vertical axes are ordered from left to right, and top to bottom, respectively. The outer index corresponds to layers and depths, with individual layers separated by the gridlines. The inner joint degree-order index is similar to Fig. 2, with the addition of (0, 0) term at the beginning of each layer block. Only the degree is annotated (See figure on next page.) Page 10 of 12 Velímský and Knopp Earth, Planets and Space (2021) 1  1  1  1  1  1  1  1  1  2  2  2  2  2  2  2  2  2  2  3  3  3  3  3  3  3  3  3  Page 11 of 12 Velímský and Knopp Earth, Planets and Space (2021) 73:4 can also be in principle used to provide the excitation signals, as was demonstrated by Honkonen et al. (2018) on forward 3-D modeling of EM fields induced by large geomagnetic storms. It remains to be seen, whether the additional complexity of this approach can help to overcome the problems of purely data-derived source models.