Three-dimensional variational assimilation with a multivariate background error covariance for the Model for Prediction Across Scales – Atmosphere with the Joint Effort for Data assimilation Integration (JEDI-MPAS 2.0.0-beta)

This paper describes the three-dimensional variational (3D-Var) data assimilation (DA) system for the Model for Prediction Across Scales – Atmosphere with the Joint Effort for Data assimilation Integration (JEDI-MPAS). Its core element is a multivariate background error covariance implemented through multiple linear variable changes, including a wind variable change from stream function and velocity potential to zonal- and meridional-wind components, a vertical linear regression representing wind–mass balance, and multiplication by a diagonal matrix of error standard deviations. The univariate spatial correlations for the “unbalanced” variables utilize the Background error on Unstructured Mesh Package (BUMP), which is one of the generic components in the JEDI framework. The variable changes and univariate correlations are modeled directly on the native MPAS unstructured mesh. BUMP provides utilities to diagnose parameters of the covariance model, such as correlation lengths, from an ensemble of forecast differences, though some manual adjustment of the parameters is necessary because of mismatches between the univariate correlation function assumed by BUMP and the correlation structure in the sample of forecast differences. The resulting multivariate covariances, as revealed by single-observation tests, are qualitatively similar to those found in previous global 3D-Var systems. Month-long cycling DA experiments using a global quasi-uniform 60 km mesh demonstrate that 3D-Var, as expected, performs somewhat worse than a pure ensemble-based covariance, while a hybrid covariance, which combines that used in 3D-Var with the ensemble covariance, significantly outperforms both 3D-Var and the pure ensemble covariance. Due to its simple workflow and minimal computational requirements, the JEDI-MPAS 3D-Var system can be useful for the research community.

1 Introduction

In the 1990s, three-dimensional variational (3D-Var) data assimilation (DA) became the algorithm of choice in operational numerical weather prediction centers (Parrish and Derber, 1992; Andersson et al., 1998; Gauthier et al., 1999; Lorenc et al., 2000), owing to its numerous advantages relative to earlier optimal-interpolation assimilation schemes. The 3D-Var system is no longer widely used operationally, both because its natural development path is to four-dimensional variational assimilation (Rabier et al., 2000; Rawlins et al., 2007) and because of the rapid development of ensemble data assimilation in the last 2 decades, including ensemble-variational techniques that employ sample covariances from a forecast ensemble within the variational framework (Lorenc, 2003; Buehner, 2005). The central components of 3D-Var systems, however, are so-called static covariance models (Bannister, 2008) that provide computationally tractable representations of complex spatial and multivariate covariances, and these remain in wide use to provide background covariances for 4D-Var and as part of hybrid techniques that consider background covariance matrices that are the sum of a static covariance and an ensemble-based covariance (Hamill and Snyder, 2000).

This paper documents 3D-Var and its associated static background covariance model for JEDI-MPAS, a data assimilation system using software infrastructure from the Joint Effort for Data assimilation Integration (JEDI; Trémolet and Auligné, 2020) and the Model for Prediction Across Scales – Atmosphere (MPAS; Skamarock et al., 2012). Two companion papers are Liu et al. (2022), which gives an overview of JEDI-MPAS and initial results from a three-dimensional ensemble-variational (3DEnVar) scheme, and Guerrette et al. (2023), which documents an ensemble of data assimilations (EDA) for JEDI-MPAS.

Our motivation for implementing 3D-Var is twofold. First, JEDI-MPAS is intended for use not only in our research, but also by the broader research community. The minimal computational cost of 3D-Var and its simple workflow make it well suited where computing is a strong constraint and when introducing new users to the system. Experience with WRFDA (Barker et al., 2012), our existing community DA system, has shown that 3D-Var is often preferred by users. Equally importantly, the static covariance model from 3D-Var can be used in hybrid ensemble-variational assimilation schemes that are known to outperform 3DEnVar alone (Wang et al., 2008; Buehner et al., 2013; Clayton et al., 2013; Kuhl et al., 2013). We show the same result for JEDI-MPAS here (see also Guerrette et al., 2023).

The formulation of the static covariance model employed here has both familiar and novel elements. We generally follow Wu et al. (2002), including (i) our choice of analysis variables, (ii) the use of linear regression from a training data set to define the approximate mass-wind balances that implicitly determine the multivariate structure of the covariances (see also Derber and Bouttier, 1999), and (iii) representing univariate correlations directly on the forecast model's grid (or mesh, in the case of MPAS). Products of vectors with the univariate spatial correlation matrices, however, are computed directly on a thinned subset of the MPAS mesh and interpolated to the full-resolution mesh using the Background error on Unstructured Mesh package (BUMP; Ménétrier, 2020). This study is the first evaluation of BUMP for use in atmospheric DA.

The outline of the paper is as follows. In the next section, we give an overview of JEDI-MPAS as configured for 3D-Var. Section 3 describes the formulation of the static background covariances and their tuning using BUMP capabilities and a training data set of forecast differences. Section 4 presents single-observation tests that illustrate the structure of the implied multivariate covariances. In Sect. 5, we summarize results from cycling DA experiments with 3D-Var and a hybrid scheme, which provide an overall evaluation of the effectiveness of the JEDI-MPAS static covariances. Section 6 concludes and offers ideas for further refinements of the static covariances.

2 JEDI-MPAS 3D-Var configuration

2.1 The forecast model

MPAS is described in detail in Skamarock et al. (2012), or see Liu et al. (2022) for a more concise summary. Briefly, MPAS integrates the nonhydrostatic equations of motion cast in a height-based, terrain-following vertical coordinate and using dry density and a modified moist potential temperature as thermodynamic variables. The equations are discretized on an unstructured mesh with the normal component of horizontal velocity defined on the edges of mesh cells and other prognostic variables defined at the cell centers. MPAS supports global and regional meshes, as well as meshes with quasi-uniform or variable resolution.

