the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
Data assimilation for the Model for Prediction Across Scales – Atmosphere with the Joint Effort for Data assimilation Integration (JEDIMPAS 1.0.0): EnVar implementation and evaluation
Chris Snyder
Jonathan J. Guerrette
ByoungJoo Jung
Junmei Ban
Steven Vahl
Yannick Trémolet
Thomas Auligné
Benjamin Ménétrier
Anna Shlyaeva
Stephen Herbener
Emily Liu
Daniel Holdaway
Benjamin T. Johnson
On 24 September 2021, JEDIMPAS 1.0.0, a new data assimilation (DA) system for the Model Prediction Across Scales – Atmosphere (MPASA) 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, JEDIMPAS capabilities include threedimensional variational (3DVar) and ensemble–variational (EnVar) schemes as well as the ensemble of DA (EDA) technique. On the observation side, one advanced feature in JEDIMPAS is the full allsky approach for satellite radiance DA with the introduction of hydrometeor analysis variables. This paper describes the formulation and implementation of EnVar for JEDIMPAS. JEDIMPAS 1.0.0 is evaluated with monthlong cycling 3DEnVar experiments with a global 30–60 km dualresolution configuration. The robustness and credible performance of JEDIMPAS are demonstrated by establishing a benchmark nonradiance DA experiment, then incrementally adding microwave radiances from three sources: Advanced Microwave Sounding UnitA (AMSUA) temperature sounding channels in clearsky scenes, AMSUA window channels in allsky scenes, and Microwave Humidity Sounder (MHS) water vapor channels in clearsky scenes. JEDIMPAS 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 bestperforming JEDIMPAS 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 coldstart 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 JEDIMPAS's great potential for both research and operations.
 Article
