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).

Fig. 1.

WRF Model domains. The inset shows terrain contours at 100-m intervals.

Fig. 1.

WRF Model domains. The inset shows terrain contours at 100-m intervals.

Fifty-one levels in the vertical direction, pwrf(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.

Fig. 2.

WRF cycling schedule.

Fig. 2.

WRF cycling schedule.

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.

Fig. 3.

Locations of convergent, nonsaturated (RH < 100% for any level) retrievals: (a) clear-sky IASI retrievals overlaid on AVHRR channel 3 background image valid for 2020 UTC 12 Apr 2018, and (b) CrIS retrievals overlaid on a VIIRS true color image valid for 2338 UTC 12 Apr 2018.

Fig. 3.

Locations of convergent, nonsaturated (RH < 100% for any level) retrievals: (a) clear-sky IASI retrievals overlaid on AVHRR channel 3 background image valid for 2020 UTC 12 Apr 2018, and (b) CrIS retrievals overlaid on a VIIRS true color image valid for 2338 UTC 12 Apr 2018.

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 yrad is related to the true continuous atmospheric fields xct by the radiative transfer operator H such that

 
yrad=H(xct)+ε0,
(1)

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

 
yrad0=H(xt)+εrad0,
(2)

where H is the forward model operator, which implements the radiative transfer equation in a finite dimensional space; xt is the discretized representation of xct, defined over a vertical pressure grid pt; and εrad0 (with covariance Rrad) is the total observation error, which includes the measurement error ε0 (which after discretization has covariance Rε) and forward model error with covariance RFM.

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.

Fig. 4.

Schematic flowchart of the data assimilation system used to feed the WRF prognostic component.

Fig. 4.

Schematic flowchart of the data assimilation system used to feed the WRF prognostic component.

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

 
J(x)=12(xxb)TB1(xxb)+12[yH(x)]TR1[yH(x)].
(3)

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 xb 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 = Rrad 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 pret, which is identical to pwrf 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 = Bret, with The assimilation module, following the same naming convention introduced in Migliorini (2012), tries to minimize Eq. (3), when y=yret are the TRs and H=Hret 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 xbretxb* and BretB*. 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 xt, given the observed radiances yrad0, is found by minimizing the cost function of the form described in Eq. (3),

 
J(x)=12(xxbret)TB1(xxbret)+12[yrad0H(x)]TRrad1[yrad0H(x)],
(4)

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 Bret is derived from B* by retaining the diagonal elements only. The elements of xbret and Bret 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 Rrad is defined as

 
Rrad=Rm+RH,
(5)

where Rm represents the errors in the measurements (instrument noise), and RH 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, ρrad=yradH(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 Rrad. Therefore it was chosen to inflate the diagonal terms of Rm in RH = RradRm, by adding the spectral residual bias ρ¯rad.

The solution after convergence becomes

 
x^MAP=xbret+(B1+H^TRrad1H^)1H^TRrad1[yradH(x^)+H^(x^xbret)],
(6)

where x^ is the state vector at convergence and H^=H/x|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 χi2=[yradH(xi)]TRrad1[yradH(xi)]. The total number of degrees of freedom for signal calculated for temperature and water vapor from Rodgers (2000),

 
ds=tr[(B1+H^TRrad1H^)1H^TRrad1H^],
(7)

gave ds = 9.65 (with dstemp = 2.26 for temperature and dswv = 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.

Fig. 5.

Example of single CrIS retrievals for temperature and dewpoint temperature in a skew T representation for 2338 UTC 12 Apr 2018. The retrieval first guess (blue) comes from the WRF forecast valid at the overpass time. The retrieval solution is in red, and the closest available soundings for (a) Hilo and (b) Lihue at 0000 UTC are in black.

Fig. 5.

Example of single CrIS retrievals for temperature and dewpoint temperature in a skew T representation for 2338 UTC 12 Apr 2018. The retrieval first guess (blue) comes from the WRF forecast valid at the overpass time. The retrieval solution is in red, and the closest available soundings for (a) Hilo and (b) Lihue at 0000 UTC are in black.

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 Lp associated with the p nonzero eigenvalues Σp of the observation error covariance Rrad:

 
Rrad=LpΣp2LpT.
(8)

The second step is to estimate the normalized transformed radiances and Jacobians, yrad and Hrad, by projecting the radiances yrad0 and Jacobians Hrad, normalized by the square root of the pm nonzero eigenvalues Σp, onto the eigenvectors Lp:

 
yrad0Σp1LpTyrad0and
(9)
 
Hrad(xt)Σp1LpTHrad(xt).
(10)

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

 
S=H^radB1/2=UrΛrVrT.
(11)

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

 
yret=UrTyrad,
(12)

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

 
Hret=UrTHrad.
(13)

The TRs yret 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 and Hret 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 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 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.

Fig. 6.

Distribution of eigenvalues (mean and standard deviation) in logarithmic scale, of the signal-to-noise matrix [Eq. (7)]:(a) the values for the 2100 UTC assimilation (135 IASI FOVs; 3303 channels per FOV) and (b) the values for the 0000 UTC assimilation (91 CrIS FOVs; 826 channels per FOV). The eigenvalues larger than 1, indicated by the horizontal dashed line, are those that actually carry significant information about the true state vector.

Fig. 6.

Distribution of eigenvalues (mean and standard deviation) in logarithmic scale, of the signal-to-noise matrix [Eq. (7)]:(a) the values for the 2100 UTC assimilation (135 IASI FOVs; 3303 channels per FOV) and (b) the values for the 0000 UTC assimilation (91 CrIS FOVs; 826 channels per FOV). The eigenvalues larger than 1, indicated by the horizontal dashed line, are those that actually carry significant information about the true state vector.

Fig. 7.

Distribution of yret (mean and standard deviation): (a) the values for the 2100 UTC assimilation (135 IASI FOVs) and (b) the values for the 0000 UTC assimilation (91 CrIS FOVs).

Fig. 7.

Distribution of yret (mean and standard deviation): (a) the values for the 2100 UTC assimilation (135 IASI FOVs) and (b) the values for the 0000 UTC assimilation (91 CrIS FOVs).

Fig. 8.

The first 10 columns of the assimilation operator Hret (mean and standard deviation) for (a),(b) temperature and (c),(d) water vapor log(q) for (left) the 2100 UTC assimilation (135 IASI FOVs) and (right) the 0000 UTC assimilation (91 CrIS FOVs).

Fig. 8.

The first 10 columns of the assimilation operator Hret (mean and standard deviation) for (a),(b) temperature and (c),(d) water vapor log(q) for (left) the 2100 UTC assimilation (135 IASI FOVs) and (right) the 0000 UTC assimilation (91 CrIS FOVs).

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,

 
δy=H(xt)H(x^MAP)H(xtx^MAP),
(14)

for values of xt within one standard deviation from x^MAP, that is, xt=x^MAP±ei where the error pattern ei=λi1/2li is the product between the square root of the ith eigenvalue and the ith eigenvector of the retrieval error covariance. Therefore, the first test is to verify that the covariance of δy is small relative to Rrad, that is, c2 = δyTRradδy is a χ2-like quantity significantly smaller than 1. As an example, Fig. 9 compares the differences δy between the forward model H(xt) and its linear approximation H(x^MAP)H(xtx^MAP), in red, with the square root of the diagonal elements of the observation error Rrad (in black) for a single FOV. For each of the selected 91 FOVs of the CrIS overpass c2 was calculated using for xt 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¯=0.0017 with a maximum value of cmax2=0.01.

Fig. 9.

CrIS linearity test for (a) longwave and (b) midwave: δy=H(xt)H(x^)H(xtx^) for one CrIS FOV near the Lihue rawinsonde station is shown in red, and the observation error diag(Rrad1/2) is shown in black.

Fig. 9.

CrIS linearity test for (a) longwave and (b) midwave: δy=H(xt)H(x^)H(xtx^) for one CrIS FOV near the Lihue rawinsonde station is shown in red, and the observation error diag(Rrad1/2) is shown in black.

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

 
yretHretx^MAP=Hretxt+εyretHret[xbret+KH(xtxbret)+Kεrad]and
(15)
 
yretHretxbret=yretHretxt+Hret(xtxbret)
(16)

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:

 
(yretHretx^MAP)(yretHretx^MAP)T¯=(ΛrΛrT+Ir)1and
(17)
 
(yretHretxbret)(yretHretxbret)T¯=I+HretBHretT=Ir+ΛrΛrT.
(18)

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. Values close to 0 indicate almost independent (informative) observations, while values close to 1 are associated with elements of yret 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.

Fig. 10.

(a),(b) The square root of the diagonal values of the covariance (yretHretx^)(yretHretx^)T¯ (in red) and the distribution of values (mean and standard deviation) for (ΛΛT + I)−1/2 from Eq. (17) (in black) and (c),(d) the square root of the diagonal values of the covariance (yretHretxb)(yretHretxb)T¯ (in red) and the distribution of values (mean and standard deviation) for (I + ΛΛT)1/2 from Eq. (17) (in black) for the (left) 2100 UTC assimilation (135 IASI FOVs) and (right) 0000 UTC assimilation (91 CrIS FOVs).

Fig. 10.

(a),(b) The square root of the diagonal values of the covariance (yretHretx^)(yretHretx^)T¯ (in red) and the distribution of values (mean and standard deviation) for (ΛΛT + I)−1/2 from Eq. (17) (in black) and (c),(d) the square root of the diagonal values of the covariance (yretHretxb)(yretHretxb)T¯ (in red) and the distribution of values (mean and standard deviation) for (I + ΛΛT)1/2 from Eq. (17) (in black) for the (left) 2100 UTC assimilation (135 IASI FOVs) and (right) 0000 UTC assimilation (91 CrIS FOVs).

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 Jan(x) (Migliorini 2012):

 
Jan(x)=12(xxb)TB*1(xxb)+12(yretHretx)T(yretHretx),
(19)

and x^MAPan can be expressed as

 
x^MAPan=xb*+Kret*(yretHretxb*),
(20)

where

 
Kret*=B*H˜retT(H˜retB*H˜retT+Ir)1,
(21)

with B* and xb* being the prior information used to constrain yret. 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* ≠ pret; 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 Bret and xbret, which means that H˜ret=HretHconvHinterp accounts for the proper vertical/horizontal interpolation through the operator Hinterp and for the proper unit conversion through the operator Hconv. By retaining the diagonal terms and by removing the off-diagonal terms of B* in the characterization of Bret, all of the cross-correlation terms have been removed, and Bret 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 in Eq. (12) as if they were separate channels from a new instrument called TR.

The implemented observation operator H˜ret computes the model equivalent (simulated data) for every observation; H˜ret consists of three components. The first component is the three-dimensional interpolator Hinterp 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˜ret is the unit converter Hconv, 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˜ret is the projection operator Hret 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(i) vector, the observation minus background variance is (1+λi2), and since this value changes from element to element in the observation vector yret and it also changes from retrieval to retrieval, it seems reasonable to use (1 + λ2)1/2 to scale the ith components of the TRs υret(i)=yret(i)/1+λi2 and the corresponding row of observation operator Hret=Hret/1+λ2. Equation (18) then can be written

 
(υretHretxb)(υretHretxb)T¯=Ir.
(22)

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 (υretHretxb) 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 and Hret, 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.

Fig. 11.

Horizontal cross section of Θ innovation (°C) for model level numbers (a),(b) 1 (≃1015 hPa) and (c),(d) 10 (≃850 hPa) for (left) 2100 UTC assimilation (IASI) and (right) 0000 UTC assimilation (CrIS).

Fig. 11.

Horizontal cross section of Θ innovation (°C) for model level numbers (a),(b) 1 (≃1015 hPa) and (c),(d) 10 (≃850 hPa) for (left) 2100 UTC assimilation (IASI) and (right) 0000 UTC assimilation (CrIS).

Fig. 12.

As in Fig. 11, but for horizontal cross section of Q innovation (g kg−1).

Fig. 12.

As in Fig. 11, but for horizontal cross section of Q innovation (g kg−1).

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.

Fig. 13.

Vertical cross sections of analysis increments: (a),(b) Θ innovation (°C) and (c),(d) Q innovation (g kg−1) for the (left) 2100 UTC assimilation (135 IASI FOVs) and (right) 0000 UTC assimilation (91 CrIS FOVs).

Fig. 13.

Vertical cross sections of analysis increments: (a),(b) Θ innovation (°C) and (c),(d) Q innovation (g kg−1) for the (left) 2100 UTC assimilation (135 IASI FOVs) and (right) 0000 UTC assimilation (91 CrIS FOVs).

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.

Fig. 14.

As in Fig. 11, but for conventional observations only.

Fig. 14.

As in Fig. 11, but for conventional observations only.

Fig. 15.

As in Fig. 12, but for conventional observations only.

Fig. 15.

As in Fig. 12, but for conventional observations only.

Fig. 16.

As in Fig. 13, but for conventional observations only.

Fig. 16.

As in Fig. 13, but for conventional observations only.

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 xt was identical to the background or the analysis profile. The differences between the radiances thus obtained and the one actually observed generate the residuals ρbg=yrad0H(xb*) and ρan=yrad0H(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.

Fig. 17.

IASI: Comparison of the (a) longwave and (b) midwave spectral residuals generated for the background (blue dots), for the analysis (red dots). The black curve represents the observation error. Averages were done over the 69 FOVs candidates for the 0000 UTC assimilation.

Fig. 17.

IASI: Comparison of the (a) longwave and (b) midwave spectral residuals generated for the background (blue dots), for the analysis (red dots). The black curve represents the observation error. Averages were done over the 69 FOVs candidates for the 0000 UTC assimilation.

Fig. 18.

As in Fig. 17, but for CrIS, with averages being done over the 91 FOVs candidates.

Fig. 18.

As in Fig. 17, but for CrIS, with averages being done over the 91 FOVs candidates.

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.

Fig. 19.

RMS of the 3-h accumulated precipitation forecast calculated for each of the eighteen 60-h forecasts for the OP (red line) and CNTR (blue line) experiments.

Fig. 19.

RMS of the 3-h accumulated precipitation forecast calculated for each of the eighteen 60-h forecasts for the OP (red line) and CNTR (blue line) experiments.

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):

 
(yretHretx^MAP)(yretHretx^MAP)T.
(A1)

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):

 
yretHretx^MAP=Hretxt+εyretHret[xbret+KH(xtxbret)+Kεrad],
(A2)

which can be further manipulated to obtain

 
yretHretx^MAP=Hret(IKH)(xtxbret)+(UTR1/2HretK)εy.
(A3)

To avoid the use of subscripts, and yret being of dimension r, in the following demonstrations U, V, and Λ will represent Ur, Vr, and Λr. Equations (A2) and (A3) then become

 
yretHretx^MAP=Hret[(IKH)(xtxbret)+(B1/2VΛ1UTR1/2K)εy].
(A4)

The second term of (A4) can be rewritten as

 
B1/2VΛ1UTR1/2K=B1/2VΛ1UTR1/2B1/2VΛ(ΛΛT+I)1UTR1/2=B1/2VΛ1[IΛΛT(ΛΛT+I)1]UTR1/2=B1/2VΛ1(ΛΛT+I)1UTR1/2.
(A5)

so that Eq. (A1) can be rewritten as

 
Hret[(IKH)B(IHTKT)+B1/2VΛ1(ΛΛT+I)2Λ1VTB1/2]HretT.
(A6)

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

 
(IKH)B(IHTKT)+B1/2VΛ1(ΛΛT+I)2Λ1VTB1/2=(BKHB)(IHTKT)+B1/2VΛ1(ΛΛT+I)2Λ1VTB1/2=BKHBBHTKT+KHBHTKT+B1/2VΛ1(ΛΛT+I)2Λ1VTB1/2=BB1/2VΛ(ΛΛT+I)I1UTR1/2R1/2UΛVTB1/2BBB1/2VΛUTR1/2R1/2U(ΛΛT+I)1ΛVTB1/2+B1/2VΛ(ΛΛT+I)1UTR1/2R1/2UΛVTB1/2BB1/2VΛUTR1/2R1/2U(ΛΛT+I)1ΛVTB1/2+B1/2VΛ1(ΛΛT+I)2Λ1VTB1/2=B2B1/2VΛ(ΛΛT+I)1ΛVTB1/2+B1/2VΛ(ΛΛT+I)1ΛΛ(ΛΛT+I)1ΛVTB1/2+B1/2VΛ1(ΛΛT+I)2Λ1VTB1/2=B1/2V(ΛΛT+I)1[(IΛΛT)(I+ΛΛT)+(ΛΛT)2+(ΛΛT)1](ΛΛT+I)1VTB1/2=B1/2V(ΛΛT+I)1(I+Λ2)(ΛΛT+I)1VTB1/2;
(A7)

thus

 
(yretHretx^MAP)(yretHretx^MAP)T=Hret[(IKH)B(IHTKT)+B1/2VΛ1(ΛΛT+I)2Λ1VTB1/2]HretT=Λ(ΛΛT+I)1[I+(ΛΛT)1](ΛΛT+I)1ΛT=(ΛΛT+I)1.
(A8)

APPENDIX B

Derivation of Eq. (18)

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

 
(yretHretxbret)(yretHretxbret)T.
(B1)

As suggested by de Haan, the demonstration of Eq. (18) starts from a manipulation of Eq. (16):

 
yretHretxbret=yretHretxt+Hret(xtxbret)=εrad+Hret(xtxbret),
(B2)

which, considering εrad to be uncorrelated with εb=(xtxbret)(xtxbret)T, implies

 
yretHretxbret=yretHretxt+Hret(xtxbret)=[εrad+Hret(xtxbret)][εrad+Hret(xtxbret)]T=εradεradT+HretBHT=I+HretBHT.
(B3)

Using the alternative expression in Eq. (24) of Migliorini et al. (2008), Hret=ΛrVrTB1/2, it is possible to write

 
(yretHretxbret)(yretHretxbret)T=I+(ΛrVTB1/2)B(ΛrVTB1/2)T=I+ΛrΛrT.
(B4)

REFERENCES

REFERENCES
Andersson
,
E.
,
J.
Pailleux
,
J.-N.
Thépaut
,
J. R.
Eyre
,
A. P.
McNally
,
G. A.
Kelly
, and
P.
Courtier
,
1994
:
Use of cloud-cleared radiances in three/four-dimensional variational data assimilation
.
Quart. J. Roy. Meteor. Soc.
,
120
,
627
653
, https://doi.org/10.1002/qj.49712051707.
Antonelli
,
P.
, and et al
,
2004
:
A principal component noise filter for high spectral resolution infrared measurements
.
J. Geophys. Res.
,
109
,
D23102
, https://doi.org/10.1029/2004jd004862.
Antonelli
,
P.
, and et al
,
2017
:
Regional retrieval processor for direct broadcast high-resolution infrared data
.
J. Appl. Meteor. Climatol.
,
56
,
1681
1705
, https://doi.org/10.1175/JAMC-D-16-0144.1.
Barker
,
D. M.
,
W.
Huang
,
Y.-R.
Guo
,
A. J.
Bourgeois
, and
Q. N.
Xiao
,
2004
:
A three-dimensional variational data assimilation system for MM5: Implementation and initial results
.
Mon. Wea. Rev.
,
132
,
897
914
, https://doi.org/10.1175/1520-0493(2004)132<0897:ATVDAS>2.0.CO;2.
Barker
,
D. M.
, and et al
,
2012
:
The Weather Research and Forecasting Model’s Community Variational/Ensemble Data Assimilation System: WRFDA
.
Bull. Amer. Meteor. Soc.
,
93
,
831
843
, https://doi.org/10.1175/BAMS-D-11-00167.1.
Bloom
,
H.
,
2001
:
The Cross-track Infrared Sounder (CrIS): A sensor for operational meteorological remote sensing. Proc. 2001 Int. Geoscience and Remote Sensing Symp., Sydney, NSW, Australia, IEEE, 1341–1343
, https://doi.org/10.1109/IGARSS.2001.976838.
Businger
,
S.
,
R.
McLaren
,
R.
Ogasawara
,
D.
Simons
, and
R. J.
Wainscoat
,
2002
:
Starcasting
.
Bull. Amer. Meteor. Soc.
,
83
,
858
871
, https://doi.org/10.1175/1520-0477(2002)083<0858:S>2.3.CO;2.
Chahine
,
M. T.
, and et al
,
2006
:
AIRS: Improving weather forecasting and providing new data on greenhouse gases
.
Bull. Amer. Meteor. Soc.
,
87
,
911
926
, https://doi.org/10.1175/BAMS-87-7-911.
Chambers
,
C. R. S.
,
2003
:
High resolution mesoscale modelling of kauai wintertime weather. M.S. thesis, Dept. of Meteorology, University of Hawai‘i at Mānoa, 104 pp
.
Cherubini
,
T.
,
S.
Businger
, and
R.
Lyman
,
2011
:
An operational perspective for modeling optical turbulence. Seeing Clearly: The Impact of Atmospheric Turbulence on the Propagation of Extraterrestrial Radiation, S. Businger and T. Cherubini, Eds., VBW Publishing, 165–182
.
Cho
,
C.
, and
D. H.
Staelin
,
2006
:
Cloud clearing of Atmospheric Infrared Sounder hyperspectral infrared radiances using stochastic methods
.
J. Geophys. Res.
,
111
,
D09S18
, https://doi.org/10.1029/2005JD006013.
de Haan
,
S.
,
S.
Barlag
,
H. K.
Baltink
,
F.
Debie
, and
H.
van der Marel
,
2004
:
Synergetic use of GPS water vapor and meteosat images for synoptic weather forecasting
.
J. Appl. Meteor.
,
43
,
514
518
, https://doi.org/10.1175/1520-0450(2004)043<0514:SUOGWV>2.0.CO;2.
Derber
,
J. C.
, and
W.-S.
Wu
,
1998
:
The use of tovs cloud-cleared radiances in the NCEP SSI analysis system
.
Mon. Wea. Rev.
,
126
,
2287
2299
, https://doi.org/10.1175/1520-0493(1998)126<2287:TUOTCC>2.0.CO;2.
Desroziers
,
G.
,
L.
Berre
,
B.
Chapnik
, and
P.
Poli
,
2005
:
Diagnosis of observation, background and analysis-error statistics in observation space
.
Quart. J. Roy. Meteor. Soc.
,
131
,
3385
3396
, https://doi.org/10.1256/qj.05.108.
Eyre
,
J. R.
,
2007
:
Progress achieved on assimilation of satellite data in numerical weather prediction over the last 30 years. Seminar on Recent Developments in the Use of Satellite Observations in Numerical Weather Prediction, Shinfield Park, Reading, United Kingdom, ECMWF
, https://www.ecmwf.int/sites/default/files/elibrary/2008/9341-progress-achieved-assimilation-satellite-data-numerical-weather-prediction-over-last-30-years.pdf.
Eyre
,
J. R.
,
2016
:
Observation bias correction schemes in data assimilation systems: A theoretical study of some of their properties
.
Quart. J. Roy. Meteor. Soc.
,
142
,
2284
2291
, https://doi.org/10.1002/qj.2819.
Eyre
,
J. R.
,
S. J.
English
, and
M.
Forsythe
,
2019
:
Assimilation of satellite data in numerical weather prediction. Part I: The early years
.
Quart. J. Roy. Meteor. Soc.
,
146
,
49
68
, https://doi.org/10.1002/qj.3654.
Gambacorta
,
A.
, and et al
,
2014
:
The NOAA Unique CrIS/ATMS Processing System (NUCAPS): Algorithm description and validation results after two years in orbit. 10th Annual Symp. on New Generation Operational Environmental Satellite Systems, Atlanta, GA, Amer. Meteor. Soc., 9.1
, https://ams.confex.com/ams/94Annual/webprogram/Paper236026.html.
Joiner
,
J.
, and
A.
Da Silva
,
1998
:
Efficient methods to assimilate remotely sensed data based on information content
.
Quart. J. Roy. Meteor. Soc.
,
124
,
1669
1694
, https://doi.org/10.1002/qj.49712454915.
Karagulian
,
F.
,
L.
Clarisse
,
C.
Clerbaux
,
A. J.
Prata
,
D.
Hurtmans
, and
P. F.
Coheur
,
2010
:
Detection of volcanic SO2, ash, and H2SO4 using the Infrared Atmospheric Sounding Interferometer (IASI)
.
J. Geophys. Res.
,
115
,
D00L02
, https://doi.org/10.1029/2009JD012786.
Kelly
,
G. A. M.
,
G. A.
Mills
, and
W. L.
Smith
,
1978
:
Impact of Nimbus-6 temperature soundings on Australian region forecasts
.
Bull. Amer. Meteor. Soc.
,
59
,
393
406
, https://doi.org/10.1175/1520-0477-59.4.393.
Klemp
,
J.
,
W. C.
Skamarock
, and
J.
Dudhia
,
2007
:
Conservative split-explicit time integration methods for the compressible non hydrostatic equations
.
Mon. Wea. Rev.
,
135
,
2897
2913
, https://doi.org/10.1175/MWR3440.1.
Lerner
,
J. A.
,
E.
Weisz
, and
G.
Kirchengast
,
2002
:
Temperature and humidity retrieval from simulated Infrared Atmospheric Sounding Interferometer (IASI) measurements
.
J. Geophys. Res.
,
107
,
4189
, https://doi.org/10.1029/2001jd900254.
Merchant
,
C.
,
S.
Saux-Picart
, and
J.
Waller
,
2020
:
Bias correction and covariance parameters for optimal estimation by exploiting matched in-situ references
.
Remote Sens. Environ.
,
237
, 111590, https://doi.org/10.1016/j.rse.2019.111590.
Migliorini
,
S.
,
2012
:
On the equivalence between radiance and retrieval assimilation
.
Mon. Wea. Rev.
,
140
,
258
265
, https://doi.org/10.1175/MWR-D-10-05047.1.
Migliorini
,
S.
,
C.
Piccolo
, and
C.
Rodgers
,
2008
:
Use of the information content in satellite measurements for an efficient interface to data assimilation
.
Mon. Wea. Rev.
,
136
,
2633
2650
, https://doi.org/10.1175/2007MWR2236.1.
Moncet
,
J.-L.
,
G.
Uymin
,
A. E.
Lipton
, and
H. E.
Snell
,
2008
:
Infrared radiance modeling by optimal spectral sampling
.
J. Atmos. Sci.
,
65
,
3917
3934
, https://doi.org/10.1175/2008JAS2711.1.
Parrish
,
D. F.
, and
J. C.
Derber
,
1992
:
The National Meteorological Center’s spectral statistical-interpolation analysis system
.
Mon. Wea. Rev.
,
120
,
1747
1763
, https://doi.org/10.1175/1520-0493(1992)120<1747:TNMCSS>2.0.CO;2.
Prates
,
C.
,
S.
Migliorini
,
L.
Stewart
, and
J.
Eyre
,
2016
:
Assimilation of transformed retrievals obtained from clear-sky IASI measurements
.
Quart. J. Roy. Meteor. Soc.
,
142
,
1697
1712
, https://doi.org/10.1002/qj.2764.
Rodgers
,
C.
,
2000
:
Inverse Methods for Atmospheric Soundings: Theory and Practice
.
World Scientific
,
238
pp.
Strow
,
L. L.
, and et al
,
2013
:
Spectral calibration and validation of the Cross-Track Infrared Sounder on the Suomi NPP satellite
.
J. Geophys. Res. Atmos.
,
118
,
12 486
12 496
, https://doi.org/10.1002/2013JD020480.
Yang
,
J.
,
Z.
Zhang
,
C.
Wei
,
F.
Lu
, and
Q.
Guo
,
2017
:
Introducing the new generation of Chinese geostationary weather satellites, Fengyun-4
.
Bull. Amer. Meteor. Soc.
,
98
,
1637
1658
, https://doi.org/10.1175/BAMS-D-16-0065.1.
Zhang
,
P.
, and et al
,
2019
:
Latest progress of the Chinese meteorological satellite program and core data processing technologies
.
Adv. Atmos. Sci.
,
36
,
1027
1045
, https://doi.org/10.1007/s00376-019-8215-x.
For information regarding reuse of this content and general copyright information, consult the AMS Copyright Policy (www.ametsoc.org/PUBSReuseLicenses).