In all experiments presented here, MPAS is configured with a global quasi-uniform mesh of 60 km resolution and 55 vertical levels up to a model top of 30 km. The physical parameterizations are those of the “mesoscale reference” suite, as listed in Table 2 of Liu et al. (2022).

2.2 The DA system

JEDI-MPAS implements various abstract classes for MPAS within the JEDI framework (Trémolet, 2020). Those abstract classes reside in object-oriented prediction systems (OOPS) and comprise all the building blocks and operations on them necessary for data assimilation algorithms. JEDI also contains generic (model independent) implementations of some building blocks, including observation operators and quality control (the Unified Forward Operator, UFO), observation storage and access (the Interface for Observation Data Access, IODA), and the background error covariance matrix (the System-Agnostic Background Error Representation, SABER).

The variational application of OOPS minimizes a cost function given in Eq. (3) of Liu et al. (2022). Denoting by x the concatenation of the analysis variables across all model levels and mesh locations, the cost function measures the simultaneous fit of x to a background state x_b, which is our best estimate of x before considering the observations, and to the observations y (also concatenated into a single vector) via an observation operator h(x) that maps a given state to the observation variables.

The minimization proceeds iteratively by linearizing the observation operator in the neighborhood of the latest iterate x_g and computing the next iterate as the minimizer of the resulting quadratic cost function,

$$\begin{array}{}\text{(1)}& \begin{array}{rl}J\left(\mathit{\delta }\mathit{x}\right)& ={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}(\mathit{\delta }\mathit{x}-\mathit{\delta }{\mathit{x}}_g{)}^T{\mathbf{B}}^{-\mathrm{1}}(\mathit{\delta }\mathit{x}-\mathit{\delta }{\mathit{x}}_g)\\ & +{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}(\mathbf{H}\mathit{\delta }\mathit{x}-\mathit{d}{)}^T{\mathbf{R}}^{-\mathrm{1}}(\mathbf{H}\mathit{\delta }\mathit{x}-\mathit{d}),\end{array}\end{array}$$

where $\mathit{\delta }\mathit{x}=\mathit{x}-{\mathit{x}}_g$ is the increment relative to x_g, $\mathit{\delta }{\mathit{x}}_g={\mathit{x}}_{\mathrm{b}}-{\mathit{x}}_g$. $\mathit{d}=\mathit{y}-h\left({\mathit{x}}_g\right)$ is the observation departure from x_g; H is the linearization of h near x_g; and B and R are the background and observational error covariance matrices, respectively. This incremental formulation is central to the architecture of OOPS and distinguishes increments, which can be operated on by B and H, from the full state, which is an argument to h. In what follows, increments will be indicated by variables preceded with δ.

All the minimization schemes for Eq. (1) and preconditioners available in OOPS involve only the application of B to increments, rather than its square root or inverse. For the single-observation test and cycling experiments shown later, we employ the B-preconditioned incremental variational application of OOPS and the Derber–Rosati Inexact Preconditioned Conjugate Gradient algorithm (Golub and Ye, 1999; Derber and Rosati, 1989).

2.3 Analysis variables and variable change

The analysis variables are the horizontal velocity (v), temperature (T), specific humidity (q), and surface pressure p_s at the MPAS cell centers, as described in Liu et al. (2022). Transformations to other variables are necessary for some observation operators and for initial conditions for MPAS forecasts. Those transformations also follow Liu et al. (2022) but with one significant improvement when computing the dry density ρ_d and potential temperature θ_d for MPAS initial conditions.

In Liu et al. (2022), ρ_d and θ_d are computed from the analyzed T, p_s, and q by assuming hydrostatic balance. Here, we instead compute increments for ρ_d and θ_d (i.e., δρ_d and δθ_d) from the increments δT, δp_s, and δq. This approach, which is implemented by linearizing the corresponding calculations of Liu et al. (2022, steps 3 and 4 of their Sect. 3.3), assumes hydrostatic balance only for the increments and not the full analyzed fields.

Assuming hydrostatic balance just for the increments is preferable because that balance is only approximate and, moreover, the discretized form of hydrostatic balance used in the variable transformation is not precisely equivalent to that implied by the discrete MPAS equations. Since the hydrostatic integral is computed from the surface upward, differences between the incremental and full-field formulations can be expected to accumulate with height. Consistently with this, JEDI-MPAS cycling experiments (not shown) using the new incremental update for ρ_d and θ_d exhibit reduced temperature bias in the stratosphere, especially near the model top.

The state (x) and increment (δx) objects in JEDI-MPAS are basically inherited from MPAS's “pool type” data structure. Thus, it is natural to choose the existing mesh decomposition and communication utilities of MPAS to handle the parallelism of the state and increment of JEDI-MPAS. The state and increment objects in JEDI-MPAS only contain their values on their own grid point without a halo region. The halo exchange (and its adjoint) is performed when needed, such as horizontal interpolation of the state or increment to the observation location and a linear variation transform containing the spatial derivatives.

3 Multivariate background error covariance

In this section, we will present how the multivariate static background error covariance is designed for JEDI-MPAS. With a series of linear variable changes, the JEDI-MPAS analysis variables are transferred into a set of variables that are independent of each other. Then, we will describe how the B statistics (or parameters) can be trained from samples. The characteristics of diagnosed B statistics at the MPAS 60 km uniform mesh will be discussed. Lastly, we will discuss what modification is made to the diagnosed B statistics to resolve the discrepancy between the assumption and actual sample data set.

3.1 Multivariate background error covariance design

