Augmented Kalman filter algorithm

Introduction

The subpackage augkf implements algorithms based on an augmented state Kalman filter (AugKF) initiated by [BGA17].

The following figure shows the Kalman filter in its simplest form:

The sequential data assimilation algorithm is composed of two kinds of operations:

Forecasts are performed every $\Delta t_f$. An ensemble of $N_e$ core states is time stepped between $t$ and $t+\Delta t_f$.
Analyses are performed every $\Delta t_a$ with $\Delta t_a=n\Delta t_f$ (analyses are performed every $n$ forecasts). A Best Linear Unbiased Estimate (BLUE) is performed using observation data.

From a technical point of view, algorithm steps take CoreState objects as inputs and return the CoreState resulting from the operations. Forecasts and analyses are handled by the AugkfAlgo object that implements the full AugKF algorithm.

More information can be found in the relevant tutorials:

The induction equation

The basic equation coupling the components of the core state at a given time is the radial induction equation at the core surface in the spherical harmonic spectral domain:

(1)\[\dot{\mathbf{b}} = \mathrm{A}(\mathbf{b})\mathbf{u} + \mathbf{e_r}\]

$\mathbf{b}$ (resp. $\mathbf{u}$ and $\dot{\mathbf{b}}$) stores the Schmidt semi-normalized spherical harmonic coefficients of the magnetic field (resp. the core flow and the secular variation) up to a truncation degree $L_b$ (resp. $L_u$ and $L_{sv}$ ). The number of stored coefficients in those vectors are respectively $N_b = L_b (L_b + 2)$, $N_u = 2L_u(L_u+2)$ and $N _{sv} = L_{sv}(L_{sv}+2)$.
$\mathrm{A}(\mathbf{b})$ is the matrix of Gaunt-Elsasser integrals [Moon79] of dimensions $N_{sv} \times N_u$ , depending on $\mathbf{b}$.
$\mathbf{e_r}$ stands for errors of representativeness (of dimension $N_{sv}$ ). This term accounts for both subgrid induction (arising due to the truncation of the fields) and magnetic diffusion.

Priors

The data assimilation steps require spatial cross-covariances that are derived from priors. The priors consist in a sampled series of spectral coefficients for the magnetic field MF, the core flow U and errors of representativeness ER. The error of representativeness (ER) can be computed from the core flow (U) and the main field (MF) using the script compute_ER.py.

The supplied priors in pygeodyn are the geodynamo free runs 0path, 50path, 71path, 100path codes but other types of priors can be used (more information in the section about `prior types`_).

AR(1) version

Forecast and AR(1) processes

We wish to fit the data with a multivariate AR-1 model of the form

\[\text{d}{\bf x} +{\sf \tilde A} {\bf x}\text{d}t = \text{d}\boldsymbol{\zeta}\]

with $\boldsymbol{\zeta}$ a Wiener process and ${\sf \tilde A}$ is the drift matrix, of dimension $P\times P$ with $P$ the dimension of the vector ${\bf x}_t={\bf x}(t)$

$\mathbf{x}$ is centered (substract average) and stands for the core flow U or the subgrid error ER or the augmented state Z ([U, ER])

Once discretized with an Euler-Maruyama scheme, it gives:

(2)\[{\bf x}_{i+1} = -{\bf x_i} {\sf A} + {\sf D}\]

with: ${\sf A} = {\Delta t \sf \tilde A} - {\bf I}$ and ${\sf D} = \sqrt{\Delta t} \mathbf{L} \mathbf{w_i}$, with $\mathbf{w_i}$ a normal draw ($\mathbf{w_i} = \mathcal{N}(0,1)$)

Thus, we seek $\mathbf{L}$ and ${\sf \tilde A}$.

Diagonal AR(1)

In the case where the geodynamo series used as priors do not allow to derive meaningful temporal statistics, the two drift matrices are simply diagonal, and controlled by a single free parameter tunable in the configuration file ($\tau_x$ for $\mathbf{x}$, U and ER only):

