An Ensemble Algorithm Based Component for Geomagnetic Data Assimilation

Geomagnetic data assimilation is one of the most recent developments in geomagnetic studies. It combines geodynamo model outputs and surface geomagnetic observations to provide more accurate estimates of the core dynamic state and provide accurate geomagnetic secular variation forecasting. To facilitate geomagnetic data assimilation studies, we develop a standalone data assimilation component for the geomagnetic community. This component is used to calculate the forecast error covariance matrices and the gain matrix from a given geodynamo solution, which can then be used for sequential geomagnetic data assimilation. This component is very flexible and can be executed independently. It can also be easily integrated with arbitrary dynamo models.


INTRODUCTION
The Earth possesses a magnetic field of internal origin (geomagnetic field).It is now widely accepted that this magnetic field is generated and maintained by convective flow in the Earth fluid outer core (geodynamo).Surface geomagnetic observations over the past several centuries, as well as paleomagnetic and acheomagnetic records, have demonstrated that the geomagnetic field varies in broad spatial and temporal scales (Jackson et al. 2000;Sabaka et al. 2004;Korte and Constable 2008;Olsen et al. 2010).Such variations are the manifestation of dynamic processes in the Earth's liquid outer core, and are therefore critical for understanding the interior, the evolution and the interactions among the different Earth components.
In addition to surface geomagnetic observations, various field and dynamic models have been utilized to understand and interpret geomagnetic variability.Among them are various dynamo models in which the convection in the outer core is described by the Navier-Stokes equation, while the equations describing magnetic induction and en-ergy conservation (when combined these are called the dynamo equations).In the dynamo models the outer core fluid properties are also defined using several non-dimensional parameters.These equations are then solved using specified boundary conditions.For more detail we refer the reader to Kono and Roberts (2002).
Numerical dynamo models have been used successfully for qualitative understanding of geomagnetic secular variation.However, due to limited knowledge of the Earth's outer core, and computational constraints, these models are still far from the "truth".This is the main reason for substantial discrepancies between the observed geomagnetic field and those from dynamo simulations (Kuang et al. 2009).Therefore, a fundamental question arises: how could one obtain the most accurate estimate of the core dynamic state with currently available information?
One answer to this question is geomagnetic data assimilation: using geomagnetic observations to constrain and improve numerical dynamo models.Generally speaking, data assimilation is a methodology in which observational data are assimilated into numerical model outputs for better estimation of the true physical states.Data assimilation has been widely used in atmospheric (e.g., Talagrand 1997) and oceanic modeling (e.g., Ghil andMalanotte-Rizzoli 1991, Sun et al. 2013).There are two types of assimilation methods: sequential and variational data assimilation.Detailed descriptions of the algorithms and history of data assimilation can be found in, e.g., Kalnay (2003).Both methods have been used in geomagnetic data assimilation studies (Fournier et al. 2007;Liu et al. 2007;Sun et al. 2007;Kuang et al. 2008Kuang et al. , 2010;;Fournier et al. 2010), with very promising results.However, extensive studies are needed to understand both the responses of numerical dynamo models to observational constraints and the consistencies of the parameters and approximations used in dynamo simulation with geomagnetic observations.Geomagnetic data assimilation is also computationally very expensive.The required computing resources could be orders of magnitude more than those for numerical dynamo simulations.There are a variety of dynamo models currently used by different research groups, thus calling for algorithm portability among these dynamo models.
One approach to solve these conflicts between research demands and computational constraints is to develop an independent assimilation component.It could then be used in the geomagnetic community.With this component community efforts can be cohesively utilized to address the fundamental questions in geomagnetic data assimilation.
In this paper we present our geomagnetic component based on an ensemble sequential assimilation algorithm.This paper is organized as follows: in section 2, the assimilation algorithm is summarized.In section 3, the basic structure and workflow of the component are described.Numerical examples are given in section 4, followed by concluding remarks in section 5.