The basic design of the JEDI-MPAS multivariate background error covariance follows that of the Gridpoint Statistical Interpolation (GSI; Wu et al., 2002) system, except in our use of BUMP Normalized Interpolated Convolution from an Adaptive Subgrid (NICAS; Ménétrier, 2020), rather than recursive filters, to model the univariate spatial correlations (see further description at the end of this section). The multivariate covariances are implemented as two linear variable changes, K₁ and K₂, applied to a block-diagonal covariance matrix

$$\begin{array}{}\text{(2)}& \mathbf{B}={\mathbf{K}}_{\mathrm{1}}{\mathbf{K}}_{\mathrm{2}}\mathbf{\Sigma }\mathbf{C}\mathbf{\Sigma }{\mathbf{K}}_{\mathrm{2}}^T{\mathbf{K}}_{\mathrm{1}}^T,\end{array}$$

where the block-diagonal covariance matrix has been written as the product of a block-diagonal correlation matrix C and a diagonal matrix Σ of standard deviations.

The linear variable changes K₁ and K₂ can be expressed in the following matrix forms.

$$\begin{array}{}\text{(3)}& {\mathbf{K}}_{\mathrm{1}}:\left[\begin{array}{c}\mathit{\delta }\mathit{v}\\ \mathit{\delta }T\\ \mathit{\delta }q\\ \mathit{\delta }{p}_{\mathrm{s}}\end{array}\right]=\left[\begin{array}{ccccc}\mathit{k}\times \mathrm{\nabla }& -\mathrm{\nabla }& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathbf{I}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathbf{I}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathbf{I}\end{array}\right]\left[\begin{array}{c}\mathit{\delta }\mathit{\psi }\\ \mathit{\delta }\mathit{\chi }\\ \mathit{\delta }T\\ \mathit{\delta }q\\ \mathit{\delta }{p}_{\mathrm{s}}\end{array}\right]\end{array}$$

and

$$\begin{array}{}\text{(4)}& {\mathbf{K}}_{\mathrm{2}}:\left[\begin{array}{c}\mathit{\delta }\mathit{\psi }\\ \mathit{\delta }\mathit{\chi }\\ \mathit{\delta }T\\ \mathit{\delta }q\\ \mathit{\delta }{p}_{\mathrm{s}}\end{array}\right]=\left[\begin{array}{ccccc}\mathbf{I}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathbf{L}& \mathbf{I}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathbf{M}& \mathrm{0}& \mathbf{I}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathbf{I}& \mathrm{0}\\ \mathbf{N}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathbf{I}\end{array}\right]\left[\begin{array}{c}\mathit{\delta }\mathit{\psi }\\ \mathit{\delta }{\mathit{\chi }}_{\mathrm{u}}\\ \mathit{\delta }{T}_{\mathrm{u}}\\ \mathit{\delta }q\\ \mathit{\delta }{p}_{\mathrm{s},\mathrm{u}}\end{array}\right]\end{array}$$

Here, K₁ computes increments for v from spatial derivatives (indicated schematically by the gradient terms in the upper-left block) of increments of stream function ψ and velocity potential χ. K₁ and the corresponding adjoint operator, ${\mathbf{K}}_{\mathrm{1}}^T$, are model-dependent JEDI components that operate on the MPAS native mesh and are implemented in Control2Analysis, a linear variable change class. Details of the calculation of v from ψ and χ are given in the Appendix.

K₂ is a linear variable change that computes δχ, δT, and δp_s from δψ and the residual or “unbalanced” increments δχ_u, δT_u, and δp_s,u (for velocity potential, temperature, and surface pressure, respectively), so called because they are by assumption independent of δψ. The relation of δχ, δT, and δp_s to δψ is based on linear regression from a training data set, following Derber and Bouttier (1999). As in Derber and Bouttier (1999), we choose to use δψ at a given mesh cell and on the full set of vertical levels as predictors for δT and δp_s at the same mesh cell and on any specific level, which makes M and N block diagonal with blocks corresponding to mesh cells and full matrices in each block. We retain δψ only on the same level as a predictor for δχ, which makes L a diagonal matrix.

Lastly, ΣCΣ^T is the covariance matrix for δψ; δq; and the unbalanced increments δχ_u, δT_u, and δp_s,u. We assume these variables are mutually independent, so C is block diagonal with blocks that give the univariate spatial (horizontal and vertical) correlation for each variable. The matrix Σ is a diagonal matrix with elements that specify the standard deviation for δψ, δχ_u, δT_u, dq, and δp_s,u.

The operations K₂, ${\mathbf{K}}_{\mathrm{2}}^T$, Σ, Σ^T, and C use the Background error on Unstructured Mesh Package (BUMP) in the System-Agnostic Background-Error Representation (SABER) repository, which is a generic component within JEDI, through the MPAS model interfaces. The BUMP Vertical BALance (VBAL) driver is used for K₂ and ${\mathbf{K}}_{\mathrm{2}}^T$. It is based on the explicit vertical covariance matrices defined for a set of latitudes and is interpolated at the model grid points' latitude. The BUMP VARiance (VAR) driver is used for Σ and Σ^T. It simply applies the pre-computed error standard deviations. The spatial correlation matrix is pre-computed from the given correlation lengths with BUMP–NICAS. Similarly to the GSI recursive filters, NICAS works in the grid-point space. However, it applies the convolution function explicitly, instead of recursively for GSI. Thus, the choice of the convolution function in NICAS is free, as long as it is positive definite. We choose a widely used fifth-order piecewise function of Gaspari and Cohn (1999), which resembles the Gaussian function but is compactly supported. To make the explicit convolution affordable for high-dimensional systems, it is actually performed on a low-resolution unstructured mesh. A linear interpolation is required to go from the unstructured mesh to the full model grid. Finally, an exact normalization factor is pre-computed and applied to ensure that the whole NICAS correlation operator is normalized (i.e., diagonal elements of the equivalent correlation matrix are “1”). Thus, the NICAS correlation matrix can be written as $\mathbf{C}=\mathbf{NS}\stackrel{\mathrm{̃}}{\mathbf{C}}{\mathbf{S}}^T{\mathbf{N}}^T$, where $\stackrel{\mathrm{̃}}{\mathbf{C}}$ is the convolution operator on the low-resolution mesh, S is the interpolation from the mesh to the full model grid, and N is the diagonal normalization operator. The low-resolution mesh density can be locally adjusted depending on the diagnosed correlation lengths (or provided by the user). Figure 1 shows a diagram for Eq. (2), with corresponding BUMP drivers and MPAS-specific linear variable change.

