## Abstract

Satellite retrievals strive to exploit the information contained in thousands of channels provided by hyperspectral sensors and show promise in providing a gain in computational efficiency over current radiance assimilation methods by transferring computationally expensive radiative transfer calculations to retrieval providers. This paper describes the implementation of a new approach based on the transformation proposed in 2008 by Migliorini et al., which reduces the impact of the a priori information in the retrievals and generates transformed retrievals (TRs) whose assimilation does not require knowledge of the hyperspectral instruments characteristics. Significantly, the results confirm both the viability of Migliorini’s approach and the possibility of assimilating data from different hyperspectral satellite sensors regardless of the instrument characteristics. The Weather Research and Forecasting (WRF) Model’s Data Assimilation (WRFDA) 3-h cycling system was tested over the central North Pacific Ocean, and the results show that the assimilation of TRs has a greater impact in the characterization of the water vapor distribution than on the temperature field. These results are consistent with the knowledge that temperature field is well constrained by the initial and boundary conditions of the Global Forecast System (GFS), whereas the water vapor distribution is less well constrained in the GFS. While some preliminary results on the comparison between the assimilation with and without TRs in the forecasting system are presented in this paper, additional work remains to explore the impact of the new assimilation approach on forecasts and will be provided in a follow-up publication.

## 1. Introduction

The aim of a data assimilation system is to create an analysis, that is, the best estimate of the atmospheric state, by blending all the available observations with a numerical model short-range forecast, accounting for the uncertainties of both the observations and the background forecast. In this regard, satellite observations of the infrared and microwave energy emitted by Earth represent a crucial source of information for the numerical weather prediction (NWP) models. In the early stage of satellite meteorology (1980s), retrievals produced operationally from broadband radiometers were used as part of the assimilation process in most NWP centers, and positive or moderate impact was found. Impact was more significant in the Southern Hemisphere than in the Northern Hemisphere because of the scarcity of conventional observations in the Southern Hemisphere (Kelly et al. 1978). However, with the improvements of the NWP systems, difficulties were encountered in using satellite soundings mainly because of the 1) limited vertical resolution of the retrieved profiles, 2) presence of a priori information in the retrievals not necessarily consistent with the background fields used in the assimilation, and 3) difficulties in properly defining the observation error covariance for the retrieved profiles. Therefore, during the late 1980s and 1990s, attention was focused on radiance assimilation in which the information in the observations is assimilated directly into the NWP systems (Andersson et al. 1994; Derber and Wu 1998) rather than retrieval assimilation.

A comprehensive overview of the theoretical aspects of infrared satellite data assimilation from retrieval to radiance assimilation can be found in Eyre (2007) and Eyre et al. (2019). Direct radiance assimilation is based on the availability of a forward model that uses the forecast background as a priori information to simulate the NWP model equivalents of the observed radiances, thus avoiding potential inconsistencies between the retrieval a priori and the assimilation background. But this comes at the cost of numerical efficiency and requires a high level of competence in dealing with different instruments. A renewed interest in assimilating satellite retrievals, especially for limited-area models, arose with the more recent launches of hyperspectral infrared sensors like the grating-based Atmospheric Infrared Sounder (AIRS; Chahine et al. 2006), the Infrared Atmospheric Sounder Interferometer (IASI; Lerner et al. 2002; Karagulian et al. 2010), the Cross-track Infrared Sounder (CrIS; Bloom 2001; Strow et al. 2013), the Hyperspectral Infrared Atmospheric Sounder (HIRAS; Zhang et al. 2019), and the Geostationary Interferometric Infrared Sounder (GIIRS) on geostationary satellite *Fengyun-4* (*FY-4*) of China (Yang et al. 2017).

Assimilation of satellite retrievals, as opposed to radiance assimilation, provides regional NWP centers running mesoscale models the opportunity to transfer complicated and computationally expensive radiative transfer calculations to the retrieval provider. In addition, with the retrievals being the highest form of compression of hyperspectral infrared observations, retrieval assimilation would allow, in clear-sky conditions, for the exploitation of the information content embedded in thousands of channels with minimal requirements for data transfer, data storage, and data manipulation. In contrast, operational direct radiance assimilation systems, even in major NWP centers, can currently handle only a subset on the order of hundreds of channels out of the thousands of channels provided by hyperspectral sensors at least until future operational use of linear combination of channels (principal components or reconstructed radiances) might overcome this limitation.

Migliorini et al. (2008) and Migliorini (2012) provide a framework to perform retrieval assimilation in an equivalent way (same information content) to direct radiance assimilation. The basic idea, which was previously introduced by Joiner and Da Silva (1998) and Rodgers (2000), is to apply a transformation to properly remove the a priori information from the retrievals before the assimilation is performed. The main goals of the transformation proposed by Migliorini are (i) to make the retrieval-forecast error cross-covariance terms negligible even when the retrieval is severely constrained by a priori information, (ii) to allow for offline radiative transfer calculations, and (iii) to reduce the number of assimilated quantities per observation to the number of effective degrees of freedom in the observation.

Inspired by the work of Migliorini, the authors of this paper generated transformed retrievals (TRs) using “Mirto” 1DVAR (Antonelli et al. 2017) and adapted a modified version of the Weather Research and Forecasting (WRF) Model 3DVAR Data Assimilation System (WRFDA; Barker et al. 2004, 2012) to ingest the TRs in an effort to exploit the information content of hyperspectral satellite data. The retrieval processor has been implemented at the Mauna Kea Weather Center (MKWC) of the University of Hawai‘i and has been running in quasi-operational mode since August 2013. The MKWC is a weather research and forecast facility funded by the astronomical observatories on Mauna Kea (Businger et al. 2002) and (Cherubini et al. 2011).

This paper describes the implementation and validation of the TR assimilation system within the WRFDA and provides additional insights regarding the retrieval process and its validation. More precisely, the paper describes how the retrievals are ingested by the WRFDA by means of a transformation that indeed reduces the influence of the a priori information used in the retrieval process, minimizes the number of variables that carry the information embedded in the original high dimensional observations, and decorrelates the errors associated with these variables. The main demonstrated advantages of such a system are (i) the computational costs of assimilating TRs are very low, which makes this methodology affordable by many; and (ii) the observation operator, unlike for direct radiance assimilation, comes with the TRs. Therefore, the assimilation module is not sensor dependent, and data from different hyperspectral sensors can be ingested by the same module.

