Articles | Volume 15, issue 20
Geosci. Model Dev., 15, 7859–7878, 2022
Geosci. Model Dev., 15, 7859–7878, 2022
Development and technical paper
26 Oct 2022
Development and technical paper | 26 Oct 2022

Data assimilation for the Model for Prediction Across Scales – Atmosphere with the Joint Effort for Data assimilation Integration (JEDI-MPAS 1.0.0): EnVar implementation and evaluation

Data assimilation for the Model for Prediction Across Scales – Atmosphere with the Joint Effort for Data assimilation Integration (JEDI-MPAS 1.0.0): EnVar implementation and evaluation
Zhiquan Liu1, Chris Snyder1, Jonathan J. Guerrette1, Byoung-Joo Jung1, Junmei Ban1, Steven Vahl1,a, Yali Wu1,b, Yannick Trémolet2, Thomas Auligné2, Benjamin Ménétrier2, Anna Shlyaeva2, Stephen Herbener2, Emily Liu2,c, Daniel Holdaway2,d, and Benjamin T. Johnson2 Zhiquan Liu et al.
  • 1National Center for Atmospheric Research, Boulder, Colorado 80301, USA
  • 2Joint Center for Satellite Data Assimilation, Boulder, Colorado 80301, USA
  • anow at: Joint Center for Satellite Data Assimilation, Boulder, USA
  • bnow at: Shenzhen Institute of Meteorological Innovation, Shenzhen, China
  • cnow at: NOAA National Centers for Environmental Prediction, College Park, USA
  • dnow at: NASA Goddard Space Flight Center, Greenbelt, USA