ALGORITHM
There are currently many dynamo models that differ in many aspects, including boundary conditions, algorithms and geophysical approximations.Regardless, these models can be symbolically written as the following initial-value problem, ( , ) ( ) where x is the state vector (an array describing the dynamo state in the outer core), F is the generalized force governing the time evolution of x, and x 0 is the initial sate at t 0 .
The numerical modeling objective is that the model output x f (called "forecast" in data assimilation) can be a good approximation of truth x t .However, it often falls short in achieving this objective for many reasons, such as large model errors, inaccurate initial guesses, turbulent dynamic processes, etc.To reduce the difference between the two, one approach is the sequential data assimilation in which the model output x f and (often partial) observations y o are combined to obtain a better approximation x a (called analysis) of truth x t : (1 ) x W x Wy The key is to find an appropriate weight W in Eq. ( 2) such that x a is the optimal estimation of x t .The analysis x a is then used as the new initial state of Eq. ( 1) to obtain a more accurate forecast in the future.If the domain of y o is different from that of x f , a mapping is then needed in Eq. ( 2).In general, the analysis Eq. ( 2) can be derived whenever both observations and forecasts are available.
The mathematical reasoning can be summarized as follows.There are two different sources of information about the truth: one is the model output x f (often called a priori information) and the other is the observation y o : Both are different from the truth, , x x e y Hx e where H is the observational operator (a mapping between the observation domain and the truth).The errors e f and e o are often assumed to have Gaussian random distributions with zero mean and variances P f and R, respectively.It is often appropriate to assume that both are independent (i.e., uncorrelated).Thus, where the error e z is also a Gaussian distribution with zero mean and the variance, To obtain the optimal estimation of truth x t from z, a linear unbiased estimator x a = Az is proposed with the transform A chosen to minimize the trace of the error covariance, where E $ ^h denotes the expectation of an ensemble.This leads to where I is the identity matrix.Thus, Using Eqs. ( 3), ( 6), ( 8), and (9), In data assimilation K is called the gain matrix, and y Hx is called the innovation.Equation ( 10b) can be also interpreted as follows: consider the sum of the distances (weighted by the error variance) from x a to the forecast x f , and from x a to the observation If K is chosen to be Eq.( 10b), then the distance J is the minimum.It can also be shown that Eq. ( 10b) is equivalent to (Kalnay 2003), Equations ( 10a) and ( 12) are used in our system.There are two extreme scenarios in Eq. ( 12): , observation replaces the forecast in the analysis; or 0 K .if R P f % , i.e., the observation is negligible in the analysis.
The observational error covariance R is determined with the measurements and the relevant analysis (e.g., geomagnetic field modeling).It is therefore an "input" to the assimilation component.For the details we refer the reader to, e.g., Olsen (2002) and Korte and Constable (2008).Evaluation of the forecast error covariance P f is the centerpiece of the algorithm.There are various means used to calculate P f , e.g., a statistical study of the innovations (Hollingsworth and Lönnberg 1986) or the lagged-forecast difference method (Parrish and Derber 1992).In our system an ensemblecovariance estimation of P f is employed that could also serve as a component for the more complex ensemble Kalman filter algorithm (Evensen 1994).This ensemble method can be carried out for any numerical dynamo model without knowing its details and can also be parallelized.
In this estimation an ensemble of N ens initial states , ,..., i N 1 2 , is first generated by adding arbitrary perturbations to the given model state vector x 0 , where the perturbations i f are generated randomly via the Monte Carlo technique and their magnitudes are only a small fraction to that of x 0 .These perturbed initial states are then integrated in time via Eq.( 1) for a specified time interval.The final states x i " , are used to evaluate P f as follows: In geodynamo modeling the state vector is comprised of spherical harmonic coefficients of the velocity, magnetic field and temperature perturbation.These harmonic coefficients also vary in the radius r and time t.Therefore P f is a function of the degree l and order m of the expansions.It also varies in radius r.From our previous experiences, the covariance between different degrees and orders are much weaker than those of the same degree and order and are therefore ignored in the evaluation.
Since the ensemble size N ens is finite and substantially smaller than that of the model output, spurious behaviors exist in the covariance of two distant points (e.g., in radius r).They could normally be eliminated by applying filters, such as the following proposed by Gaspari and Cohn (1999) where r is the given model grid location in radius, r 0 is the observational location, and c is the correlation length parameter (e.g., 10% of the outer core thickness).In Eq. ( 15), the radius is scaled by that of the the core-mantle boundary (CMB) (i.e., r = 1).Different dynamo models have different definitions of the state vector x.To make this component applicable to arbitrary dynamo models it is necessary to introduce a procedure in which the state vectors x K of an arbitrary dynamo model can be transformed back and forth with that x defined in this component: where T is the transform for the dynamo model state vectors and T y is the transform for the observational subspace.Using the definitions Eqs. ( 14) and ( 16), we have Thus, the analysis will be used for assimilation with the given dynamo model.The above formulations are for the state vectors x in the real space.If x is complex, e.g., complex spherical harmonic coefficients, one can simply replace the transpose matrices in Eqs.(10b), ( 12), and ( 14) using their complex conjugate transpose.