Figure 1 A diagram for Eq. (2). The symmetric correlation matrix, C, is wrapped with a series of forward linear (or tangent linear, TL) variable changes and corresponding adjoint (AD) variable changes.

3.2 Training the covariance model

The designed multivariate background error covariance has several parameters to be determined. These parameters are diagnosed from 366 samples of National Centers for Environmental Prediction (NCEP) Global Forecast System (GFS) 48 and 24 h forecast differences, valid at the same time of day and spanning the months of March, April, and May 2018. Here, the 24 h forecast lead time difference is chosen to remove the effect of the diurnal cycle in the perturbation samples. The GFS forecast files of 0.25° resolution on the pressure levels are interpolated to the 60 km MPAS mesh with 55 vertical levels for the following training procedures. With a recent (early June 2023) version of JEDI-MPAS source code after initial submission of this paper, we have trained the static B parameters from MPAS model's own forecast with the same methodology described here. The overall B statistics diagnosed from MPAS-based samples were similar to that from GFS-based samples reported here, except for the error standard deviations in the stratosphere, which were larger for MPAS-based samples. In the 1-month cycling experiment, this led to a reduction in temperature and wind root mean square errors (RMSEs) in 6 h forecasts in the stratosphere.

Since ΣCΣ^T depends on the statistics of δψ and δχ, we first need to calculate perturbations of ψ and χ from δv of each forecast difference. This is essentially the inverse operation of K₁ and can be expressed as solving a Poisson equation with vorticity or divergence as a source term. Because solving a Poisson equation efficiently on the unstructured grid is not an easy task, we have adopted a spherical-harmonics-based function from the NCAR Command Language (NCL, 2019) that operates on an intermediate latitude–longitude grid. We begin by interpolating δv fields to the intermediate grid and then calculate δψ and δχ through the uv2sfvpf function of NCL and interpolate back to the MPAS mesh. Note that because the definition of δχ (shown in Eq. 3) is opposite in sign to that assumed in the NCL function, multiplying δχ (from NCL) by −1 is required.

For K₂, the BUMP VBAL driver calculates the cross-variable linear regression coefficients, which are denoted as L, M, and N in Eq. (4). The vertical autocovariance matrix of δψ and the vertical cross-covariance matrices between δψ and each of δχ, δT, and δp_s are computed and averaged over latitude bands of ±10°. The desired matrices of regression coefficients are obtained in the standard way by right-multiplying the cross-variable covariance by the inverse autocovariance of the predictor variable (δψ in our design). The vertical autocovariance matrices are usually poorly conditioned, and thus direct computation of their inverses will yield noisy results in the presence of sampling error. To overcome this, we apply a pseudo-inverse, which only includes some dominant eigenmodes to calculate the inverse matrix. We have chosen the 20 leading modes (among a total of 55 modes) for the pseudo-inverse of the δψ autocovariance matrix.

For Σ, the BUMP VAR driver calculates variances for δψ, δχ_u, δT_u, δq, and δp_s,u from the samples and filters them horizontally to damp the sampling noise. The horizontal smoother is also based on NICAS, with an appropriate mean-preserving normalization factor.

The correlation matrix, C, consists of blocks that specify the univariate spatial correlation for δψ, δχ_u, δT_u, δq, and δp_s,u. The BUMP Hybrid DIAGnostic (HDIAG) driver diagnoses the horizontal and vertical correlation lengths used in modeling C parameters. HDIAG can diagnose the horizontal and vertical spatial correlations from the samples. First, it defines a low-resolution unstructured mesh. Around each mesh node, diagnostic points are randomly and isotropically drawn for different horizontal separation classes. Second, HDIAG calculates the horizontal correlation between each mesh node and its own diagnostic points from the samples, at all levels. The vertical correlation is also calculated at each mesh node, between each level and the neighboring levels. The third step is a horizontal averaging of these raw correlations, either over all the mesh nodes or over local neighborhoods. The average is binned depending on the level and the horizontal separation for the horizontal correlation and depending on the levels concerned for the vertical correlation. As a final step, HDIAG fits a Gaspari and Cohn (1999) function for each averaged correlation curve. Thus, we obtain horizontal and vertical correlation length values for each level. If the averaging and curve-fitting steps are performed over local neighborhoods, an extra interpolation step is necessary to obtain 3D fields of length scales on the model grid. These length-scale profiles or 3D fields can be stored and provided to NICAS in order to model the spatial correlation operator. In this study, the local correlation lengths were obtained from raw statistics within a 3000 km radius for a given diagnostic point.