\[\mathrm{\tilde A} = \frac{1}{\tau_{x}}\mathrm{I}_{x}\]

with $\mathrm{I}_{x}$ the identity matrix of rank $N_{x}$.

The covariance matrix of ${\bf x}$ writes: $\mathrm{P}_{xx}= \mathbb{E}\left(\mathbf{x}\mathbf{x}^T\right)$

${\sf L}$ is the lower Cholesky matrix of $\frac{\sqrt 2}{Tau_x} \mathrm{P}_{xx}$

Dense AR(1)

We only treat the case of a single geodynamo prior. To know more about the mixing of priors (71path with 100path) please refer to the scientific paper …

The filtered-version $\hat{\bf x}(t)$ of the geodynamo series ${\bf x}(t)$ relate through the smoothing window $\omega$:

\[\hat{\bf x}(t) =\int w(\tau){\bf x}(t+\tau)d\tau\]

We emphasize on the difference between:

$\Delta t$ the time step of the geodynamo series
$\Delta s$ the time duration of the smoothing window $\omega$ (configuration parameter dt_smoothing)

We write $N_\tau = \frac{\Delta s}{\Delta t}$ the size of $\omega$ and $M_\tau = \sum_{j=-N_\tau}^{N_\tau}w_j^2$, the $L_2$ norm of $\omega$

The model is discretized as:

\[\hat{\bf x}_{t_i} =\sum_{j=-N_\tau}^{N_\tau} w_j {\bf x}_{t_{i+j}}\]

We now note:

\[\hat{\bf X}_i = \hat{\bf x}_i\]

\[\Delta \hat{\bf X}_i = \hat{\bf{x}}_{i+1} - \hat{\bf{x}}_i\]

The drift matrix $\sf \tilde A$ is computed as:

\[{\sf \tilde A} = - \frac{1}{\Delta t}\left[\frac{\Delta\hat{\bf X}\hat{\bf X}^{T}}{M_\tau N_\tau} \right] \left[\frac{\hat{\bf X}\hat{\bf X}^{T}}{{M_\tau N_\tau}} \right]^{-1}\]

Then we deduce $\hat{\bf W}$:

\[\hat{\bf W} = \left(\Delta \hat{\bf X} + \Delta t {\sf \tilde A} \hat{\bf X}\right)/\sqrt{\Delta t}\]

And finally compute ${\sf S}$:

\[{\sf S} = \displaystyle\frac{1}{\displaystyle N_s} \left[ \frac{\hat{\bf W} \hat{\bf W}^{T}}{{M_\tau N_\tau}} \right]\]

with, $N_s = \frac{N_t}{N_\tau}$ the number of independant data in the smoothed geodynamo series

$\mathbf{L}$ is the lower Cholesky matrix of ${\sf S}$

Complete forecast step

The first step of the forecast is to compute $\mathbf{u}(t+\Delta t_f)$ and $\mathbf{e_r}(t+\Delta t_f)$ using Eqs. (3).

If U ER are independant:

$\mathbf{u}'$ and $\mathbf{e}'$ are centered (average zero)

\[\mathbf{u}'(t+\Delta t_f) = -{\mathbf{u}'(t)} {\sf A_u} + {\sf D_u}\]

\[\mathbf{e}'(t+\Delta t_f) = -{\mathbf{e}'(t)} {\sf A_e} + {\sf D_e}\]

Otherwise we consider the augmented state Z ([U, ER]):

$\mathbf{z}'$ is centered (average zero)

\[\mathbf{z}'(t+\Delta t_f) = -{\mathbf{z}'(t)} {\sf A_z} + {\sf D_z}\]

Then, the vector $\dot{\mathbf{b}}(t+\Delta t_f)$ is evaluated thanks to Eq. (1) by using $\mathbf{u}(t+\Delta t_f)$, $\mathbf{e_r}(t+\Delta t_f)$ and $\mathbf{b}(t)$.

\[\dot{\mathbf{b}}(t+\Delta t_f) = \mathrm{A}(\mathbf{b}(t))\mathbf{u}(t+\Delta t_f) + \mathbf{e_r}(t+\Delta t_f)\]

