Direct Sequential Simulation for spherical linear inverse problems

We present a method for obtaining efficient probabilistic solutions to geostatistical and linear inverse problems in spherical geometry. Our Spherical Direct Sequential Simulation (SDSSIM) framework combines information from possibly noisy observations, that provide either point information on the model or are related to the model by a linear averaging kernel, and statistics derived from a-priori training models. It generates realizations from marginal posterior probability distributions of model parameters that are not limited to be Gaussian. We avoid the restriction to Cartesian geometry built into many existing geostatistical simulation codes, and work instead with grids in spherical geometry relevant to problems in Earth and Space sciences. We demonstrate our scheme using a synthetic example, showing that it produces realistic posterior realizations consistent with the known solution while fitting observations within their uncertainty and reproducing the model parameter distribution and covariance statistics of a-priori training models. Secondly, we present an application to real satellite observations, estimating the posterior probability distribution for the geomagnetic field at the core-mantle boundary. Our results reproduce well-known features of the core-mantle boundary magnetic field, and also allow probabilistic investigations of the magnetic field morphology. Small-length scale features in the posterior realizations are not determined by the observations but match the covariance statistics extracted from geodynamo simulation training models. The framework presented here represents a step towards more general approaches to probabilistic inversion in spherical geometry.

Direct sequential simulation for spherical linear inverse problems years been developed to include probabilistic solutions based on Bayesian methods (e.g. Tarantola and Valette, 1982;Tarantola, 2005). In the Bayesian formulation one seeks to estimate the posterior probability density function (pdf) of the model parameters, proportional to the product of an a priori pdf and a likelihood function. Probabilistic solutions to inverse problems with non-Gaussian prior information are often obtained using sampling methods such as the Metropolis algorithm, but this becomes very expensive when working with high dimensional model spaces. Here we present an alternative approach to generating realizations of the posterior pdf for linear inverse problems in spherical geometry, based on observations related to the model by a linear averaging kernel, that can account for non-Gaussian prior distributions for the model parameters. We also include the capability of using point data of the model parameters and refer to these as direct observations.
In geostatistics, solutions based on point data (where the model parameters are known at some locations) are often obtained using kriging, i.e., interpolation through Gaussian process modelling conditioned on prior covariances, which provides the best linear unbiased predictions based on observed point data (Journel and Huijbregts, 1978;Deutsch and Journel, 1998). Such kriging schemes can be extended to systems where there is a combination of point data and observations related to the model parameters by a linear averaging kernel, and to cases where only the latter are available (Hansen et al., 2006). When the observation noise and a-priori pdfs are Gaussian, the posterior solution is provided by simple kriging through the mean and variance of a Gaussian estimate for each model parameter. Extensions of this framework, whereby the form of the posterior pdf can be obtained beyond means and covariances, is possible by sequentially simulating model parameters through sampling of local distributions conditional on prior information (Soares, 2001). Direct sequential simulation was applied to combinations of point and weighted linear average observations in Cartesian geometry in the VISIM algorithm of Hansen and Mosegaard (2008).
In spherical geometry, to the best of our knowledge, no implementation of direct sequential simulation algorithms yet exists for obtaining probabilistic solutions of linear inverse problems, although the underlying methods are well understood. Recent work by Alegría et al. (2020) introduced a method for simulating Gaussian random fields on the d-dimensional unit sphere which is computationally very efficient but does not possess the non-Gaussian capabilities of direct sequential simulation. Gneiting (2013) analyses valid positive definite correlation functions on spheres which may be used to generate the necessary spherical covariance models. Spherical harmonics are a well known method to represent continuous fields on the sphere (Wieczorek and Meschede, 2018); they also provide a means to specify isotropic prior covariance functions on a sphere which have an exact correspondence to the spherical harmonic power spectra (e.g. Moritz, 1980;Jackson, 1994;Hipkin, 2001). Building on the work of Hansen and Mosegaard (2008) in Cartesian geometry, and making use of covariance models linked to spherical harmonic spectra, we implement direct sequential simulation in spherical geometry with the aim of providing a new tool for Earth and Space science problems. While we focus here on working with non-Gaussian posterior pdfs, sequential Gaussian simulation using point data is also possible with the tool presented. We illustrate our method by obtaining a probabilistic solution for the geomagnetic field at the Earth's core-mantle boundary, taking prior information from geodynamo simulations and based on real satellite observations. Our Spherical Direct Sequential Simulation (SDSSIM) algorithm enables probabilistic solutions to this problem without assuming a priori that the model parameters are Gaussian distributed.
In section 2 we describe the linear forward problem, = , focusing on its discretization in spherical geometry. We next review the basic principles of Gaussian process based least-squares solutions to the inverse problem, sequential Gaussian simulation methods, and the theory of the direct sequential simulation method (Soares, 2001;Oz et al., 2003).
Section 3 gives a detailed description of the implementation of our SDSSIM algorithm. In section 4 we present the results of tests on both synthetic and real data, based on the geophysical problem of inferring the Earth's magnetic field at the core-mantle boundary from remote, noisy, satellite observations. Here we also demonstrate classic direct sequential simulation by using synthetic direct observations from a known simulation of the core-mantle boundary radial field. Finally we discuss the strengths and limitations of our method, along with our conclusions and some perspectives for future steps in 5.