Figure 2 Regression coefficients near 34.8° N between (a) δT and δψ $\left[\mathrm{K}\right({\mathrm{m}}^{\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}{)}^{-\mathrm{1}}]$, (b) δχ and δψ $\left[\right({\mathrm{m}}^{\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}\left)\right({\mathrm{m}}^{\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}{)}^{-\mathrm{1}}]$, and (c) δp_s and δψ $\left[\mathrm{Pa}\right({\mathrm{m}}^{\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}{)}^{-\mathrm{1}}]$. These are the nonzero elements at this mesh cell of the submatrices M, L, and N, respectively, of K₂ (see Eq. 4).

3.3 Diagnosed statistics

The regression coefficients that appear in the definition (4) of K₂, which are computed by BUMP VBAL, are shown in Fig. 2 for a location near 34.8° N latitude. Considering first Fig. 2a, the δT–δψ coefficients are largest at small separations. Their structure is dipolar, with δT at a given level positively related to δψ at nearby but higher levels, and negatively related to δψ at nearby but lower levels. The δT–δψ coefficients are generally consistent with approximate geostrophic and hydrostatic balance, which together relate ψ to the mass field and the vertical derivative of the mass field to buoyancy. The coefficient structure is different for model levels lower than 10, perhaps due to the effects of the boundary layer and terrain. Figure 2b shows δχ–δψ coefficients, which relate δχ at a given level to δψ at the same level. The balanced part of δχ depends positively on δψ near the surface, consistent with Ekman balance, under which vertical vorticity near the surface drives horizontal convergence in the boundary layer. Finally and not unexpectedly, δp_s has a positive dependence on δψ in the lower troposphere (Fig. 2c).

Figure 3 Ratio of balanced variance (i.e., that predicted by δψ) to total variance for (a) δT, (b) δχ, and (c) δp_s.

Figure 3 shows the variance that can be predicted knowing δψ normalized by the total variance, for δT, δχ, and δp_s. There are substantial variations with latitude and height. For δT, the δT–δψ regression can explain up to 70 % of the total variance in mid-latitude and high-latitude regions in the troposphere. For δχ, the regression explains up to 35 % of the total variance in the mid-latitudes near the surface (below model level 10), while for δp_s, the regression explains substantial variance everywhere except the tropics. These balanced ratios are similar to those found in Wu et al. (2002) (their Fig. 1) and Barker et al. (2012) (their Fig. 5), and their geographic variations confirm that the regressions primarily reflect dynamical balances characteristic of mid-latitudes and high latitudes.

Figure 4 Horizontally averaged standard deviations for (a) δψ, (b) δχ_u, (c) δT_u, and (d) δq. The horizontally averaged standard deviation for δp_s,u is 53.05 Pa.

The other quantities that must be estimated from the training data set are the standard deviations that form the diagonal of Σ and the fields of horizontal and vertical correlation scales for each of the variables δψ, δχ_u, δT_u, δq, and δp_s,u, which together specify the correlation matrix C. Figure 4 shows the vertical profiles of horizontally averaged standard deviations for each variable. For δψ and δχ_u, the standard deviation increases upward from the surface to a peak near the tropopause. The standard deviation of δT_u, in contrast, generally decreases upward from a peak at the surface. For δq, the profile of standard deviations has a shape similar to that for q itself, peaked in the lower troposphere and decreasing steadily with height above.

Figure 5 Horizontally averaged horizontal correlation lengths [km] for (a) δψ, (b) δχ_u, (c) δT_u, and (d) δq. The horizontally averaged length for δp_s,u is 3702.3 km.

Figure 6 Horizontally averaged vertical correlation lengths [km] for (a) δψ, (b) δχ_u, (c) δT_u, and (d) δq.

Figures 5 and 6 show the vertical profiles of horizontally averaged horizontal and vertical correlation lengths, respectively. The horizontal correlation lengths generally increase with height in the stratosphere and are nearly constant in the troposphere, though δχ_u has substantial variations in horizontal length scale throughout the profile. The horizontal correlation lengths for δψ and δχ_u are larger than those for δT_u and δq, while the horizontal correlation length for δp_s,u is roughly 3700 km, much larger than the horizontal correlation lengths for δT_u and δq near the surface. The vertical correlation lengths have minima near the surface for all variables and then increase with height toward a peak near the model top. The vertical correlation lengths for δψ and δχ_u are again larger than those for δT_u and δq.

3.4 Additional tuning

The parameters shown in the previous section are the raw statistics from the BUMP training applications. We have applied two additional modifications to these raw statistics. Without these modifications, the resulting static B performs poorly in 3D-Var and is unable to improve on 3DEnVar in hybrid applications (not shown).

First, the background error standard deviation for all variables (Σ) is scaled by a factor of $\mathrm{1}/\mathrm{3}$. While the cycling interval shown later in this study is 6 h, which is typical, the forecast differences used in the training reflect forecast-error evolution over 24 h. Thus, it is reasonable to scale the diagnosed background error standard deviation to match the error growth for a 6 h interval. Here, we choose the single scaling factor of $\mathrm{1}/\mathrm{3}$ for all variables, based on a limited set of sensitivity tests of cycling experiments with different scaling factors.

Figure 7 The isotropic correlation function for δψ on the 15th model level, based on the sample of forecast differences (blue). Also shown are the correlation function assumed in BUMP (orange) using the length scale that gives the best fit to the sample-derived correlation, as well as the tuned correlation function (green).