where, $\mathrm{A}$ is the matrix of Gaunt-Elsasser integrals [Moon79]

Finally, the vector $\mathbf{b}(t+\Delta t_f)$ is evaluated by finite difference.

\[\mathbf{b}(t+\Delta t_f) = \mathbf{b}(t) + \Delta t_f \left[\dot{\mathbf{b}}(t+\Delta t_f)\right]\]

yielding the complete CoreState at $t+\Delta t_f$.

Analysis

The analysis step takes as input the ensemble of forecast core states $\mathbf{X}^f(t_a)$, the geodynamo statistics, plus MF and SV observations at $t=t_a$ together with their uncertainties. It is performed in two steps:

First, a best linear unbiased estimate (BLUE) of an ensemble of realisations of $\mathbf{b}$ is performed from magnetic field observations $\mathbf{b}^o(t)$ and the ensemble of forecasts $\mathbf{b}^f(t)$. The two are linked through an observation operator $\mathrm{H}_{b}$ that relates $\mathbf{b}^f(t)$ to $\mathbf{b}^o(t)$ ($\mathbf{b}^o(t) = \mathrm{H}_{b}\mathbf{b}^f(t)$). The Kalman gain matrix is then given by

\[\mathrm{K}_{bb} = \mathrm{P}^f_{bb}\mathrm{H}_b^T\left(\mathrm{H}_b \mathrm{P}^f_{bb} \mathrm{H}_b^T + \mathrm{R}_{bb}\right)^{-1}\]

where $\mathrm{P}^f_{bb}$ is the correlation matrix of the forecasted magnetic field and $\mathrm{R}_{bb}$ the MF observation error matrix. We use the glasso algorithm to approximate $\mathrm{P}^f_{bb}$ [IGF23] (see remove_spurious in config).

The analysed magnetic field vector $\mathbf{b}^a(t_a)$ comes from the BLUE:

\[\mathbf{b}^a(t_a) = \mathbf{b}^f(t_a) + \mathrm{K}_{bb}\left(\mathbf{b}^o(t_a) - \mathrm{H}_b\mathbf{b}^f(t_a)\right)\]

Second, a BLUE of an ensemble of realisations of the augmented state $\mathbf{z}$, that is a concatenation of $\mathbf{u}$ and $\mathbf{e}_r$ ($\mathbf{z} = [\mathbf{u}^T\mathbf{e}_r^T]^T$), using the observations of the SV.

The observation operator $\mathrm{H}_z$ relates $\mathbf{z}^f(t)$ to $\mathbf{\dot{b}}^o(t)$. It is composed of the observation operator for SV $\mathrm{H}_{sv}$ that relates $\mathbf{\dot{b}}^f(t)$ to $\mathbf{\dot{b}}^o(t)$ and of another matrix deriving from the induction equation (1) at $t_a$:

\[\mathbf{\dot{b}^f}(t_a) = \mathrm{A}(\mathbf{b}^a(t_a))\mathbf{u}^f(t_a) + \mathbf{e_r}^f(t_a).\]

We have then $\mathrm{H}_z(t_a) = \mathrm{H}_{sv}[\mathrm{A}(\mathbf{b}^a(t_a))\mathrm{I}_e]$ where $\mathrm{I}_e$ is the identity of rank $N_e=N_{sv}$.

The Kalman gain matrix is given by

\[\mathrm{K}_{zz} = \mathrm{P}^f_{zz}\mathrm{H}_z^T\left(\mathrm{H}_z \mathrm{P}^f_{zz} \mathrm{H}_z^T + \mathrm{R}_{\dot{b}\dot{b}}\right)^{-1}\]

where $\mathrm{P}^f_{zz}$ is the correlation matrix of the forecasted augmented state and $\mathrm{R}_{\dot{b}\dot{b}}$ the SV observation error matrix. We use the glasso algorithm to approximate $\mathrm{P}^f_{zz}$ [IGF23] (see remove_spurious in config).

The analysed core flow and subgrid errors are then given by the BLUE