Correspondence: Zhiquan Liu (


On 24 September 2021, JEDI-MPAS 1.0.0, a new data assimilation (DA) system for the Model Prediction Across Scales – Atmosphere (MPAS-A) built on the software framework of the Joint Effort for Data assimilation Integration (JEDI) was publicly released for community use. Operating directly on the native MPAS unstructured mesh, JEDI-MPAS capabilities include three-dimensional variational (3DVar) and ensemble–variational (EnVar) schemes as well as the ensemble of DA (EDA) technique. On the observation side, one advanced feature in JEDI-MPAS is the full all-sky approach for satellite radiance DA with the introduction of hydrometeor analysis variables. This paper describes the formulation and implementation of EnVar for JEDI-MPAS. JEDI-MPAS 1.0.0 is evaluated with month-long cycling 3DEnVar experiments with a global 30–60 km dual-resolution configuration. The robustness and credible performance of JEDI-MPAS are demonstrated by establishing a benchmark non-radiance DA experiment, then incrementally adding microwave radiances from three sources: Advanced Microwave Sounding Unit-A (AMSU-A) temperature sounding channels in clear-sky scenes, AMSU-A window channels in all-sky scenes, and Microwave Humidity Sounder (MHS) water vapor channels in clear-sky scenes. JEDI-MPAS 3DEnVar behaves well with a substantial and significant positive impact obtained for almost all aspects of forecast verification when progressively adding more microwave radiance data. In particular, the day 5 forecast of the best-performing JEDI-MPAS experiment yields an anomaly correlation coefficient (ACC) of 0.8 for 500 hPa geopotential height, a gap of roughly a half day when compared to cold-start forecasts initialized from operational analyses of the National Centers for Environmental Prediction, whose ACC does not drop to 0.8 until a lead time of 5.5 d. This indicates JEDI-MPAS's great potential for both research and operations.

1 Introduction

Based upon probability theory, data assimilation (DA) seeks to find the optimal estimate of a model state by combining observations and prior information of the state (typically from a short-term model forecast) as well as the associated error characterizations (Tarantola2005). DA systems are an integral part of numerical weather prediction (NWP), as they provide initial conditions for numerical forecasts at operational forecast centers. Some DA systems are openly available and widely used by the community, e.g., the Data Assimilation Research Testbed (DART) (Anderson et al.2009) and the Weather Research and Forecasting model's data assimilation (WRFDA) system (Barker et al.2012) developed by the National Center for Atmospheric Research (NCAR) and the Gridpoint Statistical Interpolation (GSI) DA system (Shao et al.2016) developed at the National Centers for Environmental Prediction (NCEP). While DART implements the ensemble Kalman filter (EnKF) and variants, WRFDA and GSI mainly emphasize variational algorithms. Those DA systems are mainly programmed in Fortran.

Since 2018, the Joint Effort for Data assimilation Integration (JEDI) project (Trémolet and Auligné2020) led by the Joint Center for Satellite Data Assimilation (JCSDA) has been developing a community DA framework centered on generic, model-agnostic components, which can be interfaced to a variety of models such as those for atmosphere, ocean, land, and atmospheric chemistry. JEDI adopts object-oriented and generic programming and uses a mix of C++ and Fortran (Trémolet2020).

There are four main generic components in the JEDI framework. The Object-Oriented Prediction System (OOPS;, last access: 21 October 2022) defines the abstract representations of elements of several DA algorithms, including three- and four-dimensional variational (3D/4DVar) (e.g., Lorenc1986; Rabier et al.2007; Liu et al.2020) and 3D/4D ensemble–variational (3D/4DEnVar) (e.g., Wang et al.2013; Kleist and Ide2015) techniques with several minimization algorithms implemented. OOPS is mainly programmed in C++ using templates. The System-Agnostic Background Error Representation (SABER;, last access: 21 October 2022) repository hosts background error covariance models. All experiments discussed here employ the Background error on an Unstructured Mesh Package (BUMP) (Ménétrier2020) from SABER, which estimates and applies background error-covariance-related operators, defined on an unstructured mesh. Two observation-related components are the Unified Forward Operator (UFO;, last access: 21 October 2022) and the Interface for Observation Data Access (IODA;, last access: 21 October 2022) (Honeyager et al.2020). UFO implements observation operators for various types of observations. In addition, UFO implements other observation-space functionality, including quality control, data thinning, and bias correction. IODA provides functionality for input, output, and in-memory data access of observations. In addition to those model-agnostic components of JEDI described above, there are model-specific implementations of interfaces to those generic components that allow for the realization of model-specific DA systems within the JEDI framework (Holdaway et al.2020).

Ha et al. (2017) developed a DA system for the Model for Prediction Across Scales – Atmosphere (MPAS-A) (Skamarock et al.2012) based upon DART (i.e., MPAS-DART). MPAS-DART can only perform ensemble analysis with the ensemble Kalman filter and variants and lacks variational DA capability for deterministic analysis, which must be addressed, especially given the pervasive use of variational DA systems at operational centers. We also think that it is beneficial to implement different DA schemes for both deterministic and ensemble analyses to meet a variety of needs for research and operational applications using MPAS DA. JEDI provides such a software framework with various DA algorithms and standardizes to some extent the model interfaces to different components of JEDI. Therefore, we decided to develop a new DA system for MPAS-A based upon the JEDI framework at a similar time in 2018. In the past more than 3 years, we have developed the JEDI-MPAS DA system, which consists of MPAS interfaces to JEDI (, last access: 21 October 2022) and the four model-agnostic JEDI components. The first-release JEDI-MPAS 1.0.0 was made publicly available in September 2021. JEDI-MPAS capabilities include three-dimensional variational (3DVar) and ensemble–variational (EnVar) schemes for deterministic analyses, the ensemble of DA (EDA) technique to produce ensembles of analyses, and dual resolution, where the analysis increment and the input forecast ensemble are at a lower resolution than that of the background and analysis, thus decreasing computational requirements. Moreover, JEDI-MPAS operates directly on the native MPAS unstructured mesh and can be applied without code modifications to all MPAS meshes, whether uniform, variable resolution, global, or regional. On the observation side, one advanced feature in JEDI-MPAS is the full all-sky approach for satellite radiance DA with the introduction of hydrometeor analysis variables.

This article describes MPAS-specific aspects of JEDI-MPAS, emphasizing its implementation of EnVar, and documents JEDI-MPAS performance in global cycling experiments. The implementation and evaluation of JEDI-MPAS's EDA and 3DVar algorithms will be presented in separate articles. The next section provides a brief introduction to MPAS-A and modifications made for JEDI-MPAS. Section 3 describes the formulation and implementation of JEDI-MPAS's EnVar, followed by the experimental design in Sect. 4. The results of global cycling experiments are presented in Sect. 5, and the paper is concluded with Sect. 6.

2 MPAS-A model

MPAS-A is a nonhydrostatic model discretized horizontally on an unstructured centroidal Voronoi mesh with C-grid staggering of state variables (Skamarock et al.2012). It is run-time configurable to use either quasi-uniform or variable-resolution meshes and for either global or regional applications (Skamarock et al.2018). MPAS-A adopts a height-based terrain-following vertical coordinate ζ following Klemp (2011), in which the coordinate surfaces are progressively smoothed to remove smaller-scale terrain structure.

MPAS-A's prognostic equations employ flux variables that are coupled with density; that is, the prognostic variables in the discretized equations are

(1) X = ( U , W , Θ m , Q j , ρ ̃ d ) = ( ρ ̃ d u , ρ ̃ d w , ρ ̃ d θ m , ρ ̃ d q j , ρ ̃ d ) ,

where ρd is the density of dry air, ρ̃d=ρd(ζ/z)-1, u is the normal component of horizontal velocity at cell edges, w is the vertical velocity at cell centers, qj represents the mixing ratio of the respective water species (water vapor, cloud liquid water, cloud ice, rain, snow, graupel, hail) denoted by (qv,qc,qi,qr,qs,qg,qh), and θm=θd[1+(Rv/Rd)qv] is a modified moist potential temperature, with θd the dry potential temperature and Rd and Rv the gas constants for dry air and water vapor, respectively. The continuous prognostic equation set can symbolically be expressed as

(2) X t = F ( p , ρ d ρ m , X , ) ,

where pressure p is obtained diagnostically from the equation of state, and ρm=ρd(1+qv+qc+qi+qr+qs+qg+qh) is the density of moist air.

The fully compressible nonhydrostatic equations are solved using the time-split integration technique described by Klemp et al. (2007). To reduce the truncation error associated with the horizontal pressure gradient terms and the roundoff error associated with the right-hand-side terms in the vertical momentum equation, the governing equations are rewritten in terms of perturbation thermodynamic variables (ρd,Θm,p) relative to a hydrostatically balanced reference state that is only a function of the geometric height z.

JEDI-MPAS not only uses MPAS as the forecast model for cycling data assimilation, but also leverages and extends certain capabilities within MPAS. Crucially, JEDI-MPAS shares with MPAS the data structures for MPAS fields and accompanying utilities for manipulating those fields, such as construction of new instances, parallel decomposition across processors, and reading from or writing to files. JEDI-MPAS also depends on MPAS initialization steps that occur before a model integration, in which input variables are converted to the MPAS prognostic variables given by Eq. (1) and additional thermodynamic variables, such as p or ρd, are diagnosed as needed. In addition, JEDI-MPAS modifies MPAS routines for reading files by splitting static and time-varying fields into separate files and reducing the number of time-varying fields retained (relative to a restart file). These enhancements build on a two-stream input approach for MPAS data assimilation developed by Soyoung Ha with William Skamarock, which was expanded and implemented for JEDI-MPAS cycling experiments. All changes necessary for JEDI-MPAS 1.0.0 are in a modified version of MPAS-A (, last access: 21 October 2022) that is based upon the latest MPAS-A version 7.1.

The methodology of JEDI-MPAS DA in updating the MPAS model variables using the EnVar technique is given in the next section.

3 Methodology of JEDI-MPAS EnVar

3.1 Incremental variational analysis

The non-quadratic form of the cost function for a variational analysis is expressed as

(3) J ( x ) = 1 2 ( x - x b ) T B - 1 ( x - x b ) + 1 2 [ h ( x ) - y ] T R - 1 [ h ( x ) - y ] ,

where the column vector x represents the model analysis variables distributed in a model's horizontal grid and vertical coordinate (and at different times if it is for a four-dimensional analysis), the subscript “b” indicates a background field which is typically from a short-term model forecast, and the column vector y is the various types of observations irregularly distributed in space and time. h is the observation operator that interpolates the model fields to observation locations and transforms the analysis variables x to the observed quantities. B and R are the background and observation error covariance matrices, respectively. The superscripts “−1” and “T” denote the inverse of a matrix and the transpose of a column vector, respectively.

With the introduction of a first-guess variable xg and letting δx=x-xg and δxg=xb-xg, Eq. (3) can be rewritten in the incremental form (Courtier et al.1994) as

(4) J ( δ x ) = 1 2 ( δ x - δ x g ) T B - 1 ( δ x - δ x g ) + 1 2 ( H δ x - d ) T R - 1 ( H δ x - d ) ,

where d=y-h(xg) is the first-guess departure to observations (also called “the innovation vector”), and H is the linearized version of h in the vicinity of xg. The gradient of the cost function with respect to δx

(5) δ x J ( δ x ) = B - 1 ( δ x - δ x g ) + H T R - 1 ( H δ x - d )

shall vanish at the value δxa that minimizes the cost function, which leads to

(6) ( B - 1 + H T R - 1 H ) δ x a = B - 1 δ x g + H T R - 1 d .

Equation (6) is a linear algebra system Aδxa=b that can be solved iteratively with a minimization algorithm. The analysis of the full field can be found with xa=xg+δxa. It is common practice that the minimization iterations (also referred to as “inner loop”) start from xg=xb and δxg=0. Multiple outer loops can be used with the first-guess xg (and thus δxg) updated from the analysis of the previous outer loop. In each outer loop, d is recalculated, and H is re-linearized using the updated xg.

Several minimization algorithms are available in OOPS, and the most extensively tested one in JEDI-MPAS is the Derber–Rosati Inexact Preconditioned Conjugate Gradient (DRIPCG) algorithm (Golub and Ye1999; Derber and Rosati1989), which is documented in (last access: 21 October 2022). It is important to note that the DRIPCG algorithm introduces an additional initial vector δx^g=B-1δxg at the beginning of the inner loop so that the application of B−1 can be avoided if the initial guess is taken as δxg=0 (the case of the first outer loop) or if both δxg and δxg^ are retained from a previous outer loop. In each iteration of minimization, DRIPCG only requires the application of the background error covariance matrix B, which will be described in the next subsection.

3.2 Background error covariance model for EnVar

The modeling of the background error covariance matrix B is a key component of variational DA systems. OOPS's hybrid-EnVar follows the original form of Hamill and Snyder (2000), in which B in Eq. (4) is the weighted sum of the static background error covariance Bs and an ensemble-based covariance Be; i.e.,

(7) B = β s B s + β e L B e ,

where βs and βe are the scalar weights with βs+βe=1, and LBe denotes the Schur product (element by element) of the localization matrix L and the sample covariance matrix Be of the ensemble. Note that L is a correlation matrix, with diagonal elements being 1 and off-diagonal elements smaller than 1 that reduce to zero for a certain distance between two model grid points. Localization ameliorates the detrimental effects of sampling noise in Be by reducing spurious correlations between state variables separated by large distances and by increasing the rank of the localized sample covariance relative to Be (Hamill et al.2001; Houtekamer and Mitchell2001).

The modeling of the static covariance Bs for JEDI-MPAS will be presented in a separate paper. This article only describes the pure EnVar algorithm, in which the weight of the static covariance is zero with βs=0 and βe=1. In each iteration of DRIPCG minimization, Brk is calculated, where rk=b-Aδxk is the residual vector at the kth iteration, and A and b are defined by Eq. (6). In the case of pure EnVar with B=LBe, Brk is actually evaluated in OOPS in the form

(8) ( L B e ) r = m = 1 N e x m L ( x m r ) ,

where Ne is the ensemble size, and xm is (Ne-1)-1/2 times the deviation of the mth member from the ensemble mean. That is, (LBe)r is given by taking the Schur product of each xm with r, passing the result through L, taking the Schur product of the result again with xm, and finally summing over m.

To see that Eq. (8) is correct, write (LBe)ij=Lijm=1Neximxjm and so the ith element of the vector (LBe)r

(9) [ ( L B e ) r ] i = j = 1 N x m = 1 N e L i j x i m x j m r j = m = 1 N e x i m j = 1 N x L i j x j m r j = m = 1 N e x i m [ L ( x m r ) ] i ,

where Nx is the dimension of the state vector. This gives Eq. (8). It is worth mentioning that OOPS does not take the α control variable form of the hybrid-EnVar cost function used in other DA systems (e.g., Lorenc2003; Wang2010; Wang et al.2013; Schwartz et al.2015), which do not perform the explicit calculation of Eq. (8).

Since the dimension of the localization matrix L is Nx×Nx, in-memory storage of the full matrix L is prohibitive, and the multiplication of L with a vector of size Nx can be computationally expensive for a high-resolution setting of the model. JEDI-MPAS 1.0.0 adopts the Normalized Interpolated Convolution from an Adaptive Subgrid (NICAS) method to model L such that

(10) L = NSCS T N T ,

where the correlation matrix C is defined on a coarse grid (i.e., subgrid), S is an interpolation operator to transform C on the coarse grid to L on the full grid, and the diagonal normalization matrix N ensures that resulting L has diagonal entries equal to 1. The local resolution of the subgrid in NICAS adapts to the correlation length-scale to reduce computational cost while retaining a reasonable approximation of the underlying continuous correlation function. Therefore, the convolution of L with a vector becomes a sequence of computationally less expensive operations with reduced memory usage compared to the use of the full matrix L. NICAS is available in the BUMP part of the SABER code repository.

BUMP/NICAS is implemented based upon an unstructured mesh that is naturally suited for JEDI-MPAS (also built on an unstructured mesh). This is a key feature that enables JEDI-MPAS to work seamlessly for global or regional domains with a quasi-uniform or variable-resolution MPAS mesh. More details about BUMP and NICAS are given by Ménétrier (2020) (, last access: 21 October 2022).

To further reduce the computational cost, a dual-resolution capability is implemented in JEDI-MPAS's EnVar, in which the background and analysis (thus first-guess departure calculation) are at the same resolution as that for the model forecast, but the analysis increment and ensemble input (thus inner loop minimization) are at a lower resolution. While this study evaluated 3DEnVar, JEDI-MPAS can be run with the 4DEnVar mode, which assimilates observations at more appropriate times by having the background and analysis files as well as the ensemble input at multiple times within the assimilation time window.

3.3 Analysis variables and variable transforms

One important aspect in designing a DA system is the choice of analysis variables (i.e., x). JEDI-MPAS uses as analysis variables the cell-center values of temperature (T), specific humidity (q), zonal and meridional components of horizontal velocity (v), and surface pressure (ps). Optionally, mixing ratios of five hydrometeors (cloud liquid water, cloud ice, rain, snow, and graupel) can also be included as analysis variables when assimilating cloud-/precipitation-affected observations. We chose these analysis variables because of our experience with modeling background-error covariances for similar variables in WRFDA and because use of T rather than θ will reduce nonlinearity in the observation operators for many conventional observations and for satellite radiances.

After obtaining the analysis increment for (T, q, v, ps) by iteratively solving Eq. (6), several additional steps are followed to arrive at a full-field analysis for variables closer to the MPAS prognostic variables:

  1. Transform the cell-center v increments to increments of u, the normal component of horizontal velocity at cell edges, by averaging the corresponding component of the v increment from the cells adjacent a given edge. Obtain the full-field analysis of u by adding increments to the first guess.

  2. Get the full-field analyses of T, q, and ps (and optionally hydrometeors) by adding analysis increments to the first guess. Calculate the water vapor mixing ratio analysis with qv=q/(1-q).

  3. Diagnose the full-field analyses of 3D pressure p and dry-air density ρd using T, qv, and ps by applying the hypsometric equation p(z2)=p(z1)exp(-gRdz1z2dzTv) and the equation of state of moist air ρd=p/[(1+qv+)RdTv] with virtual temperature Tv=(1+0.608qv)T. Note that the hypsometric equation, which implies hydrostatic balance, is applied upward layer by layer, starting from the ground with the analysis field of ps.

  4. Update dry potential temperature θd=T(p0p)Rd/cp, with cp the specific heat of dry air at constant pressure and p0=1000 hPa a reference pressure.

Sensitivity experiments showed that the hydrostatic balance constraint in Step 3 is important for JEDI-MPAS to cycle stably by allowing surface pressure observations to constrain upper-air fields.

3.4 Observation handling and all-sky radiance DA

In JEDI-MPAS, observation I/O and in-memory storage are handled through IODA, which is one of model-agnostic components of JEDI and reads specific netCDF4 (based on hdf5) format observation input files and outputs DA feedback files. UFO, another model-agnostic component of JEDI, consists of nonlinear observation operators and associated tangent linear and adjoint (TL/AD) (and/or Jacobian) operators. At the time of this study, quality control (QC) procedures and radiance bias correction algorithms within the UFO are still under development, and thus the GSI DA system (Shao et al.2016) is used as a data preprocessor to perform QC of observations and radiance bias correction. GSI observation feedback files in the so-called gsi-ncdiag format are converted into the IODA-v2 format and then assimilated into JEDI-MPAS. Therefore, mainly two QC filters are applied in JEDI-MPAS experiments (see Sect. 4). One uses QC flags generated by GSI, and another is the “background check” filter, which discards observations with first-guess departures larger than 3 times the observation error standard deviation. For surface pressure (ps) data, additional QC is applied to reject observations with differences between the model terrain height and station elevation larger than 200 m. Terrain correction is also applied to ps data following Ingleby (2014).

In recent years, leading operational NWP centers assimilate increasingly more satellite radiance data under cloudy situations, with most successes in using all-sky microwave radiances (e.g., Geer et al.2017, 2018; Zhu et al.2016, 2019; Migliorini and Candy2019). JEDI-MPAS allows the assimilation of all-sky radiance data with the introduction of hydrometeor analysis variables, cloudy radiance observation operator via the Community Radiative Transfer Model (CRTM) (Liu and Weng2006), and situation-dependent all-sky observation error models (Geer and Bauer2011). While all-sky radiance DA in JEDI-MPAS can be applied to data from both microwave and infrared sensors, this study evaluates only microwave all-sky radiance DA from Advanced Microwave Sounding Unit-A (AMSU-A) window channels 1–4 and 15 over water. The AMSU-A all-sky radiance error model follows Zhu et al. (2016), in which the observation error is a piecewise-linear ramp function (see Sect. 4) of the model–observation-averaged cloud liquid water path (cobs+cfg)/2 retrieved from the AMSU-A channel 1's and 2's brightness temperatures. Zhu et al. (2016) assimilated AMSU-A window channels only for non-precipitating clouds due to the lack of precipitation variables in the microphysics scheme used in their study. In contrast, with the WSM6 microphysics scheme (Hong and Lim2006) used in this study and all five hydrometeor variables included in the analysis variables, JEDI-MPAS assimilates AMSU-A window channel radiances in both non-precipitating and precipitating cloudy situations with a full all-sky approach. Note that microwave radiances are assimilated in terms of brightness temperature in units of kelvin. For the sake of simplicity, we will interchangeably use radiance and brightness temperature throughout the article.

4 Experimental design

We evaluate JEDI-MPAS 3DEnVar's performance with month-long (15 April–14 May 2018) cycling experiments. Although this is not an observational impact study, incrementally adding satellite radiance data and evaluating their impact is informative about the implementation correctness and the robustness and performance of the system. In this regard, four 6-hourly cycling dual-resolution 3DEnVar experiments are conducted with a global quasi-uniform 30 km mesh for high-resolution analysis/background and a 60 km mesh for low-resolution increment and ensemble input. The first analysis is at 00:00 UTC, 15 April 2018, with the background being a 6 h forecast initialized from the NCEP's Global Forecast System (GFS) analysis at 18:00 UTC, 14 April. From the second analysis through the final analysis at 18:00 UTC, 14 May, the background is a 6 h forecast initialized from the previous cycle's JEDI-MPAS analysis. The first-guess departure is computed using the 30 km background, and the minimization is performed on the 60 km mesh. Each DA analysis uses two outer loops, each with 60 inner iterations. The 20-member ensemble input on the 60 km mesh is generated from 6 h MPAS model forecasts initialized from the NCEP's 20-member Global Ensemble Forecast System (GEFS) ensemble analysis (Zhou et al.2017).

Table 1Observations assimilated in the four cycling 3DEnVar experiments.

Download Print Version | Download XLSX

The four experiments differ from each other by assimilating different sets of observations, as listed in Table 1. The conv experiment assimilates temperature, specific humidity, and zonal and meridional wind components from radiosondes and aircraft; GPS refractivity; atmospheric motion vectors (AMVs) retrieved from geostationary and polar-orbiting satellite sensors; and surface pressure (ps). In addition to the observations used in conv, the clrama experiment assimilates clear-sky AMSU-A radiances (channels 5–9) from six satellites (NOAA-15/18/19, Metop-A/B, Aqua; omitting platforms and channels that had failed before the experimental period). The cldama experiment further adds over-water AMSU-A window channels (channels 1–4 and 15) all-sky radiances from five satellites. Note that Aqua AMSU-A channels 1 and 2, which are needed for the all-sky observation error model (see Sect. 3.4), are not available, and thus Aqua AMSU-A window channels are excluded in cldama. We use the fourth experiment, clrmhs, to evaluate the additional benefit of assimilating three water vapor channels from the Microwave Humidity Sounder (MHS) on board four satellites (NOAA-18/19, Metop-A/B) under clear-sky conditions. Data thinning is applied to AMSU-A and MHS radiance data with a 145 km thinning mesh. All four experiments use a ±3 h time window.

Figure 1All-sky observation error model for NOAA-19 AMSU-A channel 15. The magenta curve is the statistics of the standard deviation of O minus B (left y axis) as a function of the model–observation-averaged cloud liquid water path (x axis). The cyan curve gives the number of observations (right y axis) used in each bin's statistics (0.01 kg m−2 bin interval). The dashed black curve is the fitted piecewise-linear ramp function of the magenta curve, serving as the all-sky observation error model. Variables in the title are those of Eq. (11) in Sect. 4.


All-sky observation error statistics were computed separately for five AMSU-A sensors using the observations minus the CRTM-simulated brightness temperatures from the 6 h forecast background fields (i.e., O minus B) of the clrama experiment over the month-long period. Figure 1 gives an example of the all-sky observation error model for NOAA-19 AMSU-A channel 15 at 89 GHz. The blue curve is the statistics of the standard deviation of O minus B as a function of the model–observation-averaged cloud liquid water path c=(cobs+cfg)/2. The red curve gives the number of observations (the right y axis) used in each bin's statistics (0.01 kg m−2 bin interval). The dashed black curve is the fitted piecewise-linear ramp function of the blue curve in the form (Geer and Bauer2011)

(11) σ o = σ clr , if  c c clr σ clr + ( σ cld - σ clr ) c - c clr c cld - c clr , if  c clr < c c cld σ cld , if  c > c cld ,

which serves as the observation error model for this particular channel with cclr=0.03, ccld=0.24, σclr=6.33, and σcld=19.24. The four numbers for each of channels 1–4 and 15 of five AMSU-A sensors are provided in the supplemental configuration file of JEDI-MPAS (jedi_mpas.yaml in the Supplement), which specifically describes the configuration for the clrmhs experiment.

BUMP has the ability to diagnose horizontal and vertical localization length scales with geographical variation. For a relatively small ensemble size, however, BUMP-diagnosed localization length scales may not be optimal. Therefore, we used the fixed localization scale values in this study, which is not uncommon for ensemble DA. After some sensitivity testing, we use 1200 and 6 km for the horizontal and vertical full widths of Gaspari and Cohn (1999) localization function (Eq. 4.10), respectively. That function reduces the spatial correlation to zero when the distance between two locations is greater than the full width.

The MPAS-A model is configured with a 30 km quasi-uniform mesh and 55 vertical levels with a 30 km model top. The “mesoscale reference” physics suite is used with the parameterization schemes listed in Table 2. A time step of 180 s is applied for the model time integration. We used the same “mesoscale reference” physics suite to generate the 6 h MPAS ensemble forecasts on the 60 km mesh that comprise the sample ensemble covariances.

(Tiedtke1989; Zhang and Wang2017)(Hong and Lim2006)(Chen and Dudhia2001)(Hong et al.2006)(Jiménez et al.2012)(Iacono et al.2008)(Xu and Randall1996)(Choi and Hong2015)

Table 2Parameterization schemes in the “mesoscale reference” physics suite.

Download Print Version | Download XLSX

5 Evaluation results

5.1 Short-term forecasts

Here we assess the DA cycling performance of the four experiments by examining the first-guess departure (i.e., observation minus background or OMB) statistics, similar to the diagnostic commonly monitored at operational NWP centers (e.g.,, last access: 21 October 2022). Additionally, we present error statistics of 6 h forecast backgrounds with respect to the NCEP's GFS analysis to more clearly assess the geographic and vertical variations of the cycling performance.

Figure 2Time series (00:00 UTC, 15 April–18:00 UTC, 14 May 2018) of (a, b) mean and (c, d) rms of OMB for AMV (a, c) zonal and (b, d) meridional wind components (units: m s−1) from the four experiments. The background is a 6 h forecast. Statistics are computed over the globe and averaged for all vertical levels.


5.1.1 Observation space verification

Figure 2 displays time series of mean and root mean square (rms) values of OMB for AMV zonal and meridional wind components from the four experiments. Overall the mean and rms give the consistent performance indication for the four experiments. For zonal wind, the conv experiment has the largest bias and root mean square error (RMSE); adding clear-sky AMSU-A radiances leads to substantial error reduction in terms of both bias and RMSE; further adding all-sky AMSU-A window-channels' radiances additionally reduces the RMSE. The red curves (clrmhs) are basically overlaid with the green curves (cldama), indicating a neutral impact on 6 h wind forecasts by adding clear-sky MHS water vapor radiances to all-sky AMSU-A radiances. The statistics of OMB for AMV meridional wind component indicate similar relative experiment performances overall, except for a slightly positive impact for clrmhs versus cldama.

Figure 3Vertical distributions of relative rms reduction of OMB for radiosondes' (a) virtual temperature, (b) specific humidity, (c) zonal wind, and (d) meridional wind. Statistics are computed over the period from 00:00 UTC, 18 April, to 18:00 UTC, 14 May 2018, for the three radiance DA experiments against the conv experiment. A negative percentage indicates the positive impact from radiance DA experiments relative to the non-radiance DA experiment.


Figure 3 compares the 6 h forecast fits to radiosonde observations from the three radiance DA experiments with those of the conv experiment. Shown in Fig. 3 are relative differences between rms of OMB from the three radiance DA experiments and that from the conv experiment, given by

(12) ratio = 100 × ( rms EXP - rms conv ) / rms conv ,

where “EXP” can be clrama, cldama, or clrmhs. Negative ratios indicate positive impacts from EXP relative to conv, and vice versa. Different magnitudes of ratio among the three radiance experiments give their relative performance. A 95 % confidence interval is shown at each pressure level as an error bar, which is computed using bootstrap resampling (Hamill1999), with a resample size of 10 000. The rms difference between EXP and conv at each cycle is treated as an independent sample in the bootstrap. Assimilating clear-sky temperature sounding channels from clrama has small impacts on humidity (blue curve in Fig. 3b), whereas large positive impacts on moisture are obtained by adding all-sky AMSU-A and clear-sky MHS radiances. The impact on virtual temperature (Tv) is mostly from AMSU-A clear-sky temperature sounding channels, and the additional benefit on Tv from all-sky AMSU-A and clear-sky MHS radiances is relatively small (Fig. 3a). The reason we show Tv here instead of T is that the majority of radiosonde temperature observations in GSI's ncdiag files are expressed in terms of Tv. The rms decreases up to 5 % for zonal wind (Fig. 3c) and 2 % for meridional wind (Fig. 3d) from clrama, likely caused by the mass–wind balance implied in the ensemble background error covariance. The extra benefit on wind fields from AMSU-A all-sky window channel radiances is confined to 200–800 hPa. Additional positive impact on winds between 300 and 650 hPa is achieved by adding MHS clear-sky water vapor channels' radiances. Each of the three radiance DA experiments exhibit, to a different extent, improvement with respect to the conv experiment for all four observed radiosonde variables at almost all pressure levels.

Figure 4Similar to Fig. 3 but for latitudinal variation of relative rms reduction of OMB for NOAA-19 AMSU-A ch. 1–3 and 15 for clrmhs versus cldama.


Among the four experiments, both cldama and clrmhs assimilate AMSU-A all-sky radiances from the window channels (ch. 1–4 and 15). It is expected a priori that assimilating MHS water vapor channels' radiances in clrmhs will improve background fitting to those AMSU-A all-sky radiances. That is indeed the case in Fig. 4; clrmhs achieves statistically significant improvement up to 2.5 % at tropical latitudes for AMSU-A channel 1 (Fig. 4a), while positive impacts from MHS radiances are smaller and less significant for other AMSU-A window channels (Fig. 4b–d). This is because AMSU-A channel 1 (23.8 GHz) is more sensitive to water vapor and clouds than other AMSU-A window channels.

Figure 5Time series (00:00 UTC, 15 April–18:00 UTC, 14 May 2018) of RMSE of cycling background with respect to GFS analyses for global (a) temperature (units: kelvin) and (b) specific humidity (units: g kg−1) at 2 m, (c) zonal wind component (units: m s−1) at 10 m, and (d) surface pressure (units: Pa) from the four experiments.


Figure 6Vertical distributions of relative RMSE reduction (with respect to GFS analyses) of cycling background for global (a) temperature, (b) water vapor mixing ratio, (c) zonal wind, and (d) meridional wind (cldama and clrmhs versus clrama). Statistics are for the period from 00:00 UTC, 18 April, to 18:00 UTC, 14 May 2018.


Figure 7Similar to Fig. 6 but for latitudinal variations of RMSE reduction for (a) surface pressure, (b) water vapor mixing ratio of model level 21 to 30, and the whole-column (c) zonal wind and (d) meridional wind.


5.1.2 Model space verification

Similar to Figs. 2–4, Figs. 5–7 evaluate the 6 h forecast backgrounds in model space relative to GFS analyses. Figure 5 shows the time series of global RMSE of the four experiments for four surface variables. Note that the background for the first analysis at 00:00 UTC, 15 April 2018, is a 6 h forecast initialized from the GFS analysis. GFS utilizes a hybrid-4DEnVar DA technique via the GSI system (Kleist and Ide2015) and includes assimilation of many more satellite observations than those used in our own cycling experiments. Smaller RMSEs for the first 2–3 d come from the memory of the superior performance of the operational analysis. Cycling performance quickly converges to a stable level for the three radiance DA experiments due to the good global coverage of observations. Gradually increasing errors are clearly seen in the conv experiment for 2 m temperature (T2), 10 m zonal wind (u10), and ps, which is caused by sparse coverage of non-radiance data over the Southern Hemisphere. The impacts of AMSU-A clear-sky and all-sky radiances and MHS clear-sky radiances are consistent with those shown in upper-air observation space of Fig. 3.

Similar to Fig. 3, Fig. 6 presents the vertical distribution of the global percentage RMSE difference of cldama and clrmhs relative to clrama. Figure 6 yields similar conclusions as Fig. 3, but the magnitudes of improvement in model space (Fig. 6) are larger than those in observation space (Fig. 3) because the model space verification covers the entire globe. The most drastic improvement is for water vapor mixing ratio, with nearly 20 % RMSE decrease at model level 15 (∼2.2 km or ∼700 hPa) by assimilating AMSU-A all-sky radiances and >20 % improvement at model level 30 (∼8.8 km or ∼330 hPa) when adding both AMSU-A all-sky radiances and MHS clear-sky radiances (Fig. 6b). One noticeable issue is the negative impact for temperature near the model top (Fig. 6a), which is believed to be related to an increasing cold bias during the cycling (not shown). A biased model top makes radiance DA less effective. Note that the MPAS-A model currently lacks processes for representing O3 chemistry and that a climatological O3 field is used, leading to inaccurate radiation calculations attributed to O3. Work is underway to add a parameterized O3 chemistry scheme (McCormack et al.2006) and include O3 as a prognostic variable in MPAS-A.

Figure 7 is similar to Fig. 6 but depicts the latitudinal variation of relative RMSE reduction for several variables. For ps, the greatest impact is between 25–50 S (Fig. 7a), and large RMSE reductions are also seen in North and South Pole regions for the clrmhs experiment (Fig. 7a). The large RMSE reduction from clrmhs in the Arctic is apparently associated with a reduction of a ∼1 hPa low pressure bias in the conv, clrama, and cldama experiments (not shown). For the other three variables, the greatest impact is within tropical regions (Fig. 7b–d). In addition to the large benefit on moisture from the addition of MHS water vapor channels' radiances (Fig. 7b), the impact of MHS radiances on the meridional wind component (Fig. 7d) is also substantial and significant, possibly indicating an improved representation of equatorial convergence–divergence patterns associated with convective activity.

These observation-space and model-space evaluations of short-term forecasts demonstrate that the JEDI-MPAS 3DEnVar scheme behaves well for the assimilation of clear-sky and all-sky microwave radiance data. More accurate background fields are obtained by adding more and more radiance data.

5.2 The 1–10 d forecasts

This subsection examines the 1–10 d forecast performance relative to GFS analyses, independent radiance data from the Advanced Himawari Imager (AHI), and a global precipitation product. A total of 27 10 d forecasts initialized each day at 00:00 UTC from 18 April to 14 May 2018 are used to produce the error statistics. The first 3 d (15–17 April) is considered the spin-up period for each experiment to converge from the GFS climatology to JEDI-MPAS's own cycling climatology and thus are not included in the forecast evaluation. The three radiance DA experiments significantly outperform the conv experiment. Therefore, we only evaluate the cldama and clrmhs experiments, from here considering the clrama experiment as the benchmark.

Figure 8Relative RMSE reductions (with respect to GFS analyses) as a function of forecast lead time for global (a) surface pressure, (b) specific humidity at 2 m, (c) zonal wind at 10 m, (d) temperature and (e) water vapor mixing ratio below model level 30, and (f) the whole-column zonal wind (cldama and clrmhs versus clrama). Statistics are aggregated over all 27 forecasts.


5.2.1 Model space verification

Figure 8 displays relative RMSE reductions (with respect to GFS analyses) as a function of forecast lead time for three surface and three upper-air variables. Due to a temperature bias at the model top, as described in Sect. 5.1 in relation to Fig. 6a, the evaluation on temperature is limited to below model level 30. We also exclude the moisture assessment above model level 30, where the air is very dry. In general, adding AMSU-A all-sky radiances produces more accurate forecasts for all variables compared to the AMSU-A clear-sky radiance experiment, and extra MHS clear-sky radiances further reduce forecast errors. Percentage improvements and associated statistical significance decrease with forecast range but can last up to 10 d, especially for the clrmhs experiment. Note that forecast errors increase with forecast lead time, and thus absolute error reduction could be increased or flat at certain forecast lead times, even with decreasing percentage improvement. The greatest improvement from clrmhs is for upper-air moisture with up to 10 % and >1 % RMSE reduction for day 1 and day 10, respectively (Fig. 8e). The improvement for moisture remains statistically significant for the whole 10 d forecast range, while that is not the case for other variables. This result is consistent with AMSU-A window channels' and the three MHS channels' large sensitivity to water vapor (and clouds for AMSU-A all-sky radiances). Large impacts are also seen on Q2m (Fig. 8b), though to a lesser extent than is seen for upper-air moisture. Conversely, the impact is the smallest on temperature (Fig. 8d), with only ∼1 % and 3 % improvement for the day 1 forecast from cldama and clrmhs, respectively. The magnitudes of improvement on ps and the zonal wind are comparable to each other, gradually decreasing from an initial ∼5 % rms reduction for the day 1 forecast.

Figure 9Similar to Fig. 8 but for water vapor mixing ratio below model level 30 and zonal wind over the three regions (30–90 N, 30 S–30 N, and 30–90 S).


Figure 9 shows the impact on moisture and zonal wind over three different regions. The largest impact is over the tropical region for both variables, which is consistent with the short-term forecast impact as shown in Fig. 7b and c. The smallest impact is between 30 and 90 N, where conventional observations have good coverage. One interesting thing is that the impact over the tropics exhibits a clear decreasing trend with forecast range, but this decreasing trend is not always present poleward of 30. This might be related to the different levels of convective weather activity between low and middle/high latitudes.

Figure 10Similar to Fig. 8 but for clrmhs versus clrdma of four upper-air variables.


Figure 10 more clearly depicts the relative impact of clrmhs with respect to cldama for four upper-air variables. Again, the greatest and most significant improvement (up to 7 d) achieved by adding MHS water vapor channels' radiances is on moisture with 2.5 % RMSE reduction for day 1 and still retaining ∼1 % for day 7 (Fig. 10b). The magnitude of MHS radiance impact on temperature (Fig. 10a) is comparable with that of moisture but less statistically significant. The MHS impact on the wind field is smaller, but the percentage of improvement remains nearly unchanged, with the forecast lead time for up to day 8 for zonal wind (Fig. 10c) and day 6 for meridional wind (Fig. 10d). Similar to the impact on short-term forecasts shown in Fig. 7, the MHS radiance impact on meridional wind (∼1.5 %) is 50 % greater than that for the zonal wind (∼1 %) up to day 6.

5.2.2 Verification against AHI radiances

JEDI-MPAS also has the capability of simulating and assimilating clear-sky and cloudy infrared (IR) radiances such as those from the geostationary satellite sensor Advanced Himawari Imager (AHI). CRTM's cloudy IR radiance simulation capability with JEDI-MPAS is used to evaluate 1–10 d forecasts using independent (i.e., not assimilated) AHI full-disk radiance data centered at 140.7 E. Using the JEDI-MPAS HofX3D application, and in conjunction with UFO, CRTM takes the forecast temperature, moisture, and hydrometeor profiles at AHI locations to compute modeled AHI radiances. Then we compare the simulated and observed AHI brightness temperatures offline. Before being ingested into JEDI-MPAS, the raw AHI IR channels' data at 2 km resolution undergo cloud detection at each pixel following Wu et al. (2020) and are preprocessed into super-observations averaged over 15×15 pixels. The cloud detection and superobbing are carried out offline using WRFDA. These super-observations better match the model's 30 km mesh, reducing representativeness errors. A by-product of cloud detection on the raw AHI pixels and superobbing is the cloud fraction for each super-observation, which allows for statistical evaluation separately over clear, partly cloudy, and overcast observations.

Figure 11Similar to Fig. 8 but for RMSE in AHI's six channels' brightness temperature space.


Figure 11 shows the relative RMSE reduction as a function of forecast range in AHI radiance space for all superobbed pixels of six channels. Among them, channels 8–10 are water vapor (WV) channels sensing the middle-to-upper tropospheric moisture, and bands 13–15 are window channels sensitive to surface properties and clouds. The patterns of the RMSE reduction across different channels are very similar and overall consistent with those in the model space variables shown in Figs. 8 and 9. The impact of cldama and clrmhs for the higher peaking channel 8 (350–400 hPa) is larger than that for the lower peaking channel 10 (600–700 hPa), which is consistent with the vertical distribution of humidity RMSE reduction shown in Fig. 6b. Improved fitting to the three window channels (relative to clrama) is most likely evidence of improved cloud forecasts, with significance at least up to day 4. Figure 12 further splits error statistics into clear, partly cloudy, and fully cloudy pixels for the lower peaking WV channel 10. Over the verification period, clear, partly cloudy, and fully cloudy pixels account for 28 %, 33 %, and 39 % of observations, respectively. The relative impact from cldama and clrmhs decreases with more cloudiness. The absolute RMSE reduction, however, is comparable for the three categories due to increased RMSE with more cloudiness.

Figure 12Similar to Fig. 11 but for AHI channel 10's error statistics over (a) clear, (b) partly cloudy, and (c) fully cloudy pixels, respectively.


Figure 13(a) Gilbert skill scores (GSSs) of 24 h accumulated precipitation forecast verified against CMORPH observations for a threshold of 10 mm. (b) Relative difference of GSSs between cldama/clrmhs and clrama, also for a threshold of 10 mm.


5.2.3 Precipitation forecast verification

The 24 h accumulated precipitation forecasts from the four experiments are verified against the Climate Prediction Center's morphing technique (CMORPH) 0.25 precipitation product (Xie et al.2017). Figure 13a gives the Gilbert skill scores (GSSs) (also known as equitable threat scores, ETSs) (Hogan et al.2010) computed for the observed and predicted 24 h accumulated rainfall ≥10 mm using the Model Evaluation Tool (MET) (Brown et al.2021). The GSS ranges from -1/3 to 1, with ≤0 corresponding to no skill and 1 to perfect skill. An improvement of 5 %–8 % is gained up to day 7 precipitation forecast from AMSU-A temperature-sounding radiance DA when compared to the non-radiance DA experiment. A further 5 %–12 % improvement is achieved by adding AMSU-A all-sky radiances and MHS radiances for up to day 6 (Fig. 13b). From precipitation forecast verification here and other aspects of forecast verification shown earlier, JEDI-MPAS's AMSU-A all-sky DA exhibits a larger positive impact than that at operational centers (e.g., Zhu et al.2016; Migliorini and Candy2019). This might be partly related to the use of the full all-sky approach in JEDI-MPAS, whereas Zhu et al. (2016) and Migliorini and Candy (2019) assimilated AMSU-A all-sky radiances only under non-precipitating conditions. Also, obtaining additional benefit above a well-tuned operational system with abundant satellite observations, as in Zhu et al. (2016) and Migliorini and Candy (2019), is much harder than for the benchmark experimental settings used in this study.

Figure 14Global anomaly correlation coefficients of 500 hPa geopotential height as a function of forecast lead time for the four cycling experiments together with the cold-start forecasts initialized from GFS analyses (blue).


5.2.4 Anomaly correlation coefficient for 500 hPa geopotential height

Figure 14 provides global anomaly correlation coefficients (ACCs) of 500 hPa geopotential height as a function of forecast range for the four JEDI-MPAS experiments together with the cold-start forecast experiment using GFS analyses as initial conditions (light blue, marked as cold). The climatology derived from the 1980–2010 NCEP/NCAR reanalysis products (Kalnay et al.1996) is used in the ACC calculation. The largest improvement on ACC is from AMSU-A temperature sounding channels, gaining a half day of skill (from 4.5 to 5 d) with respect to the non-radiance experiment for 0.8 ACC level. Further ACC improvement by adding AMSU-A all-sky window channels and MHS clear-sky WV channels is visible but much smaller than that gained by AMSU-A clear-sky temperature-sensitive channels. This could be understandable with the tight link between temperature and geopotential height. The cold-start forecasts from GFS analyses reach 5.5 d for 0.8 ACC level and outperforms JEDI-MPAS's three radiance DA experiments by half a day. Given that the operational GFS analysis is produced with a more advanced hybrid-4DEnVar algorithm (Kleist and Ide2015) at a higher resolution and with ensemble covariances formed by 80-member ensemble input, plus the assimilation of many more satellite observations, this half-day gap on ACC is not surprising with a less advanced pure 3DEnVar algorithm, a smaller ensemble size, and fewer observations assimilated in JEDI-MPAS experiments.

6 Conclusions and future perspectives

A new DA system for the MPAS-A model has been developed based upon the JEDI software framework and was publicly released as JEDI-MPAS 1.0.0 for community use. This article describes the EnVar's formulation and technical implementation of JEDI-MPAS and evaluates its robustness and performance through four global cycling 3DEnVar experiments. JEDI-MPAS is robust, with stable cycling over a month-long period, and behaves well, with positive impacts obtained when adding more and more microwave radiance data for most aspects of assessment in both observation and model space. One advanced feature developed with JEDI-MPAS is all-sky radiance DA capability along with the introduction of hydrometeor analysis variables, which is first applied to the assimilation of all-sky radiances from AMSU-A's five window channels. Large impacts are achieved from AMSU-A all-sky DA, especially for moisture, clouds, and precipitation forecasts. Credible forecast skill is seen from the assimilation of a subset of conventional and satellite observations used in operational NWP centers, with a half-day skill gap relative to forecasts initialized from GFS analyses for the 0.8 ACC level of 500 hPa geopotential height. It is expected that the forecast performance can be further improved to be closer to operational skill using a more advanced DA algorithm such as 4DEnVar, a larger ensemble size, more satellite data, a higher model resolution, and fine tuning of DA settings, as well as by improving the MPAS model. JEDI-MPAS can also perform deterministic analyses using 3DVar and ensemble analyses using the EDA technique. The implementation and evaluation of 3DVar and EDA will be described in separate articles. JEDI-MPAS in regional mode and with a variable resolution will also be evaluated with higher resolutions.

The general JEDI software framework and more specific JEDI-MPAS DA system are still undergoing active development since the JEDI-MPAS 1.0.0 release. One new development specific to JEDI-MPAS is to change the hydrostatic balance constraint in the full fields to the increment fields, which can still allow surface pressure observations to constrain the upper-air variables but is better suited for the non-hydrostatic MPAS-A model. Initial tests of this new linear hydrostatic balance constraint show improvement, especially near the model top. Other ongoing work on the model side is to introduce a parameterized O3 chemistry scheme (McCormack et al.2006) and include O3 as a prognostic variable in MPAS-A, which should reduce temperature bias in the stratosphere and also allow for the assimilation of satellite-based O3 observations. Improving stratosphere processes will allow the use of a high model top, which will enable the assimilation of higher-peaking radiance data. All-sky radiance DA for both microwave and infrared sensors continues to be refined with JEDI-MPAS. One important feature missing in the released JEDI/UFO but available in the “develop” branch of the UFO is the variational bias correction (VarBC) for satellite radiance data (Auligné et al.2007; Dee and Uppala2009; Zhu et al.2013), which is undergoing active testing and evaluation with JEDI-MPAS at the time of writing. The quality control procedures of several types of observations implemented in the UFO are also being tested with assimilation of raw NCEP-bufr converted observation input files, with the goal of eliminating reliance on GSI for quality control. Moreover, the computational efficiency of JEDI is being improved. We intend to periodically release new and improved versions of JEDI-MPAS in upcoming years.

Code and data availability

JEDI-MPAS 1.0.0 is publicly released on GitHub, accessible from (last access: 21 October 2022). It is also available on Zenodo at (Joint Center for Satellite Data Assimilation and National Center for Atmospheric Research2021). Global Forecast System analysis data are downloaded from NCAR Research Data Archive (last access: 21 October 2022; National Centers For Environmental Prediction/National Weather Service/NOAA/U.S. Department Of Commerce2015). Global Ensemble Forecast System ensemble analysis data are downloaded from (last access: 21 October 2022). Conventional and satellite observations assimilated are downloaded from (last access: 21 October 2022; National Centers For Environmental Prediction/National Weather Service/NOAA/U.S. Department Of Commerce2008) and (last access: 21 October 2022; National Centers For Environmental Prediction/National Weather Service/NOAA/U.S. Department Of Commerce2009). CMORPH precipitation product and AHI radiance data used in the forecast verification are downloaded from (last access: 21 October 2022; Xie et al.2019) and (last access: 21 October 2022), respectively.


The supplement related to this article is available online at:

Author contributions

ZL designed and conducted all experiments and wrote the paper. All co-authors contributed to the development of JEDI-MPAS and the article.

Competing interests

The contact author has declared that none of the authors has any competing interests.


Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.


The National Center for Atmospheric Research is sponsored by the National Science Foundation of the United States. We thank Craig Schwartz for reviewing an initial version of the manuscript and making valuable suggestions that improved the manuscript. We also thank three anonymous reviewers for valuable comments and suggestions.

Financial support

This research has been supported by the United States Air Force (grant no. NA21OAR4310383).

Review statement

This paper was edited by Po-Lun Ma and reviewed by three anonymous referees.


Anderson, J., Hoar, T., Raeder, K., Liu, H., Collins, N., Torn, R., and Avellano, A.: The Data Assimilation Research Testbed: A Community Facility, B. Am. Meteorol. Soc., 90, 1283–1296,, 2009. a

Auligné, T., McNally, A. P., and Dee, D. P.: Adaptive bias correction for satellite data in a numerical weather prediction system, Q. J. Roy. Meteor. Soc., 133, 631–642,, 2007. a

Barker, D., Huang, X.-Y., Liu, Z., Auligné, T., Zhang, X., Rugg, S., Ajjaji, R., Bourgeois, A., Bray, J., Chen, Y., Demirtas, M., Guo, Y.-R., Henderson, T., Huang, W., Lin, H.-C., Michalakes, J., Rizvi, S., and Zhang, X.: The Weather Research and Forecasting Model's Community Variational/Ensemble Data Assimilation System: WRFDA, B. Am. Meteorol. Soc., 93, 831–843,, 2012. a

Brown, B., Jensen, T., Gotway, J. H., Bullock, R., Gilleland, E., Fowler, T., Newman, K., Adriaansen, D., Blank, L., Burek, T., Harrold, M., Hertneky, T., Kalb, C., Kucera, P., Nance, L., Opatz, J., Vigh, J., and Wolff, J.: The Model Evaluation Tools (MET): More than a Decade of Community-Supported Forecast Verification, B. Am. Meteorol. Soc., 102, E782–E807,, 2021. a

Chen, F. and Dudhia, J.: Coupling an Advanced Land Surface–Hydrology Model with the Penn State–NCAR MM5 Modeling System. Part II: Preliminary Model Validation, Mon. Weather Rev., 129, 587–604,<0587:caalsh>;2, 2001. a

Choi, H.-J. and Hong, S.-Y.: An updated subgrid orographic parameterization for global atmospheric forecast models, J. Geophys. Res.-Atmos., 120, 12445–12457,, 2015. a

Courtier, P., Thépaut, J.-N., and Hollingsworth, A.: A strategy for operational implementation of 4D-Var, using an incremental approach, Q. J. Roy. Meteor. Soc., 120, 1367–1387,, 1994. a

Dee, D. P. and Uppala, S.: Variational bias correction of satellite radiance data in the ERA-Interim reanalysis, Q. J. Roy. Meteor. Soc., 135, 1830–1841,, 2009. a

Derber, J. and Rosati, A.: A Global Oceanic Data Assimilation System, J. Phys. Oceanogr., 19, 1333–1347,<1333:agodas>;2, 1989. a

Gaspari, G. and Cohn, S. E.: Construction of correlation functions in two and three dimensions, Q. J. Roy. Meteor. Soc., 125, 723–757,, 1999. a

Geer, A. J. and Bauer, P.: Observation errors in all-sky data assimilation, Q. J. Roy. Meteor. Soc., 137, 2024–2037,, 2011. a, b

Geer, A. J., Baordo, F., Bormann, N., Chambon, P., English, S. J., Kazumori, M., Lawrence, H., Lean, P., Lonitz, K., and Lupu, C.: The growing impact of satellite observations sensitive to humidity, cloud and precipitation, Q. J. Roy. Meteor. Soc., 143, 3189–3206,, 2017. a

Geer, A. J., Lonitz, K., Weston, P., Kazumori, M., Okamoto, K., Zhu, Y., Liu, E. H., Collard, A., Bell, W., Migliorini, S., Chambon, P., Fourrié, N., Kim, M.-J., Köpken-Watts, C., and Schraff, C.: All-sky satellite data assimilation at operational weather forecasting centres, Q. J. Roy. Meteor. Soc., 144, 1191–1217,, 2018. a

Golub, G. H. and Ye, Q.: Inexact Preconditioned Conjugate Gradient Method with Inner-Outer Iteration, SIAM J. Sci. Comput., 21, 1305–1320,, 1999. a

Ha, S., Snyder, C., Skamarock, W. C., Anderson, J., and Collins, N.: Ensemble Kalman Filter Data Assimilation for the Model for Prediction Across Scales (MPAS), Mon. Weather Rev., 145, 4673–4692,, 2017. a

Hamill, T. M.: Hypothesis Tests for Evaluating Numerical Precipitation Forecasts, Weather Forecast., 14, 155–167,<0155:htfenp>;2, 1999. a

Hamill, T. M. and Snyder, C.: A Hybrid Ensemble Kalman Filter – 3D Variational Analysis Scheme, Mon. Weather Rev., 128, 2905–2919,<2905:ahekfv>;2, 2000. a

Hamill, T. M., Whitaker, J. S., and Snyder, C.: Distance-Dependent Filtering of Background Error Covariance Estimates in an Ensemble Kalman Filter, Mon. Weather Rev., 129, 2776–2790,<2776:ddfobe>;2, 2001. a

Hogan, R. J., Ferro, C. A. T., Jolliffe, I. T., and Stephenson, D. B.: Equitability Revisited: Why the “Equitable Threat Score” Is Not Equitable, Weather Forecast., 25, 710–726,, 2010. a

Holdaway, D., Vernieres, G., Wlasak, M., and King, S.: Status of Model Interfacing in the Joint Effort for Data Assimilation Integration (JEDI), JCSDA Quarterly Newsletter, 66, 15–24,, 2020. a

Honeyager, R., Herbener, S., Zhang, X., Shlyaeva, A., and Trémolet, Y.: Observations in the Joint Effort for Data Assimilation Integration (JEDI) – Unified Forward Operator (UFO) and Interface for Observation Data Access (IODA), JCSDA Quarterly Newsletter, 66, 24–33,, 2020. a

Hong, S.-Y. and Lim, J.-O. J.: The WRF Single-Moment 6-Class Microphysics Scheme (WSM6), J. Korean Meteor. Soc., 42, 129–151, 2006. a, b

Hong, S.-Y., Noh, Y., and Dudhia, J.: A New Vertical Diffusion Package with an Explicit Treatment of Entrainment Processes, Mon. Weather Rev., 134, 2318–2341,, 2006. a

Houtekamer, P. L. and Mitchell, H. L.: A Sequential Ensemble Kalman Filter for Atmospheric Data Assimilation, Mon. Weather Rev., 129, 123–137,<0123:asekff>;2, 2001. a

Iacono, M. J., Delamere, J. S., Mlawer, E. J., Shephard, M. W., Clough, S. A., and Collins, W. D.: Radiative forcing by long-lived greenhouse gases: Calculations with the AER radiative transfer models, J. Geophys. Res., 113, D13103,, 2008. a

Ingleby, B.: Global assimilation of air temperature, humidity, wind and pressure from surface stations, Q. J. Roy. Meteor. Soc., 141, 504–517,, 2014. a

Jiménez, P. A., Dudhia, J., González-Rouco, J. F., Navarro, J., Montávez, J. P., and García-Bustamante, E.: A Revised Scheme for the WRF Surface Layer Formulation, Mon. Weather Rev., 140, 898–918,, 2012. a

Joint Center for Satellite Data Assimilation and National Center for Atmospheric Research: JEDI-MPAS Data Assimilation System v1.0.0, Zenodo [code],, 2021. a

Kalnay, E., Kanamitsu, M., Kistler, R., Collins, W., Deaven, D., Gandin, L., Iredell, M., Saha, S., White, G., Woollen, J., Zhu, Y., Chelliah, M., Ebisuzaki, W., Higgins, W., Janowiak, J., Mo, K. C., Ropelewski, C., Wang, J., Leetmaa, A., Reynolds, R., Jenne, R. and Joseph, D.: The NCEP/NCAR 40-year reanalysis project, B. Am. Meteorol. Soc., 77, 437–472,<0437:TNYRP>2.0.CO;2, 1996. a

Kleist, D. T. and Ide, K.: An OSSE-Based Evaluation of Hybrid Variational – Ensemble Data Assimilation for the NCEP GFS. Part II: 4DEnVar and Hybrid Variants, Mon. Weather Rev., 143, 452–470,, 2015. a, b, c

Klemp, J. B.: A Terrain-Following Coordinate with Smoothed Coordinate Surfaces, Mon. Weather Rev., 139, 2163–2169,, 2011. a

Klemp, J. B., Skamarock, W. C., and Dudhia, J.: Conservative Split-Explicit Time Integration Methods for the Compressible Nonhydrostatic Equations, Mon. Weather Rev., 135, 2897–2913,, 2007. a

Liu, Q. and Weng, F.: Advanced Doubling – Adding Method for Radiative Transfer in Planetary Atmospheres, J. Atmos. Sci., 63, 3459–3465,, 2006. a

Liu, Z., Ban, J., Hong, J.-S., and Kuo, Y.-H.: Multi-resolution incremental 4D-Var for WRF: Implementation and application at convective scale, Q. J. Roy. Meteor. Soc., 146, 3661–3674,, 2020. a

Lorenc, A. C.: Analysis methods for numerical weather prediction, Q. J. Roy. Meteor. Soc., 112, 1177–1194,, 1986. a

Lorenc, A. C.: The potential of the ensemble Kalman filter for NWP – a comparison with 4D-Var, Q. J. Roy. Meteor. Soc., 129, 3183–3203,, 2003. a

McCormack, J. P., Eckermann, S. D., Siskind, D. E., and McGee, T. J.: CHEM2D-OPP: A new linearized gas-phase ozone photochemistry parameterization for high-altitude NWP and climate models, Atmos. Chem. Phys., 6, 4943–4972,, 2006. a, b

Ménétrier, B.: Normalized Interpolated Convolution from an Adaptive Subgrid documentation, (last access: 21 October 2022), 2020. a, b

Migliorini, S. and Candy, B.: All-sky satellite data assimilation of microwave temperature sounding channels at the Met Office, Q. J. Roy. Meteor. Soc., 145, 867–883,, 2019. a, b, c, d

National Centers For Environmental Prediction/National Weather Service/NOAA/U.S. Department Of Commerce: NCEP ADP Global Upper Air and Surface Weather Observations (PREPBUFR format),, 2008. a

National Centers For Environmental Prediction/National Weather Service/NOAA/U.S. Department Of Commerce: NCEP GDAS Satellite Data 2004-continuing,, 2009. a

National Centers For Environmental Prediction/National Weather Service/NOAA/U.S. Department Of Commerce: NCEP GFS 0.25 Degree Global Forecast Grids Historical Archive,, 2015. a

Rabier, F., Järvinen, H., Klinker, E., Mahfouf, J.-F., and Simmons, A.: The ECMWF operational implementation of four-dimensional variational assimilation. I: Experimental results with simplified physics, Q. J. Roy. Meteor. Soc., 126, 1143–1170,, 2007. a

Schwartz, C. S., Liu, Z., and Huang, X.-Y.: Sensitivity of Limited-Area Hybrid Variational-Ensemble Analyses and Forecasts to Ensemble Perturbation Resolution, Mon. Weather Rev., 143, 3454–3477,, 2015. a

Shao, H., Derber, J., Huang, X.-Y., Hu, M., Newman, K., Stark, D., Lueken, M., Zhou, C., Nance, L., Kuo, Y.-H., and Brown, B.: Bridging Research to Operations Transitions: Status and Plans of Community GSI, B. Am. Meteorol. Soc., 97, 1427–1440,, 2016. a, b

Skamarock, W. C., Klemp, J. B., Duda, M. G., Fowler, L. D., Park, S.-H., and Ringler, T. D.: A Multiscale Nonhydrostatic Atmospheric Model Using Centroidal Voronoi Tesselations and C-Grid Staggering, Mon. Weather Rev., 140, 3090–3105,, 2012. a, b

Skamarock, W. C., Duda, M. G., Ha, S., and Park, S.-H.: Limited-Area Atmospheric Modeling Using an Unstructured Mesh, Mon. Weather Rev., 146, 3445–3460,, 2018. a

Tarantola, A.: Inverse Problem Theory and Methods for Model Parameter Estimation, Society for Industrial and Applied Mathematics,, 2005. a

Tiedtke, M.: A Comprehensive Mass Flux Scheme for Cumulus Parameterization in Large-Scale Models, Mon. Weather Rev., 117, 1779–1800,<1779:acmfsf>;2, 1989. a

Trémolet, Y.: Joint Effort for Data Assimilation Integration (JEDI) Design and Structure, JCSDA Quarterly Newsletter, 66, 6–14,, 2020. a

Trémolet, Y. and Auligné, T.: The Joint Effort for Data Assimilation Integration (JEDI), JCSDA Quarterly Newsletter, 66, 1–5,, 2020.  a

Wang, X.: Incorporating Ensemble Covariance in the Gridpoint Statistical Interpolation Variational Minimization: A Mathematical Framework, Mon. Weather Rev., 138, 2990–2995,, 2010. a

Wang, X., Parrish, D., Kleist, D., and Whitaker, J.: GSI 3DVar-Based Ensemble – Variational Hybrid Data Assimilation for NCEP Global Forecast System: Single-Resolution Experiments, Mon. Weather Rev., 141, 4098–4117,, 2013. a, b

Wu, Y., Liu, Z., and Li, D.: Improving forecasts of a record-breaking rainstorm in Guangzhou by assimilating every 10-min AHI radiances with WRF 4DVAR, Atmos. Res., 239, 104912,, 2020. a

Xie, P., Joyce, R., Wu, S., Yoo, S.-H., Yarosh, Y., Sun, F., and Lin, R.: Reprocessed, Bias-Corrected CMORPH Global High-Resolution Precipitation Estimates from 1998, J. Hydrometeorol., 18, 1617–1641,, 2017. a

Xie, P., Joyce, R., Wu, S., Yoo, S.-H., Yarosh, Y., Sun, F., Lin, R., and NOAA CDR Program: NOAA Climate Data Record (CDR) of CPC Morphing Technique (CMORPH) High Resolution Global Precipitation Estimates, Version 1,, 2019. a

Xu, K.-M. and Randall, D. A.: A Semiempirical Cloudiness Parameterization for Use in Climate Models, J. Atmos. Sci., 53, 3084–3102,<3084:ascpfu>;2, 1996. a

Zhang, C. and Wang, Y.: Projected Future Changes of Tropical Cyclone Activity over the Western North and South Pacific in a 20-km-Mesh Regional Climate Model, J. Climate, 30, 5923–5941,, 2017. a

Zhou, X., Zhu, Y., Hou, D., Luo, Y., Peng, J., and Wobus, R.: Performance of the new NCEP Global Ensemble Forecast System in a parallel experiment, Weather Forecast., 32, 1989–2004,, 2017. a

Zhu, Y., Derber, J., Collard, A., Dee, D., Treadon, R., Gayno, G., and Jung, J. A.: Enhanced radiance bias correction in the National Centers for Environmental Prediction's Gridpoint Statistical Interpolation data assimilation system, Q. J. Roy. Meteor. Soc., 140, 1479–1492,, 2013. a

Zhu, Y., Liu, E., Mahajan, R., Thomas, C., Groff, D., Delst, P. V., Collard, A., Kleist, D., Treadon, R., and Derber, J. C.: All-Sky Microwave Radiance Assimilation in NCEP's GSI Analysis System, Mon. Weather Rev., 144, 4709–4735,, 2016. a, b, c, d, e, f

Zhu, Y., Gayno, G., Purser, R. J., Su, X., and Yang, R.: Expansion of the All-Sky Radiance Assimilation to ATMS at NCEP, Mon. Weather Rev., 147, 2603–2620,, 2019. a

Short summary
JEDI-MPAS 1.0.0, a new data assimilation (DA) system for the MPAS model, was publicly released for community use. This article describes JEDI-MPAS's implementation of the ensemble–variational DA technique and demonstrates its robustness and credible performance by incrementally adding three types of microwave radiances (clear-sky AMSU-A, all-sky AMSU-A, clear-sky MHS) to a non-radiance DA experiment. We intend to periodically release new and improved versions of JEDI-MPAS in upcoming years.