We also reduce the diagnosed horizontal correlation lengths for δψ and δχ_u by half. Figure 7 shows the raw horizontal correlation function for δψ on model level 15 together with the best-fit correlation function, which is diagnosed by BUMP by adjusting the length scale of the fifth-order, compactly supported function from Gaspari and Cohn (1999). There is a clear discrepancy between the sample correlation function and that assumed in BUMP – the raw correlation decreases more rapidly for small separations and has larger correlations at large separations. The implied velocity variance in the modeled covariance (Eq. 2) is proportional to the second derivative at the origin of the δψ (and δχ) correlation (Lorenc, 1981; Daley, 1985). That is, δψ correlations that more strongly peak at the origin will produce larger velocity variance even if the δψ variance is fixed. Thus, the modeled covariances greatly underestimate the velocity variance relative to the statistics of the original training data. Reducing the horizontal correlation length for δψ and δχ_u by a factor of 2 increases the second derivative of their correlation, and therefore the velocity variance, by a factor of 4, leading to a better fit to the velocity variance in the training data. Ideally, the mismatch between the assumed correlation structure and that of the training data would be addressed by a more flexible correlation model in BUMP. A capability for this is now available, and we hope to report on its use in the future. In other systems, the necessary flexibility has been achieved using sums of recursive filters with different length scales (Wu et al., 2002; Kleist et al., 2009) or modeling the correlations in spectral space (Parrish and Derber, 1992; Lorenc et al., 2000).

4 Single-observation tests

To explore the structure of diagnosed and tuned multivariate B, two sets of a single-observation test were performed; one for assimilating a single zonal-wind observation with 1 m s⁻¹ innovation and 1 m s⁻¹ observation-error standard deviation and the other for assimilating a single temperature observation with 1 K innovation and 1 K observation-error standard deviation. Both single observations are placed at (40.4113° N, 38.68° W) and at roughly 800 hPa, a location that corresponds to one of the MPAS cell centers and vertical level 15.

Figure 8 Analysis increments for the (left column) zonal component of v [m s⁻¹], (center column) meridional component of v [m s⁻¹], and (right column) T [K] on model level (upper row) 20, (middle row) 15, and (lower row) 10, from a single temperature observation with 1 K innovation and 1 K observation-error standard deviation, located at (40.41° N, 38.68° W) on model level 15 with a marker ×.

Figure 8 shows the analysis increments from the single temperature observation, for the zonal and meridional components of v (left and center columns, respectively) and for T (right column), on the three different vertical levels (10, 15, and 20, shown in the bottom, middle, and top rows). The T increments have a horizontally isotropic structure with maximum values at level 15, where the observation is located. The wind increments are, to a first approximation, linked to the T increment through the thermal-wind relation: cyclonic circulation is introduced on level 10, below the maximum temperature increment, and an anti-cyclonic circulation appears above, on level 20. This reflects the approximate geostrophic and hydrostatic relations captured by K₂ and is consistent with Parrish and Derber (1992) (their Fig. 2), which uses the linear balance equation between mass and momentum variables.

Figure 9 Same as Fig. 8, except from a single zonal-wind observation with 1 m s⁻¹ innovation and 1 m s⁻¹ observation-error standard deviation, located at (40.41° N, 38.68° W) on model level 15 with a marker ×.

Similarly, Fig. 9 shows the analysis increments from the single zonal-wind observation. The positive zonal-wind increment is maximized at the observation location on the 15th model level. To the north and south of the observation location, negative zonal-wind increments are introduced, which – together with the increments of meridional winds – represent a cyclonic circulation to the north of the observation and an anti-cyclonic circulation to the south. Temperature increments are negative below the cyclonic circulation (i.e., on level 10) and positive above (on level 20) and vice versa for the anti-cyclonic circulation. The structure of the increments again reflects thermal-wind balance, and in this case is consistent with Wu et al. (2002) (their Fig. 9) or Kleist et al. (2009) (their Fig. 3).

Figure 10 Analysis increments for p_s [Pa] (left) from a single temperature observation with 1 K innovation and 1 K observation-error standard deviation and (right) from a single zonal-wind observation with 1 m s⁻¹ innovation and 1 m s⁻¹ observation-error standard deviation. The observation location is marked with ×.

Figure 10 shows the surface pressure increments from two single-observation tests. When the single temperature observation is assimilated, cyclonic circulation is introduced into the lower troposphere. The negative surface pressure increment is approximately geostrophically related to this cyclonic circulation. When the single zonal-wind observation is assimilated, the zonal-wind increments extend throughout the troposphere, including to the surface. The dipole of positive and negative surface pressure increments, south and north, respectively, of the observation location, is geostrophically related to the increment of the surface wind.

5 Cycling experiments

5.1 Experimental design and assimilated observations

For further evaluation of the multivariate static B for JEDI-MPAS, three sets of 1-month (15 April–14 May 2018) cycling experiments were performed on NCAR's high-performance computing system, Cheyenne. As a reference experiment, the “3DEnVar” experiment was performed using the pure ensemble covariances, as in Liu et al. (2022). At each cycle, a 20-member ensemble of 6 h MPAS forecasts was made using initial conditions from the Global Ensemble Forecast System (GEFS; Zhou et al., 2017). Covariance localization was applied to the ensemble covariances via BUMP's generic localization scheme, using globally constant localization scales of 1200 km horizontally and 6 km vertically. To demonstrate the static covariances, the “3D-Var” experiment was performed with the static B formulated and tuned as described in the preceding sections. Lastly, the “Hybrid-3DEnVar” experiment was performed using a hybrid covariance given by a weighted sum of static and ensemble B (Hamill and Snyder, 2000). Here, we choose a weight of 0.5 for each component, similarly to Wang et al. (2013), Clayton et al. (2013), and Kuhl et al. (2013). This final experiment evaluates the effectiveness of our static B for hybrid applications. In all three experiments, the same global MPAS quasi-uniform 60 km mesh is used both for analysis and background fields and for analysis increment. For the minimization, two outer loops are used, with 60 inner iterations for each outer loop.