\[\mathbf{z}^a(t_a) = \mathbf{z}^f(t_a) + \mathrm{K}_{zz}\left(\mathbf{\dot{b}}^o(t_a) - \mathrm{H}_z\mathbf{z}^f(t_a)\right)\]

Having now $\mathbf{b}^a(t_a)$, $\mathbf{u}^a(t_a)$, $\mathbf{e_r}^a(t_a)$, all is left is to derive $\mathbf{\dot{b}}^a(t_a)$ from the induction equation (1):

\[\mathbf{\dot{b}}^a(t_a) = \mathrm{A}(\mathbf{b}^a(t_a))\mathbf{u}^a(t_a) + \mathbf{e_r}^a(t_a).\]

Complete analysis step

The complete analysis step is therefore:

\[ \begin{align}\begin{aligned}\mathbf{b}^a(t_a) &= \mathbf{b}^f(t_a) + \mathrm{K}_{bb}\left(\mathbf{b}^o(t_a) - \mathrm{H}_b\mathbf{b}^f(t_a)\right)\\\mathbf{z}^a(t_a) &= \mathbf{z}^f(t_a) + \mathrm{K}_{zz}\left(\mathbf{\dot{b}}^o(t_a) - \mathrm{H}_z\mathbf{z}^f(t_a)\right)\\\mathbf{\dot{b}}^a(t_a) &= \mathrm{A}(\mathbf{b}^a(t_a))\mathbf{u}^a(t_a) + \mathbf{e_r}^a(t_a)\end{aligned}\end{align} \]

yielding the complete CoreState at $t_a$.

AR(3) version

Forecast and AR(3) processes

We wish to fit the data with a multivariate AR-3 model of the form

\[\text{d}{\bf x}''(t) + {\sf \tilde A} \text{d}{\bf x}'(t) + {\sf \tilde B} \text{d}{\bf x}(t) + {\sf \tilde C}{\bf x}(t) = {\sf \tilde D} \text{d}\boldsymbol{\zeta}\]

with $\boldsymbol{\zeta}$ a Wiener process and ${\sf \tilde A}$, ${\sf \tilde B}$, ${\sf \tilde C}$ are drift matrices, of dimension $P\times P$ with $P$ the dimension of the vector ${\bf x}_t={\bf x}(t)$

$\mathbf{x}$ is centered (substract average) and stands for the core flow U or the subgrid error ER or the augmented state Z ([U, ER])

We can write the AR(3) process as an AR(1) process:

\[\begin{split}\text{d}{\bf z}(t) = - \begin{bmatrix} 0&-I&0\\ 0&0&-I\\ {\sf \tilde C}&{\sf \tilde B}&{\sf \tilde A} \end{bmatrix} {\bf z}(t) \text{d}t + \begin{bmatrix} 0&0&0\\ 0&0&0\\ 0&0&{\sf \tilde D} \end{bmatrix} \begin{bmatrix} {\bf 0}\\ {\bf 0}\\ {\sf \tilde D} \end{bmatrix}\end{split}\]

with, ${\bf z}(t) = \begin{bmatrix} {\bf x}(t)\\ {\bf x}'(t)\\ {\bf x}''(t) \end{bmatrix}$

A first order differenciation gives

\[\begin{split}{\bf z}_i = \begin{bmatrix} {\bf x}_i\\ \frac{{\bf x}_{i+1} - {\bf x}_i}{\Delta t}\\ \frac{{\bf x}_{i+2} - 2 {\bf x}_{i+1} + {\bf x}_i}{\Delta t^2} \end{bmatrix}, {\Delta}{\bf z}_i = \begin{bmatrix} {\bf x}_{i+1} - {\bf x}_i\\ \frac{{\bf x}_{i+2} - 2 {\bf x}_{i+1} + {\bf x}_i}{\Delta t}\\ \frac{{\bf x}_{i+3} - 3 {\bf x}_{i+2} + 3{\bf x}_{i+1} - {\bf x}_i}{\Delta t^2} \end{bmatrix}\end{split}\]

Once discretized with an Euler-Maruyama scheme, it gives:

(3)\[{\bf x}_{i+3} = - {\bf x_{i+2}} {\sf A} - {\bf x_{i+1}} {\sf B} - {\bf x_i} {\sf C} + {\sf D}\]

with:

${\sf A} = {\Delta t \sf \tilde A} - 3 {\bf I}$ ${\sf B} = {\Delta t^2 \sf \tilde B} - 2 {\Delta t \sf \tilde A} + 3 {\bf I}$ ${\sf C} = {\Delta t^3 \sf \tilde C} - {\Delta t^2 \sf \tilde B} + {\Delta t \sf \tilde A} - {\bf I}$ ${\sf D} = \Delta t^{5/2} \mathbf{L} \mathbf{w_i}$, with $\mathbf{w_i}$ a normal draw ($\mathbf{w_i} = \mathcal{N}(0,1)$)

Thus, we seek $\mathbf{L}$, ${\sf \tilde A}$, ${\sf \tilde B}$ and ${\sf \tilde C}$.

The computation of these matrices is similar to the Dense AR(1) process.

Complete forecast step

The first step of the forecast is to compute $\mathbf{u}(t+\Delta t_f)$ and $\mathbf{e_r}(t+\Delta t_f)$ using Eqs. (3).

If U ER are independant:

$\mathbf{u}'$ and $\mathbf{e}'$ are centered (average zero)

\[\mathbf{u}'(t+\Delta t_f) = -{\mathbf{u}'(t)} {\sf A_u} -{\mathbf{u}'(t-\Delta t_f)} {\sf B_u} -{\mathbf{u}'(t- 2 \Delta t_f)} {\sf C_u} + {\sf D_u}\]