COMPONENT STRUCTURE AND WORKFLOW
In this component the state vector includes the velocity field, the magnetic field and the density perturbation: , , , , where P v and T v are the poloidal and toroidal velocity scalars; P B and T B are the poloidal and toroidal magnetic field scalars; and H is the density perturbation.Each scalar in Eq. ( 20) is described by a spherical harmonic expansion.For example, where C.C. stands for the complex conjugate (and L is the truncation order of the expansion).Surface geomagnetic observations (with maximal degree L obs ) are downward continued to the radial position r 0 .Thus ( ) ( , ) . .

P r r P r r H H P P H P H P H H
For example, using Eq. ( 14), we have  24) and ( 25) stands for the complex conjugate transpose.
The component structure is described in Fig. 1: it is comprised of three independent sections, each with its own functions: the gain matrix calculation (the right column in Fig. 1), the transformation between the state vectors (the center column in Fig. 1), and the dynamo simulation (the left column in Fig. 1).
The dynamo section is normally imported from a chosen dynamo model.This part does not need any reengineering effort.The transformation section needs input from users to define the transformation matrices (T and T -1 ) for the chosen dynamo model.The gain matrix section includes all "intrinsic" procedures and functions.Users do not need to know the details for their applications.By introducing such a structure this component is simple to implement and easy to learn.
The workflow of this component can be summarized as follows: first, select an initial state vector (a solution of a given dynamo model).Then, transform the state vector into the standard representation for the component.Next, add randomly generated perturbations to the transformed state vector.The modified state vector will then be inversely transformed and sent back to the dynamo model for time integration.The final output from the dynamo simulation will be sent back to the component for calculating the forecast error covariance P f .This process will be repeated until N ens final states are successfully obtained.The covariance P f and the gain matrix K are calculated and exported for geomagnetic data assimilation.