The observation files are converted from GSI's ncdiag files, which contain the observation location, observation value, observation error, GSI's quality control, and satellite bias correction information. The observation quality control basically follows the GSI's quality control (called “PREQC”), and the background innovation check is added, which filters out the observation when the absolute value of observation departure is larger than 3 times the given observation error. In all three experiments, the surface pressure, radiosondes (wind, temperature or virtual temperature, and specific humidity), aircraft (wind, temperature, specific humidity), atmospheric motion vectors, Global Navigation Satellite System Radio Occultation (GNSS RO) refractivity, and clear-sky Advanced Microwave Sounding Unit-A (AMSU-A) radiances from NOAA-15, NOAA-18, and NOAA-19; MetOp-A and MetOp-B; and Aqua satellites are assimilated. The AMSU-A radiances are bias-corrected from GSI's information and pre-thinned with a 145 km mesh.

Figure 11 Time series (00:00 UTC, 15 April, to 18:00 UTC, 14 May 2018) of background RMSEs for p_s verified with GFS analysis over (a) northern extratropical (30–90° N; NXTro), (b) tropical (30° S–30° N; Tro), and (c) southern extratropical (30–90° S; SXTro) regions.

5.2 Results

Figure 11 shows the time series of background RMSEs for surface pressure during the cycling period. The RMSEs are calculated with respect to GFS analysis at the valid time as a reference. For surface pressure, 3D-Var gives somewhat smaller RMSEs compared to the 3DEnVar experiment, except for the southern extratropical region. Hybrid-3DEnVar gives smaller RMSEs compared to 3D-Var. Figure 12 shows the vertical profiles of relative RMSE changes for background fields during the cycling period, with the RMSE of 3DEnVar as reference. The confidence intervals with 95 % significant level are also shown as error bars, from a bootstrap resampling method with a resampling size of 10 000. Compared to 3DEnVar, the 3D-Var experiment shows some degradation in the troposphere and some improvement in the stratosphere in general. Hybrid-3DEnVar shows overall improvement over both 3DEnVar and 3D-Var experiments and throughout levels in both the troposphere and the stratosphere.

Figure 12 Vertical profiles of relative background RMSE changes (with respect to GFS analysis) for (a) T, (b) Q_v, (c) zonal wind, and (d) meridional wind, compared to the 3DEnVar experiment. Statistics are aggregated for the period from 00:00 UTC, 18 April, to 18:00 UTC, 14 May 2018, with 95 % confidence intervals.

Figure 13 Vertical distribution of relative rms changes of observation minus background, or first-guess departure, for the (a) virtual temperature T_v, (b) specific humidity (c) zonal wind, and (d) meridional wind of radiosonde observations. Statistics are aggregated for the period from 00:00 UTC, 18 April, to 18:00 UTC, 14 May 2018, with 95 % confidence intervals.

Figure 13 shows the observation space verification for radiosondes. The relative changes of the root mean square (rms) first-guess departure (OMB, observation minus background) are mostly consistent with the model space verification in Fig. 12. Compared to the RMSEs of 3DEnVar, the RMSEs of Hybrid-3DEnVar are significantly improved, except for temperature observation. The RMSEs of 3D-Var are degraded by ∼5 %. In the observation space verification for AMSU-A radiance observations (not shown), which assimilates the channels sensitive to the atmospheric temperature profile, Hybrid-3DEnVar shows neutral to slightly improved impact in the RMSEs for channels 5 and 6. For channels 7, 8, and 9, both 3D-Var and Hybrid-3DEnVar show significant improvement over 3DEnVar. Larger improvement is shown over both high latitudes. This is consistent with the large temperature RMSE reduction in the model space verification (Fig. 12a).

Figure 14 Relative RMSE changes as a function of forecast lead time for (a–d) near-surface and (e–h) three-dimensional variables, compared to the 3DEnVar experiment. Shown are (a) T at 2 m, (b) q at 2 m, (c) zonal wind at 10 m, (d) p_s, (e) T, (f) Q_v, (g) zonal wind, and (h) meridional wind. Statistics are aggregated over 27 extended forecasts with 95 % confidence intervals.

Additional 10 d extended forecasts were conducted at each 00:00 UTC initialization time to evaluate the impact of analysis on the longer forecast lengths. The changes in RMSE for 3D-Var and Hybrid-3DEnVar relative to 3DEnVar are shown in Fig. 14. At short forecast lead times, the relative RMSEs look similar to the relative RMSEs for the 6 h background forecasts shown in Figs. 11 and 12. It is notable that the larger water vapor mixing ratio (Q_v) RMSE for 3D-Var lasts until 6 d forecasts (Fig. 14f). This might be because the moisture variable is univariate in the current B design (Sect. 3.1), while the moisture analysis in 3DEnVar or Hybrid-3DEnVar can be done through the multivariate ensemble covariances. The benefit of hybrid background error covariance can be found for up to 5 d lead times for surface pressure, temperature, and zonal- and meridional-wind fields. The benefit of hybrid covariance is only kept for ∼2 d lead times for humidity fields.

6 Conclusions

This study has described the multivariate static background error covariances for JEDI-MPAS 3D-Var. Similarly to Liu et al. (2022), JEDI-MPAS 3D-Var utilizes generic JEDI components through interfaces that are specific to MPAS.

The formulation of the JEDI-MPAS static B generally follows Wu et al. (2002) but with the novel use of BUMP for multiple elements of the covariance model. Two linear variable transforms represent the multivariate relationship. One transform is a variable change from stream function ψ and velocity potential χ to v, which directly operates on the MPAS native mesh. The other transform, which uses the BUMP driver VBAL from JEDI's SABER repository, is based on linear regression over vertical columns of increments in other variables against increments in ψ. The full multivariate covariances are then given by these linear transforms (and their adjoints) applied to a univariate covariance. The univariate correlation matrix employs BUMP–NICAS, which efficiently computes the three-dimensional convolution of an input vector with a specified spatial correlation function on an optimally subsampled mesh and then interpolates back to the full MPAS mesh.