\[\mathbf{e}'(t+\Delta t_f) = -{\mathbf{e}'(t)} {\sf A_e} -{\mathbf{e}'(t-\Delta t_f)} {\sf B_e} -{\mathbf{e}'(t- 2 \Delta t_f)} {\sf C_e} + {\sf D_e}\]

Otherwise we consider the augmented state Z ([U, ER]):

$\mathbf{z}'$ is centered (average zero)

\[\mathbf{z}'(t+\Delta t_f) = -{\mathbf{z}'(t)} {\sf A_z} -{\mathbf{z}'(t-\Delta t_f)} {\sf B_z} -{\mathbf{z}'(t- 2 \Delta t_f)} {\sf C_z} + {\sf D_z}\]

Then, the vector $\dot{\mathbf{b}}(t+\Delta t_f)$ is evaluated thanks to Eq. (1) by using $\mathbf{u}(t+\Delta t_f)$, $\mathbf{e_r}(t+\Delta t_f)$ and $\mathbf{b}(t)$.

\[\dot{\mathbf{b}}(t+\Delta t_f) = \mathrm{A}(\mathbf{b}(t))\mathbf{u}(t+\Delta t_f) + \mathbf{e_r}(t+\Delta t_f)\]

where, $\mathrm{A}$ is the matrix of Gaunt-Elsasser integrals [Moon79]

Finally, the vector $\mathbf{b}(t+\Delta t_f)$ is evaluated by finite difference.

\[\mathbf{b}(t+\Delta t_f) = \mathbf{b}(t) + \Delta t_f \left[\dot{\mathbf{b}}(t+\Delta t_f)\right]\]

yielding the complete CoreState at $t+\Delta t_f$.

Analysis

In the following, we note $ta^- = ta - dta$ and $ta^+ = ta + dta$

The analysis step takes as input the ensemble of forecast core states $\mathbf{X}^f$ between $t = ta^-$ and $t = ta^+$, the geodynamo statistics, plus MF and SV observations at $t = ta^-$, $t=t_a$ and $t = ta^+$ together with their uncertainties.

It is performed in two steps:

MF analysis

First, we perform a linear regression to fit every realisation of the forecast core states from $t=ta^-$ to $t=ta^+$ using a 2nd order Taylor expansion.