(2597 KB) 
Supplement
(408 KB)  BibTeX
 EndNote
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 shortterm model forecast) as well as the associated error characterizations (Tarantola, 2005). 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, modelagnostic components, which can be interfaced to a variety of models such as those for atmosphere, ocean, land, and atmospheric chemistry. JEDI adopts objectoriented and generic programming and uses a mix of C++ and Fortran (Trémolet, 2020).
There are four main generic components in the JEDI framework. The ObjectOriented Prediction System (OOPS; https://github.com/JCSDA/oops, last access: 21 October 2022) defines the abstract representations of elements of several DA algorithms, including three and fourdimensional variational (3D/4DVar) (e.g., Lorenc, 1986; Rabier et al., 2007; Liu et al., 2020) and 3D/4D ensemble–variational (3D/4DEnVar) (e.g., Wang et al., 2013; Kleist and Ide, 2015) techniques with several minimization algorithms implemented. OOPS is mainly programmed in C++ using templates. The SystemAgnostic Background Error Representation (SABER; https://github.com/JCSDA/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étrier, 2020) from SABER, which estimates and applies background errorcovariancerelated operators, defined on an unstructured mesh. Two observationrelated components are the Unified Forward Operator (UFO; https://github.com/JCSDA/ufo, last access: 21 October 2022) and the Interface for Observation Data Access (IODA; https://github.com/JCSDA/ioda, last access: 21 October 2022) (Honeyager et al., 2020). UFO implements observation operators for various types of observations. In addition, UFO implements other observationspace functionality, including quality control, data thinning, and bias correction. IODA provides functionality for input, output, and inmemory data access of observations. In addition to those modelagnostic components of JEDI described above, there are modelspecific implementations of interfaces to those generic components that allow for the realization of modelspecific 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 (MPASA) (Skamarock et al., 2012) based upon DART (i.e., MPASDART). MPASDART 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 MPASA based upon the JEDI framework at a similar time in 2018. In the past more than 3 years, we have developed the JEDIMPAS DA system, which consists of MPAS interfaces to JEDI (https://github.com/JCSDA/mpasjedi, last access: 21 October 2022) and the four modelagnostic JEDI components. The firstrelease JEDIMPAS 1.0.0 was made publicly available in September 2021. JEDIMPAS capabilities include threedimensional 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, JEDIMPAS 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 JEDIMPAS is the full allsky approach for satellite radiance DA with the introduction of hydrometeor analysis variables.
This article describes MPASspecific aspects of JEDIMPAS, emphasizing its implementation of EnVar, and documents JEDIMPAS performance in global cycling experiments. The implementation and evaluation of JEDIMPAS's EDA and 3DVar algorithms will be presented in separate articles. The next section provides a brief introduction to MPASA and modifications made for JEDIMPAS. Section 3 describes the formulation and implementation of JEDIMPAS'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.
MPASA is a nonhydrostatic model discretized horizontally on an unstructured centroidal Voronoi mesh with Cgrid staggering of state variables (Skamarock et al., 2012). It is runtime configurable to use either quasiuniform or variableresolution meshes and for either global or regional applications (Skamarock et al., 2018). MPASA adopts a heightbased terrainfollowing vertical coordinate ζ following Klemp (2011), in which the coordinate surfaces are progressively smoothed to remove smallerscale terrain structure.
MPASA's prognostic equations employ flux variables that are coupled with density; that is, the prognostic variables in the discretized equations are
where ρ_{d} is the density of dry air, ${\stackrel{\mathrm{\u0303}}{\mathit{\rho}}}_{\mathrm{d}}={\mathit{\rho}}_{\mathrm{d}}(\partial \mathit{\zeta}/\partial z{)}^{\mathrm{1}}$, u is the normal component of horizontal velocity at cell edges, w is the vertical velocity at cell centers, q_{j} represents the mixing ratio of the respective water species (water vapor, cloud liquid water, cloud ice, rain, snow, graupel, hail) denoted by (${q}_{\mathrm{v}},{q}_{\mathrm{c}},{q}_{\mathrm{i}},{q}_{\mathrm{r}},{q}_{\mathrm{s}},{q}_{\mathrm{g}},{q}_{\mathrm{h}})$, and ${\mathit{\theta}}_{\mathrm{m}}={\mathit{\theta}}_{\mathrm{d}}[\mathrm{1}+({R}_{\mathrm{v}}/{R}_{\mathrm{d}}\left){q}_{\mathrm{v}}\right]$ is a modified moist potential temperature, with θ_{d} the dry potential temperature and R_{d} and R_{v} the gas constants for dry air and water vapor, respectively. The continuous prognostic equation set can symbolically be expressed as
where pressure p is obtained diagnostically from the equation of state, and ${\mathit{\rho}}_{\mathrm{m}}={\mathit{\rho}}_{\mathrm{d}}(\mathrm{1}+{q}_{\mathrm{v}}+{q}_{\mathrm{c}}+{q}_{\mathrm{i}}+{q}_{\mathrm{r}}+{q}_{\mathrm{s}}+{q}_{\mathrm{g}}+{q}_{\mathrm{h}})$ is the density of moist air.
The fully compressible nonhydrostatic equations are solved using the timesplit 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 righthandside terms in the vertical momentum equation, the governing equations are rewritten in terms of perturbation thermodynamic variables $({\mathit{\rho}}_{\mathrm{d}}^{\prime},\phantom{\rule{0.125em}{0ex}}\phantom{\rule{0.125em}{0ex}}{\mathrm{\Theta}}_{\mathrm{m}}^{\prime},\phantom{\rule{0.125em}{0ex}}{p}^{\prime})$ relative to a hydrostatically balanced reference state that is only a function of the geometric height z.
JEDIMPAS not only uses MPAS as the forecast model for cycling data assimilation, but also leverages and extends certain capabilities within MPAS. Crucially, JEDIMPAS 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. JEDIMPAS 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, JEDIMPAS modifies MPAS routines for reading files by splitting static and timevarying fields into separate files and reducing the number of timevarying fields retained (relative to a restart file). These enhancements build on a twostream input approach for MPAS data assimilation developed by Soyoung Ha with William Skamarock, which was expanded and implemented for JEDIMPAS cycling experiments. All changes necessary for JEDIMPAS 1.0.0 are in a modified version of MPASA (https://github.com/JCSDA/MPASModel, last access: 21 October 2022) that is based upon the latest MPASA version 7.1.
The methodology of JEDIMPAS DA in updating the MPAS model variables using the EnVar technique is given in the next section.
3.1 Incremental variational analysis
The nonquadratic form of the cost function for a variational analysis is expressed as
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 fourdimensional analysis), the subscript “b” indicates a background field which is typically from a shortterm 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 firstguess variable x_{g} and letting $\mathit{\delta}\mathit{x}=\mathit{x}{\mathit{x}}_{\mathrm{g}}$ and $\mathit{\delta}{\mathit{x}}_{\mathrm{g}}={\mathit{x}}_{\mathrm{b}}{\mathit{x}}_{\mathrm{g}}$, Eq. (3) can be rewritten in the incremental form (Courtier et al., 1994) as
where $\mathit{d}=\mathit{y}h\left({\mathit{x}}_{\mathrm{g}}\right)$ is the firstguess departure to observations (also called “the innovation vector”), and H is the linearized version of h in the vicinity of x_{g}. The gradient of the cost function with respect to δx
shall vanish at the value δx_{a} that minimizes the cost function, which leads to
Equation (6) is a linear algebra system Aδx_{a}=b that can be solved iteratively with a minimization algorithm. The analysis of the full field can be found with ${\mathit{x}}_{a}={\mathit{x}}_{\mathrm{g}}+\mathit{\delta}{\mathit{x}}_{a}$. It is common practice that the minimization iterations (also referred to as “inner loop”) start from x_{g}=x_{b} and δx_{g}=0. Multiple outer loops can be used with the firstguess x_{g} (and thus δx_{g}) updated from the analysis of the previous outer loop. In each outer loop, d is recalculated, and H is relinearized using the updated x_{g}.
Several minimization algorithms are available in OOPS, and the most extensively tested one in JEDIMPAS is the Derber–Rosati Inexact Preconditioned Conjugate Gradient (DRIPCG) algorithm (Golub and Ye, 1999; Derber and Rosati, 1989), which is documented in https://jointcenterforsatellitedataassimilationjedidocs.readthedocshosted.com/en/latest/inside/jedicomponents/oops/algorithmic_details/solvers.html (last access: 21 October 2022). It is important to note that the DRIPCG algorithm introduces an additional initial vector ${\widehat{\mathit{\delta}\mathit{x}}}_{\mathrm{g}}={\mathbf{B}}^{\mathrm{1}}\mathit{\delta}{\mathit{x}}_{\mathrm{g}}$ at the beginning of the inner loop so that the application of B^{−1} can be avoided if the initial guess is taken as δx_{g}=0 (the case of the first outer loop) or if both δx_{g} and $\widehat{\mathit{\delta}{\mathit{x}}_{\mathrm{g}}}$ 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 hybridEnVar 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 B_{s} and an ensemblebased covariance B_{e}; i.e.,
where β_{s} and β_{e} are the scalar weights with ${\mathit{\beta}}_{\mathrm{s}}+{\mathit{\beta}}_{\mathrm{e}}=\mathrm{1}$, and L∘B_{e} denotes the Schur product (element by element) of the localization matrix L and the sample covariance matrix B_{e} of the ensemble. Note that L is a correlation matrix, with diagonal elements being 1 and offdiagonal 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 B_{e} by reducing spurious correlations between state variables separated by large distances and by increasing the rank of the localized sample covariance relative to B_{e} (Hamill et al., 2001; Houtekamer and Mitchell, 2001).
The modeling of the static covariance B_{s} for JEDIMPAS 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, Br_{k} is calculated, where ${\mathit{r}}_{k}=\mathit{b}\mathbf{A}\mathit{\delta}{\mathit{x}}_{k}$ is the residual vector at the kth iteration, and A and b are defined by Eq. (6). In the case of pure EnVar with $\mathbf{B}=\mathbf{L}\circ {\mathbf{B}}_{\mathrm{e}}$, Br_{k} is actually evaluated in OOPS in the form
where N_{e} is the ensemble size, and ${\mathit{x}}_{\mathrm{m}}^{\prime}$ is $({N}_{\mathrm{e}}\mathrm{1}{)}^{\mathrm{1}/\mathrm{2}}$ times the deviation of the mth member from the ensemble mean. That is, (L∘B_{e})r is given by taking the Schur product of each ${\mathit{x}}_{\mathrm{m}}^{\prime}$ with r, passing the result through L, taking the Schur product of the result again with ${\mathit{x}}_{\mathrm{m}}^{\prime}$, and finally summing over m.
To see that Eq. (8) is correct, write $(\mathbf{L}\circ {\mathbf{B}}_{\mathrm{e}}{)}_{ij}={L}_{ij}{\sum}_{m=\mathrm{1}}^{{N}_{\mathrm{e}}}{x}_{im}^{\prime}{x}_{jm}^{\prime}$ and so the ith element of the vector (L∘B_{e})r
where N_{x} 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 hybridEnVar cost function used in other DA systems (e.g., Lorenc, 2003; Wang, 2010; 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 N_{x}×N_{x}, inmemory storage of the full matrix L is prohibitive, and the multiplication of L with a vector of size N_{x} can be computationally expensive for a highresolution setting of the model. JEDIMPAS 1.0.0 adopts the Normalized Interpolated Convolution from an Adaptive Subgrid (NICAS) method to model L such that
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 lengthscale 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 JEDIMPAS (also built on an unstructured mesh). This is a key feature that enables JEDIMPAS to work seamlessly for global or regional domains with a quasiuniform or variableresolution MPAS mesh. More details about BUMP and NICAS are given by Ménétrier (2020) (https://github.com/benjaminmenetrier/nicas_doc/blob/master/nicas_doc.pdf, last access: 21 October 2022).
To further reduce the computational cost, a dualresolution capability is implemented in JEDIMPAS's EnVar, in which the background and analysis (thus firstguess 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, JEDIMPAS 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). JEDIMPAS uses as analysis variables the cellcenter values of temperature (T), specific humidity (q), zonal and meridional components of horizontal velocity (v), and surface pressure (p_{s}). 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/precipitationaffected observations. We chose these analysis variables because of our experience with modeling backgrounderror 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, p_{s}) by iteratively solving Eq. (6), several additional steps are followed to arrive at a fullfield analysis for variables closer to the MPAS prognostic variables:

Transform the cellcenter 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 fullfield analysis of u by adding increments to the first guess.

Get the fullfield analyses of T, q, and p_{s} (and optionally hydrometeors) by adding analysis increments to the first guess. Calculate the water vapor mixing ratio analysis with ${q}_{\mathrm{v}}=q/(\mathrm{1}q)$.

Diagnose the fullfield analyses of 3D pressure p and dryair density ρ_{d} using T, q_{v}, and p_{s} by applying the hypsometric equation $p\left({z}_{\mathrm{2}}\right)=p\left({z}_{\mathrm{1}}\right)\mathrm{exp}(\frac{g}{{R}_{\mathrm{d}}}{\int}_{{z}_{\mathrm{1}}}^{{z}_{\mathrm{2}}}\frac{\mathrm{d}z}{{T}_{\mathrm{v}}})$ and the equation of state of moist air ${\mathit{\rho}}_{\mathrm{d}}=p/\left[\right(\mathrm{1}+{q}_{\mathrm{v}}+\mathrm{\dots}\left){R}_{\mathrm{d}}{T}_{\mathrm{v}}\right]$ with virtual temperature ${T}_{\mathrm{v}}=(\mathrm{1}+\mathrm{0.608}{q}_{\mathrm{v}})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 p_{s}.

Update dry potential temperature ${\mathit{\theta}}_{\mathrm{d}}=T(\frac{{p}_{\mathrm{0}}}{p}{)}^{{R}_{\mathrm{d}}/{c}_{p}}$, with c_{p} the specific heat of dry air at constant pressure and p_{0}=1000 hPa a reference pressure.
Sensitivity experiments showed that the hydrostatic balance constraint in Step 3 is important for JEDIMPAS to cycle stably by allowing surface pressure observations to constrain upperair fields.
3.4 Observation handling and allsky radiance DA
In JEDIMPAS, observation I/O and inmemory storage are handled through IODA, which is one of modelagnostic components of JEDI and reads specific netCDF4 (based on hdf5) format observation input files and outputs DA feedback files. UFO, another modelagnostic 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 socalled gsincdiag format are converted into the IODAv2 format and then assimilated into JEDIMPAS. Therefore, mainly two QC filters are applied in JEDIMPAS experiments (see Sect. 4). One uses QC flags generated by GSI, and another is the “background check” filter, which discards observations with firstguess departures larger than 3 times the observation error standard deviation. For surface pressure (p_{s}) 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 p_{s} 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 allsky microwave radiances (e.g., Geer et al., 2017, 2018; Zhu et al., 2016, 2019; Migliorini and Candy, 2019). JEDIMPAS allows the assimilation of allsky radiance data with the introduction of hydrometeor analysis variables, cloudy radiance observation operator via the Community Radiative Transfer Model (CRTM) (Liu and Weng, 2006), and situationdependent allsky observation error models (Geer and Bauer, 2011). While allsky radiance DA in JEDIMPAS can be applied to data from both microwave and infrared sensors, this study evaluates only microwave allsky radiance DA from Advanced Microwave Sounding UnitA (AMSUA) window channels 1–4 and 15 over water. The AMSUA allsky radiance error model follows Zhu et al. (2016), in which the observation error is a piecewiselinear ramp function (see Sect. 4) of the model–observationaveraged cloud liquid water path $({c}_{\mathrm{obs}}+{c}_{\mathrm{fg}})/\mathrm{2}$ retrieved from the AMSUA channel 1's and 2's brightness temperatures. Zhu et al. (2016) assimilated AMSUA window channels only for nonprecipitating 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 Lim, 2006) used in this study and all five hydrometeor variables included in the analysis variables, JEDIMPAS assimilates AMSUA window channel radiances in both nonprecipitating and precipitating cloudy situations with a full allsky 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.
We evaluate JEDIMPAS 3DEnVar's performance with monthlong (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 6hourly cycling dualresolution 3DEnVar experiments are conducted with a global quasiuniform 30 km mesh for highresolution analysis/background and a 60 km mesh for lowresolution 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 JEDIMPAS analysis. The firstguess 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 20member ensemble input on the 60 km mesh is generated from 6 h MPAS model forecasts initialized from the NCEP's 20member Global Ensemble Forecast System (GEFS) ensemble analysis (Zhou et al., 2017).
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 polarorbiting satellite sensors; and surface pressure (p_{s}). In addition to the observations used in conv, the clrama experiment assimilates clearsky AMSUA radiances (channels 5–9) from six satellites (NOAA15/18/19, MetopA/B, Aqua; omitting platforms and channels that had failed before the experimental period). The cldama experiment further adds overwater AMSUA window channels (channels 1–4 and 15) allsky radiances from five satellites. Note that Aqua AMSUA channels 1 and 2, which are needed for the allsky observation error model (see Sect. 3.4), are not available, and thus Aqua AMSUA 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 (NOAA18/19, MetopA/B) under clearsky conditions. Data thinning is applied to AMSUA and MHS radiance data with a 145 km thinning mesh. All four experiments use a ±3 h time window.
Allsky observation error statistics were computed separately for five AMSUA sensors using the observations minus the CRTMsimulated brightness temperatures from the 6 h forecast background fields (i.e., O minus B) of the clrama experiment over the monthlong period. Figure 1 gives an example of the allsky observation error model for NOAA19 AMSUA 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–observationaveraged cloud liquid water path $\stackrel{\mathrm{\u203e}}{c}=({c}_{\mathrm{obs}}+{c}_{\mathrm{fg}})/\mathrm{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 piecewiselinear ramp function of the blue curve in the form (Geer and Bauer, 2011)
which serves as the observation error model for this particular channel with ${\stackrel{\mathrm{\u203e}}{c}}_{\mathrm{clr}}=\mathrm{0.03}$, ${\stackrel{\mathrm{\u203e}}{c}}_{\mathrm{cld}}=\mathrm{0.24}$,
σ_{clr}=6.33, and σ_{cld}=19.24.
The four numbers for each of channels 1–4 and 15 of five AMSUA sensors are provided in the supplemental configuration file of JEDIMPAS (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, BUMPdiagnosed 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 MPASA model is configured with a 30 km quasiuniform 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.
(Tiedtke, 1989; Zhang and Wang, 2017)(Hong and Lim, 2006)(Chen and Dudhia, 2001)(Hong et al., 2006)(Jiménez et al., 2012)(Iacono et al., 2008)(Xu and Randall, 1996)(Choi and Hong, 2015)5.1 Shortterm forecasts
Here we assess the DA cycling performance of the four experiments by examining the firstguess departure (i.e., observation minus background or OMB) statistics, similar to the diagnostic commonly monitored at operational NWP centers (e.g., https://www.ecmwf.int/en/forecasts/charts/obstat/, 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.
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 clearsky AMSUA radiances leads to substantial error reduction in terms of both bias and RMSE; further adding allsky AMSUA windowchannels' 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 clearsky MHS water vapor radiances to allsky AMSUA 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 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
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 (Hamill, 1999), 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 clearsky 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 allsky AMSUA and clearsky MHS radiances. The impact on virtual temperature (T_{v}) is mostly from AMSUA clearsky temperature sounding channels, and the additional benefit on T_{v} from allsky AMSUA and clearsky MHS radiances is relatively small (Fig. 3a). The reason we show T_{v} here instead of T is that the majority of radiosonde temperature observations in GSI's ncdiag files are expressed in terms of T_{v}. 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 AMSUA allsky window channel radiances is confined to 200–800 hPa. Additional positive impact on winds between 300 and 650 hPa is achieved by adding MHS clearsky 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.
Among the four experiments, both cldama and clrmhs assimilate AMSUA allsky 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 AMSUA allsky radiances. That is indeed the case in Fig. 4; clrmhs achieves statistically significant improvement up to 2.5 % at tropical latitudes for AMSUA channel 1 (Fig. 4a), while positive impacts from MHS radiances are smaller and less significant for other AMSUA window channels (Fig. 4b–d). This is because AMSUA channel 1 (23.8 GHz) is more sensitive to water vapor and clouds than other AMSUA window channels.
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 hybrid4DEnVar DA technique via the GSI system (Kleist and Ide, 2015) 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 (T_{2}), 10 m zonal wind (u_{10}), and p_{s}, which is caused by sparse coverage of nonradiance data over the Southern Hemisphere. The impacts of AMSUA clearsky and allsky radiances and MHS clearsky radiances are consistent with those shown in upperair 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 AMSUA allsky radiances and >20 % improvement at model level 30 (∼8.8 km or ∼330 hPa) when adding both AMSUA allsky radiances and MHS clearsky 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 MPASA model currently lacks processes for representing O_{3} chemistry and that a climatological O_{3} field is used, leading to inaccurate radiation calculations attributed to O_{3}. Work is underway to add a parameterized O_{3} chemistry scheme (McCormack et al., 2006) and include O_{3} as a prognostic variable in MPASA.
Figure 7 is similar to Fig. 6 but depicts the latitudinal variation of relative RMSE reduction for several variables. For p_{s}, 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 observationspace and modelspace evaluations of shortterm forecasts demonstrate that the JEDIMPAS 3DEnVar scheme behaves well for the assimilation of clearsky and allsky 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 spinup period for each experiment to converge from the GFS climatology to JEDIMPAS'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.
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 upperair 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 AMSUA allsky radiances produces more accurate forecasts for all variables compared to the AMSUA clearsky radiance experiment, and extra MHS clearsky 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 upperair 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 AMSUA window channels' and the three MHS channels' large sensitivity to water vapor (and clouds for AMSUA allsky radiances). Large impacts are also seen on Q2m (Fig. 8b), though to a lesser extent than is seen for upperair 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 p_{s} and the zonal wind are comparable to each other, gradually decreasing from an initial ∼5 % rms reduction for the day 1 forecast.
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 shortterm 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 10 more clearly depicts the relative impact of clrmhs with respect to cldama for four upperair 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 shortterm 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
JEDIMPAS also has the capability of simulating and assimilating clearsky 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 JEDIMPAS is used to evaluate 1–10 d forecasts using independent (i.e., not assimilated) AHI fulldisk radiance data centered at 140.7^{∘} E. Using the JEDIMPAS 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 JEDIMPAS, 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 superobservations averaged over 15×15 pixels. The cloud detection and superobbing are carried out offline using WRFDA. These superobservations better match the model's 30 km mesh, reducing representativeness errors. A byproduct of cloud detection on the raw AHI pixels and superobbing is the cloud fraction for each superobservation, which allows for statistical evaluation separately over clear, partly cloudy, and overcast observations.
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 middletoupper 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.
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 $\mathrm{1}/\mathrm{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 AMSUA temperaturesounding radiance DA when compared to the nonradiance DA experiment. A further 5 %–12 % improvement is achieved by adding AMSUA allsky radiances and MHS radiances for up to day 6 (Fig. 13b). From precipitation forecast verification here and other aspects of forecast verification shown earlier, JEDIMPAS's AMSUA allsky DA exhibits a larger positive impact than that at operational centers (e.g., Zhu et al., 2016; Migliorini and Candy, 2019). This might be partly related to the use of the full allsky approach in JEDIMPAS, whereas Zhu et al. (2016) and Migliorini and Candy (2019) assimilated AMSUA allsky radiances only under nonprecipitating conditions. Also, obtaining additional benefit above a welltuned 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.
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 JEDIMPAS experiments together with the coldstart 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 AMSUA temperature sounding channels, gaining a half day of skill (from 4.5 to 5 d) with respect to the nonradiance experiment for 0.8 ACC level. Further ACC improvement by adding AMSUA allsky window channels and MHS clearsky WV channels is visible but much smaller than that gained by AMSUA clearsky temperaturesensitive channels. This could be understandable with the tight link between temperature and geopotential height. The coldstart forecasts from GFS analyses reach 5.5 d for 0.8 ACC level and outperforms JEDIMPAS's three radiance DA experiments by half a day. Given that the operational GFS analysis is produced with a more advanced hybrid4DEnVar algorithm (Kleist and Ide, 2015) at a higher resolution and with ensemble covariances formed by 80member ensemble input, plus the assimilation of many more satellite observations, this halfday gap on ACC is not surprising with a less advanced pure 3DEnVar algorithm, a smaller ensemble size, and fewer observations assimilated in JEDIMPAS experiments.
A new DA system for the MPASA model has been developed based upon the JEDI software framework and was publicly released as JEDIMPAS 1.0.0 for community use. This article describes the EnVar's formulation and technical implementation of JEDIMPAS and evaluates its robustness and performance through four global cycling 3DEnVar experiments. JEDIMPAS is robust, with stable cycling over a monthlong 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 JEDIMPAS is allsky radiance DA capability along with the introduction of hydrometeor analysis variables, which is first applied to the assimilation of allsky radiances from AMSUA's five window channels. Large impacts are achieved from AMSUA allsky 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 halfday 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. JEDIMPAS 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. JEDIMPAS in regional mode and with a variable resolution will also be evaluated with higher resolutions.
The general JEDI software framework and more specific JEDIMPAS DA system are still undergoing active development since the JEDIMPAS 1.0.0 release. One new development specific to JEDIMPAS 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 upperair variables but is better suited for the nonhydrostatic MPASA 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 O_{3} chemistry scheme (McCormack et al., 2006) and include O_{3} as a prognostic variable in MPASA, which should reduce temperature bias in the stratosphere and also allow for the assimilation of satellitebased O_{3} observations. Improving stratosphere processes will allow the use of a high model top, which will enable the assimilation of higherpeaking radiance data. Allsky radiance DA for both microwave and infrared sensors continues to be refined with JEDIMPAS. 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 Uppala, 2009; Zhu et al., 2013), which is undergoing active testing and evaluation with JEDIMPAS 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 NCEPbufr 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 JEDIMPAS in upcoming years.
JEDIMPAS 1.0.0 is publicly released on GitHub, accessible from https://www.jcsda.org/jedimpas (last access: 21 October 2022). It is also available on Zenodo at https://doi.org/10.5281/zenodo.6784274 (Joint Center for Satellite Data Assimilation and National Center for Atmospheric Research, 2021). Global Forecast System analysis data are downloaded from NCAR Research Data Archive https://rda.ucar.edu/datasets/ds084.1/ (last access: 21 October 2022; National Centers For Environmental Prediction/National Weather Service/NOAA/U.S. Department Of Commerce, 2015). Global Ensemble Forecast System ensemble analysis data are downloaded from https://www.ncei.noaa.gov/products/weatherclimatemodels/globalensembleforecast (last access: 21 October 2022). Conventional and satellite observations assimilated are downloaded from https://rda.ucar.edu/datasets/ds337.0/ (last access: 21 October 2022; National Centers For Environmental Prediction/National Weather Service/NOAA/U.S. Department Of Commerce, 2008) and https://rda.ucar.edu/datasets/ds735.0/ (last access: 21 October 2022; National Centers For Environmental Prediction/National Weather Service/NOAA/U.S. Department Of Commerce, 2009). CMORPH precipitation product and AHI radiance data used in the forecast verification are downloaded from https://www.ncei.noaa.gov/products/climatedatarecords/precipitationcmorph/ (last access: 21 October 2022; Xie et al., 2019) and https://registry.opendata.aws/noaahimawari/ (last access: 21 October 2022), respectively.
The supplement related to this article is available online at: https://doi.org/10.5194/gmd1578592022supplement.
ZL designed and conducted all experiments and wrote the paper. All coauthors contributed to the development of JEDIMPAS and the article.
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.
This research has been supported by the United States Air Force (grant no. NA21OAR4310383).
This paper was edited by PoLun 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, https://doi.org/10.1175/2009bams2618.1, 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, https://doi.org/10.1002/qj.56, 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, https://doi.org/10.1175/bamsd1100167.1, 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 CommunitySupported Forecast Verification, B. Am. Meteorol. Soc., 102, E782–E807, https://doi.org/10.1175/bamsd190093.1, 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, https://doi.org/10.1175/15200493(2001)129<0587:caalsh>2.0.co;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, https://doi.org/10.1002/2015jd024230, 2015. a
Courtier, P., Thépaut, J.N., and Hollingsworth, A.: A strategy for operational implementation of 4DVar, using an incremental approach, Q. J. Roy. Meteor. Soc., 120, 1367–1387, https://doi.org/10.1002/qj.49712051912, 1994. a
Dee, D. P. and Uppala, S.: Variational bias correction of satellite radiance data in the ERAInterim reanalysis, Q. J. Roy. Meteor. Soc., 135, 1830–1841, https://doi.org/10.1002/qj.493, 2009. a
Derber, J. and Rosati, A.: A Global Oceanic Data Assimilation System, J. Phys. Oceanogr., 19, 1333–1347, https://doi.org/10.1175/15200485(1989)019<1333:agodas>2.0.co;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, https://doi.org/10.1002/qj.49712555417, 1999. a
Geer, A. J. and Bauer, P.: Observation errors in allsky data assimilation, Q. J. Roy. Meteor. Soc., 137, 2024–2037, https://doi.org/10.1002/qj.830, 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, https://doi.org/10.1002/qj.3172, 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öpkenWatts, C., and Schraff, C.: Allsky satellite data assimilation at operational weather forecasting centres, Q. J. Roy. Meteor. Soc., 144, 1191–1217, https://doi.org/10.1002/qj.3202, 2018. a
Golub, G. H. and Ye, Q.: Inexact Preconditioned Conjugate Gradient Method with InnerOuter Iteration, SIAM J. Sci. Comput., 21, 1305–1320, https://doi.org/10.1137/s1064827597323415, 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, https://doi.org/10.1175/mwrd170145.1, 2017. a
Hamill, T. M.: Hypothesis Tests for Evaluating Numerical Precipitation Forecasts, Weather Forecast., 14, 155–167, https://doi.org/10.1175/15200434(1999)014<0155:htfenp>2.0.co;2, 1999. a
Hamill, T. M. and Snyder, C.: A Hybrid Ensemble Kalman Filter – 3D Variational Analysis Scheme, Mon. Weather Rev., 128, 2905–2919, https://doi.org/10.1175/15200493(2000)128<2905:ahekfv>2.0.co;2, 2000. a
Hamill, T. M., Whitaker, J. S., and Snyder, C.: DistanceDependent Filtering of Background Error Covariance Estimates in an Ensemble Kalman Filter, Mon. Weather Rev., 129, 2776–2790, https://doi.org/10.1175/15200493(2001)129<2776:ddfobe>2.0.co;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, https://doi.org/10.1175/2009waf2222350.1, 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, https://doi.org/10.25923/RB190Q26, 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, https://doi.org/10.25923/RB190Q26, 2020. a
Hong, S.Y. and Lim, J.O. J.: The WRF SingleMoment 6Class 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, https://doi.org/10.1175/mwr3199.1, 2006. a
Houtekamer, P. L. and Mitchell, H. L.: A Sequential Ensemble Kalman Filter for Atmospheric Data Assimilation, Mon. Weather Rev., 129, 123–137, https://doi.org/10.1175/15200493(2001)129<0123:asekff>2.0.co;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 longlived greenhouse gases: Calculations with the AER radiative transfer models, J. Geophys. Res., 113, D13103, https://doi.org/10.1029/2008jd009944, 2008. a
Ingleby, B.: Global assimilation of air temperature, humidity, wind and pressure from surface stations, Q. J. Roy. Meteor. Soc., 141, 504–517, https://doi.org/10.1002/qj.2372, 2014. a
Jiménez, P. A., Dudhia, J., GonzálezRouco, J. F., Navarro, J., Montávez, J. P., and GarcíaBustamante, E.: A Revised Scheme for the WRF Surface Layer Formulation, Mon. Weather Rev., 140, 898–918, https://doi.org/10.1175/mwrd1100056.1, 2012. a
Joint Center for Satellite Data Assimilation and National Center for Atmospheric Research: JEDIMPAS Data Assimilation System v1.0.0, Zenodo [code], https://doi.org/10.5281/ZENODO.6784274, 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 40year reanalysis project, B. Am. Meteorol. Soc., 77, 437–472, https://doi.org/10.1175/15200477(1996)077<0437:TNYRP>2.0.CO;2, 1996. a
Kleist, D. T. and Ide, K.: An OSSEBased Evaluation of Hybrid Variational – Ensemble Data Assimilation for the NCEP GFS. Part II: 4DEnVar and Hybrid Variants, Mon. Weather Rev., 143, 452–470, https://doi.org/10.1175/mwrd1300350.1, 2015. a, b, c
Klemp, J. B.: A TerrainFollowing Coordinate with Smoothed Coordinate Surfaces, Mon. Weather Rev., 139, 2163–2169, https://doi.org/10.1175/mwrd1005046.1, 2011. a
Klemp, J. B., Skamarock, W. C., and Dudhia, J.: Conservative SplitExplicit Time Integration Methods for the Compressible Nonhydrostatic Equations, Mon. Weather Rev., 135, 2897–2913, https://doi.org/10.1175/mwr3440.1, 2007. a
Liu, Q. and Weng, F.: Advanced Doubling – Adding Method for Radiative Transfer in Planetary Atmospheres, J. Atmos. Sci., 63, 3459–3465, https://doi.org/10.1175/jas3808.1, 2006. a
Liu, Z., Ban, J., Hong, J.S., and Kuo, Y.H.: Multiresolution incremental 4DVar for WRF: Implementation and application at convective scale, Q. J. Roy. Meteor. Soc., 146, 3661–3674, https://doi.org/10.1002/qj.3865, 2020. a
Lorenc, A. C.: Analysis methods for numerical weather prediction, Q. J. Roy. Meteor. Soc., 112, 1177–1194, https://doi.org/10.1002/qj.49711247414, 1986. a
Lorenc, A. C.: The potential of the ensemble Kalman filter for NWP – a comparison with 4DVar, Q. J. Roy. Meteor. Soc., 129, 3183–3203, https://doi.org/10.1256/qj.02.132, 2003. a
McCormack, J. P., Eckermann, S. D., Siskind, D. E., and McGee, T. J.: CHEM2DOPP: A new linearized gasphase ozone photochemistry parameterization for highaltitude NWP and climate models, Atmos. Chem. Phys., 6, 4943–4972, https://doi.org/10.5194/acp649432006, 2006. a, b
Ménétrier, B.: Normalized Interpolated Convolution from an Adaptive Subgrid documentation, https://github.com/benjaminmenetrier/nicas_doc/blob/master/nicas_doc.pdf (last access: 21 October 2022), 2020. a, b
Migliorini, S. and Candy, B.: Allsky satellite data assimilation of microwave temperature sounding channels at the Met Office, Q. J. Roy. Meteor. Soc., 145, 867–883, https://doi.org/10.1002/qj.3470, 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), https://doi.org/10.5065/Z83FN512, 2008. a
National Centers For Environmental Prediction/National Weather Service/NOAA/U.S. Department Of Commerce: NCEP GDAS Satellite Data 2004continuing, https://doi.org/10.5065/DWYZQ852, 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, https://doi.org/10.5065/D65D8PWK, 2015. a
Rabier, F., Järvinen, H., Klinker, E., Mahfouf, J.F., and Simmons, A.: The ECMWF operational implementation of fourdimensional variational assimilation. I: Experimental results with simplified physics, Q. J. Roy. Meteor. Soc., 126, 1143–1170, https://doi.org/10.1002/qj.49712656415, 2007. a
Schwartz, C. S., Liu, Z., and Huang, X.Y.: Sensitivity of LimitedArea Hybrid VariationalEnsemble Analyses and Forecasts to Ensemble Perturbation Resolution, Mon. Weather Rev., 143, 3454–3477, https://doi.org/10.1175/mwrd1400259.1, 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, https://doi.org/10.1175/bamsd1300245.1, 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 CGrid Staggering, Mon. Weather Rev., 140, 3090–3105, https://doi.org/10.1175/mwrd1100215.1, 2012. a, b
Skamarock, W. C., Duda, M. G., Ha, S., and Park, S.H.: LimitedArea Atmospheric Modeling Using an Unstructured Mesh, Mon. Weather Rev., 146, 3445–3460, https://doi.org/10.1175/mwrd180155.1, 2018. a
Tarantola, A.: Inverse Problem Theory and Methods for Model Parameter Estimation, Society for Industrial and Applied Mathematics, https://doi.org/10.1137/1.9780898717921, 2005. a
Tiedtke, M.: A Comprehensive Mass Flux Scheme for Cumulus Parameterization in LargeScale Models, Mon. Weather Rev., 117, 1779–1800, https://doi.org/10.1175/15200493(1989)117<1779:acmfsf>2.0.co;2, 1989. a
Trémolet, Y.: Joint Effort for Data Assimilation Integration (JEDI) Design and Structure, JCSDA Quarterly Newsletter, 66, 6–14, https://doi.org/10.25923/RB190Q26, 2020. a
Trémolet, Y. and Auligné, T.: The Joint Effort for Data Assimilation Integration (JEDI), JCSDA Quarterly Newsletter, 66, 1–5, https://doi.org/10.25923/RB190Q26, 2020. a
Wang, X.: Incorporating Ensemble Covariance in the Gridpoint Statistical Interpolation Variational Minimization: A Mathematical Framework, Mon. Weather Rev., 138, 2990–2995, https://doi.org/10.1175/2010mwr3245.1, 2010. a
Wang, X., Parrish, D., Kleist, D., and Whitaker, J.: GSI 3DVarBased Ensemble – Variational Hybrid Data Assimilation for NCEP Global Forecast System: SingleResolution Experiments, Mon. Weather Rev., 141, 4098–4117, https://doi.org/10.1175/mwrd1200141.1, 2013. a, b
Wu, Y., Liu, Z., and Li, D.: Improving forecasts of a recordbreaking rainstorm in Guangzhou by assimilating every 10min AHI radiances with WRF 4DVAR, Atmos. Res., 239, 104912, https://doi.org/10.1016/j.atmosres.2020.104912, 2020. a
Xie, P., Joyce, R., Wu, S., Yoo, S.H., Yarosh, Y., Sun, F., and Lin, R.: Reprocessed, BiasCorrected CMORPH Global HighResolution Precipitation Estimates from 1998, J. Hydrometeorol., 18, 1617–1641, https://doi.org/10.1175/jhmd160168.1, 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, https://doi.org/10.25921/W9VAQ159, 2019. a
Xu, K.M. and Randall, D. A.: A Semiempirical Cloudiness Parameterization for Use in Climate Models, J. Atmos. Sci., 53, 3084–3102, https://doi.org/10.1175/15200469(1996)053<3084:ascpfu>2.0.co;2, 1996. a
Zhang, C. and Wang, Y.: Projected Future Changes of Tropical Cyclone Activity over the Western North and South Pacific in a 20kmMesh Regional Climate Model, J. Climate, 30, 5923–5941, https://doi.org/10.1175/jclid160597.1, 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, https://doi.org/10.1175/WAFD170023.1, 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, https://doi.org/10.1002/qj.2233, 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.: AllSky Microwave Radiance Assimilation in NCEP's GSI Analysis System, Mon. Weather Rev., 144, 4709–4735, https://doi.org/10.1175/mwrd150445.1, 2016. a, b, c, d, e, f
Zhu, Y., Gayno, G., Purser, R. J., Su, X., and Yang, R.: Expansion of the AllSky Radiance Assimilation to ATMS at NCEP, Mon. Weather Rev., 147, 2603–2620, https://doi.org/10.1175/mwrd180228.1, 2019. a