For the experiments presented here, we estimated various parameters in the covariance model from a training data set of 366 differences between 48 and 24 h forecasts on the 60 km MPAS mesh. In general, the regression coefficients capture the linear, approximately geostrophic balance between mass and momentum variables that holds outside the tropics. While the error standard deviations for ψ and χ_u become larger at higher vertical levels, up to a peak near the tropopause, the error standard deviations for T_u and q are larger in the middle to lower troposphere. The horizontal and vertical correlation lengths for errors in ψ and χ_u are, in general, larger than those for errors in T_u and q. We also made two modifications to the parameters objectively estimated by BUMP. The error standard deviation is scaled by a factor of $\mathrm{1}/\mathrm{3}$ to match the 24 h time differences (i.e., 48 h forecast and 24 forecast pairs) in the training samples to the typical 6 h DA cycling interval. In addition, the horizontal correlation lengths for increments of ψ and χ_u are reduced by half to compensate for the discrepancy between the raw correlation structure from the training data set and the correlation function assumed in BUMP, which has much less curvature at small separations.

We evaluated the JEDI-MPAS static B in cycling data assimilation experiments extending over a month and assimilating observations from a significant fraction of the global observing network, including GNSS RO, AMSU-A, and conventional observations. The 3D-Var system using this static B generally performs close to, but worse than, EnVar using pure ensemble covariances as in Liu et al. (2022), while using a hybrid background covariance that is a weighted sum of the static B and ensemble covariances improves significantly on both 3D-Var and EnVar. Neither of these results is novel, as numerous studies have shown the advantage of EnVar over 3D-Var and of the hybrid algorithm over EnVar, but they do demonstrate clearly the effectiveness of the static B.

The static background covariances presented here are an initial implementation, with plenty of room for further refinements. Two extensions that are already underway are training the covariance model based on an ensemble from JEDI-MPAS, such as that provided by the EDA of Guerrette et al. (2023), and including hydrometeor increments, which will be especially important for all-sky assimilation of radiances. There are also BUMP capabilities that we have yet to exercise, including more general correlation functions that should remove the need for manual retuning of correlation lengths diagnosed by BUMP and joint estimation of hybridization and localization coefficients (Ménétrier and Auligné, 2015). In the current B design, the moisture variable is univariate, which may limit the performance of moisture analysis even with Hybrid-3DEnVar configuration. To this end, we will explore more sophisticated moisture variables, such as pseudo-relative humidity (Dee and da Silva, 2003) or normalized relative humidity (Hólm et al., 2002).

Appendix A: Diagnosing velocity from the stream function and velocity potential in MPAS

To compute horizontal velocity v from ψ and χ on the MPAS mesh, we rely on the fact that the irrotational component of velocity is given by (minus) the gradient of χ, while the nondivergent component of velocity is given by the cross-product of the vertical unit vector and the gradient of ψ. The edge-normal component of velocity u on the MPAS mesh is oriented parallel to the segment connecting the centers of adjacent cells and, naturally, normal to the edge itself. The natural finite difference relation is then

$$\begin{array}{}\text{(A1)}& u=-{\mathit{\delta }}_{\mathrm{c}}\mathit{\chi }-{\mathit{\delta }}_{\mathrm{v}}\mathit{\psi },\end{array}$$

where, in an abuse of our previous notation, δ_c is the difference operator between the centers of the mesh cells adjoining the edge and δ_v is the difference operator between the cell vertices at either end of the edge.

Implementing Eq. (A1) is straightforward given the MPAS mesh conventions and the mesh information in MPAS initialization files (MPAS Mesh Specification, 2015). The difference operators are defined as

$$\begin{array}{}\text{(A2)}& \begin{array}{rl}{\mathit{\delta }}_{\mathrm{c}}\mathit{\chi }& ={\displaystyle \frac{{\mathit{\chi }}_{\mathrm{c},\mathrm{2}}-{\mathit{\chi }}_{\mathrm{c},\mathrm{1}}}{\mathrm{\Delta }c}},\\ {\mathit{\delta }}_{\mathrm{v}}\mathit{\psi }& ={\displaystyle \frac{{\mathit{\psi }}_{\mathrm{v},\mathrm{2}}-{\mathit{\psi }}_{\mathrm{v},\mathrm{1}}}{\mathrm{\Delta }v}},\end{array}\end{array}$$

where Δc is the distance between cell centers sharing the edge and Δv is the distance between vertices on the edge. The ordering of the cells and vertices is such that this formula will give the correct signs for the velocities on the MPAS mesh. For a given edge in the files, the lengths Δc and Δv are in the variables dcEdge(edge) and dvEdge(edge), respectively. The cells sharing an edge are cellsOnEdge(2,edge) and cellsOnEdge(1,edge), and the vertices are verticesOnEdge(2,edge) and verticesOnEdge(1,edge).

Although Eqs. (A1) and (A2) use the full variables (i.e., ψ and χ), they are also applicable to the incremental variables (i.e., δψ and δχ) because of their linear form. The computation begins with δψ and δχ at the cell centers. Values of δψ at the cell vertices are computed by interpolating δψ from the centers of the three cells containing each vertex, before applying Eq. (A1). MPAS provides a utility that employs radial basis functions to reconstruct the vector wind at a cell center from δu on the edges of the cell (following Bonaventura et al., 2011), which provides the final δv. Because these three steps involve a different subset of the MPAS mesh (i.e., cells, vertices, and edges), a halo exchange routine from MPAS is required and called between each step. We have also implemented the adjoint operator of the halo exchange, which is needed when applying ${\mathbf{K}}_{\mathrm{1}}^T$.