We obtain an ensemble $\mathbf{B}_{t^{a}}^{f}=\left[\mathbf{b}_{t^{a}}^{f }, \mathbf{b}_{t^{a}}^{f^{\prime}}, \mathbf{b}_{t^{a}}^{f^{\prime \prime}}\right]$

We write: $\mathbf{Y}_{t^{a}}^{M F}=\mathrm{H}_{B} \mathbf{B}_{t^{a}}+\mathbf{E}_{t^{a}}^{M F}$

where $\mathbf{E}_{t^{a}}^{M F}=\left[\mathbf{e}_{t^{a-}}^{M F T}, \mathbf{e}_{t^{a}}^{M F T}, \mathbf{e}_{t^{a+}}^{M F T}\right]^{T}, \mathbf{Y}_{t^{a}}^{M F}=\left[\mathbf{y}_{t^{a-}}^{M F T}, \mathbf{y}_{t^{a}}^{M F T}, \mathbf{y}_{t^{a+}}^{M F T}\right]^{T}$, and

$\mathrm{H}_{B}=\left[\begin{array}{ccc} \mathrm{H}_{b} & -\delta t^{a} \mathrm{H}_{b} & \frac{1}{2} \delta t^{a 2} \mathrm{H}_{b} \\ \mathrm{H}_{b} & 0 & 0 \\ \mathrm{H}_{b} & +\delta t^{a} \mathrm{H}_{b} & \frac{1}{2} \delta t^{a 2} \mathrm{H}_{b} \end{array}\right]$

From an ensemble of $\mathbf{B}_{t^{a}}^{f}$ we compute empirically the model state forecast covariance matrix $\mathrm{P}_{B B}^{f}\left(t^{a}\right)=\mathbb{E}\left(\delta \mathbf{B}_{t^{a}}^{f} \delta \mathbf{B}_{t^{a}}^{f T}\right)$. We use the glasso algorithm to approximate $\mathrm{P}_{B B}^{f}\left(t^{a}\right)$ [IGF23] (see remove_spurious in config).

The Kalman gain matrix writes:

\[\mathrm{K}_{B}=\mathrm{P}_{B B}^{f} \mathrm{H}_{B}^{T}\left[\mathrm{H}_{B} \mathrm{P}_{B B}^{f} \mathrm{H}_{B}^{T}+\tilde{\mathrm{R}}^{M F}\right]^{-1}\]

to be used for the analysis step for each ensemble member:

\[\mathbf{B}_{t^{a}}^{a}=\mathbf{B}_{t^{a}}^{f}+\mathrm{K}_{B}\left(\mathbf{Y}_{t^{a}}^{M F}-\mathrm{H}_{B} \mathbf{B}_{t^{a}}^{f}\right)\]

The analyzed MF at the three epochs can then be deduced from $\mathbf{B}_{t^{a}}$ using a Taylor expansion:

\[\begin{split}\left[\begin{array}{l} \mathbf{b}_{t^{a-}}^{a} \\ \mathbf{b}_{t^{a}}^{a} \\ \mathbf{b}_{t^{a+}}^{a} \end{array}\right]=\left[\begin{array}{ccc} \mathbf{I}_{b} & -\delta t^{a} \mathbf{I}_{b} & \frac{1}{2} \delta t^{a 2} \mathbf{I}_{b} \\ \mathbf{I}_{b} & 0 & 0 \\ \mathbf{I}_{b} & +\delta t^{a} \mathbf{I}_{b} & \frac{1}{2} \delta t^{a 2} \mathbf{I}_{b} \end{array}\right] \mathbf{B}_{t^{a}}^{a}\end{split}\]

$\mathrm{I}_{b}$ is the identity matrix of the size of the vector $\mathbf{b}$.

SV analysis

First, we perform a linear regression to fit every realisation of the forecast core states from $t=ta^-$ to $t=ta^+$ using a 2nd order Taylor expansion.

We obtain an ensemble $\mathbf{Z}_{t^{a}}^{f}=\left[\mathbf{z}_{t^{a}}^{f }, \mathbf{z}_{t^{a}}^{f^{\prime}}, \mathbf{z}_{t^{a}}^{f^{\prime \prime}}\right]$