This paper also anticipates how the designed system will become an integral part of a Rapid Update Cycle (RUC) system able to run every three hours and ingest, in near–real time, the TRs from IASI and CrIS data over the central Pacific, an area that is nearly devoid of conventional in situ observations. Thus, a local data assimilation system, which has to greatly rely on satellite radiance data, will be limited by computational costs and by the timeliness of radiance data availability in a format suitable for ingestion. Tests to replace the current MKWC modeling system with an operational implementation of this system are under way.

This paper is structured as follows: section 2 provides the working framework for the implementation of the designed system in terms of model descriptions and experimental cases under consideration; section 3 provides a description of the theoretical background, briefly reviewing the retrieval process and Migliorini’s approach, while focusing on the mathematical framework behind the assimilation module implemented within the WRFDA; validation strategies and results are presented in section 4; and conclusions and future work are outlined in section 5. Assessment of the impact of transformed retrieval assimilation on forecasting accuracy of the RUC system is left to a subsequent paper.

## 2. The underlying working framework

A coincidence of advantageous factors has made possible the implementation of the assimilation system presented in this paper: (i) the paucity of in situ conventional atmospheric observations in the north-central Pacific area; (ii) the availability of a direct broadcast system operated by the Space Science Engineering Center (SSEC) in collaboration with the National Weather Service and the University of Hawai‘i at Mānoa, which provides nearly real-time (about 15 min from satellite overpass time) high-resolution infrared data from the CrIS sensors on *Suomi NPP*/*NOAA-20* and the IASI sensors on *MetOp-A/B*; (iii) the availability of operational WRF forecasts at the MKWC, which provides the a priori knowledge of the atmospheric state needed by the retrieval system; and (iv) the availability of near-real-time TRs, generated by the Mirto inversion package/regional processor (Antonelli et al. 2017).

The MKWC routinely runs the WRF Model system (Klemp et al. (2007), http://www.wrf-model.org). Since Antonelli et al. (2017), the operational WRF modeling system has been updated in terms of both software and configuration, and the latest model configuration encompasses 2 two-way nested domains, with horizontal grids spacing of 4500 and 900 m, centered over the north-central Pacific area and the island of Hawaii, respectively. A one-way 300 m nested domain centered over the summit of Mauna Kea is also part of the modeling configuration and it is dedicated to provide high-spatial-resolution predictions of weather and observing conditions for the community of astronomers that operate there (Fig. 1).

Fifty-one levels in the vertical direction, *p*_{wrf}(*i* = 1, …, 51), are used for the 4500- and 900-m domains. The vertical spacing is on the order of tens of meters for the levels nearest the ground and gradually increases with height. The vertical resolution is doubled for the 300-m domain. The model top is fixed at 40 hPa, which corresponds to a height of ~22 km above the ground level. Although a higher model top would be desirable for satellite data assimilation, the MKWC chosen model top is retained because it minimizes the number of vertical levels used (and therefore computational costs) while maintaining a vertical resolution of at most ~600 m for domain 1 and 2 and ~300 m for domain 3; a constraint required by the broader MKWC applications, geared to improve optical turbulence forecasting. The WRF modeling system can run in cold-start mode, which uses as initial and boundary conditions the National Centers for Environmental Prediction (NCEP) Global Forecast System (GFS) analyses and forecast, or in cycling/hot-start mode, where each forecast cycle is initialized using the forecast from a previous cycle as background, creating a custom analysis with all the observations available within the operational MKWC data flow (meteorological surface, maritime radiosondes, aircraft observations, and satellite-derived winds), which is hereafter referred to as MK-OBS. An important limit of the current system is the lack of microwave observations in the assimilation process. Inclusion of this kind of data into Mirto or as direct radiance assimilation is envisioned for the near future. In practice, the forecast cycle can work with different cycling frequencies. The MKWC routine operation uses a 6-h updating frequency and refreshes the cycle with a cold-start run every 3 days. Whether in cold or cycling mode, each forecast cycle produces a 60-h-long forecast.

To assess the performance of the assimilation system presented in this work, a variation of the MKWC routine configuration is used: only the coarse 4.5-km domain is retained and the WRF modeling system is run in cycling mode with a 3-h cycling frequency (Fig. 2), which allows timely assimilation of both sensors given the overpasses times of *Suomi NPP* and *NOAA-20* and *MetOp-A* and *MetOp-B*. The polar orbiter *Suomi NPP* and *NOAA-20*, carrying CrIS, overpass the central Pacific area between 1030 and 1230 UTC and 2230 and 0030 UTC, whereas the two polar orbiters, *MetOp-A* and *MetOp-B*, carrying IASI, overpass the area between 0730 and 0930 UTC and 1930 and 2130 UTC, respectively. Because of this overpass schedule, hyperspectral data can be assimilated at 0000, 0900, 1200, and 2100 UTC. Choosing a ±1.5-h assimilation window and depending on the overpass times of *MetOp-A* and *MetOp-B*, data could also be available at times for the 0600 and 1800 UTC cycles.

Level-1 data from CrIS and IASI are fed to the Mirto processor, which produces 1DVAR physical retrievals of temperature and relative humidity, instability indices, and TRs, along with their observation operator, all in near–real time. TRs and their corresponding observation operators are the quantities used in the assimilation process within the modified WRFDA system, as explained in next section. Figure 3 provides examples of the distribution of clear-sky fields of view (FOVs) associated with successful Mirto retrievals. These distributions were obtained by patching together the adjacent overpasses that overlap the domain in Fig. 1 and that occur during the 3-h windows centered at 2100 and 0000 UTC for IASI and CrIS data, respectively. The aggregated data for each overpass are quality controlled and then thinned to 90 km, which means that assimilated observations have been subselected to ensure a minimum distance of 90 km between them. Quality control is performed by excluding 1) FOVs with cloud contamination over 5%, 2) FOVs whose retrievals show relative humidity ≥100%, and 3) FOVs that did not reach convergence. The thresholds used in the quality control (QC) have been set according to the authors’ experience. For point 2, the IASI cloud mask is estimated from the Advanced Very High Resolution Radiometer (AVHRR), and it is distributed in the level-1C products, whereas the CrIS cloud mask is estimated from Visible Infrared Imaging Radiometer Suite (VIIRS) radiance at run time via a collocation software package developed by AdaptiveMeteo S.R.L.

Preliminary results from a case study of a rainfall forecast over Kauai in which TRs are assimilated is presented and discussed in section 4. The full investigation of this case as well as a more in-depth analysis of the impact of the new assimilation approach on the forecast's accuracy will be material for two subsequent papers.

## 3. Assimilation system design and underlying theory

In this section an overview of the theory underlying the assimilation design is provided. The Earth-emitted infrared spectrum at the top of the atmosphere *y*_{rad} is related to the true continuous atmospheric fields $xct$ by the radiative transfer operator $H$ such that

where *ε*^{0} is the observation error with variance $R$.

Let us denote with $yrad0$ a vector representing the radiances observed by a satellite instrument; Eq. (1) can then be discretized in the form

where *H* is the forward model operator, which implements the radiative transfer equation in a finite dimensional space; **x**^{t} is the discretized representation of $xct$, defined over a vertical pressure grid *p*_{t}; and $\epsilon rad0$ (with covariance $R$_{rad}) is the total observation error, which includes the measurement error *ε*^{0} (which after discretization has covariance $R$_{ε}) and forward model error with covariance $R$_{FM}.

Under these circumstances, given a real satellite hyperspectral observation $yrad0$, the retrieval assimilation problem is reduced to 1) finding the maximum a posteriori (MAP) estimate of $xct$ (inverse problem solution) for each observation, 2) applying the transformation of the MAP solution into a small number of statistically independent quantities that can be directly ingested by an assimilation system, and 3) assimilating the transformed MAP solution with WRFDA.

The overall system implemented at the MKWC consists of two modules: the inversion module (section 3a), responsible for steps 1 and 2, that is, for mapping the observed radiances into vertical profiles of atmospheric parameters and for applying Migliorini’s transformation and the assimilation module (section 3b), responsible for step 3, designed to combine the information embedded in the TRs with the information on atmospheric state (background) available prior the observations.

Both modules implement Bayesian approaches to solve noisy inverse problems in which prior understanding or expectation about vertical atmospheric temperature and water vapor profiles are updated in the light of the information made available by new observations. The first module updates a forecast valid at satellite overpass time using the information provided by hyperspectral IR observations (inversion module), while the second module updates the initial conditions, or background, of a forecast model, using the output of the first module along with MK-OBS. Figure 4 details the structure of the system.

The modules are both designed to minimize cost functions of similar form:

The first term on the right-hand side of Eq. (3) can be thought as the volume in the state space, where the solution is likely to be found with a certain confidence, before the observation is made, and it is given by the square difference of the atmospheric state **x** to the a priori or background state **x**_{b} weighted by the a priori covariance $B$. The second term on the right-hand side of Eq. (3) represents the region in which the states consistent with the observations could be found, and it is given by the distance between the measurements **y** and the retrieved state mapped into measurement space using the observation operator *H* weighted by the total error covariance matrix $R$.

The inversion module tries to minimize Eq. (3), when the observation $y=yrad0$ is a vector representing the radiances observed by a satellite instrument, *H* is the forward model [Eq. (2)] that maps the atmospheric state into radiances, and $R$ = $R$_{rad} is the uncertainty on the observations due to instrument noise and forward model errors. The a priori $xbret$ is the best characterization of the atmospheric state available prior the satellite overpass. In this study the atmospheric temperature and water vapor components of $xbret$ are from the WRF RUC forecast valid at the time of the satellite overpass and are defined over the vertical pressure grid **p**_{ret}, which is identical to **p**_{wrf} for the first 51 levels and has equally spaced in log(*p*) levels from the top of the model to the top of the atmosphere at 0.01 hPa. The uncertainties related to the a priori state $xbret$ are described by the background covariance error $B$ = $B$_{ret}, with The assimilation module, following the same naming convention introduced in Migliorini (2012), tries to minimize Eq. (3), when $y=yret\u2032$ are the TRs and $H=Hret\u2032$ is the observation operator, which maps the atmospheric space into the transformed retrieval space [Eq. (19)]. The a priori $xb=xb*$ is the best characterization of the atmospheric state available prior to the assimilation time. In this study the atmospheric temperature and water vapor are from the most recent WRF forecast validating the assimilation time. The uncertainties related to the a priori state $xb*$ are described by the background covariance error $B$ = $B$*, with $xbret\u2260xb*$ and $B$_{ret} ≠ $B$*. This choice represents the most general case to mimic the situation where the TRs are generated from an independent data provider and are based on an a priori potentially incompatible with the assimilation background. In radiance assimilation, the final analysis is prior dependent. According to Migliorini’s work, provided that the indicated conditions on the linearity and on the covariances are valid, radiance and retrieval assimilation are equivalent, even when the TRs are generated from an a priori independent from the model background. Therefore, even if the dependence of the analysis from the background prior becomes more complicated for the TR assimilation, the equivalence holds true, and the contamination of the analysis from the retrieval prior is mitigated. It is worth stressing how important it is to ensure that TRs are generated using a *good* background forecast and background error covariance, that is, consistent as much as possible with those used in operational meteorological centers to meet the above mentioned conditions.

Section 3a provides a description of the inversion module (radiance inversion and retrieval transformation), and section 3b describes the assimilation module. Finally, section 3c describes the forecasting experiment designed to test the assimilation system.

### a. Inversion module

The inversion module is based on Mirto, a 1DVAR system that allows for the retrieval of atmospheric temperature, water vapor mixing ratio, carbon dioxide, ozone, and surface temperature and emissivity from calibrated radiances observed by high-spectral-resolution infrared instruments (Antonelli et al. 2017). Besides inverting the radiances into physically based retrievals, the inversion module is also in charge of applying Migliorini’s transformation to generate physically based TRs.

#### 1) From radiances to physical retrievals

Mapping of observed radiances into nonobservable atmospheric variables is achieved through Bayesian inversion and, at this stage, is intended for clear-sky cases only. From a theoretical point of view retrievals of temperature and water vapor can be also done in cloudy conditions. In practice, in infrared remote sensing, the inversion process requires special techniques such as cloud clearing for partially cloudy conditions (Cho and Staelin 2006; Gambacorta et al. 2014), and adaptive surface level, for overcast conditions. Neither of these techniques are currently implemented in Mirto (Antonelli et al. 2017), which is the system used to perform the inversion. The package uses the Atmospheric and Environmental Research, Inc. (AER), Optimal Spectral Sampling (OSS) radiative transfer model (Moncet et al. 2008) for the computation of simulated radiances and Jacobians [Eq. (2)]. At the time of writing, the inversion package Mirto works with IASI (3303 channels), and CrIS (826 channels at reduced resolution) observations. In Mirto the MAP solution for **x**^{t}, given the observed radiances $yrad0$, is found by minimizing the cost function of the form described in Eq. (3),

and it is obtained through an iterative process, expressed in terms of departure from $xbret$.

For the chosen case study, in the retrieval module, the $xbret$ components for atmospheric temperature and water vapor comes from the RUC system, that is the 3-h forecast from the run started at 1800 and 2100 UTC 12 April 2018 for the IASI (2024 and 2139 UTC) and CrIS (at 2338 UTC) overpasses, respectively, while $B$_{ret} is derived from $B$* by retaining the diagonal elements only. The elements of $xbret$ and $B$_{ret} relative to carbon dioxide, ozone, surface temperature, and surface emissivity were generated using the method described in Antonelli et al. (2017).

The total error covariance matrix used in Mirto $R$_{rad} is defined as

where $R$_{m} represents the errors in the measurements (instrument noise), and $R$_{H} represents the errors in the observation operator (forward model). This last term represents the errors associated with the radiative transfer model like the errors in the spectroscopical parameters (line position and half-width) and the errors due to the atmospheric/surface variables not specifically accounted for in the radiative transfer calculations.

The total error covariance matrix as defined in Eq. (5) was estimated from retrieval residuals, $\rho rad=yrad\u2212H\u2061(x^ret)$, obtained from a large dataset of clear-sky retrievals that reached convergence. Considering that the diagonal of the spectral residuals covariance is an estimator of the instrument noise, for the first term on the right side of Eq. (5), cov(*ρ*_{rad}) was calculated to represent a lower limit for $R$_{rad}. Therefore it was chosen to inflate the diagonal terms of $R$_{m} in $R$_{H} = $R$_{rad} − $R$_{m}, by adding the spectral residual bias $\rho \xafrad$.

The solution after convergence becomes

where $x^$ is the state vector at convergence and $H^=\u2202H/\u2202x|x=x^MAP$ is the Jacobian of the observation operator *H* at convergence [Rodgers (2000), their Eq. (5.9)]. The forecast profile is used both as the prior profile $xbret$ in the retrieval and also as the first guess for the iterative retrieval (as is normal practice in 1DVAR, except for the cloud variables). In principle, the choice of the first guess should not affect the value of the solution, unless there are multiple minima, but it affects the speed of convergence, while the prior $xbret$ clearly affects the value of the solution through Eq. (6).

Retrieved profiles shown in Fig. 5 were obtained for a CrIS FOVs located ~13 km from Hilo (Fig. 5a) and ~40 km from Lihue (Fig. 5b). The retrieval near Hilo achieved convergence in five iterations, with a final *χ*^{2} = 1.54 where $\chi i2=\u2061[yrad\u2212H\u2061(xi)]TRrad\u22121\u2061[yrad\u2212H\u2061(xi)]$. The total number of degrees of freedom for signal calculated for temperature and water vapor from Rodgers (2000),

gave *d*_{s} = 9.65 (with *d*_{stemp} = 2.26 for temperature and *d*_{swv} = 7.39 for water vapor). This is the number of independent quantities in the selected CrIS observed spectrum useful to characterize the vertical profiles of temperature and water vapor.

#### 2) From physical retrieval to transformed retrievals

The inversion module also applies the linear transformation proposed by Migliorini (Migliorini et al. 2008; Migliorini 2012). Under the assumption that the observation operator is linear within the error bounds of the retrievals, this transformation can be applied to the retrievals to 1) mitigate the dependence of the retrieval from the a priori by isolating the portion of the profile effectively constrained by the observations, 2) compress as much as possible the information content embedded in the observations, and 3) map the retrievals into an orthogonal space where the associated errors are uncorrelated. The transformation goes through several steps to map the retrieval solution $x^MAP$ in Eq. (6) into the TRs. The first step is to calculate the *p* nonsingular eigenvectors **L**_{p} associated with the *p* nonzero eigenvalues **Σ**_{p} of the observation error covariance $R$_{rad}:

The second step is to estimate the normalized transformed radiances and Jacobians, $yrad\u2032$ and $Hrad\u2032$, by projecting the radiances $yrad0$ and Jacobians $H$_{rad}, normalized by the square root of the *p* ≤ *m* nonzero eigenvalues **Σ**_{p}, onto the eigenvectors **L**_{p}:

The third step is to calculate and diagonalize the signal-to-noise matrix:

The last step is then to project the normalized transformed radiances and Jacobians, $yrad\u2032$ and $Hrad\u2032$ on the *r* eigenvectors of the signal-to-noise matrix associated with eigenvalues *λ*_{i=1,…,r} ≥ 1 that provide an orthonormal base of dimension *r* ≤ *n* onto which the projection of normalized TRs and the normalized Jacobians produce the actual components to be assimilated, that is, the TRs:

and the forward model *H* in Eq. (2) into the transformed retrieval observation operator:

The TRs $yret\u2032$ represent the fraction of the *true* state “readable” by the instrument (along with its indetermination *ε*_{ret}). They are the candidate observations for data assimilation, as they do not introduce significant bias in the analysis due to the a priori information used in the inversion problem that is potentially inconsistent with the analysis background. Thus, the TRs are compressed to retain the atmospheric signal embedded in the original observed radiances, while filtering out most of the random component of the observation error (Antonelli et al. 2004).

For the data used in this study, Fig. 6 shows the mean and the standard deviation of the eigenvalues of the signal-to-noise matrix [Eq. (11)], whereas Figs. 7 and 8 show $yret\u2032$ and $Hret\u2032$ for temperature in kelvins and water vapor in terms of ln(*q*) with *q* in kilograms per kilogram. In particular, the first figure shows mean and the standard deviation of $yret\u2032$ for the 135 IASI FOVs (Fig. 7a), and the 91 CrIS FOVs (Fig. 7b) selected respectively for the 2100 and 0000 UTC assimilations, while Fig. 8 shows the mean and the standard deviation of $Hret\u2032$ for the same IASI (Fig. 8a) and CrIS (Fig. 8b) FOVs. The 10 functions showed in each panel of Fig. 8 represent the first 10 columns of the observation operator corresponding to the eigenvalues of the signal-to-noise matrix [Eq. (11)] larger than one as shown in Fig. 6. Differences between the IASI and CrIS functions for the two overpasses taken into consideration are visible in the figures; any speculation on the possible reasons for such differences is deferred to an analysis that will be done on a statistically significant dataset.

To properly and efficiently assimilate hyperspectral observations in the form of TRs and to obtain an unbiased analysis, the forward model operator *H* in Eq. (2) needs to be linear within the error bounds of the retrievals. This requirement is checked to ensure the applicability of the approach suggested by Migliorini. Whether a forward model can be considered linear is addressed by evaluating the difference between the forward model itself and its linearization about the retrieval solution $x^MAP$,

for values of **x**^{t} within one standard deviation from $x^MAP$, that is, $xt=x^MAP\xb1ei$ where the error pattern $ei=\lambda i1/2li$ is the product between the square root of the *i*th eigenvalue and the *i*th eigenvector of the retrieval error covariance. Therefore, the first test is to verify that the covariance of *δ***y** is small relative to $R$_{rad}, that is, *c*^{2} = *δ***y**^{T}$R$_{rad}*δ***y** is a *χ*^{2}-like quantity significantly smaller than 1. As an example, Fig. 9 compares the differences *δ***y** between the forward model *H*(**x**^{t}) and its linear approximation $H\u2061(x^MAP)\u2212H\u2061(xt\u2212x^MAP)$, in red, with the square root of the diagonal elements of the observation error $R$_{rad} (in black) for a single FOV. For each of the selected 91 FOVs of the CrIS overpass *c*^{2} was calculated using for **x**^{t} 10 error patterns all with the same sign so to get in one of the two extreme cases of error, and we found a mean value $c2\xaf=0.0017$ with a maximum value of $cmax2=0.01$.

A second test proposed by de Haan aims to detect possible anomalies in the transformation of the retrieval solutions by using the differences

whose covariance should have expected values that are dependent on the eigenvalues of the signal-to-noise matrix [Eq. (11)] as shown in appendix A:

The expected values in Eq. (17) lies between 0 and 1. The smallest values occur for the largest eigenvalues and are related to the information content of $yret\u2032$. Values close to 0 indicate almost independent (informative) observations, while values close to 1 are associated with elements of $yret\u2032$ that do not contain additional information.

Figures 10a and 10b show that for all the eigenvectors associated with eigenvalues greater than one, the observed values are close to their expected values, for both IASI (Fig. 10a) and CrIS (Fig. 10b), while for higher-order eigenvectors, for IASI, there is a departure from the observed and expected values that needs further investigation.

The expected values in Eq. (18) are always larger than 1, with the maximum occurring for the largest eigenvalues. Figures 10c and 10d show that also in this case for all the eigenvectors associated with eigenvalues greater than one, the observed values are close to their expected values for both IASI (Fig. 10c) and CrIS (Fig. 10d). In addition to demonstrating the correct functioning of the inversion module, this test could play an important role for the quality control in the assimilation of the TRs (section 3b).

### b. Assimilation of transformed retrievals with WRFDA

The assimilation problem consists in estimating the analysis $x^MAPan$ that minimizes the cost function *J*^{an}(**x**) (Migliorini 2012):

and $x^MAPan$ can be expressed as

where

with $B$* and $xb*$ being the prior information used to constrain $yret\u2032$. Specifically, $xb*$ is obtained from the WRF operational forecast from the 9- and 12-h runs started at 1200 UTC 12 April 2018, and it is defined on a pressure vertical grid **p*** ≠ **p**_{ret}; the background error covariance $B$* has been estimated following the National Meteorological Center (NMC) method (Parrish and Derber, 1992): over a month of WRF output from both the 0000 and 1200 UTC cycles were used to calculate the differences between each WRF 24-h forecast and the WRF 12-h forecast valid at the same time from the following cycle and then to generating their covariance.

For Eq. (21) to be implemented, $B$* and $xb*$ are expected to have the same order and units of $B$_{ret} and $xbret$, which means that $H\u02dcret\u2032=Hret\u2032HconvHinterp$ accounts for the proper vertical/horizontal interpolation through the operator $H$_{interp} and for the proper unit conversion through the operator $H$_{conv}. By retaining the diagonal terms and by removing the off-diagonal terms of $B$* in the characterization of $B$_{ret}, all of the cross-correlation terms have been removed, and $B$_{ret} becomes less constraining than $B$* as suggested by Migliorini et al. (2008) and Prates et al. (2016).

In the WRFDA implementation, the TRs are assimilated by an adapted version of the original satellite radiance module. The WRFDA radiance assimilation module has been expanded to include a new instrument called TRs. In WRFDA, each sensor is identified by a triplet used in the radiative transfer model [either RTTOV or Community Radiative Transfer Model (CRTM)], which includes the platform number, satellite identification number, and the sensor number. TR is uniquely identified by its new sensor number. Each eigenvector of observation operator in Eq. (13) is considered a channel and thus the WRFDA assimilates the observations $yret\u2032$ in Eq. (12) as if they were separate channels from a new instrument called TR.

The implemented observation operator $H\u02dc\u2032ret$ computes the model equivalent (simulated data) for every observation; $H\u02dcret\u2032$ consists of three components. The first component is the three-dimensional interpolator $H$_{interp} that interpolates the model background temperature and humidity horizontally to the observation location (using latitude and longitude from the observation FOV), and vertically to the number of pressure levels used by Mirto. For the vertical interpolation, piecewise functions were adopted using the model vertical levels at the observation location (e.g., de Haan et al. 2004). The second component of $H\u02dc\u2032ret$ is the unit converter $H$_{conv}, which converts the background units for water vapor from grams per kilogram to the Mirto internal units of log(*q*) with *q* in kilograms per kilogram. The third component of $H\u02dcret\u2032$ is the projection operator $Hret\u2032$ that maps the interpolated temperature and humidity profiles into the TR space.

We emphasize that, in the current configuration, neither Mirto nor the WRFDA module for TR assimilation implement any bias correction scheme. Since it has been demonstrated that the bias correction is very relevant for an assimilation experiments to be successful (Eyre 2016), a variational bias correction scheme (Desroziers et al. 2005; Eyre 2016) is under development for the WRFDA module. The opportunity to implement a variational bias correction (Merchant et al. 2020) for Mirto, and hence to redefine the observation error covariance matrix, will be evaluated in a follow-up study.

When no cloud mask is available, or in the presence of not-detected clouds, to screen out cloud-contaminated retrievals, the implemented WRFDA module for the TR observations has also the option to performs a QC based on Eq. (22). Since for every element *i* = 1, …, *r* of the $yret\u2032\u2061(i)$ vector, the observation minus background variance is $\u2061(1+\lambda i2)$, and since this value changes from element to element in the observation vector $yret\u2032$ and it also changes from retrieval to retrieval, it seems reasonable to use (1 + *λ*^{2})^{1/2} to scale the *i*th components of the TRs $\upsilon ret\u2032\u2061(i)=yret\u2032\u2061(i)/1+\lambda i2$ and the corresponding row of observation operator $Hret\u2032=Hret\u2032/1+\lambda 2$. Equation (18) then can be written

Because of this normalization, the difference between the observation and background becomes more independent from the element position in the array and from the individual retrieval. The QC would reject any normalized observation whose difference $\u2061(\upsilon ret\u2032\u2212Hret\u2032xb)$ is greater than 3 (i.e., outside the 3*σ* confidence interval).

### c. Forecasting experiment

To test the TR assimilation within the MKWC 3-h cycling system, the following experiment was designed for a record breaking heavy rainfall event that occurred in April 2018 over the western Hawaiian Islands of Oahu and Kauai. The underlying meteorological conditions on 12 April include a deep tropical upper-tropospheric trough (TUTT) developing to the northwest of the State of Hawaii, bringing increasingly cold air aloft with it. At the surface, a low-level ridge to the northeast of Hawaii produced easterly trade winds that advected a continuous supply of deep moisture into the area. During subsequent days the TUTT deepened, further destabilizing the air mass, and slowly shifted westward. On 15 April 2018, the U.S. rainfall record for a 24-h period was broken, with ~50 in. (127 cm) rainfall recorded on the north coast of Kauai. This historic event provides a uniquely important opportunity to investigate the potential of assimilating hyperspectral TRs to create a more complete picture of the state of the atmosphere in the hours before the flooding began. IASI and CrIS TRs have been produced by the Mirto system, using as a priori the default MKWC WRF 12/6 h forecast validating the satellites’ overpass times, for each MetOp and NPP overpass available in the time frame from 2100 UTC 12 April to 0000 UTC 15 April. MK-OBS for this time frame have also been considered for assimilation. The 3-h WRF RUC system was started on 2100 UTC 12 April with an assimilation on the 9-h WRF forecast initialized on 1200 UTC 12 April (cold run), and it ran up to 0000 UTC 15 April for a total of 18 assimilation processes and relative forecasts (OP experiment). Each assimilation used as background the previous cycle 3-h forecast. The cycling system so designed is expected to enhance the cumulative impact of assimilating observations through the various cycles. A second WRF RUC system was also run as a control experiment (CNTR) in which assimilation of the TRs was turned off. Results from these experiments are presented in section 4. The Kauai 2018 case study and the impact of this new assimilation approach on the forecast accuracy of the modeling system will be carefully investigated and detailed in a subsequent paper.

## 4. Results

One of the strengths of the assimilation system described in this paper is its independence from the characteristics of the hyperspectral instruments. The two assimilation scenarios presented in this section used data from two different instruments, IASI and CrIS. All the instrument related variables (spectral resolution, instrument noise, viewing geometry, etc.) are embedded in the 1DVAR inversion (Mirto). Therefore, the TRs derived from two different hyperspectral sensors can be treated in exactly the same way by the assimilation system as all the information regarding the sensors and the observation error covariances are embedded in $yret\u2032$ and $Hret\u2032$, both made available by the retrieval provider (Mirto). The assimilation was done following Migliorini et al.’s (2008) recipe with the covariance of the observation error used for the TRs being the identity matrix. The magnitude of the temperature and water vapor analysis increments, defined as the difference between the analysis $x^MAPan$ and the background $xb*$, is comparable for the two instruments, with the increments being, in both cases, smaller near the surface (Figs. 11 and 12a,b ) than middle model levels (Figs. 11 and 12c,d). This is to be expected given the vertical structure of the observation operator (Fig. 8): smaller near the surface than in the middle atmospheric levels. The results shown imply that most of the atmospheric signal contributed by the hyperspectral observations comes from the region between the upper part of the boundary layer around 800 and 250 hPa.

The temperature field is in general well predicted over the central North Pacific region, likely due to its smooth spatial variability that is well constrained by the initial and boundary conditions, while the water vapor field is less well constrained. The results show that, unlike the analysis increment for potential temperature, the ones for water vapor do reflect a higher granularity of the highly variable humidity field. Therefore, a larger impact is expected from the assimilation of the TRs in the water vapor analyses.

The analysis increments in potential temperature at 2100 UTC (Fig. 13a) are smaller in magnitude, compared to the 0000 assimilation (Fig. 13b), while for water vapor (Fig. 13c) they are comparable to the 0000 UTC ones (Fig. 13d). Figure 13b shows a layer of negative temperature analysis increments between 100 and 300 hPa, neutral down to 400 hPa, and then positive from 400 hPa down to 850 hPa for the CrIS assimilation. As for humidity, Figs. 13c and 13d show larger analysis increments both negative and positive in the 950–700-hPa layers, that is, a layer around the trade wind inversion, which is where the satellite would be expected to help the most.

As a reference to compare the analysis increments obtained with the assimilation of TRs, a second run was carried out assimilating only MK-OBS for the 0000 UTC run assimilation. Figures 14 and 15 show the temperature and water vapor analysis increments for the MK-OBS runs for 2100 and 0000 UTC so that they can be directly compared to Figs. 11 and 12. In addition, Fig. 16 shows the analysis increments produced by MK-OBS only on the same vertical cross section of Fig. 13. From the assimilation of MK-OBS, it clearly emerges that the dominant contribution to analysis increments is due to radiosondes. From the direct comparison, the temperature analysis increments obtained in the two assimilations are consistent in size and sign for both the assimilations at 2100 and 0000 UTC. The analysis increments related to the assimilation of TRs, however, have a significantly wider spatial coverage and much higher representativeness, especially in terms of water vapor distribution. This is likely due to the paucity of water vapor observations available for the MK-OBS over the central Pacific area. If compared to an analysis done with satellite radiances, the representativeness would likely be comparable.

To further assess the correct functioning of the assimilation module, it is necessary to show that the analysis increments actually bring the background fields closer to the observations. To measure the distance from the atmospheric fields and the observations, the background and the analyses are mapped into radiance space through the application of the forward model *H* as in Eq. (2). For each FOV used in the assimilation, a state vector **x** is created using the vertical profiles of temperatures and water vapor interpolated from the background or from the analyses and borrowing the components for carbon dioxide, ozone, and surface parameters (temperatures and emissivity) from physical (untransformed) retrievals. The state vectors $xb*$ and $x^MAPan$ are then fed to the forward model *H* to create the synthetic spectra, that is, the radiances that would be observed by the instrument if the atmospheric part of **x**^{t} was identical to the background or the analysis profile. The differences between the radiances thus obtained and the one actually observed generate the residuals $\rho bg=yrad0\u2212H\u2061(xb*)$ and $\rho an=yrad0\u2212H\u2061(x^MAPan)$. Figures 17 and 18 show the root-mean-square (RMS) of the residuals for the 2100 (IASI) and 0000 (CrIS) UTC scenarios calculated over the 135 and 91 FOV used in the assimilation. The RMS of *ρ*_{bg} is represented by the blue points, and it is larger than the RMS of *ρ*_{an}, red points, while the observation error used as reference is represented by the dotted line in black. The reduction in the spectral residual derived from the analysis suggests that the assimilation of the TRs brings the simulated radiances closer to the ones actually observed by the satellite instruments as it is supposed to do.

Figure19 shows the comparison between the OP and CNTR experiments, in terms of RMS of the 3-h accumulated precipitation predicted by the 18 forecasts completed during the full cycle. The observed (true) 3-h precipitation used in the RMS calculation comes from the Mount Wai‘ale‘ale rain gauge (WLLH1). While a precipitation analysis on the scale of the WRF domain is not available at the time of writing, the WLLH1 rain gauge is chosen among all available because it is the one showing the early onset of this precipitation event and reflects the large-scale signal (Chambers 2003), and also because its observations are independent because they have not been used in the process of assimilation and creation of the analysis. The large-scale signal is reflected by the model prediction, given the chosen horizontal grid resolution of 4.5 km. The other available rain gauges show instead the local effects of the strong orographic enhancement that, along with the large moisture content of the air mass, contributed to the record-breaking event recorded at Waipa Garden, in the northern part of Kauai. Predicting with high accuracy the distribution of this precipitation event is beyond the scope of this paper and would require a much higher horizontal grid resolution; therefore, it is left to a separate study. The RMS comparison suggests an improvement in the forecast accuracy when TRs are assimilated versus when they are not. The improvement is more evident during the second day of assimilation both because of the expected cumulative positive impact and because the magnitude of the rainfall was increasing then. Further studies to asses the efficacy of this methodology and to characterize the full impact on the forecast accuracy are needed. A statistical analysis over a longer period of assimilation to assess the impact of the TR assimilation is the subject of a second experiment currently in progress and of a subsequent paper.

## 5. Conclusions

The launch of satellites carrying hyperspectral sensors has created renewed interest in assimilating satellite retrievals, especially for limited-area weather models. Satellite retrievals that attempt to exploit the information content embedded in thousands of channels provided by these very high-resolution sensors represent an opportunity for a formidable gain in the assimilation computational efficiency over current radiance assimilation methods by off-loading the radiative transfer calculation to the retrieval provider. An obstacle to successful retrieval assimilation has historically been the bias introduced by the retrieval prior. However, Rodgers (2000) and Migliorini et al. (2008) provide a theoretical framework to mitigate possible inconsistencies between retrieval and assimilation a priori information and identify the circumstances under which the radiance and retrieval assimilations are equivalent. In this paper, the Migliorini approach was implemented and tested in WRF 3DVAR, using data provided by the 1DVAR inversion package Mirto (Antonelli et al. 2017). The resulting assimilation system is integral to a WRF RUC system that nominally runs every three hours and allows for consistent assimilation of IASI and CrIS data over the central Pacific area surrounding Hawaii. At the MKWC, Mirto takes on the role of the data provider and TRs are produced in house; however, TRs provided by an external source could also be ingested without any change to the assimilation system. Another positive aspect of the assimilation system detailed here is its independence from the characteristics of the hyperspectral instruments. The analysis results presented in this paper are derived by assimilating data from two different instruments, IASI and CrIS. All the instrument related variables (spectral coverage, spectral resolution, instrument noise, viewing geometry, etc.) are embedded in the 1DVAR inversion (Mirto) and in the Migliorini transformation. Therefore, the TRs derived from the two different hyperspectral systems have been treated in exactly the same way by the assimilation system, that is, the assimilation module did not need to know whether the TRs were generated from IASI from CrIS. Significantly, the results of the data assimilation experiments described in this paper, compared to the case where only MK-OBS were assimilated, confirm the implementability of Migliorini’s approach and open the way to dedicated studies to assess the validity of the theory.

The results indicate that a larger impact of the assimilation of hyperspectral data is expected from changes in the water vapor distribution rather than in the temperature field because of the higher spatial variability of the water vapor concentrations and the independence of the water vapor field from temperature and pressure. Therefore, beyond demonstrating the practicality of the Migliorini approach for hyperspectral data assimilation, the potential for improving the regional weather prediction system currently implemented at the MKWC was also demonstrated. Preliminary results comparing forecasts from a cycling system with and without TR assimilation have shown a decrease in RMS of the 3-h accumulated precipitation forecast. Significant additional work remains to fully explore the impact of the new assimilation approach on the accuracy of forecasts and to further assess the validity of Migliorini’s theoretical work. A subsequent paper that is in preparation will detail the impact of the assimilation system on the accuracy of the WRF forecast of precipitation for the record-breaking Kauai flood event using a high-resolution nested grid. In addition, a longer-running experiment to demonstrate the impact of the assimilation of TRs more generally is also in the works.

## Acknowledgments

We thank the SSEC and the MKWC for supporting the implementation of the regional retrieval system based on Mirto; Dr. Graziano Giuliani for his contribution in creation of Mirto; Drs. Thomas Auligne and Francois Vendenberge for their contribution in creating the prototype of the TR assimilation module in WRFDA; Drs. Stefano Piani and Livio Bernardini for their continuous and highly qualified support to the development of both Mirto and the WRFDA assimilation module; Ms. Nancy Hulbirt for her support with the paper editing; Atmospheric and Environmental Research (AER) for the use of the Optimal Spectral Sampling forward model; the MTG-IRS program of EUMETSAT for supporting years of investigation on Migliorini’s approach; the UCAR Unidata Program for the dissemination of real-time GFS products; the University of Wyoming for their soundings archive; the MADIS program of NOAA’s Office of Oceanic and Atmospheric Research Earth System Research Laboratory/Global Systems Division for the real-time dissemination of meteorological data. We also acknowledge high-performance computing support from Cheyenne (doi:10.5065/D6RX99HX) provided by NCAR’s Computational and Information Systems Laboratory, sponsored by the National Science Foundation. Last, we desire to thank the reviewers for their precious contribution to this paper. This research was supported by the Office of Naval Research under ONR Award N00014-18-1-2166.

### APPENDIX A

#### Derivation of Eq. (17)

This appendix provides the description of the steps needed to calculate the quantity in Eq. (17):

As suggested by de Haan, the demonstration of Eq. (17) starts from a manipulation of Eq. (15), which is obtained from Eqs. (16), (22), (24), and (26) of Migliorini et al. (2008):

which can be further manipulated to obtain

To avoid the use of subscripts, and $yret\u2032$ being of dimension *r*, in the following demonstrations $U$, $V$, and **Λ** will represent $U$_{r}, $V$_{r}, and **Λ**_{r}. Equations (A2) and (A3) then become

The second term of (A4) can be rewritten as

so that Eq. (A1) can be rewritten as

With a focus on the products inside the most external parentheses,

thus

### APPENDIX B

#### Derivation of Eq. (18)

This appendix provides the description of the steps needed to calculate the quantity in Eq. (18):

which, considering $\epsilon rad\u2032$ to be uncorrelated with $\epsilon b=(xt\u2212xbret)\u2061(xt\u2212xbret)T$, implies

Using the alternative expression in Eq. (24) of Migliorini et al. (2008), $Hret\u2032=\Lambda rVrTB\u22121/2$, it is possible to write

## REFERENCES

*Quart. J. Roy. Meteor. Soc.*

*J. Geophys. Res.*

*J. Appl. Meteor. Climatol.*

*Mon. Wea. Rev.*

*Bull. Amer. Meteor. Soc.*

*Bull. Amer. Meteor. Soc.*

*Bull. Amer. Meteor. Soc.*

*Seeing Clearly: The Impact of Atmospheric Turbulence on the Propagation of Extraterrestrial Radiation*, S. Businger and T. Cherubini, Eds., VBW Publishing, 165–182

*J. Geophys. Res.*

*J. Appl. Meteor.*

*Mon. Wea. Rev.*

*Quart. J. Roy. Meteor. Soc.*

*Seminar on Recent Developments in the Use of Satellite Observations in Numerical Weather Prediction*, Shinfield Park, Reading, United Kingdom, ECMWF

*Quart. J. Roy. Meteor. Soc.*

*Quart. J. Roy. Meteor. Soc.*

*10th Annual Symp. on New Generation Operational Environmental Satellite Systems*, Atlanta, GA, Amer. Meteor. Soc., 9.1

*Quart. J. Roy. Meteor. Soc.*

_{2}, ash, and H

_{2}SO

_{4}using the Infrared Atmospheric Sounding Interferometer (IASI)

*J. Geophys. Res.*

*Nimbus-6*temperature soundings on Australian region forecasts

*Bull. Amer. Meteor. Soc.*

*Mon. Wea. Rev.*

*J. Geophys. Res.*

*Remote Sens. Environ.*

*Mon. Wea. Rev.*

*Mon. Wea. Rev.*

*J. Atmos. Sci.*

*Mon. Wea. Rev.*

*Quart. J. Roy. Meteor. Soc.*

*Inverse Methods for Atmospheric Soundings: Theory and Practice*

*J. Geophys. Res. Atmos.*

*Bull. Amer. Meteor. Soc.*

*Adv. Atmos. Sci.*

Proc. 2001 Int. Geoscience and Remote Sensing Symp., Sydney, NSW, Australia, IEEE, 1341–1343