Incremental analysis update (IAU) in the Model for Prediction Across Scales coupled with the Joint Effort for Data assimilation Integration (MPAS–JEDI 2.0.0)

In a cycling system where data assimilation (DA) and model simulation are executed consecutively, the model forecasts initialized from the analysis (or data assimilation) can be systematically affected by dynamic imbalances generated during the analysis process. The high-frequency noise arising from the imbalances in the initial conditions can impose constraints on computational stability and efficiency during subsequent model simulations and can potentially become the low-frequency waves of physical significance. To mitigate these initial imbalances, the incremental analysis update (IAU) has long been utilized in the cycling context. This study introduces our recent implementation of the IAU in the Model for Prediction Across Scales – Atmospheric (MPAS-A) coupled with the Joint Effort for Data assimilation Integration (JEDI) through the cycling system called MPAS-Workflow. During the integration of the compressible nonhydrostatic equations in MPAS-A, analysis increments are distributed over a predefined time window (e.g., 6 h) as fractional forcing at each time step. In a real case study with the assimilation of all conventional and satellite radiance observations every 6 h for 1 month, starting from mid-April 2018, model forecasts with the IAU show that the initial noise illustrated by surface pressure tendency becomes well constrained throughout the forecast lead times, enhancing the system reliability. The month-long cycling with the assimilation of real observations demonstrates the successful implementation of the IAU capability in the MPAS–JEDI cycling system. Along with the comparison between the forecasts with and without the IAU, several aspects regarding the implementation in MPAS–JEDI are discussed. Corresponding updates have been incorporated into the MPAS-A model (originally based on version 7.1), which is now publicly available in MPAS–JEDI and MPAS-Workflow version 2.0.0.

1 Introduction

Data assimilation (DA) is a mathematical or statistical procedure that incorporates observations, unevenly distributed in time and space, into adjacent grids in the model forecast (or background) using relative weights based on the error statistics of the forecasts and observations. It is not required to account for dynamical or physical balances across model grids or variables, nor does it ensure the conservation of mass, momentum, or energy. Consequently, the initial balance of the atmospheric flow can be disrupted by data assimilation when the initial state is replaced by the analysis state. Such imbalances can introduce artificial high-frequency noise, which is amplified and propagated throughout the model simulation.

In the cycling system that alternates between analysis and forecast, noise can continuously accumulate through cycles, which can degrade numerical stability and efficiency (with the time step smaller than 6Δx for a nominal grid spacing of Δx). If the noise is not properly controlled and its nonlinear interactions with lower-frequency modes of physical interest are triggered, subsequent forecasts can be contaminated. If the forecast error growth is accelerated at each cycle, it not only limits the predictability of the atmospheric flow (Hohenegger and Schär, 2007) but can eventually cause the model simulation to crash.

To mitigate abrupt changes (or shocks) originating from inconsistencies or imbalances in the initial state, the incremental analysis update (IAU) method was introduced by Bloom et al. (1996) and has been widely used in leading operational or research centers (Polavarapu et al., 2004; Buehner et al., 2015; Lorenc et al., 2015; Lei and Whitaker, 2016; Ha et al., 2017). In this method, instead of initializing a numerical prediction model with the analysis state, analysis increments (e.g., analysis minus background) are distributed over a certain time window as fractional forcing applied to each time step during the model integration. Compared to other conventional continuous data assimilation approaches such as nudging, the IAU has an attractive time-filtering feature, which can suppress spurious noise caused by data assimilation, providing a more consistent and balanced initial state before the next assimilation cycle.

The Model for Prediction Across Scales – Atmosphere (MPAS-A or MPAS, hereafter) utilizes centroidal Voronoi tessellations for its horizontal meshes and terrain-following height as a vertical coordinate and employs a fully compressible nonhydrostatic model with the capability of high-resolution forecasting over a local area for both global and regional applications (Skamarock et al., 2012, 2018). Ha et al. (2017) first implemented the IAU in MPAS version 4 for a global cycling system with the ensemble Kalman filter (EnKF) in the Data Assimilation Research Testbed (DART; Anderson et al., 2009) and demonstrated that it effectively suppressed spurious high-frequency noise, leading to better forecast skills. Since then, the model has been significantly upgraded for various new features, including acoustic filtering (Klemp et al., 2018), which effectively damps out noise in the acoustic waves using horizontally explicit, vertically implicit (HEVI) and split-explicit time integration schemes. The unstructured MPAS mesh can be configured with horizontally variable-resolution meshes that vary smoothly between low- and high-resolution regions, but the model integration time step is limited by the finest grid spacing of the mesh and uniformly applied across the domain to maintain the stability in the fine-mesh region.