We write: $\Delta \mathbf{Y}_{t^{a}}^{S V}=\mathrm{H}_{Z}\left(\mathbf{B}_{t^{a}}\right) \mathbf{Z}_{t^{a}}+\mathbf{E}_{t^{a}}^{S V}$ where $\Delta \mathbf{y}_{t}^{S V}$ is a SV observation vector at epoch $t$ once corrected from the SV prediction from the average dynamo state.

where $\mathbf{E}_{t^{a}}^{S V}=\left[\mathbf{e}_{t^{a-}}^{S V T}, \mathbf{e}_{t^{a}}^{S V T}, \mathbf{e}_{t^{a+}}^{S V T}\right]^{T}, \Delta \mathbf{Y}_{t^{a}}^{S V}=\left[\Delta \mathbf{y}_{t^{a-}}^{S V T}, \Delta \mathbf{y}_{t^{a}}^{S V T}, \Delta \mathbf{y}_{t^{a+}}^{S V T}\right]^{T}$, and

\[\begin{split}\mathrm{H}_{Z}=\left[\begin{array}{ccc} \mathrm{H}_{z}\left(\mathbf{b}_{t^{a-}}^{a}\right) & -\delta t^{a} \mathrm{H}_{z}\left(\mathbf{b}_{t^{a-}}^{a}\right) & \frac{1}{2} \delta t^{a 2} \mathrm{H}_{z}\left(\mathbf{b}_{t^{a-}}^{a}\right) \\ \mathrm{H}_{z}\left(\mathbf{b}_{t^{a}}^{a}\right) & 0 & 0 \\ \mathrm{H}_{z}\left(\mathbf{b}_{t^{a+}}^{a}\right) & +\delta t^{a} \mathrm{H}_{z}\left(\mathbf{b}_{t^{a+}}^{a}\right) & \frac{1}{2} \delta t^{a} \mathrm{H}_{z}\left(\mathbf{b}_{t^{a+}}^{a}\right) \end{array}\right]\end{split}\]

with,

\[\mathbf{H}_{z}\left(\mathbf{b}_{t}^{a}\right)=\left[\begin{array}{ll} \mathbf{H}_{b} \tilde{\mathrm{A}}\left(\mathbf{b}_{t}^{a}\right) & \mathbf{H}_{b} \end{array}\right],\]

where, $\mathrm{A}$ is the matrix of Gaunt-Elsasser integrals [Moon79]

From an ensemble of $\mathbf{Z}_{t^{a}}^{f}$ we compute empirically the model state forecast covariance matrix $\mathrm{P}_{Z Z}^{f}\left(t^{a}\right)=\mathbb{E}\left(\delta \mathbf{Z}_{t^{a}}^{f} \delta \mathbf{Z}_{t^{a}}^{f T}\right)$. We use the glasso algorithm to approximate $\mathrm{P}_{Z Z}^{f}\left(t^{a}\right)$ [IGF23] (see remove_spurious in config).

The Kalman gain matrix writes:

\[\mathrm{K}_{Z}=\mathrm{P}_{Z Z}^{f} \mathrm{H}_{Z}^{T}\left[\mathrm{H}_{Z} \mathrm{P}_{Z Z}^{f} \mathrm{H}_{Z}^{T}+\tilde{\mathrm{R}}^{S V}\right]^{-1}\]

used for the analysis step for each ensemble member:

\[\mathbf{Z}_{t^{a}}^{a}=\mathbf{Z}_{t^{a}}^{f}+\mathrm{K}_{Z}\left(\Delta \mathbf{Y}_{t^{a}}^{S V}-\mathrm{H}_{Z} \mathbf{Z}_{t^{a}}^{f}\right)\]

Complete analysis step

The complete analysis step is therefore:

\[ \begin{align}\begin{aligned}\mathbf{B}_{t^{a}}^{a} &= \mathbf{B}_{t^{a}}^{f}+\mathrm{K}_{B}\left(\mathbf{Y}_{t^{a}}^{M F}-\mathrm{H}_{B} \mathbf{B}_{t^{a}}^{f}\right)\\\mathbf{Z}_{t^{a}}^{a} &= \mathbf{Z}_{t^{a}}^{f}+\mathrm{K}_{Z}\left(\Delta \mathbf{Y}_{t^{a}}^{S V}-\mathrm{H}_{Z} \mathbf{Z}_{t^{a}}^{f}\right)\\\mathbf{\dot{b}}_{t^{a}}^{a} &= \mathrm{A}(\mathbf{b}_{t^{a}}^{a})\mathbf{u}_{t^{a}}^{a} + \mathbf{e_r}_{t^{a}}^{a}\end{aligned}\end{align} \]

where, $\mathrm{A}$ is the matrix of Gaunt-Elsasser integrals [Moon79]

yielding the complete CoreState at $t_a$.

Restart of a forecast for an AR-3 model

Using an AR model of order higher than 3 requires to know, at each initialization of the forecast, the model state at 3 successive epochs.

During the analysis we have computed: $\mathbf{Z}_{t^{a}}^{a}=\left[\mathbf{z}_{t^{a}}^{a}, \mathbf{z}_{t^{a}}^{a^{\prime}}, \mathbf{z}_{t^{a}}^{a^{\prime \prime}}\right]$

So we restart from a Taylor expansion:

\[\begin{split}\left\{\begin{array}{rl} \mathbf{z}_{t^{a}} &= \mathbf{z}_{t^{a}}^{a}\\ \mathbf{z}_{t^{a}-\Delta t^{f}} & =\mathbf{z}_{t^{a}}^{a}-\Delta t^{f} \mathbf{z}_{t^{a}}^{a^{\prime}}+\frac{(\Delta t^{f})^{2}}{2}\mathbf{z}_{t^{a}}^{a^{\prime \prime}}\\ \mathbf{z}_{t^{a}-2 \Delta t^{f}} & =\mathbf{z}_{t^{a}}^{a}-2 \Delta t^{f} \mathbf{z}_{t^{a}}^{a^{\prime}}+\frac{(2 \Delta t^{f})^{2}}{2}\mathbf{z}_{t^{a}}^{a^{\prime \prime}} \end{array} .\right.\end{split}\]

References

[BGA17] (1,2)

Barrois, O., Gillet, N. & Aubert, J. Contributions to the geomagnetic secular variation from a reanalysis of core surface dynamics, Geophysical Journal International, 211, 50–68, 2017. doi:10.1093/gji/ggx280

[Moon79] (1,2,3,4,5)

Moon, W. Numerical evaluation of geomagnetic dynamo integrals (Elsasser and Adams-Gaunt integrals), Computer Physics Communications, 16, 267–271, 1979. doi:10.1016/0010-4655(79)90092-4

For more details on the above steps, we refer to [BGA17] and to:

[BHF18]

Barrois, O., Hammer, M. D., Finlay, C. C., Martin, Y. & Gillet, N. Assimilation of ground and satellite magnetic measurements: inference of core surface magnetic and velocity field changes, Geophysical Journal International, 2018. doi:10.1093/gji/ggy297

[BGA19]

Barrois, O. et al. Erratum: ‘Contributions to the geomagnetic secular variation from a reanalysis of core surface dynamics’ and ‘Assimilation of ground and satellite magnetic measurements: inference of core surface magnetic and velocity field changes’, Geophysical Journal International, 216, 2106–2113, 2019. doi:10.1093/gji/ggy471

[GHA19]

Gillet, N., Huder, L. & Aubert, J. A reduced stochastic model of core surface dynamics based on geodynamo simulations, Geophysical Journal International, 219, 522–539, 2019. doi:10.1093/gji/ggz313

[IGF23] (1,2,3,4)

Istas, M., Gillet, N., Finlay, C. C., Hammer, M. D., & Huder, L. Transient core surface dynamics from ground and satellite geomagnetic data, Geophysical Journal International, 233, 1890–1915, 2023. doi:10.1093/gji/ggad039