APPLICATION EXAMPLE
This component can be easily integrated with arbitrary dynamo models.To demonstrate this we consider the MoSST core dynamics model series developed by Kuang and Bloxham (1999), Kuang and Chao (2003), Jiang and Kuang (2008).The state vector x K defined in MoSST is slightly dif- ferent from that used in this component: the radial variation of the velocity scalars P v and T v are expressed by Chebyshev polynomial expansions in MoSST, while they are defined on the radial grid points in the component.Therefore, the transform matrix T is simply the Chebyshev transform.
In this application the entire MoSST dynamo model (including its source codes and executable) is in one directory (hereafter called the "dynamo directory"), while the component is in a separate directory (call the "component directory").The component directory includes the three subdirectories for the dynamo, interface and gain matrix sections, respectively.The first step is to copy the necessary source codes for defining the dynamo state vector and parameters from the dynamo directory into the component directory.Users then modify the interface section accordingly, including the dynamo directory information.They also need to define the parameters for the component, e.g., the ensemble size N ens , the correlation length c defined in Eq. ( 15), and the perturbation amplitude.The component execution script includes the MoSST execution for dynamo simulation.The simulation process ends when an ensemble of N ens successful final states is obtained.It should be pointed out that dynamo simulation results from a given initial state are stored in a designated subdirectory (thus N ens subdirectories will be created in the process).These outputs are used for calculating the gain matrix K (part of the component output) and will be used for future analysis if needed.
The magnitude of the random perturbations is small: it is set as 1% relative to the spectral coefficients of the original dynamo solutions.In Fig. 2 is the coefficient ( ) b r 2 1 of the initial state defined in Eq. ( 21). Figure 3 shows an example of random perturbation.The original random perturbation is applied to the 2 nd order radial derivative of ( ) b r 2 1 (Fig. 3b): perturbations are added at every four grid point; those in between are set via linear interpolation.The final perturbation f [see Eq. ( 13)] is shown in Fig. 3a.
In our application, the forecast error covariance P f is simply of the form defined in Eqs. ( 24) and ( 25).The filter Eq. (15) used in our example here is defined with c = 0.1D (D is the thickness of the outer core).As an example, we show in Fig. 4a

DISCUSSION
This paper described our geomagnetic data assimilation component with an ensemble-based algorithm Eq. ( 14).This component is stand-alone.It can be implemented independent of dynamo models and the final result can be transformed easily for the designated dynamo model via Eqs.( 17), (18), and (19).Given an initial dynamo state and a dynamo model (either serial or parallel), this data assimilation component can generate the gain matrix K automatically without any interactions from users.Its implementation is very simple: only the transform section is customer-designed.The gain matrix calculation component is very userfriendly: its details are not needed for applications.The dynamo section is imported directly from the selected dynamo model.However, since the definitions of the dynamo state vectors x and x K are often very simple, transform matrix T construction is not very complicated.The component is easy to maintain.Only the transform section needs to be updated if any change is made to the dynamo model.
This component is efficient for parallelization: its scalability can be achieved either through the original dynamo model (since the dynamo executable will be called directly by the component), or, if the original dynamo model is not scalable, by evenly distributing individual dynamo runs.
Simple modification to the component can be achieved by selecting different values for the ensemble size N ens and the filtering length scale c.However, more complicated modifications to the algorithm will need substantial reprogramming effort.For example, if the observation data y o are different, e.g., the core flow inverted from the observed geomagnetic secular variation is included as part of the observations, the observation operator H in Eq. ( 12) will be different.Thus, the gain matrix of Eq. ( 12) will be more complicated.Another example is that this component is not very effective when the time variation of the gain matrix K becomes important.In this case it will be optimal to integrate it directly to the dynamo model.This is beyond the scope of this paper.This component can be improved fur-ther as more users implement it for their geomagnetic data assimilation studies.the stand-alone Fortran modules are stored in the subdirectory named "modules_EnK".
The component application can be summarized as follows: (1) copy all files to a separate subdirectory (we strongly recommend users for this procedure for better maintenance); (2) compile the code; (3) modify the parameter files according to the dynamo package of the users; (4) execute the component and check the output.The last step may be repeated if the desired ensemble members cannot be reached in a single execution.
. The superscript H in Eqs. ( corresponding gain matrix elements with R = 0 are shown in Figs.6 and 7, respectively.

Fig. 1 .Fig. 3 . 1 0
Fig. 1.Sketch of the component structure: in the left column are the dynamo modules and the initial state imported by users; in the right column are the main modules of the component for computing K; in the center column are the interface modules to compute the transform T.

Fig. 5 .
Fig. 5. Similar to Fig. 4, but for the error covariance ( , ) P r r l jb 1 0 .Fig. 6.Similar to Fig. 4, but for the gain matrix K corresponding to ( , ) P r r l bb 1 0 .