The linear forward problem in spherical geometry
In spherical geometry, given observations, ( ), and model parameters on a spherical surface, ( ), related through a forward kernel operator ( , ) we consider a linear forward problem of the form shown in equation (1).
where = sin ′ ′ ′ , with = ( , , ) indicating locations of the observations, and = ( ′ , ′ , ′ ) indicating locations of the model parameters on a spherical surface of radius ′ . This system describes observations that are related to the model parameters by a linear averaging kernel. The integral equation in (1) may be approximated numerically via quadrature rules; here we use a Gauss-Legendre quadrature scheme appropriate for spherical geometry (e.g. Atkinson, 1982;Wieczorek and Meschede, 2018), in which the integration is carried out on a (2 − 1) × grid, where is the number of latitudinal nodes. cos ′ are then the Gauss-Legendre nodes on the interval [−1, 1], with corresponding integration weights, . ′ is chosen such that the points are equally spaced with separation ∕( − 1∕2) on the interval [0, 2 [. For model parameters on a sphere distributed according to Gauss-Legendre quadrature rules, this numerical integration is exact for polynomials of degrees less than 2 (Atkinson, 1982). The integral may then be discretized according to equation (2).
where = (2 − 1) × is the number of model parameters. For a series of observations, = 1 ( 1 ), … , ( ), … , ( ) , with a vector of model parameters, = 1 ( 1 ), … , ( ), … , ( ) ; we absorb the constant, −1∕2 , and integration weights into the elements of a matrix (size × ) such that Other grids on the sphere could alternatively be used, along with suitable quadrature weights. We adopted the Gauss-Legendre grid for simplicity and due to the ease of transforming to a spherical harmonic representation. Any linear forward problem in spherical geometry may then be written in the familiar form = Here we are concerned with the inverse problem of how best to estimate , a vector of parameter values on a spherical surface grid, given noisy observed data linearly related to the model, along with any prior information on the model parameters.

Equivalent least-squares solution to the linear inverse problem
A simple solution to the above inverse problem exists if we are able to assume the spherical surface model parameters can be represented by a Gaussian probability density function (pdf) with a priori mean 0 and covariance , while the observations, , represent realizations of Gaussian random variables with data error covariance (Tarantola, 2005). The least-squares solution is then also a Gaussian pdf with mean and covariancê This solution is identical to the solution of a simple kriging system (Hansen et al., 2006). Before presenting spherical direct sequential simulation (section 2.4) as an alternative solution method which avoids these often restrictive Gaussian assumptions, we first briefly describe the method of sequential Gaussian simulation on the sphere.

Sequential Gaussian Simulation on the sphere
The method of sequential simulation (e.g. Deutsch and Journel, 1998;Hansen and Mosegaard, 2008)

involves inferring
Gaussian posterior realizationŝ of the random variables , from the observations . For a joint distribution of random variables, , conditioned on a set of known observations, , the variate cumulative distribution function where denotes probability. Sequential simulation involves drawing an variate sample based on (8) making use of the product rule of probability, such that each probability term on the right-hand side is sampled in succession (Deutsch and Journel, 1998). Realizations are thus obtained in a series of sequential steps, gradually increasing the conditioning, beginning with the observations, .
For a Gaussian random field, drawing samples satisfying (8) equates to drawing from the Gaussian pdf,  ( , 2 ), where and 2 are the kriging mean and variance found by solving the kriging system (e.g. Journel and Huijbregts, 1978;Deutsch and Journel, 1998;Hansen and Mosegaard, 2008), which in our notation is where is a member of the available conditional variables in each step (observations, point data, and previously simulated model parameters),̂ is the model parameter currently being simulated, and are known as the kriging weights which determine the desired Gaussian pdf,  ( , 2 ). ( , ) are the a priori covariances between conditional variables including any measurement error covariance, and ( , ) are the covariances between conditional variables and the target model parameter. Solving equation (9) for the kriging weights , the kriging mean and variance are = ⋅ ( − 0̄ ) + 0 wherē = 1, … , 1 of length (10) where 0 and 2 0 are a-priori estimates of the mean and variance of the model parameters.
Covariances between observation pairs and observation/model parameter pairs can be obtained from the forward relation between the observation and model parameters defined in (2) and (4), the a-priori model parameter covariance measure, {}, and an example a-priori model, 0 ( ). With the covariance of an observation pair defined by (2), and again absorbing the constant and integration weights into ( , ), we have and All required covariances are thus available given an a priori model parameter covariance based on prior information regarding the random variable on the spherical surface and knowledge of the forward problem. The kriging system (9) can therefore be expanded as follows, in order to explicitly show the contributing parts of the covariance matrix is a matrix of observation to observation covariances computed using (12), holds observation data error covariances. and the vector contain covariances between previously simulated model parameters and between previously simulated model parameters and the target model parameter, both are obtained directly from the a-priori model covariance. and are covariances from observations to previously simulated model parameters, and to the target model parameter respectively, as given by (13).
Solving the kriging system sequentially using the above covariances, results in a sequential Gaussian simulation model realization  ( , 2 ). Example model realizations are drawn from the posterior distribution by visiting model parameters in a random order for each realization.
Such sequential Gaussian simulation schemes are well-known and widely used in geostatistics. However many physical processes in Earth and Space physics are fundamentally nonlinear which results in non-Gaussian statistics for the model parameters . In the next section we extend the above treatment to permit non-Gaussian model parameter distributions, based on the method of direct sequential simulation (Journel, 1994;Tran et al., 2001;Oz et al., 2003).
This approach ensures the linear relationship (4) between and is preserved, which is not the case for sequential Gaussian simulations after transforming and/or to Gaussian variables. In sequential Gaussian simulation such a transformation destroys the linear relationship so the necessary covariance matrices cannot be expressed simply as a function of and .

Direct Sequential Simulation for non-Gaussian fields on a sphere
The system described in section 2.3 allows one to sequentially simulate model parameters on a spherical surface given observations, leading to a realization of a Gaussian random field which fits the observations to within measurement error, and as far as this fit allows, reproduces the mean and covariance of the a priori information. We now go further and simulate non-Gaussian random fields using direct sequential simulation with histogram reproduction for a given training model, following the methods proposed by Journel (1994), Tran et al. (2001), and Oz et al. (2003).
A normal-score transform of a training model to variables, , that follow a standard Gaussian distribution (i.e. zero mean, variance one), and the associated back-transformation to original values may be performed through (15) and (16).
Where −1 is a standard Gaussian quantile function with cdf, , is the training model cdf with quantile function −1 , and 0 are the training model values. This transformation describes the connection between random variables with a Gaussian distribution, and non-Gaussian distributions defined by the training model. It allows one to generate a collection of non-Gaussian cdf's through which the sampling in the sequential simulation steps described in equation (8) can occur. Generating a non-Gaussian cdf is achieved by substituting the standard Gaussian representation of the training model, , for a Gaussian distribution, , with mean, , and variance 2 . This can be achieved as follows through inverse transform sampling using a vector, , of uniformly spaced quantiles between zero and one, which divides the Gaussian distribution into intervals of equal probability. The ranges of the mean and variance should cover approximately [−3.5, 3.5] and [0, 2] respectively, to fully utilize the training model in characterizing conditional distributions (Oz et al., 2003). A variance range of [0,2] is used in order to alleviate cases where the obtained kriging variances lie outside the domain of the generated local conditional distributions as shown by Deutsch et al. (2001).
Performing transformation with as in (16) results in a discrete vectorized quantile function, , conditional on the training model and with length describing a distribution with mean, , and variance, 2 , from which a sample, , can be drawn using a uniform distribution, , discretized in intervals Solving the kriging system yields an estimated mean and variance of the local Gaussian distributions and a distribution is then assigned from based on the mean, , and variance, 2 , closest to the kriging mean and variance. In this step a distance measure must be used and our implementation is described in section 3.3.2. We refer to the chosen distributions collectively as the local distributions. However, reproduction of the training model is only ensured if the applied local distribution has mean and variance equal to the kriging mean and variance (Journel, 1994). This further requires that the local distributions are scaled to have exactly the kriging mean and variance. For a value sampled from one of the local distributions, this is achieved bŷ wherê is the final simulated model parameter value that makes up the vector of model parameterŝ = ̂ 1 , … ,̂ , … ,̂ in a given realization. A probabilistic solution is achieved by collecting model parameter realizations in the matrix,̂ , and computing the sample covariance as follows Direct sequential simulation for spherical linear inverse problemŝ wherê is an length column vector of the model parameter means. This procedure ensures that simulations represent samples from the posterior probability density function of the model parameters based on the mean, variance, covariance structure, and histogram provided by the training model, while honoring the data Oz et al., 2003).

Implementation
We have implemented the methods described in section 2 as a Python repository called Spherical Direct Sequential Simulation, which is available on Github at github.com/mikkelotzen/spherical_direct_sequential_simul ation. The implementation includes five modules. (i) The geometry of the problem and the forward operator, relating the observations to the model parameters, (ii) the prior information, (iii) the measured observations, (iv) the simulation itself, and (v) the posterior pdf output. Figure 1 gives an overview of these modules; their content is described in more detail below. The propagation of information is shown with arrows.

Geometry and the forward problem
In SDSSIM the spherical polar coordinates for the surface points to be modelled and observations are stored in arrays of length and respectively, with the radius of the points to be modelled on the spherical surface of interest being a constant, ′ . Equation (1) defines the spherical linear inverse problem by connecting random variables on a spherical surface to observations. This spatial structure is illustrated with an example in figure 2, where a unit radius Gauss-Legendre quadrature (GLQ) grid is displayed among equal area distributed observations with varying radii.
The discretized forward operator, , which defines the connection between the observations and spherical surface in equations (3) -(4), is an array of size ( , ), and its construction is problem dependent. In section 4 an example of a forward operator is presented.

Observations
We store observations, ( ), in the vector array of length , with one value for each observation coordinate = ( , , ). These values are noisy and linearly related to the desired model parameters being simulated as described by the forward problem in (1)

Posterior pdf Estimation statistics
Spherical Direct Sequential Simulation

Prior information
Geometry and the forward problem Results Figure 1: The implementation of Spherical Direct Sequential Simulation (SDSSIM) shown as a flowchart. We consider our implementation as five distinct modules. Geometry and forward operator to set up the spherical inverse problem, prior information as required to solve the inverse problem, observations conditioning the solution, SDSSIM performing inversion by solving system equations, and results as the final output.
covariance indices. In section 4 we consider a case using the magnitude of the radial component of Earth's magnetic field at satellite altitude and a case using synthetic direct observations at the core-mantle boundary.

Prior information
Prior information is contained in the training model, 0 , an array of length . The values in the training model provides a distribution, with a priori mean, 0 , and variance, 2 0 . The training model is also used as conditioning for a range of possible local distributions and further prior information in the form of a covariance model. Included in Figure 2: Illustration of unit Gauss-Legendre quadrature grid and observations distributed according to an equal area grid with varying radii in spherical space. It is also possible to work with other grids, including approximate equal area grids, depending on the application. Grids on spherical surfaces are typical of Earth observation and geophysics problems utilizing satellites or global ground networks for data collection.
the implementation is the possibility of semi-variogram modelling. In our implementation this allows for estimating a covariance model from the training model based on the assumption that it is second-order stationary and isotropic. It is also possible to use an a priori power spectrum derived from the model to specify the covariance model.

Covariance model
For an isotropic field one can compute the spherical harmonic power spectrum and use this to define the covariance model, and hence the covariance matrix (e.g. Jackson, 1994;Hipkin, 2001, in the geomagnetic framework). Given any such covariance matrix, , the data to data and data to model parameter covariances are computed through (12) and (13) respectively. The observation to observation data error covariance, , will later be added to ; this is a diagonal array of the data error covariance level.

Local distributions conditional on prior
Local distributions conditional on the a priori training model used to sample model parameters are implemented as a lookup-table (LUT) based on equations (15)-(18). The procedure is shown as pseudo-code in algorithm 1. The first input is the training model and the second is the number of local distribution quantiles, . The number of quantiles is chosen by the user and should at most equal the size of the training model, as it controls the level of detail expressed in the conditional distributions, which cannot exceed the level of detail in the histogram on which they are based. Two further inputs are the discretization levels, and , the ranges of the mean and standard deviation used in generating Gaussian distributions as shown in equation (17). The discretization level of the mean and standard deviation range determines the number of local distributions in the LUT. From these inputs, a uniformly spaced array, , in the range zero to one is generated containing equally spaced values. and the range of mean and standard deviation values determined by and are then used iteratively in equation (17) and (18)  Having generated , we require a measure for finding the nearest local distribution given a kriging mean and variance, and 2 , such that a simulated model parameter (20) can be computed. This is achieved using an array of measures () + = 1 + = 1

Spherical Direct Sequential Simulation
Expanding on the outline of SDSSIM given in the algorithm flowchart of figure 1, the algorithm contains the following steps.
1. Determine a random path through the model parameters on the spherical surface.
2. At each location in the random path solve equation (9) with the appropriate covariance matrices based on all available observations and previously simulated values. The kriging mean, , and variance, 2 , are then determined through (10) and (11). 3. The nearest local distribution in is found through (22). This provides the local distribution closest to the kriging mean and variance.
4. Draw a sample from this nearest local distribution. (20)

A case study from geophysics: Core-mantle boundary magnetic field estimation
We now demonstrate SDSSIM on an example spherical linear inverse problem, estimating the radial component of Earth's magnetic field on the approximately spherical core-mantle boundary (CMB) from globally-distributed satellite magnetic observations (e.g. Langel, 1987;Bloxham et al., 1989;Gubbins, 2004;Finlay, 2020).
First we test our method on a synthetic case considering two distinct scenarios where different types of observations are generated from a known source, referred to below as the synthetic truth. In subsection 4.1.1 we consider synthetic satellite observations and in subsection 4.1.2 we consider direct observations of some of the model parameters (i.e. of the synthetic truth radial field at the CMB with added noise). In all cases, for prior information we use as traning models an ensemble of 487 instances of the core-mantle boundary field (up to spherical harmonic degree 30) from a numerical model of the magnetic field generating dynamo process in Earth's outer core (Aubert et al., 2017). These dynamo fields contain highly localized field structures that lead to a more Laplacian than Gaussian distribution of the radial field at the CMB. The synthetic truth is chosen to be another snapshot from the dynamo model, not included in the training set.
Having validated the method in this test case we go on to use real satellite magnetic field observations from Swarm Alpha, one satellite from ESA's Swarm constellation mission (e.g. Friis-Christensen et al., 2008). Swarm data are freely available and were downloaded through the Swarm Virtual Research Environment 1 .
The geomagnetic forward problem is of the form (1) with a forward operator based on the Green's function describing a potential field solution to Laplace's equation in spherical geometry for internal source Neumann boundary conditions (e.g. Gubbins and Roberts, 1983;Hammer and Finlay, 2019). For simplicity we consider only observations of the radial component of the field; the forward operator linking the radial geomagnetic field at an observation location, = ( , , ), to the radial field at a source location on the spherical core-mantle boundary , = ( ′ , ′ , ′ ), is then where ℎ = ′ , = √ 2 + ′ 2 − 2 ′ cos Υ with cos Υ = cos cos ′ + sin sin ′ cos( − ′ ) We represent the radial magnetic field in physical space at the CMB at a radius of ′ = 3480.0km, on a Gauss Legendre quadrature grid with = 31 latitudinal nodes and thus have = 1891 model parameters. This allows accurate transformation to a spherical harmonic representation up to degree = 30.
When only satellite observations are available we use the training ensemble of dynamo model realizations to specify an a-priori covariance function for the CMB radial magnetic field model. Assuming isotropy and stationarity over the spherical surface the covariance function for the radial magnetic field may be written (Jackson, 1994;Hipkin, 2001) where ( ′ ) = ( + 1) ′ 2 +4 ∑ =0 ( ) 2 + (ℎ ) 2 is the Lowes spherical harmonic power spectrum (Lowes, 1966) at ′ , with and ℎ Schmidt quasi-normalized spherical harmonic coefficients of degree and order , is the reference radius of the spherical harmonics, and are Legendre polynomials of degree . We used (24) to compute an a-priori model covariance matrix between all grid points on the CMB that is consistent with the isotropic second order statistics given by the power spectra of the realizations in the dynamo model training ensemble. The mean of these defines our a-priori model covariance matrix, . The training models in the dynamo ensemble are provided to spherical harmonic degree = 30, whereas equation (24) involves a summation to infinity. Generating a covariance matrix using (24) based on truncation to degree 30 can thus result in non-positive-definite covariance matrices due to the truncation. To avoid this problem we implement a function that gradually tapers the spectra to zero beyond degree 30, using = 0.5 −5 + 0.5 −2 . The a-priori local distribution for the model parameters is obtained by concatenating histograms from the training ensemble of numerical dynamo simulations of the CMB radial field. Figure   3(b) shows the training ensemble histograms and (c) shows the ensemble of Lowes power spectra used to generate the a priori covariance model.
For the test based on direct observations (i.e. on sampled values of the radial magnetic field at the CMB with noise added), the a-priori covariance is instead specified using an exponential semi-variogram model estimated from the available direct observations and a-priori local distributions are generated from their distribution shape.

Synthetic satellite observations
For this test we used = 2773 synthetic observations located at satellite altitude at positions sampled along real Swarm Alpha orbits taken from April-June 2018. To the synthetic radial field computed on the satellite orbits (based on the synthetic truth model) we add zero mean random Gaussian noise with std.dev. of 2nT. In the simulation we characterize this with a diagonal data error covariance matrix of (2nT) 2 .
where P is here the marginal posterior distribution and Q is a Gaussian distribution with equal mean and variance sam-

Synthetic direct observations of the CMB magnetic field
Validation tests based on synthetic direct observations of the model parameters are common in geostatistical studies.
We report here briefly the results of such a test in order to demonstrate that our method can also be used within a more  conventional geostatistical setup . To obtain direct observations, we sampled 27% of the synthetic truth radial magnetic field at the CMB and added zero mean random Gaussian noise with a std.dev. of 2nT. 1000 posterior realizations are generated using the covariance models and a-priori histograms described above. In figure 6 we display the resulting posterior mean and std. deviation maps of the radial magnetic field at the CMB. The large scale structures are well reproduced for areas containing direct observations, whereas areas without direct observations lack the structures present in the synthetic truth (bottom of figure 4). The posterior std. deviation clearly aligns with the locations of the sampled direct observations and show that these areas are more informed. The histogram and semi-variogram reproduction is successful (not shown). This test shows that SDSSIM is capable of performing classical direct sequential simulation using only direct observations of the model.

Probabilistic inversion of real satellite magnetic observations
We now move on to the more interesting case of inferring the radial magnetic field at the CMB using real satellite data. We use data from Swarm Alpha sampled at 5 minute intervals, during dark times, over one year from 01-11-2018 to 01-11-2019, applying standard data selection criteria to remove observations collected during periods of high solar-driven field disturbances (Kauristie et al., 2017;Finlay, 2020).

This resulted in
= 4884 observations at altitudes between 432km and 452km. In order to isolate the core field, we removed the LCS-1  model of the lithospheric magnetic field and the magnetospheric field and secular variation (changes over time of the core field) predicted by the CHAOS-7.2 model . The resulting radial field observations are presented in figure 7. For simplicity we assume the data error to be independent, Gaussian, with zero mean, and a standard deviation of 6nT at all latitudes. CHAOS-7 is truncated as is conventional at degree 13 while the other models are visualized after truncating at degree 30. Our results show slightly higher power and more detailed structures, particularly the maximum of the marginal posterior. There are similarities to structures seen in other studies which have attempted to infer the core field above degree 13 (e.g. Baerenzung et al., 2020;Aubert, 2020). In particular the strong flux patch in the equatorial Atlantic is split into two as also seen by Baerenzung et al. (2020) while e.g. above and slightly to the west of this patch, small scale patches are present which were also seen in Aubert (2020). The difference seen in the maximum of the marginal posterior map compared to the equivalent LSQ and posterior mean solutions indicate the presence of non-Gaussian features in the posterior realizations.
Although the differences between the maps of the maximum of the marginal posterior and the posterior mean are minor, and less than the differences between either of them and traditional spherical harmonic-based models such as CHAOS-7, there are nevertheless some interesting features. We note that the maximum of the marginal posterior shows higher amplitude features in the regions under the South Atlantic south-west of Africa; such features are important for understanding recent changes in the South Atlantic weak field anomaly at Earth's surface (Finlay, 2020). Higher amplitude features are seen in the central Pacific region, around latitude 15 degrees South, longitude 100 degrees West.
Finally there is a noticeable East-West elongation of a positive flux feature South East of Madagascar around latitude 30 degrees South and 60 degrees East. Overall the map of the maximum of the marginal posterior shows generally sharper features than the map of the posterior mean or that from CHAOS-7.
An example of a probabilistic investigation of CMB radial field structures is shown in Figure 11. This presents histograms of the integrated radial magnetic field inside the cylinder tangent to the Earth's solid inner core (Livermore et al., 2017), separated into normal and reversed polarities, in the north and south hemispheres. This analysis demonstrates the northern hemisphere has with high probability more reversed magnetic flux and weaker normal flux, a result of importance in geodynamo studies that was difficult to quantify with conventional field models.

Discussion and conclusions
In the case studies presented we found that the posterior mean is close to the equivalent least-squares solution. Why then go to all the trouble of generating posterior realizations? A key point here is that marginal posterior distributions at particular locations can still be non-Gaussian (see e.g. Fig. 5) and hence are not necessarily well described by the posterior mean and variance. The importance of this has previously been highlighted in the Cartesian case (Hansen and Mosegaard, 2008) and will doubtless also prove crucial for some applications in spherical geometry, particularly when the prior model distributions are strongly non-Gaussian.
In the presented applications we have routinely transformed from the simulated grids in physical space to spherical harmonic representations. This was found to be useful for comparisons with existing geomagnetic field models and for visualization, but care is needed with this procedure. The transform from the grid in physical space to spherical harmonics is only exact for real square-integrable functions and when the spherical harmonic degree of the underlying function is limited to the level of chosen Gauss-Legendre quadrature grid (e.g. Wieczorek and Meschede, 2018). We observe a smoothing/loss of power on the grid scale compared to the originally simulated values in some of the results presented here. This is acceptable if one wishes to compare models only up to some specific spherical harmonic degree, but for applications with covariance functions that allow discontinuities between neighbouring grid points, it is recommended to work instead with the simulated grids in physical space. There may be important advantages to working directly in the physical domain because the a priori covariance information can then be allowed to vary with position. For geophysical problems involving Earth's lithosphere and upper mantle it may for instance be important to allow different covariance models for positions in the continents versus oceans, or to use locally defined information based on auxiliary variables such as geological composition or features. Allowing spatial variations in the a-priori covariance models is an obvious next step for the framework presented here. The presented geomagnetic application was somewhat limited in the sense it involved only sources at one depth and ignored any time dependence. The extension to sources at multiple depths (ideally with independently specified prior information) can be achieved by superposing the sources and visiting the model parameters at all depths during the sequential simulation. Similarly the model grids could be extended to a sequence of times in order to account for time-dependent source processes, provided the necessary time-dependent covariance matrices are specified (see e.g. Gillet et al., 2013;Ropp et al., 2020;Baerenzung et al., 2020).
A major limitation of the present implementation of the SDSSIM algorithm is the use of two-point statistics (covariances) for describing the a-priori conditional relationship between the model parameters. Given the complexity of natural phenomena on the sphere, these are not capable of fully capturing all the essential details. In order to move beyond this limitation, similar algorithms in Cartesian geometry have utilized multiple-point statistics (Strebelle, 2002;Gravey and Mariethoz, 2020). Use of multiple-point statistics in spherical geometry is not yet well developed, but would certainly be of interest for improving on the results obtained here.
The SDSSIM scheme is in principle applicable to a wide variety of problems involving linear inversion or interpolation on a sphere. For example, possible applications could involve meteorological data such as the case presented by Jun and Stein (2008) where non-stationary covariance models are used to analyse global ozone levels or Jeong et al. (2017) where isotropic and non-stationary covariance models are used with global surface temperature data. Extensions to 3D using grids at different radii and radial covariance functions is also possible, e.g. for inversion problems in seismology (Meschede and Romanowicz, 2015) or gravity (Save et al., 2016). The success in such applications will rely on the availability of suitable prior information, for example in the form of training images or covariance fuctions. In such cases SDSSIM may allow for an improved exploitation of prior information and probabilistic descriptions of models in spherical geometry.

Funding sources
This study was funded by the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (Grant agreement No. 772561).

Computer Code Availability
• Name of code: spherical_direct_sequential_simulation