Through the collaborative work with the Joint Effort for Data assimilation Integration (JEDI) team, an interface between MPAS and JEDI has been recently developed for a new community data assimilation system based on the variational approach (Liu et al., 2022). MPAS–JEDI (or JEDI–MPAS) shares with other geophysical models the object-oriented framework for generic components for data assimilation such as the Interface for Observation Data Access (IODA; Honeyager et al., 2020), the Unified Forward Operator (UFO; Honeyager et al., 2020), the background error covariance models through the System Agnostic Background Error Representation (SABER; http://data.jcsda.org/doxygen/Release/saber/1.2.0/index.html, last access: 20 December 2023), and several minimization algorithms.

Since all of the analysis variables in MPAS–JEDI are either derived or diagnostic variables in the MPAS model, their analysis increments need to be converted back to the model's prognostic variables after the minimization so that the analysis updates can be reflected in the model forecast. The variable transformations of mass fields are carried out based on the reconstructed pressure coordinate derived from surface pressure, assuming hydrostatic balance and an ideal gas state. The recent version of MPAS–JEDI is updated to transform analysis increments to the increments of the model's prognostic variables (instead of their full fields), as stated in Guerrette et al. (2023b). This approach can reduce the errors arising from various assumptions or approximations made during variable transformations but cannot avoid introducing imbalances between the prognostic fields and across the meshes due to the local nature of the observed quantity and the spatial localization applied to ensemble background error covariances, which can lead to potentially destructive effects on the quality of solutions from the compressible nonhydrostatic model solutions.

This article reports on our new implementation of the IAU feature in MPAS–JEDI using the cycling system called MPAS-Workflow (https://github.com/NCAR/MPAS-Workflow; last access: 2 August 2023). For mathematical completeness, the IAU module in the MPAS model code (first released in v5.0) has been updated with several corrections. The main goal of the IAU is to stabilize the cycling system and maintain numerical efficiency by reducing imbalances in the initial state. But it will increase the forecast lead times because the IAU forcing is incorporated into the model prognostic equations at every time step within a time window centered on the analysis time, followed by free forecasts. For 6 h cycling with a 6 h IAU time window, the total model integration time becomes 9 h, meaning the computation is increased by at least 30 %. As a technical paper, default namelist parameters for the IAU are listed, as defined in MPAS/namelist.atmosphere. This study does not discuss comprehensive characteristics of the complicated MPAS–JEDI or MPAS-Workflow systems but focuses on the implementation of the IAU, ensuring its reliability and performance through the forecast verification in surface pressure. Details of the implementation are described in Sect. 2, followed by the cycling system using MPAS-Workflow in Sect. 3. A real case study is presented in Sect. 4 and then concluded in Sect. 5.

2 The IAU implementation in MPAS–JEDI

The MPAS-A model uses a height-based coordinate following Klemp (2011), where terrain influences are progressively smoothed out toward the model top. The geometric height of coordinate surfaces (z) is defined by the combination of the nominal height (ignoring terrain) of coordinate surfaces (ζ) and terrain height (h_s) on a two-dimensional (x,z) model domain:

$$\begin{array}{}\text{(1)}& z=\mathit{\zeta }+A\left(\mathit{\zeta }\right)h_{\mathrm{s}}(\mathrm{x},\mathit{\zeta }),\end{array}$$

where $A\left(\mathit{\zeta }\right)=\mathrm{1}-\mathit{\zeta }/z_{\mathrm{t}}$, indicating a weight for the terrain influence on the coordinate surfaces, and z_t is the height of the model top. With $\mathrm{0}\le A\le \mathrm{1}-\mathit{\zeta }/z_{\mathrm{t}}$, an increasing amount of smoothing is applied at higher levels of ζ (e.g., decreasing A(ζ)) such that the vertical coordinate transitions from terrain-following at the surface (e.g., $z=h_{\mathrm{s}}(\mathrm{x},\mathit{\zeta })$ for ζ=0) toward constant-height surfaces aloft. While rectangular grids are usually represented as (x,y), the horizontal unstructured mesh is defined solely by individual grid cells (in a random order). For simplicity, we denote the horizontal grid cell dimension as x, in addition to the model level (z), to specify the two-dimensional model dimension (x,z).

In the MPAS model, the nonhydrostatic compressible equation is formulated for conserved quantities (mass, momentum, and moisture) represented in the flux form using the split-explicit time integration techniques introduced in Klemp et al. (2007). As described in Skamarock et al. (2012), defining a dry-air density adjusted by a terrain transformation ($\stackrel{\mathrm{̃}}{{\mathit{\rho }}_{\mathrm{d}}}={\mathit{\rho }}_{\mathrm{d}}/{\mathit{\zeta }}_z$, where ρ_d is dry-air density and $\mathrm{\nabla }\mathit{\zeta }=({\mathit{\zeta }}_x,{\mathit{\zeta }}_z)$), the flux variables X become

$$\begin{array}{}\text{(2)}& X=({\mathrm{\Theta }}_{\mathrm{m}},V_{\mathrm{e}},W,Q_j)=\stackrel{\mathrm{̃}}{{\mathit{\rho }}_{\mathrm{d}}}\cdot ({\mathit{\theta }}_{\mathrm{m}},u_{\mathrm{e}},w,q_j),\end{array}$$

where ${\mathit{\theta }}_{\mathrm{m}}=\mathit{\theta }(\mathrm{1}+(R_{\mathrm{v}}/R_{\mathrm{d}}\left)q_{\mathrm{v}}\right)$ is a potential temperature (θ) modified by water vapor mixing ratio (q_v). R_v and R_d are gas constants for water vapor and dry air, respectively. q_j represents the mixing ratio of each hydrometeor – water vapor (q_v), cloud liquid water (q_c), cloud ice (q_i), rain (q_r), snow (q_s), graupel (q_g), and hail (q_h). The horizontal momentum is predicted in terms of the wind speed normal to cell edges (u_e), while w stands for the vertical velocity. Coupled with $\stackrel{\mathrm{̃}}{{\mathit{\rho }}_{\mathrm{d}}}$, all the flux-form variables in the prognostic equations now include terrain effects through transformation of the vertical coordinate. Note that all the variables associated with the vertical coordinate surfaces (z, ζ, and ∇ζ) are constant fields.

The nonhydrostatic equations are integrated by updating total tendencies computed from each component of the modeling process. In the default configuration without the IAU, the model is simply integrated from the initial condition updated with analysis variables at the initial time. However, if the IAU is activated in the forecast (e.g., config_IAU_option = “on” in namelist.atmosphere), instead of changing the initial state, we compute the analysis increments by subtracting the background forecast from the analysis and divide them by the total number of time steps within the IAU window for a three-dimensional IAU (3DIAU). In this context, the background forecast (or first guess) is the forecast valid at the initial (or analysis) time from the previous cycle. While the analysis file is produced by MPAS–JEDI (and is employed as a new initial condition for the model run without the IAU), a separate analysis increment file at the analysis time (“AmB.nc”) is created for the IAU forcing in MPAS-Workflow and provided to the model through another data stream called “iau”. During the time integration, the total tendencies are adjusted by the IAU forcing (${\left.\frac{\partial X}{\partial t}\right|}_{\mathrm{amb}}$) at every time step as below.

$$\begin{array}{}\text{(3)}& {\left.{\displaystyle \frac{\partial X}{\partial t}}\right|}_{\mathrm{total}}={\left.{\displaystyle \frac{\partial X}{\partial t}}\right|}_{\mathrm{dyn}}+{\left.{\displaystyle \frac{\partial X}{\partial t}}\right|}_{\mathrm{phys}}+{\left.{\displaystyle \frac{\partial X}{\partial t}}\right|}_{\mathrm{amb}},\end{array}$$

where the first and the second term on the right-hand side represent the total tendencies from the dynamics and physics schemes, respectively.

In MPAS–JEDI, analysis variables are defined as temperature (T), specific humidity (s), surface pressure (P_s), and zonal and meridional wind components (u and v, respectively) at cell centers by default. Except for surface pressure, all of them are two-dimensional (2-D) variables in cells and levels (x,z). Note that in the unstructured mesh, 1-D indicates a horizontal plane, while 2-D includes the vertical dimension (similar to a traditional 3-D Cartesian coordinate). As defined in Eq. (2), the MPAS model predicts 2-D $\stackrel{\mathrm{̃}}{{\mathit{\rho }}_{\mathrm{d}}}$, θ_m, u_e, w, and hydrometeors in mixing ratios. In other words, none of the analysis variables are prognostic in the model, meaning that once their increments ($\mathit{\delta }v=v^{\mathrm{a}}-v^{\mathrm{b}}$, where “a” and “b” stand for the analysis and background for the variable v, respectively) are computed through minimization, they should be transformed to the prognostic variables for model integration. To reduce potential errors resulting from approximations such as hydrostatic balance and the equation of state, variable transformations are only applied to the increments (δv) rather than the full analysis fields (v^a), keeping prior states (v^b) as nonhydrostatic forecasts from the previous cycle. The edge wind speed (u_e) is updated by the increments in the horizontal wind components at cell centers (e.g., du and dv). Mass variable transformations begin with the increments in the approximated water vapor mixing ratio (δq_v), converted from the increments in specific humidity (δs) and its prior state (s). To linearize the equation $q_{\mathrm{v}}=s/(\mathrm{1}-s)$ for small increments δq_v and δs, the first derivative of q_v with respect to s is written as

$$\begin{array}{}\text{(4)}& {\displaystyle }{\displaystyle }{q}_{\mathrm{v}}^{\prime }\left(s\right)={\displaystyle \frac{\mathrm{d}}{\mathrm{d}s}}\left({\displaystyle \frac{s}{\mathrm{1}-s}}\right)={\displaystyle \frac{\mathrm{1}}{(\mathrm{1}-s{)}^{\mathrm{2}}}}.\end{array}$$

Using the first-order Taylor series expansion, we can linearize q_v for the increments δq_v and δs as

$$\begin{array}{}\text{(5)}& q_{\mathrm{v}}+\mathit{\delta }{q}_{\mathrm{v}}\approx q_{\mathrm{v}}\left(s\right)+q_{\mathrm{v}}^{\prime }\left(s\right)\mathit{\delta }s\approx {\displaystyle \frac{s}{\mathrm{1}-s}}+{\displaystyle \frac{\mathrm{1}}{(\mathrm{1}-s{)}^{\mathrm{2}}}}\mathit{\delta }s.\end{array}$$

Then the increments in the pressure field are derived based on the changes in virtual temperature ($T_{\mathrm{v}}=T(\mathrm{1}+\mathrm{0.608}{q}_{\mathrm{v}})$). Assuming hydrostatic balance, the two-dimensional pressure field (P(x,z)) is integrated upward from the surface based on the analyzed surface pressure ($P_{\mathrm{s}}^{\mathrm{a}}=P_{\mathrm{s}}^{\mathrm{b}}+\mathit{\delta }{P}_{\mathrm{s}}$). Dry-air density (ρ_d) is updated based on the approximated equation of state (${\mathit{\rho }}_{\mathrm{d}}=P/\left[R_{\mathrm{d}}{T}_{\mathrm{v}}\right(\mathrm{1}+q_{\mathrm{v}}\left)\right]$), and potential temperature is derived from the relationship $\mathit{\theta }=T(\frac{P_{\mathrm{0}}}{P}{)}^{R_{\mathrm{d}}/C_p}$, where P₀=1000 hPa is a reference pressure and C_p is the specific heat at constant pressure. Details are summarized in Appendix B.

While the nonhydrostatic MPAS-A model solves a prognostic equation for ${\mathrm{\Theta }}_{\mathrm{m}}=\stackrel{\mathrm{̃}}{{\mathit{\rho }}_{\mathrm{d}}}\cdot {\mathit{\theta }}_{\mathrm{m}}$ and diagnoses the full pressure ($p=p_{\mathrm{0}}(R_{\mathrm{d}}{\mathit{\zeta }}_z{\mathrm{\Theta }}_{\mathrm{m}}/p_{\mathrm{0}}{)}^{\mathit{\gamma }}$, where $\mathit{\gamma }=C_p/C_{\mathrm{v}}$), the analysis updates in MPAS–JEDI primarily rely on surface pressure (in both the analysis and background) and several approximations such as hydrostatic balance and the simplified moist conversion (e.g., ignoring all hydrometeors except for q_v).

As the dry density multiplied by the vertical coordinate Jacobian ($\stackrel{\mathrm{̃}}{{\mathit{\rho }}_{\mathrm{d}}}$) is also updated during the analysis, once the prognostic variables are updated through the variable transformations, they should be recoupled with the updated $\stackrel{\mathrm{̃}}{{\mathit{\rho }}_{\mathrm{d}}}$, as on the right-hand side of Eq. (2), to compute the tendencies for the IAU forcing (${\left.\frac{\partial X}{\partial t}\right|}_{\mathrm{amb}}$). When the model is configured without the IAU, the recoupling step is carried out as part of the MPAS initialization process, but in the case of the IAU, this is taken care of within the IAU module in the MPAS model. The tendency for θ_m in the flux form, for instance, can be expressed in the following partial equation:

$$\begin{array}{}\text{(6)}& {\displaystyle \frac{\partial X}{\partial t}}={\displaystyle \frac{\partial }{\partial t}}(\stackrel{\mathrm{̃}}{{\mathit{\rho }}_{\mathrm{d}}}\cdot {\mathit{\theta }}_{\mathrm{m}})={\displaystyle \frac{\partial }{\partial t}}(\stackrel{\mathrm{̃}}{{\mathit{\rho }}_{\mathrm{d}}}\cdot \mathit{\theta }(\mathrm{1}+{\displaystyle \frac{R_{\mathrm{v}}}{R_{\mathrm{d}}}}{q}_{\mathrm{v}}\left)\right).\end{array}$$

Based on the differential rule, Eq. (6) can be rewritten as

$$\begin{array}{ll}{\displaystyle \frac{\partial }{\partial t}}\left(\stackrel{\mathrm{̃}}{{\mathit{\rho }}_{\mathrm{d}}}\cdot \mathit{\theta }\left(\mathrm{1}+{\displaystyle \frac{R_{\mathrm{v}}}{R_{\mathrm{d}}}}{q}_{\mathrm{v}}\right)\right)& {\displaystyle }=\mathit{\theta }\left(\mathrm{1}+{\displaystyle \frac{R_{\mathrm{v}}}{R_{\mathrm{d}}}}{q}_{\mathrm{v}}\right){\displaystyle \frac{\partial \stackrel{\mathrm{̃}}{{\mathit{\rho }}_{\mathrm{d}}}}{\partial t}}\\ {\displaystyle }& {\displaystyle }+\stackrel{\mathrm{̃}}{{\mathit{\rho }}_{\mathrm{d}}}\left(\mathrm{1}+{\displaystyle \frac{R_{\mathrm{v}}}{R_{\mathrm{d}}}}{q}_{\mathrm{v}}\right){\displaystyle \frac{\partial \mathit{\theta }}{\partial t}}\\ \text{(7)}& {\displaystyle }& {\displaystyle }+\stackrel{\mathrm{̃}}{{\mathit{\rho }}_{\mathrm{d}}}\mathit{\theta }{\displaystyle \frac{\partial }{\partial t}}\left(\mathrm{1}+{\displaystyle \frac{R_{\mathrm{v}}}{R_{\mathrm{d}}}}{q}_{\mathrm{v}}\right){\displaystyle }& {\displaystyle }=\left(\mathrm{1}+{\displaystyle \frac{R_{\mathrm{v}}}{R_{\mathrm{d}}}}{q}_{\mathrm{v}}\right)\left(\mathit{\theta }{\displaystyle \frac{\partial \stackrel{\mathrm{̃}}{{\mathit{\rho }}_{\mathrm{d}}}}{\partial t}}+\stackrel{\mathrm{̃}}{{\mathit{\rho }}_{\mathrm{d}}}{\displaystyle \frac{\partial \mathit{\theta }}{\partial t}}\right)\\ \text{(8)}& {\displaystyle }& {\displaystyle }+\stackrel{\mathrm{̃}}{{\mathit{\rho }}_{\mathrm{d}}}\mathit{\theta }{\displaystyle \frac{R_{\mathrm{v}}}{R_{\mathrm{d}}}}{\displaystyle \frac{\partial q_{\mathrm{v}}}{\partial t}}{\displaystyle }& {\displaystyle }=\left(\mathrm{1}+{\displaystyle \frac{R_{\mathrm{v}}}{R_{\mathrm{d}}}}{q}_{\mathrm{v}}\right){\displaystyle \frac{\partial \left(\stackrel{\mathrm{̃}}{{\mathit{\rho }}_{\mathrm{d}}}\mathit{\theta }\right)}{\partial t}}\\ \text{(9)}& {\displaystyle }& {\displaystyle }+\stackrel{\mathrm{̃}}{{\mathit{\rho }}_{\mathrm{d}}}\mathit{\theta }{\displaystyle \frac{R_{\mathrm{v}}}{R_{\mathrm{d}}}}{\displaystyle \frac{\partial q_{\mathrm{v}}}{\partial t}}.\end{array}$$

The rightmost term in Eq. (9) contains $\stackrel{\mathrm{̃}}{{\mathit{\rho }}_{\mathrm{d}}}\frac{\partial q_{\mathrm{v}}}{\partial t}$, which can be derived from the tendency equation for the water vapor mixing ratio, as below.

$$\begin{array}{}\text{(10)}& {\displaystyle \frac{\partial }{\partial t}}\left(\stackrel{\mathrm{̃}}{{\mathit{\rho }}_{\mathrm{d}}}{q}_{\mathrm{v}}\right)=\stackrel{\mathrm{̃}}{{\mathit{\rho }}_{\mathrm{d}}}{\displaystyle \frac{\partial q_{\mathrm{v}}}{\partial t}}+q_{\mathrm{v}}{\displaystyle \frac{\partial \stackrel{\mathrm{̃}}{{\mathit{\rho }}_{\mathrm{d}}}}{\partial t}}\text{(11)}& {\displaystyle }\stackrel{\mathrm{̃}}{{\mathit{\rho }}_{\mathrm{d}}}{\displaystyle \frac{\partial q_{\mathrm{v}}}{\partial t}}={\displaystyle \frac{\partial \left(\stackrel{\mathrm{̃}}{{\mathit{\rho }}_{\mathrm{d}}}{q}_{\mathrm{v}}\right)}{\partial t}}-q_{\mathrm{v}}{\displaystyle \frac{\partial \stackrel{\mathrm{̃}}{{\mathit{\rho }}_{\mathrm{d}}}}{\partial t}}\end{array}$$

By incorporating Eq. (11) into Eq. (9), Eq. (6) can be rewritten as

$$\begin{array}{}\text{(12)}& \begin{array}{rl}{\displaystyle \frac{\partial }{\partial t}}(\stackrel{\mathrm{̃}}{{\mathit{\rho }}_{\mathrm{d}}}\cdot {\mathit{\theta }}_{\mathrm{m}})& ={\displaystyle \frac{\partial }{\partial t}}\left(\stackrel{\mathrm{̃}}{{\mathit{\rho }}_{\mathrm{d}}}\cdot \mathit{\theta }\left(\mathrm{1}+{\displaystyle \frac{R_{\mathrm{v}}}{R_{\mathrm{d}}}}{q}_{\mathrm{v}}\right)\right)\\ & =\left(\mathrm{1}+{\displaystyle \frac{R_{\mathrm{v}}}{R_{\mathrm{d}}}}{q}_{\mathrm{v}}\right){\displaystyle \frac{\partial \left(\stackrel{\mathrm{̃}}{{\mathit{\rho }}_{\mathrm{d}}}\mathit{\theta }\right)}{\partial t}}\\ & +\mathit{\theta }{\displaystyle \frac{R_{\mathrm{v}}}{R_{\mathrm{d}}}}\left({\displaystyle \frac{\partial \left(\stackrel{\mathrm{̃}}{{\mathit{\rho }}_{\mathrm{d}}}{q}_{\mathrm{v}}\right)}{\partial t}}-q_{\mathrm{v}}{\displaystyle \frac{\partial \stackrel{\mathrm{̃}}{{\mathit{\rho }}_{\mathrm{d}}}}{\partial t}}\right).\end{array}\end{array}$$

In the original IAU implementation in MPAS-A (version 7.1 and prior), the last term ($-\mathit{\theta }\frac{R_{\mathrm{v}}}{R_{\mathrm{d}}}{q}_{\mathrm{v}}\frac{\partial \stackrel{\mathrm{̃}}{{\mathit{\rho }}_{\mathrm{d}}}}{\partial t}$) was omitted but is now added back in the latest release of MPAS–JEDI. The missing term was a mathematical error which needed to be fixed for accuracy. Our experiment on the coarse mesh with a nominal resolution of 120 km (as described in the following section) is not suitable for quantifying the error magnitude caused by this term. But there might be cases where the absence of this term can introduce sizable forecast errors in both thermodynamic and dynamic tendencies, especially in high-resolution simulations of strong pressure gradients.

The tendency for edge wind (u_e) can be described in the same manner, with the same correction applied.

$$\begin{array}{}\text{(13)}& {\displaystyle \frac{\partial }{\partial t}}(\stackrel{\mathrm{̃}}{{\mathit{\rho }}_{\mathrm{d}}}\cdot u_{\mathrm{e}})=\stackrel{\mathrm{̃}}{{\mathit{\rho }}_{\mathrm{d}}}{\displaystyle \frac{\partial u_{\mathrm{e}}}{\partial t}}+u_{\mathrm{e}}{\displaystyle \frac{\partial \stackrel{\mathrm{̃}}{{\mathit{\rho }}_{\mathrm{d}}}}{\partial t}}\end{array}$$

Here u_e is coupled with dry-air density ($\stackrel{\mathrm{̃}}{{\mathit{\rho }}_{\mathrm{d}}}$) at the edges, rather than the cell centers. In the horizontally unstructured MPAS mesh, the center of each grid cell (mostly in a hexagonal shape) serves as the center of mass (centroidal), and cell edges bisect the lines connecting the two cell centers sharing the edges. Thus, the corresponding dry-air density at the edge is defined as the mean of $\stackrel{\mathrm{̃}}{{\mathit{\rho }}_{\mathrm{d}}}$ at the two cell centers that share the edge. (For the detailed description on the mesh characteristics, users should refer to Fig. 1 in Skamarock et al., 2012.) Note that all the model's prognostic variables except the edge wind (u_e) are defined at the cell centers, coupled with $\stackrel{\mathrm{̃}}{{\mathit{\rho }}_{\mathrm{d}}}$ at the cell centers. While w is defined at the cell center of the horizontal mesh, it is C grid staggered in the vertical.

In the IAU module, analysis increments (Δx) in the so-called “uncoupled” variables are read from the iau stream

$$\begin{array}{}\text{(14)}& \mathrm{\Delta }x=x^{\mathrm{a}}-x^{\mathrm{b}}\text{, where }x=({\mathit{\rho }}_{\mathrm{d}},\mathit{\theta },q_{\mathrm{v}}),\end{array}$$

so that the IAU tendencies (${\left.\frac{\partial X}{\partial t}\right|}_{\mathrm{amb}}$) are computed for the coupled variables (X) as in Eqs. (12) and (13). They are then multiplied by a predefined weighting function (ω_k) to be applied every time step (Δt; config_dt) over the IAU time window (Δτ; config_IAU_window_length_s).

$$\begin{array}{}\text{(15)}& {\displaystyle }{\displaystyle }\underset{-\frac{\mathrm{\Delta }\mathit{\tau }}{\mathrm{2}}}{\overset{\frac{\mathrm{\Delta }\mathit{\tau }}{\mathrm{2}}}{\int }}{\left.{\displaystyle \frac{\partial X}{\partial t}}\right|}_{\mathrm{amb}}\mathrm{d}t=\sum _{k=\mathrm{1}}^n{\mathit{\omega }}_k\cdot {\displaystyle \frac{\mathrm{\Delta }X}{\mathrm{\Delta }t}}=\sum _{k=\mathrm{1}}^n{\mathit{\omega }}_k\cdot {\displaystyle \frac{\mathrm{\Delta }(\stackrel{\mathrm{̃}}{{\mathit{\rho }}_{\mathrm{d}}}\cdot x)}{\mathrm{\Delta }t}},\end{array}$$

where $\mathrm{\Delta }X=X^{\mathrm{a}}-X^{\mathrm{b}}$, $n=\mathrm{\Delta }\mathit{\tau }/\mathrm{\Delta }t$, and ${\mathit{\omega }}_k=\mathrm{1}/n$ for each time step k. In the default configuration for 3DIAU, the IAU forcing computed from analysis increments is evenly distributed across time steps for the 6 h IAU window, followed by a 3 h free forecast to the next analysis time, making a 9 h background forecast at each analysis cycle (Fig. 1). Although a simple 3DIAU is currently implemented with constant forcing, it could be extended for a four-dimensional IAU (4DIAU) with varying weights over the IAU time window.

Figure 1 A diagram for cycling with the IAU. Δτ is the IAU time window (config_IAU_window_length_s), and Δt is an integration time step (config_dt) defined in namelist.atmosphere.

3 Cycling data assimilation with MPAS-Workflow

Our open-source MPAS-Workflow was introduced in Guerrette et al. (2023b) and has only been tested on Cheyenne, one of the National Center for Atmospheric Research (NCAR) high-performance computers (HPCs), but we describe some technical details here since it is at the heart of the cycling system for MPAS–JEDI and it has gone through major updates for the IAU feature. MPAS-Workflow uses the Cylc general purpose workflow engine (v7.8.3) (Oliver et al., 2019, https://cylc.github.io/; last access: 29 August 2023), and it is designed for end-to-end processes of the MPAS–JEDI cycling system. It controls all the parameters for running the MPAS model and data assimilation, with high flexibility for a number of different configurations. It defines various lists of variables such as those for the analysis, background, and observation types to assimilate. Furthermore, it is equipped with a Python-based post-processing package (https://github.com/JCSDA/mpas-jedi/tree/2.0.0/graphics; last access: 2 August 2023), including diagnostics and plotting utilities.

Once the IAU is activated in MPAS-Workflow, the model initial and run times are automatically adjusted to −3 and 9 h, respectively, for a 6 h IAU time window, as depicted in Fig. 1. Also, the analysis increment file (AmB.nc) is created at the analysis time (e.g., t=0), while the background forecast valid at −3 h (instead of 0 h) is employed from the previous cycle. Due to the availability of the first-guess file, the IAU option is activated only from the second cycle, and the model output interval is changed from 6 to 3 h for 6 h cycling.

MPAS–JEDI employs its own customized version of the MPAS-A model, using a two-stream I/O (input/output) approach by default to run DA cycling more efficiently. The two-stream (or split) I/O approach was originally developed for DA cycling in a restart mode to avoid writing time-invariant fields in every restart file while ensuring the model forecasts are reproducible. In the restart mode, the MPAS-A model produces a restart file with about 230 variables, among which only ∼70 variables vary with time, while the rest of them are time-invariant. Note that these static fields are not the same as those in the static file for the MPAS initialization because the MPAS static file only contains horizontal fields, with no vertical dimension or coordinate. This is because the MPAS initialization consists of two consecutive steps – the first step for constructing the 1-D horizontal mesh in the static file and the second for producing 2-D fields based on the terrain-following height vertical coordinate in the initial-condition file. By splitting the restart file for the variables between time-variant and time-invariant, we can cycle with a much smaller file containing time-variant fields only, which is approximately $\mathrm{1}/\mathrm{6}$ of the original size of the restart file. We refer to the new I/O stream as “da_state” in the MPAS-A model, as the file serves as both input and output for DA cycling and can be updated by the analysis process. In fact, this I/O stream can be used in a restart mode regardless of DA or DA cycling. One caveat of the split I/O approach is that the reproducibility might depend on the model configuration (and potentially the version), meaning that the variable list of da_state is not always applicable to or guaranteed for all different namelist options for MPAS-A. We had initially developed this new I/O stream in the MPAS-A model based on version 6.1 and ensured its bit-for-bit reproducibility. However, the MPAS–JEDI cycling system is run in a cold-start mode, initializing all the physics tendencies at the analysis time, so it only uses the da_state stream as a shorter version of the model input and output files, not to replace a restart file. At the time of writing, this new I/O stream is available only in MPAS–JEDI, but it will be merged into the official version of the MPAS-A model (https://github.com/MPAS-Dev/MPAS-Model, last access: 21 May 2024) in the future.

The variational approach essentially linearizes the model and constructs a static background (or forecast) error covariance to find an analysis solution closest to observations iteratively (e.g., through a minimization process). Although the static background error covariance only represents the climatological information (e.g., with no temporal variations), it is a key component of variational data assimilation algorithms, modeling the relationships between control variables through physical transformation or balance operators as well as spatial auto-correlations of each control variable to determine how to propagate the observed information across model grids and variables (Descombes et al., 2015). In this study, we use a pure ensemble-variational (EnVar) approach with zero static error covariance. For the minimization process, we employ an incremental approach (e.g., minimizing the cost function for increments) (Courtier et al., 1994). Following Liu et al. (2022), the ensemble background error covariance composed of 20-member 6 h MPAS forecasts initialized from the National Centers for Environmental Prediction (NCEP) 20-member Global Ensemble Forecast System (GEFS) ensemble analysis at each cycle. Due to the small ensemble size, we also apply the distance-based correlation function by Gaspari and Cohn (1999) using 1200 and 6 km as full-width radii for horizontal and vertical covariance localization, respectively.

In the JEDI system, the Unified Forward Operator (UFO) not only provides observation operators (named “HofX”) to compute innovations (e.g., differences between observations and the corresponding forecasts, O − F) for all different observation types but also handles all the data quality control (QC) filters, the data range or coverage area, data manipulation (such as data thinning or averaging onto model levels), bias correction, and the specification of observation errors (including observation uncertainties and error correlations, if applicable). Due to various filters (or QCs) that can be applied to observations, the UFO produces observation errors before and after the QCs, named “ObsError” and “EffectiveErrors”, respectively. While ObsError indicates the initial observation error values from the input observation file, the actual observation errors applied to data assimilation can be found in EffectiveErrors. When multiple filters are applied to observations or observation operators (such as data thinning or variable transformation), users can specify the order of filters through the YAML configuration files (YAML is a human-readable data format; https://yaml.org, last access: 27 December 2023). The default quality control (called “PreQC”) rejects observations when the data QC flag is larger than 3, as provided by input observation files (such as GSI-ncdiag files used in this study). In the assimilation of surface observations, for example, they are also discarded if innovations exceed 3 times the standard deviation of the observation error (σ_o) or if the station height differs from the model's surface elevation by more than 200 m. In this study, the height difference threshold is reduced to 100 m and the surface pressure terrain height correction (called “SfcPCorrected”) uses the so-called Weather Research and Forecasting Data Assimilation (WRFDA) method rather than the default UK Met Office (UKMO) method (Ingleby, 2013). In our initial test (not shown), the difference between different height correction methods was insignificant, but it would be worth revisiting in the future, especially for regions with significant orography. While the common functions or modules are located under UFO in JEDI, most of these options are controlled through YAML configurations in MPAS-Workflow. An example for surface data assimilation is provided in Appendix A.

4 Experiments

After the new implementation of the IAU in MPAS v7, global analysis and forecast cycling was conducted over a global 120 km quasi-uniform mesh every 6 h for 1 month, starting from 15 April 2018, using MPAS-Workflow for the hybrid 3DEnVar in the MPAS–JEDI system. During the cycling, all the conventional observations, satellite winds, and clear-sky microwave radiances from six Advanced Microwave Sounding Unit-A (AMSU-A) sensors aboard NOAA-15, NOAA-18, NOAA-19, Aqua, Metop-A, and Metop-B were assimilated together, using diagonal observation error covariances and a pure ensemble background error covariance (computed from GEFS), like in Guerrette et al. (2023b).

In the model simulation, a “mesoscale_reference” physics suite is used that includes WSM6 (WRF single-moment 6-class) microphysics, new Tiedtke cumulus, YSU PBL (Yonsei University planetary boundary layer), YSU gravity wave drag over orography (GWDO), RRTMG SW and LW (Rapid Radiative Transfer Model for general circulation models shortwave and longwave), and Noah LSM (land surface model) variables. Ozone climatology is activated, and radiation effective radii for cloud water (q_c), cloud ice (q_i), and snow (q_s) are computed in the microphysics scheme (e.g., config_microp_re = true).

Figure 2 Time series of the globally averaged absolute value of the surface pressure tendency ($|\mathrm{d}{P}_{\mathrm{s}}/\mathrm{d}t|$) in the forecast from the analysis at 00:00 UTC on 1 May 2018 over the 120 km quasi-uniform mesh.

The 2-month-long cycling experiments are conducted with and without the IAU, named IAU and CTRL (control), respectively. Figure 2 illustrates a comparison of the absolute value of the surface pressure tendency ($\left|\frac{\mathrm{d}{P}_{\mathrm{s}}}{\mathrm{d}t}\right|$) as an area-weighted global mean, in the background forecasts from the analysis (shown at 0 on the x axis), valid at 00:00 UTC on 1 May 2018. Even after a 2-week spinup period, it shows that the initial noise arising from the analysis increments is very high in the control run (“CTRL” as a dashed line), but the IAU effectively suppresses such noise throughout the 9 h forecast. The asymptotic value is ∼0.013 Pa s⁻¹ for the 9 h forecast, which is comparable to ∼0.01 Pa s⁻¹ for the 6 h forecast, as presented in Lynch and Huang (1992) in their digital filter initialization (DFI) study. In Fig. 3, the horizontal distribution reveals that noise from the initial state in CTRL is widespread across the globe from the first time step, whereas it almost disappears with IAU. It is noted that, in the cycling with the IAU, the forecast starts from −3 h, leading to a different initial time for model integration compared to the control run (which starts at the analysis time, i.e., 0 h). However, regardless of the actual initial time, our focus here is on comparing the deviation of the first time step from the initial state.

Figure 3 Horizontal distribution of $\mathrm{d}{P}_{\mathrm{s}}/\mathrm{d}t$ simulated in the first time step from the initial conditions in (a) the control run without the IAU and (b) the IAU run.

Figure 4 Time series of the difference in the total number of surface pressure observations assimilated in CTRL and IAU (red) and the total number of the observations available at each cycle (gray).

In DA cycling, it is common to monitor the total number of observations assimilated at each cycle. A time series with a declining trend might be indicative of the analysis with poor quality, rejecting more observations with cycles. As shown in Fig. 4, positive differences (IAU − CTRL; red) indicate that slightly more observations are assimilated with the IAU. The difference is small compared to the total number of observations (as a gray line with the right y axis), but the IAU tends to assimilate more observations in most cycles, another sign of our successful implementation.

Figure 5 Time series of (a) root mean square (rms) and (b) mean bias errors in surface pressure (Pa) (O − B) in both the CTRL and IAU runs during the cycles.

To examine the performance of cycling DA, Fig. 5 compares observation minus the background (O − B) in (a) root mean square (rms) errors and (b) mean errors over the globe for each cycle. Note that O − B or (O − F) is computed for 6 h background forecasts in CTRL but for 9 h forecasts in the IAU run; both are at the same validation time. The rms errors in IAU (dots on the red line) are slightly yet consistently smaller than those in CTRL (+ symbols on the black line) throughout the cycling period. In terms of mean errors, however, the month-long average in surface pressure is consistent at 20 Pa for both experiments. The time series illustrates that both rms and bias errors initially start with large error magnitudes, but after about 1 week of cycling, the errors tend to stabilize. While the global-mean rms error in CTRL, averaged over all the cycles following the 1-week spinup period, is ∼1.2 hPa, IAU leads to an improvement in the forecast error by approximately 3 % (e.g., (CTRL − IAU) $/$ CTRL $\times \mathrm{100}=-\mathrm{3}$ %).

Figure 6 Percent difference of rms (O − B) in IAU with respect to the one in CTRL for radiosonde verification in (a) zonal wind (m s⁻¹) and (b) temperature (K).

In the sounding verification over the globe, the percentage difference of rms errors in IAU with respect to the one in CTRL for the entire cycling period also shows slight but systematic improvements in the background forecasts, as depicted in Fig. 6. The impact on zonal wind is almost neutral, but temperature forecasts near the surface (e.g., 1000 hPa) are improved by 3.6 %, similar to the improvement in surface pressure. The moisture verification for the 6 h forecast is provided in the Supplement.

Figure 7 The rms errors in 5 d forecasts in surface pressure at different latitudes relative to the GFS analysis are displayed (a) for CTRL and (b) as the rms differences from CTRL in percentages in IAU with statistical significance at the 95 % confidence level indicated by light circles.

We also run 5 d forecasts from the 00:00 UTC analysis every day and compute forecast errors with respect to the Global Forecast System (GFS) analysis. Figure 7 displays rms errors in CTRL for 5 d forecasts on the x axis across latitudes (on the y axis) on the right, while the relative differences in rms errors in IAU (%) are depicted on the right (b), with red indicating degradation and blue implying forecast improvements over the baseline (CTRL). Statistical significance is denoted by light circles at the 95 % confidence level. Compared to the GFS analysis, MPAS forecasts in CTRL exhibit the largest (or the fastest) error growth in the Southern Hemisphere. Forecasts in the IAU run, on the other hand, tend to reduce errors in the tropics while increasing errors near the North Pole region. This aligns with the findings of Ha et al. (2017), where forecast errors were significantly reduced over the tropics in a variable-resolution mesh, including both resolution-transition and high-resolution parts. Because the IAU is implemented on the model's unstructured mesh (which is in a random order), it is not associated with particular geographic locations or mesh configurations. Given its time-filtering feature, the IAU might be more effective in simulating low-frequency modes dominant over the tropics. It is also noted that the impact of the IAU may be nonlinearly intertwined with model errors in data-sparse regions, such as the poles. However, model errors are not accounted for in the hybrid 3DEnVar system used in this study. Additional area-specific features in the verification are provided in the Supplement. As the GFS analysis also suffers from its own errors (Bhargava et al., 2018), the forecast verification against the analysis is not intended for a thorough investigation of the model performance. It rather serves to introduce post-processing capabilities in both model and observation spaces provided by MPAS-Workflow.

The overall results from the cycling experiments are promising, showing the reliability throughout the month-long period. We can examine the impact of the IAU option through more extensive diagnostics against other observation types and variables and through the evaluation of longer forecasts in the future. As a proof of concept, only 3DIAU on a 120 km global mesh was tested here, but it was implemented in a way to make it extensible with different weighting functions or even to 4DIAU in the MPAS model. Also, it is applicable to variable-resolution meshes in case one would want to examine the impact of the IAU over the area with mesh refinement.

5 Conclusions

This study introduces the incremental analysis update (IAU) implemented in the MPAS–JEDI cycling system operated by MPAS-Workflow. Through a real case study for 1 month, starting from mid-April 2018, assimilating all conventional, aircraft, and satellite radiance observations, we demonstrate that the IAU is successfully implemented on the model's unstructured mesh, effectively suppressing the artificial noise produced by initial imbalances during the analysis process. Although the current implementation is a simple three-dimensional IAU (3DIAU) with the same fractional forcing applied to each time step, there are several aspects that might be worth pointing out in regards to our development effort in MPAS–JEDI. (i) Computational stability and efficiency might be critical to any numerical weather prediction (NWP) models, but special attention was taken regarding the MPAS-A model which solves the compressible nonhydrostatic equations employing an unstructured mesh based on centroidal Voronoi tessellations. In the model integration, a time step is set based on the smallest grid spacing of a given unstructured mesh and is then uniformly applied across the entire mesh (e.g., regardless of the nominal grid spacing of individual grid cells). To ensure numerical stability even for the unstructured mesh applications, various filtering techniques are carefully designed and applied to the model numerics (Klemp et al., 2007, 2018). Aside from the modeling approach, the IAU is considered another efficient way to control high-frequency oscillations produced by the analysis procedure so that energy does not accumulate in the acoustic or unbalanced gravity wave modes through cycles due to the initialization. (ii) The MPAS-A model treats prognostic variables in the flux form, meaning that the variables are coupled with dry-air density. Because the density is also updated as part of the analysis variables, analysis increments in the IAU forcing term should also be computed in the tendency form of each analysis variable recoupled with the updated density. During our implementation, the IAU module is corrected to properly represent changes in the density from the analysis. (iii) Analysis increments in MPAS–JEDI are basically computed from the linearized version of the hydrostatic balance equation that is vertically integrated from the analyzed surface pressure. Horizontally, MPAS–JEDI updates all the analysis variables in the model's native (e.g., unstructured) mesh, but vertically the MPAS-A height coordinate is mapped to the hydrostatic pressure coordinate for the increments based on the surface pressure analysis. A global 120 km mesh employed for our cycling experiments demonstrates our technical implementation of the IAU capability in MPAS–JEDI. As we move towards convective-scale data assimilation with unconventional observations such as radars or lidars, however, it would be worth revisiting the influence of the IAU on either variable- or higher-resolution mesh, in conjunction with various DA techniques available in MPAS–JEDI.

Appendix A: Cycling with MPAS-Workflow

MPAS–JEDI is an interface between the MPAS-A model and the JEDI data assimilation system, including all the model-specific components such as variable transformation and HofX, which computes analysis increments in the MPAS variables. To cycle MPAS–JEDI over a period of time in a synthetic way, MPAS-Workflow controls the entire data stream as well as all the configurations for data assimilation and model forecasts. As the IAU changes the input/output file stream and the model configuration, they are all accounted for in MPAS-Workflow (mostly through various YAML configurations).

MPAS-Workflow offers high flexibility for a number of applications such as data assimilation with 3DVar, pure 3/4DEnVar, hybrid 3DEnVar with dual resolution, ensemble data assimilation, and more recently the local ensemble transform Kalman filter (LETKF) algorithm in MPAS–JEDI. In addition, it can generate observations in the IODA format, the GFS analyses in MPAS initial-condition format, and free deterministic forecasts from the GFS analyses using specific Cylc suites (no cycling). Observations and GFS analyses can be obtained from the NCAR Research Data Archive (RDA) or the NCEP FTP server. With all of these capabilities, it has been tested for near-real-time cycling runs using the 3DVar algorithm.

At the time of writing, MPAS-Workflow does not build either MPAS or JEDI, which should be built separately after downloading the source codes from https://github.com/JCSDA/mpas-bundle/tree/release/2.0.0 (last access: 2 August 2023); mpas-bundle is built using cmake, a set of CMake macros provided by the European Centre for Medium-Range Weather Forecasts (ECMWF), along with their libraries. For the installation guide for those tools, readers can refer to https://jointcenterforsatellitedataassimilation-jedi-docs.readthedocs-hosted.com/en/latest/inside/developer_tools/cmake.html (last access: 2 August 2023). Once built, the path should be specified in initialize/framework/Build.py in MPAS-Workflow.

To run a cycling experiment with the IAU, all users need to do is edit a single YAML file. Each section controls the specific configuration to build up the YAML file of the MPAS–JEDI application (i.e., pure 3DEnVar) and other components of the workflow to construct the Cylc suite that will orchestrate the cycling experiment. For the IAU, a new logical parameter has to be added in a new line as “IAU: True”. Here is the configuration employed this study.

<MPAS-Workflow/scenarios/3denvar_OIE120km_ WarmStart_IAU.yaml>

workflow:

first cycle point: 20180414T18

final cycle point: 20180510T00

experiment:

suffix: “_IAU”

observations:

resource: PANDACArchive

members:

n: 1

model:

outerMesh: 120 km

innerMesh: 120 km

ensembleMesh: 120 km

firstbackground:

resource: “PANDAC.GFS”

externalanalyses:

resource: “GFS.PANDAC”

variational:

DAType: 3denvar

ensemble:

forecasts:

resource: “PANDAC.GEFS”

forecast:

IAU: True

The YAML file used to run MPAS–JEDI with 3DEnVar is built by adding the observers snippets to the application sections (see MPAS-Workflow/config/jedi/applications/3denvar.yaml). Here we also provide two more sample YAML configurations used in this study for assimilating surface observations.

• (1)

  MPAS-Workflow/config/jedi/ObsPlugs/da/filters/sfc.yaml (The following is for the data QC flag.)

  obs filters:

  • –

    filter: PreQC

    maxvalue: 3 # Maximum data QC flag

  • –

    filter: Difference Check

    reference: MetaData/stationElevation value: GeoVaLs/surface_altitude

    threshold: 100.0 # default: 200.

  • –

    filter: Background Check

    threshold: 3.0

• (2)

  MPAS-Workflow/config/jedi/ObsPlugs/da/base/sfc.yaml (The following are options for an observation operator for surface observations.)

  • –

    obs space:

    name: SfcPCorrected

    _obsdatain: &ObsDataIn

    engine:

    type: H5File

    obsfile: InDBDir/sfc_obs_thisValidDate.h5

    _obsdataout: &ObsDataOut

    engine:

    type: H5File

    obsfile: OutDBDirMemberDir/obsPrefix_sfcObsOutSuffix.h5

    obsdatain: *ObsDataIn

    ObsDataOut

    simulated variables: [stationPressure]

    obs error: *ObsErrorDiagonal

    obs operator:

    name: SfcPCorrected

    da_psfc_scheme: WRFDA # default: UKMO

    linear obs operator:

    name: Identity

    observation alias file: obsop_name_map.yaml

Appendix B: Linearized equations for incremental variable transformations

As described in Sect. 2, MPAS–JEDI first computes the analysis through the iterative minimization procedure and then converts the increments in the analysis variables to the model's prognostic fields. Based on Eq. (5), the increments in the water vapor mixing ratio are computed as $\mathit{\delta }{q}_{\mathrm{v}}^k=\mathit{\delta }{s}^k/(\mathrm{1}-s^k{)}^{\mathrm{2}}$ at each model level k. Then, the increments in virtual temperature ($\mathit{\delta }{T}_{\mathrm{v}}^k$) and pressure (δP^k) are derived using the first derivative as follows:

$$\begin{array}{}\text{(B1)}& {\displaystyle }\mathit{\delta }{T}_{\mathrm{v}}^k=\mathit{\delta }{T}^k(\mathrm{1}+\mathrm{0.608}{q}_{\mathrm{v}}^k)+\mathrm{0.608}\cdot \mathit{\delta }{q}_{\mathrm{v}}^k\cdot T^k,\text{(B2)}& {\displaystyle }\begin{array}{rl}\mathit{\delta }{P}^k& =\mathit{\delta }{P}^{k-\mathrm{1}}\exp \mathit{\{}-{\displaystyle \frac{g(z^k-z^{k-\mathrm{1}})}{R_{\mathrm{d}}{T}_{\mathrm{v}}^k}}\mathit{\}}\\ & +P^{k-\mathrm{1}}\exp \mathit{\{}-{\displaystyle \frac{g(z^k-z^{k-\mathrm{1}})}{R_{\mathrm{d}}{T}_{\mathrm{v}}^k}}\mathit{\}}{\displaystyle \frac{g(z^k-z^{k-\mathrm{1}})}{R_{\mathrm{d}}{T_{\mathrm{v}}^k}^{\mathrm{2}}}}\mathit{\delta }{T}_{\mathrm{v}}^k,\end{array}\end{array}$$

where R_d is the gas constant for dry air. With P⁰=P_s and z⁰=z_s (e.g., terrain height) at the lowest level (k=1), the equations are applied to each model level upward from the surface.

The increments in dry-air density (δρ_d) and potential temperature (δθ) are also derived at each level through the linearized formulas below. (As all the variables are computed at the same level of k, we omit the superscript k.)

$$\begin{array}{ll}{\displaystyle }{\mathit{\rho }}_{\mathrm{d}}+\mathit{\delta }{\mathit{\rho }}_{\mathrm{d}}& {\displaystyle }\approx {\mathit{\rho }}_{\mathrm{d}}+{\displaystyle \frac{\partial {\mathit{\rho }}_{\mathrm{d}}}{\partial P}}\cdot \mathit{\delta }P+{\displaystyle \frac{\partial {\mathit{\rho }}_{\mathrm{d}}}{\partial T_{\mathrm{v}}}}\cdot \mathit{\delta }{T}_{\mathrm{v}}+{\displaystyle \frac{\partial {\mathit{\rho }}_{\mathrm{d}}}{\partial q_{\mathrm{v}}}}\cdot \mathit{\delta }{q}_{\mathrm{v}}\\ {\displaystyle }& {\displaystyle }\approx {\displaystyle \frac{P}{R_{\mathrm{d}}{T}_{\mathrm{v}}(\mathrm{1}+q_{\mathrm{v}})}}+{\displaystyle \frac{\mathrm{1}}{R_{\mathrm{d}}{T}_{\mathrm{v}}(\mathrm{1}+q_{\mathrm{v}})}}\cdot \mathit{\delta }P\\ {\displaystyle }& {\displaystyle }-{\displaystyle \frac{P}{R_{\mathrm{d}}{T_{\mathrm{v}}}^{\mathrm{2}}(\mathrm{1}+q_{\mathrm{v}})}}\cdot \mathit{\delta }{T}_{\mathrm{v}}-{\displaystyle \frac{P}{R_{\mathrm{d}}{T}_{\mathrm{v}}(\mathrm{1}+q_{\mathrm{v}}{)}^{\mathrm{2}}}}\cdot \mathit{\delta }{q}_{\mathrm{v}}\\ {\displaystyle }\mathit{\theta }+\mathit{\delta }\mathit{\theta }& {\displaystyle }\approx \mathit{\theta }+{\displaystyle \frac{\partial \mathit{\theta }}{\partial T}}\cdot \mathit{\delta }T+{\displaystyle \frac{\partial \mathit{\theta }}{\partial P}}\cdot \mathit{\delta }P\\ {\displaystyle }& {\displaystyle }\approx T{\left({\displaystyle \frac{P_{\mathrm{0}}}{P}}\right)}^{\frac{R_{\mathrm{d}}}{C_p}}+{\left({\displaystyle \frac{P_{\mathrm{0}}}{P}}\right)}^{\frac{R_{\mathrm{d}}}{C_p}}\cdot \mathit{\delta }T\\ {\displaystyle }& {\displaystyle }-{\displaystyle \frac{R_{\mathrm{d}}}{C_p}}{\displaystyle \frac{T}{P}}{\left({\displaystyle \frac{P_{\mathrm{0}}}{P}}\right)}^{\frac{R_{\mathrm{d}}}{C_p}}\cdot \mathit{\delta }P\end{array}$$

Here, the first term on the right-hand side (RHS) is replaced with the background, and the full fields (e.g., variables without δ in front of them) represent the prior states as well.