# The capability of 4D-Var systems to assimilate cloud-affected satellite infrared radiances

Frederic Chevallier, Philippe Lopez, Adrian Tompkins, Marta Janisková, Emmanuel Moreau

► To cite this version:

Frederic Chevallier, Philippe Lopez, Adrian Tompkins, Marta Janisková, Emmanuel Moreau. The capability of 4D-Var systems to assimilate cloud-affected satellite infrared radiances. Quarterly Journal of the Royal Meteorological Society, 2004, 130 (598), pp.917-932. 10.1256/QJ.03.113. hal-02947749

HAL Id: hal-02947749

https://hal.science/hal-02947749v1

Submitted on 23 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# The capability of 4D-Var systems to assimilate cloud-affected satellite infrared radiances

By FRÉDÉRIC CHEVALLIER*, PHILIPPE LOPEZ, ADRIAN M. TOMPKINS, MARTA JANISKOVÁ and EMMANUEL MOREAU

*European Centre for Medium-Range Weather Forecasts, UK*

*(Received 30 June 2003; revised 27 October 2003)*

## SUMMARY

Four-dimensional variational (4D-Var) assimilation schemes assume the linearity of their forward model in the vicinity of prior information and usually do not properly handle variables that have finer temporal and spatial scales in the real world than in the forward model. Hence cloud-affected satellite infrared radiances are discarded from numerical weather-prediction 4D-Var systems, despite the critical need of observations within the cloudy regions. This paper suggests the reappraisal of that choice, subject to achieving improvements in the numerical simulation of cloudiness.

A new observation operator, that computes cloud-affected infrared radiances from 4D-Var control variables, namely atmospheric temperature, humidity, ozone, surface temperature and surface pressure, is presented. The vertical distributions of cloud cover and of cloud condensate are diagnosed in the operator itself. The goal of this paper is to assess the feasibility of using it to assimilate cloud-affected infrared radiances, such as those from the narrow-band Advanced Infrared Sounder on-board the Aqua platform or those from the broad-band Meteosat Visible and Infrared Imager. It is shown that there is a potential benefit in assimilating directly in 4D-Var some of the upper-tropospheric channels at 4.5, 6.3 and 14.3 µm in the presence of clouds, for instance the 6.3 µm channel on board all the geostationary satellites. The approach is illustrated with one-dimensional variational retrievals collocated with radiosonde observations.

KEYWORDS: Clouds Data assimilation Satellite observations

## 1. INTRODUCTION

The improvement of weather forecast skill in recent years owes much to the development of Bayesian estimation techniques for atmospheric data assimilation. In particular, an increasing number of numerical-weather-prediction (NWP) centres opt for three- and four-dimensional variational assimilation systems (respectively, 3D-Var and 4D-Var) to perform their atmospheric analyses (e.g. Rabier et al. 2000; Lorenc et al. 2000). 3D- and 4D-Var systems estimate the atmospheric variables from a background state (usually based on an earlier forecast) and from available observations. They perform best when the background and observation-error statistics are Gaussian, unbiased and perfectly known, and when the forward model, which relates the background variables to the observed variables, is linear. Linearity is actually only needed in the vicinity of the background for perturbations whose magnitude is of the same order as the background errors. Yet observations for which significant nonlinearities affect the forward model are discarded from 3D- and 4D-Var systems. Consequently, infrared satellite radiances are currently not assimilated in the presence of clouds, even though they would inform the NWP systems about regions of the atmosphere that strongly influence the forecasts (McNally 2002).

At the European Centre for Medium-Range Weather Forecasts (ECMWF) an observation operator that computes cloud-affected radiances from some of the ECMWF 4D-Var control variables (temperature, humidity and ozone profiles, surface temperature and surface pressure) has been developed for data assimilation. The operator diagnoses the vertical distributions of cloud cover and of cloud condensate by taking both large-scale and convective processes into consideration. The goal of this paper is to assess the

<footnote>* Corresponding author, present address: LSCE CE L'Orme des Merisiers, 91191 Gif sur Yvette, Cedex, France. e-mail: frederic.chevallier@cea.fr</footnote>

<footnote>© Royal Meteorological Society, 2004.</footnote>
possibility of its use for assimilating cloud-affected radiances within a 4D-Var system. The assessment is based on the examination of the accuracy and of the linearity of this new observation operator for the simulation of the narrow-band channels from the Advanced InfraRed Sounder (AIRS) on-board the Aqua platform and of the broad-band channels from the Meteosat Visible and InfraRed Imager (MVIRI).

The plan of the paper is as follows. Section 2 outlines the 4D-Var method and discusses the implications for the selection of observations. Section 3 describes the data and the observation operator. Accuracy and linearity of the observation operator are estimated in sections 4 and 5 respectively. A one-dimensional variational (1D-Var) scheme is used in section 6 to illustrate the previous results. Concluding discussion follows in section 7.

## 2. THEORETICAL FRAMEWORK OF 4D-VAR

### (a) General formalism

4D-Var systems are a practical formulation of Bayesian estimation theory for the particular case of a linear problem with unbiased Gaussian errors. They seek a model trajectory, $\mathbf{x}(t)$, that is statistically consistent with the information provided by the observations, $\mathbf{y}$, available during a given time period $[t_0, t_n]$ and that provided by an a-priori model state, $\mathbf{x}^b$, called the background state. The background, $\mathbf{x}^b$, is usually based on a short-range forecast. The model trajectory, $\mathbf{x}(t)$, is completely defined by the initial state, $\mathbf{x}_0$, at time $t_0$.

The balance of the model initial state, $\mathbf{x}_0$, between the observations, $\mathbf{y}$, and the background, $\mathbf{x}^b$, is measured by an objective cost-function defined as follows:

$$J(\mathbf{x}_0) = \frac{1}{2}(\mathbf{x}_0 - \mathbf{x}_0^b)^T \mathbf{B}^{-1} (\mathbf{x}_0 - \mathbf{x}_0^b) + \frac{1}{2} \sum_i (\psi_i[\mathbf{x}_0] - \mathbf{y}_i)^T \mathbf{R}_i^{-1} (\psi_i[\mathbf{x}_0] - \mathbf{y}_i), \quad (1)$$

where, at any time $t_i$, $\mathbf{y}_i$ is the vector of observations, $\psi_i$ is the forward operator providing the equivalent of the observations from the model initial state $\mathbf{x}_0$, $\mathbf{R}_i$ is the observation-error covariance matrix and $\mathbf{B}$ is the background-error covariance matrix of the state $\mathbf{x}^b$. Superscripts $-1$ and $T$ denote, respectively, the matrix inverse and transpose. Subscript $i$ denotes the time index.

$\psi_i$ is the combination of a dynamic *model operator*, usually called $\mathbf{M}_i$, and of a static *observation operator*, usually called $\mathbf{H}_i$ (Ide et al. 1997). Observation errors described by $\mathbf{R}_i$ are assumed to be uncorrelated from one time step $i$ to another. It is important to understand that they are defined with respect to the forward model $\psi_i$, which is assumed to be perfect. Consequently, $\mathbf{R}_i$ actually contains the errors of $\psi_i$ in addition to the measurement errors.

The control vector, $\mathbf{x}_0$, includes the prognostic variables to be initialized in the forecast model: vorticity, divergence, temperature, humidity, ozone, surface temperature and surface pressure in the case of the ECMWF 4D-Var. Surface temperature is actually defined at satellite observation points only, since soil variables are specifically analysed in a separate optimal interpolation system. In theory, cloud cover and cloud condensate could be 4D-Var control variables as well, but the definition of their background errors and of their coupling with temperature and humidity background errors would be particularly challenging, for instance when cloud cover is zero or unity.

The minimization uses a descent algorithm which requires several computations of the gradient of $J$ with respect to the initial state $\mathbf{x}_0$. Given the dimension of the state
vector, the adjoint technique is used to provide an efficient estimate of $\nabla J$:

$$ \nabla J(\mathbf{x}_0) = \mathbf{B}^{-1}(\mathbf{x}_0 - \mathbf{x}_0^b) + \sum_i \Psi_i^T \mathbf{R}_i^{-1} (\psi_i[\mathbf{x}_0] - \mathbf{y}_i), \quad (2) $$

where $\Psi_i^T$ is the adjoint of the forward operator (i.e. the derivative or Jacobian matrix).

Additionally, an incremental formulation of the variational system is used at ECMWF (Courtier et al. 1994): departures $\psi_i[\mathbf{x}_0] - \mathbf{y}_i$ are computed at resolution 40 km, the analysis increments are computed at the lower resolution of 125 km. The system copes with weak nonlinearities through an ‘inner-loop/outer-loop algorithm’, where the cost-function is minimized in the inner loop based on a linearization of $\psi_i$ in the vicinity of the initial state which is updated in the outer loop. A simplified version of the derivatives of $\psi_i$ are used both for the linearization (tangent-linear model $\Psi_i$) and for the adjoint model $\Psi_i^T$. In the following, cycle 25r5 of the ECMWF forecasting system, which became operational in March 2003, is used.

### (b) Implications for the selection of cloud-affected observations

The variational formulation of the inverse problem outlined above, together with currently available computer power, allows the operational handling of large numbers of control variables (about five million currently at ECMWF) and of observations (about 1.5 million per 12-hour analysis cycle). It would provide statistically optimal analyses if the error statistics $\mathbf{B}$ and $\mathbf{R}_i$ were unbiased, Gaussian and perfectly known, and if the problem was linear. For instance, significant nonlinearities may exist in a 4D-Var system, but they degrade the realism of the corresponding analyses and tend to limit the impact in the subsequent forecasts to short ranges. As a consequence, attempts are made to bring the 4D-Var systems as close as possible to optimality by the removal of biases (e.g. Harris and Kelly 2001) by choosing Gaussian-error control variables (e.g. Dee and da Silva 2003), by a careful estimation of the error statistics (e.g. Derber and Bouttier 1999), by the improvement of the parametrizations of the forward operator, and by avoiding observations for which $\psi_i$ is significantly nonlinear with respect to the analysis increments $\mathbf{x}_0 - \mathbf{x}_0^b$.

To account for cloud processes in such a framework is obviously a challenge. Indeed, fine-scale atmospheric processes significantly impact the cloud fields and result in significant nonlinearities at the spatial and temporal scales of the NWP models. Further, they make cloud parametrizations particularly difficult to formulate. Conversely, cloud observations also contain large-scale information, through the dynamics, which allows a realistic representation of cloud systems in NWP (e.g. Chevallier and Kelly 2002). In the next sections the accuracy of the recently developed observation operator, and its linearity with respect to perturbations of the order of the expected size of analysis increments, is estimated in order to select cloud-affected channels for assimilation within 4D-Var.

## 3. DATA AND MODEL

### (a) The satellite data

The present study exploits the observations from two space-borne instruments, MVIRI and AIRS.

MVIRI has been flown on board the Meteosat 1 to 7 geostationary satellites, operated by the European organization for the exploitation of meteorological satellites (EUMETSAT). It includes three spectrally-broad measurement channels: two infrared
ones (from 10.5 to 12.5 µm and from 5.7 to 7.1 µm, respectively) and one visible (0.45 to 1 µm). An image of the earth disc as seen from the geostationary orbit is generated at 30 min intervals in each channel. It consists of 2500 × 2500 pixels for the infrared channels (5 × 5 km² resolution at the subsatellite point) and 5000 × 5000 pixels (2.5 × 2.5 km²) for the visible one. One Meteosat satellite is normally stationed close to 0° longitude but, depending on operational choice, an additional one may be positioned above the Atlantic or the Indian Ocean.

For the present study, the radiances from the two infrared Meteosat channels are extracted from the EUMETSAT Radiance Data from Clouds (RDC) product in the form of mean brightness temperatures in 16 × 16 pixel quadrants. This averaging results in a resolution (80 km at the subsatellite point) similar to the ECMWF analysis (currently about 125 km). It is useful to note that data from a related EUMETSAT product, the Clear Sky Radiances (CSR), are operationally assimilated at ECMWF (Köpken *et al.* 2004). Those radiances result from averages for the pixels not affected by clouds in each quadrant. Only data from the quadrants where more than 70% of the pixels are judged clear are actually assimilated. Consistently, quadrants will be considered cloudy in the following when more than 30% of the pixels are not classified as clear.

The AIRS instrument, operated by the National Aeronautics and Space Administration (NASA), provides significantly different measurements of the infrared spectrum. On board the Aqua sun-synchronous polar orbiter, it observes nearly all points of the globe twice a day, moving northward across the equator at about 01:30 pm local time. It samples the infrared spectrum between 3.7 and 15.4 µm with 2378 channels. Additionally four channels are located in the visible (from 0.4 to 1.0 µm). Horizontal resolution reaches 13.5 km and 2.3 km at nadir for the infrared and the visible channels, respectively. No attempt is made to average the data. As a starting point, a subset of 324 channels for one satellite spot in eighteen has been made operationally available to ECMWF by the National Environment Satellite Data and Information Service (NESDIS). Data from cloud-free channels are currently assimilated in a research mode, with a cloud detection method described by McNally and Watts (2003).

AIRS observations for wave numbers below 2000 cm⁻¹ and MVRI ones are bias-corrected using a constant offset in each channel in order to account for possible erroneous knowledge of the instrument characteristics. The offset is estimated independently from cloud-free departure statistics (Köpken and McNally 2003, personal communication).

### (b) Observation operator

(i) *Diagnostic and prognostic models.* A prognostic model implies a scheme that computes the tendencies $\partial\eta/\partial t$ of some cloud quantity $\eta$ with respect to time $t$, therefore retaining cloud information from previous time steps of the integration. A diagnostic model, alternatively, diagnoses the state of $\eta$ at time $t$ from other variables, thus precluding a memory of cloud variables, and thus implying that cloud mass is not necessarily conserved.

The ECMWF forecast model includes a prognostic cloud scheme (Tiedtke 1989, 1993) which could, in principle, be used in the 4D-Var physics for the assimilation of cloud information, provided that cloud variables are added in the 4D-Var control vector. Such a strategy would pose the acute problem of defining background-error statistics for the new variables, as mentioned in section 2(a). As a consequence, a diagnostic approach has been preferred. A model has been developed and is described hereafter. It uses only the existing 4D-Var control variables as input and is kept relatively simple
so that thresholds and strong nonlinearities do not make the 4D-Var minimization stop before reaching the absolute minimum of the cost function.

The present study does not cover the full 4D-Var forward model. It focuses mainly on its static part, or observation operator **H** (see section 2(a)), that computes satellite radiances from the relevant 4D-Var control variables, namely atmospheric temperature, humidity, ozone, surface temperature and surface pressure. Indeed, the dynamic part, or model operator **M**, computes the time evolution of all the control variables and is not specific to the assimilation of radiances in the presence of clouds. The reader is referred to Janisková *et al.* (2002) for the impact of the introduction of cloud–radiation interaction within the model operator.

(ii) Convection model. The convection model of the observation operator is based on the ECMWF operational mass-flux scheme, initiated by Tiedtke (1989), but it uses a simplified algorithm. All types of convection (shallow, mid-level, and deep) are treated in a similar way. In particular, the link between the model control variables and the subgrid-scale convective quantities (the so-called closure assumption) is expressed through a single formulation that depends on the release of convective available potential energy in time. The equations that describe the vertical evolution of the updraught mass flux $M_u$ and of the updraught thermodynamic variables $\Phi_u$ are uncoupled:

$$ \frac{\partial M_u}{\partial z} = -(\epsilon - \delta)M_u \quad (3) $$

$$ \frac{\partial \Phi_u}{\partial z} = -\epsilon(\Phi_u - \bar{\Phi}), \quad (4) $$

where $\bar{\Phi}$ denotes field values in the environment, and $\epsilon$ and $\delta$ are the entrainment and detrainment fractional rates, respectively. The bulk convective updraught is assumed to originate from the surface only if its initial vertical velocity, as calculated from the surface heat fluxes, is positive. The departures of the updraught from the environment are also assumed to be dependent on the surface heat fluxes. If convection cannot be initiated from the surface, the convective ascent may originate from higher levels provided relative humidity exceeds 80%. In this case, the initial vertical velocity of the bulk updraught is set equal to 1 m s⁻¹. Regardless of whether the updraught originates from the surface or higher up, the vertical evolution of its kinetic energy is computed following Simpson and Wiggert (1969), and involves the entrainment of environmental air into the updraught. Convection is assumed to be active only if the updraught vertical velocity remains positive at cloud base. Simplified calculations of downdraughts and convective momentum transport, based on the operational scheme, are also included in the modified parametrization. The calculation of the precipitation formation rate from the updraught cloud condensate specific content is inspired from the work of Tiedtke (1989). Further details of the convection scheme have been given by Moreau *et al.* (2003).

(iii) Large-scale cloud model. As described by Tompkins and Janisková (2003), the stratiform part of the diagnostic cloud scheme assumes that the subgrid-scale fluctuations of total water are uniformly distributed. Two parameters only are, therefore, needed to describe the corresponding probability density function (PDF) in any model layer i: $RH_i^c$ and $\kappa_i$. $RH_i^c$ is the critical relative humidity above which clouds are allowed to form, and $\kappa_i$ sets the total water variance or, equivalently, the distribution
width. Stratiform cloud cover in layer *i* is consequently expressed as:

$$C_i^{\text{strat}} = 1 - \sqrt{\frac{1 - R H_i}{1 - R H_i^c - \kappa_i (R H_i - R H_i^c)}}, \quad (5)$$

where $R H_i$ is the relative humidity in layer *i*. $R H_i^c$ and $\kappa_i$ are empirically defined in each layer once for all. Initially the $R H_i^c$'s were estimated using $\kappa_i = 0.7$ so that the diagnostic model fitted the prognostic scheme (Tiedtke 1993) well. Values were obtained for $R H_i^c$ of 0.55 in the middle troposphere up to about 0.9 at the tropopause and at the surface. The coefficients $\kappa_i$ were then adjusted to agree with top-of-the-atmosphere radiation observations. In addition, precipitation fluxes were estimated in a manner close to the prognostic scheme, in order to pave the way for the assimilation of rain-affected microwave radiances, but this aspect is not exploited here.

(iv) *The radiance model.* The radiative-transfer model to compute satellite radiances has already been documented within previous studies (Chevallier *et al.* 2001, 2002; Chevallier and Kelly 2002), and only its salient features are recalled here. The scheme applies in the infrared as in the microwave spectrum, but does not take scattering into account. Cloud absorption is computed in each model layer using parametrizations similar to those used in the ECMWF forecast model for the broad-band radiative fluxes. Cloud layers are assumed to overlap according to the maximum-random scheme of Raïsänen (1998). The Radiative Transfer for the TIROS* Operational Vertical Sounder (RTTOV) (Eyre 1991; Saunders *et al.* 2002) computes the gaseous absorption. The model has been used since 2000 for the routine monitoring of the clouds modelled by the ECMWF operational forecasting system, and in the 40-year reanalysis programme.

# 4. ACCURACY OF THE OBSERVATION OPERATOR

## (a) MVIRI

The accuracy of the observation operator is illustrated for the Meteosat-7 data for 12 UTC 30 November 2002. The atmospheric circulation is very similar to that described by Chevallier and Kelly (2002) for December 2000. Various cloud systems are present in the 11 µm full-resolution image (Fig. 1): several frontal systems in mid-latitudes, and tropical convection that is as active over the oceans as over the land. Inputs to the observation operator are extracted from the 12-hour forecast initialized from the 00 UTC analysis. Cloud variables are diagnosed by the moist physics described above, but can also be simply extracted from the forecast archive. In the latter case, the cloud variables are computed by the prognostic cloud scheme. The same radiation model is used in both cases. Although the forecast resolution reaches 40 km, a reduced grid of 80 km is used here for better compatibility with the 4D-Var system.

Statistics of the model versus the observations are presented in Table 1 for the two infrared channels in the Meteosat cloudy quadrants. From this simple test, the diagnostic and the prognostic physics appear to perform similarly with respect to observations, despite significant standard deviations (9.6 K at 11 µm) of the differences between each other. This result is particularly important for 4D-Var applications, where the diagnostic model would be used in place of the prognostic one (see section 3(b)(i)). The largest departures from the observations occur around 11 µm. Indeed in this atmospheric spectral window little gaseous absorption takes place, so that any cloud above the boundary layer has a strong impact on this channel. Comparatively, the 6.3 µm channel

<footnote>* Television and Infrared Observation Satellite.</footnote>
<figure><img src="image_3.png" /><figcaption>Figure 1. Observed Meteosat-7 11 µm image at 12 UTC 30 November 2002.</figcaption></figure>

<table><caption>TABLE 1. MEAN AND STANDARD DEVIATION (Std) OF THE OBSERVATIONS (O) AND OF THE DIFFERENCES BETWEEN THE TWO MODEL BRIGHTNESS TEMPERATURES (D-P) OR BETWEEN MODELLED AND OBSERVED BRIGHTNESS TEMPERATURES (D-O) AND (P-O)</caption><thead><tr><th rowspan="2">MVIRI channel<br>(µm)</th><th colspan="2">O</th><th colspan="2">D-P</th><th colspan="2">D-O</th><th colspan="2">P-O</th></tr><tr><th>Mean<br>(K)</th><th>Std<br>(K)</th><th>Mean<br>(K)</th><th>Std<br>(K)</th><th>Mean<br>(K)</th><th>Std<br>(K)</th><th>Mean<br>(K)</th><th>Std<br>(K)</th></tr></thead><tbody><tr><td>6.3</td><td>234.1</td><td>7.9</td><td>-0.2</td><td>2.6</td><td>-1.9</td><td>5.8</td><td>-2.0</td><td>5.7</td></tr><tr><td>11</td><td>260.9</td><td>18.4</td><td>-1.5</td><td>9.6</td><td>8.2</td><td>18.0</td><td>9.7</td><td>16.5</td></tr></tbody></table>

The model cloud variables are either prognostic (P) or diagnosed (D). Statistics are presented for the two MVIRI infrared channels and are restricted to the Meteosat quadrants where less than 70% of the pixels are judged clear (10770 quadrants out of 22718).
<figure><img src="image_4.png" /><figcaption>Figure 2. Probability density function (PDF) of the departures of the diagnostic-model values (D) from the MVIRI observations (O) for (a) 6.3 µm and (b) 11 µm brightness temperatures in the Meteosat-7 cloudy quadrants at 12 UTC 30 November 2002. The Gaussian distributions with the same means and standard deviations are also reported on the graphs.</figcaption></figure>

is very sensitive to upper-tropospheric humidity (Schmetz and Turpeinen 1988). It is, therefore, hardly affected by low-level clouds and less affected by upper- and mid-level clouds. Table 1 also indicates the existence of biases between the model and the observations that should be removed before those observations are assimilated.

The fraction of observed variance explained by the model is defined as $f = [\text{Var}\{y\} - \text{Var}\{H(x) - y\}]/ \text{Var}\{y\}$, where $\text{Var}\{y\}$ is the variance of the observations and $\text{Var}\{H(x) - y\}$ is the variance of the departures ($f = 1$ if the model is perfect and $f = 0$ if the model does not correlate with the observations). $f$ is equal to 0.48 and 0.20 at 6.3 and 11 µm, respectively, for the prognostic model and 0.46 and 0.04, respectively, for the diagnostic model. This means that both models have similar significant skill in predicting 80 km resolution 6.3 µm cloud optical thickness, but have nearly no significant skill at 11 µm, at least on average in the Meteosat disk. The deficiencies at 11 µm can be caused by the cloud scheme, by the radiation model, or by the temperature and humidity fields. One may note that the prognostic scheme is used in dynamic mode and is combined with other parametrizations in the forecast model to predict the temperature and humidity variations.

It would be interesting to investigate the PDF of the observation-operator errors, because they would dominate the observation-error covariance **R** (the time index *i* has been dropped and model operator is ignored) and could consequently be assumed to be Gaussian within 4D-Var (see section 2(a)). Two methods could be used to infer **R**. The first one would directly measure it with accurate validation data, including both the inputs and the outputs of the observation operator. Such data, for instance from the Atmospheric Radiation Measurement program (ARM, see Stokes and Schwartz 1994), do not exist in sufficient numbers to build reliable error statistics. The second method would infer **R** indirectly from the PDF of the departures of the model from the observations, which include both observation and background errors, and from accurate estimation of the background-error statistics **B**. Rough estimates of **B** exist and are used below, but do not reach a sufficient accuracy for this task. Therefore, the PDF of the departures in the cloudy quadrants is used here as a poor surrogate. They are shown in Fig. 2 for the cloudy quadrants, together with the Gaussian distributions that have the same mean and standard deviation. There is no indication in the figure that the assumption of Gaussian shape for the observation errors is a significant issue in comparison with the accuracy and with the linearity properties.
<figure><img src="image_5.png" /><figcaption>Figure 3. Bias and standard deviation (Std) of the differences between the model and the cloud-affected observed AIRS brightness temperatures on 30 November 2002. The model uses the prognostic cloud scheme. Above 2000 cm⁻¹ daytime pixels are discarded. For comparison, the filter functions of the infrared MVIRI channel cover the wave numbers between 1351 and 1818 cm⁻¹ around 6.3 µm and those between 770 and 976 cm⁻¹ around 11 µm.</figcaption></figure>

### (b) AIRS

The model evaluation for the 324 available AIRS channels is less straightforward than for Meteosat because the AIRS spectrometer provides only a few observations for a specific time. Model fields from several time steps are needed in order to accumulate a significant amount of collocated data. As a consequence, AIRS data are being passively monitored in the forecasting system (with ranges from 3 to 15 hours and at resolution 40 km), which takes the model data at observation time. Figure 3 presents the corresponding global statistics of the differences between the prognostic model and the observation for the cloud-affected AIRS channels on 30 November 2002. Other periods have been investigated, and very similar results to those of Fig. 3 have been obtained. For comparison, the wave numbers covered by the MVIRI infrared filter functions are indicated in the legend of the figure. In addition, Fig. 4 displays the fraction of the cloud-affected observed variance explained by the model.

For technical reasons, the diagnostic model cannot be used yet for passive monitoring in the forecasting system. As a consequence, diagnostic and prognostic brightness temperatures are compared independently with the real AIRS data, using the Meteosat-7 6.3 µm cloud mask for 12 UTC 30 November 2002. For each cloud-affected quadrant, equivalent AIRS brightness temperatures are computed using the diagnostic and the prognostic scheme. Rather than using a constant zenith angle, the Meteosat-7 angle is used. Corresponding statistics are presented in Fig. 5.

From Fig. 3, it is obvious that the model statistics are the best in the channels least affected by clouds. Biases are mainly positive, showing that the model underestimates the cloud radiative forcing, consistent with the Meteosat validation and with previous studies (e.g. Chevallier et al. 2001). A different behaviour occurs for the near-infrared channels above 2500 cm⁻¹. Although they are window channels, the bias reduces with
<figure><img src="image_6.png" /><figcaption>Figure 4. Fraction of the observation variance explained by the model for the cloud-affected AIRS data on 30 November 2002. The model uses the prognostic cloud scheme.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 5. Bias and standard deviation of the differences between the model diagnostic and prognostic AIRS brightness temperatures. Model data correspond to the Meteosat-7 disk at 12 UTC 30 November 2002. Clear points are removed using the Meteosat-7 cloud detection. In contrast to Fig. 3, the stratospheric channels are not removed.</figcaption></figure>
increasing wave number and finally changes sign. This is likely to be caused by the absence of cloud scattering and/or of cloud reflection in the radiation model. Large negative values occur around $2300 \text{ cm}^{-1}$ for cloud-affected and for clear channels (McNally 2003, personal communication) and are being investigated. The fraction of variance explained by the model is much higher than for the Meteosat comparison (mentioned in section 4(a)) despite the discrepancy between the observation and the model resolutions (respectively 13 km at nadir and 40 km), which did not happen with the Meteosat RDC. This apparent contradiction actually reflects the different orbits of the two instruments. Indeed mid-latitudes are better modelled than the tropics, in particular over land, which favours the AIRS-based statistics. For instance, if the latitudes within 30° of the equator are removed from the Meteosat statistics, the fraction of the variance explained by the prognostic model reaches 0.66 and 0.89 at 6.3 and 11 µm, respectively.

Since the diagnostic model has been tuned to radiation observations, the biases between diagnostic and prognostic brightness temperatures nearly cancel the biases between prognostic ones and observations (Fig. 5). Standard deviations are slightly smaller between the two models than with observations, but are here much smaller than the observation ‘random’ variations (not shown).

## 5. LINEARITY OF THE OBSERVATION OPERATOR

The linearity assumption is tested here for perturbations $\delta\mathbf{x} = \mathbf{x} - \mathbf{x}^b$ that are of the order of magnitude expected in 4D-Var, i.e. comparable to the background errors. Consequently, the perturbations are defined based on the principal components of the ECMWF operational background-error matrix **B** (Rabier *et al.* 1998; Derber and Bouttier 1999). Temperature errors vary with latitude, and humidity-error statistics are a function of relative humidity. Temperature and humidity errors are uncoupled. One $\delta\mathbf{x}$ is then a Gaussian perturbation applied to all principal components at once. This ensures that **B** is the covariance matrix of the perturbations.

The choice is made here to use the atmospheric profiles within the Meteosat-7 disk at 12 UTC 30 November 2002 as a dataset sampling very diverse atmospheric conditions. For each cloud-affected quadrant, the correlation between the tangent-linear perturbations $\mathbf{H}\delta\mathbf{x}$ (where **H** is the adjoint of the observation operator) of the AIRS brightness temperatures and the nonlinear ones $\mathbf{H}[\mathbf{x}^b + \delta\mathbf{x}] - \mathbf{H}[\mathbf{x}^b]$ is computed using an ensemble of 100 perturbations. The zenith angle is set to that of Meteosat-7. The PDF of the correlations is shown for each channel in Fig. 6(a). Stratospheric channels are easily identified because they correspond to the narrow PDFs close to unity, around wave numbers $500 \text{ cm}^{-1}$ and $2300 \text{ cm}^{-1}$. Channels that both have a sensitivity in the troposphere and systematically correspond to high correlations (e.g. above 0.85) can be found in the H$_2$O $\nu_2$ band (around $1500 \text{ cm}^{-1}$) and next to the stratospheric channels (in the lower-wave-number part of the CO$_2$ $\nu_2$ band (around $700 \text{ cm}^{-1}$) and around $2250 \text{ cm}^{-1}$) only. Those few channels sound the upper troposphere only and are less affected by clouds. Channels with sensitivities lower down in the troposphere show highly nonlinear behaviour. Reducing the humidity perturbations by a factor of two only slightly increases the correlations (not shown). This indicates that improvements in the quality of the background in forthcoming years are not likely to change the status of those channels with respect to linearity. Results for MVIRI, with high correlations at 6.3 µm and low at 11 µm, do not bring additional information and are not reported here.
<figure><img src="image_8.png" /><figcaption>Figure 6. (a) The bi-dimensional histogram of the correlation between linear and nonlinear brightness temperature perturbations for each one of the 324 AIRS channel subset. The input temperature and humidity perturbations follow the statistics of the ECMWF background error. Negative correlations are not represented. (b) The correspondence between channel index and wave number.</figcaption></figure>

## 6. APPLICATION IN ONE DIMENSION

A 1D-Var scheme is used to illustrate the previous findings about the accuracy and the linearity of the infrared satellite radiances. The principle of 1D-Var is similar to that of 4D-Var, but the control vector **x** in Eq. (1) represents only a single column and there is no time dimension. In the present case, the background **x**₀ᵇ comes from the ECMWF short-range forecasts, the observations **y** are either the MVIRI 6.3 µm channel or 35 near-linear AIRS channels at 6.3 and 14.3 µm, as described below. Consistent with the present discussion, the forward operator is linearized around the background
<figure><img src="image_9.png" /></figure>

Figure 7. Root-mean-square (RMS) difference between the model relative humidity and collocated radiosonde measurements for November 2002 and February 2003 in the northern hemisphere mid-latitudes. Relative humidity is defined between 0 and 1. The model profile is either the background or the 1D-Var retrieval. Assimilated observations are either (a) 35 AIRS channels, or (b) the MVIRI 6.3 µm channel. For AIRS, the data are processed when cloudiness is detected in at least 22 channels among the 35. Cloudy MVIRI quadrants only are used in (b). In both cases, statistics include about 250 cases in the upper troposphere and about 1500 in the lower troposphere.

state $x_0^b$ during the 1D-Var minimization. Background-error statistics **B** for the control variables (temperature and humidity) are the ones used in the ECMWF operational model. Observation-error standard deviations in **R** are defined conservatively from the departure statistics presented in Table 1 for MVIRI and Fig. 3 for AIRS, as in section 4(a). Error correlations of 0.8 are arbitrarily specified between AIRS channels. Observations are bias-corrected with respect to the background by removing the biases shown in Table 1 for Meteosat and Fig. 3 and 5 for AIRS.

The 35 near-linear AIRS channels are selected among the subset of 324 from the following three quantitative criteria: the cloud impact on the brightness temperature (estimated from the model simulations) must be more than 0.5 K on average, the correlations between linear and nonlinear increments must exceed 0.85 (from Fig. 6), and the standard deviations of the differences between diagnosed and observed brightness temperatures (computed from the numbers in Figs. 3 and 5 and simply assuming a lack of correlation between one another) must be below 6 K. Further, 4.5 µm AIRS channels are not used here because of the solar radiation, but could be used during night-time. Among the 35, 13 channels are located about 14.3 µm and 22 channels are located about 6.3 µm. At 6.3 µm water-vapour absorption impedes cloud absorption and accurate linear channels can be found with lower weighting functions (i.e. which peak as low as about 400 hPa) than at 14.3 µm.

The 1D-Var is applied to satellite radiances for November 2002 and February 2003, which have been collocated with 00 and 12 UTC operational radiosondes. AIRS radiances are processed only when cloudiness is detected in 22 channels at least among the 35 and MVIRI is used in cloudy quadrants only. Results are presented for collocations in the northern hemisphere mid-latitudes (between 30°N and 70°N) and for 1.5 hours × 80 km time-space windows. They actually mainly represent Europe because of the orbital characteristics of the satellites and of the locations of the operational radiosondes. Inputs and outputs of the 1D-Var are compared with the collocated radiosondes in terms of relative humidity. For atmospheric temperatures below 243 K, only Vaisala RS90 radiosondes are used (Nash 2002).
1D-Var relative-humidity increments reach a maximum at 350 and 400 hPa for MVIRI and AIRS, respectively, with root-mean-square values of 0.017 and 0.10, respectively (not shown). Smaller increments are obtained from MVIRI since a single channel is used. Figure 7 indicates that the 1D-Var method reduces the differences between the model upper-tropospheric relative humidity and the radiosondes measurements by up to 0.005 for MVIRI and by up to 0.015 for AIRS. Interestingly, the difference reduction is about the same when cloud-free observations are processed with clear-sky R statistics and no cloud processes in the observation operator H (not shown). Obviously, no perfect fit between the model and the observations can be achieved because of significant measurement and collocation errors. Most importantly, the high temporal frequency of the MVIRI images is not exploited by the 1D-Var, but it impacts the 4D-Var wind analysis significantly, as was demonstrated by the 4D-Var assimilation of MVIRI in the clear-sky quadrants (Köpken *et al.* 2004).

7. CONCLUSION

An observation operator has been developed that computes cloud-affected satellite brightness temperatures from some of the ECMWF 4D-Var control variables: temperature, humidity and ozone profiles, surface temperature and surface pressure. It comprises a diagnostic cloud scheme with a representation of large-scale and convective processes and a radiation model. In order to evaluate the capability of 4D-Var systems to handle satellite infrared observations in the presence of clouds, its accuracy and its linearity have been assessed. Results have been illustrated within a 1D-Var framework.

A first important result of the present study concerns the diagnostic cloud scheme. It is expected not to perform as well as the prognostic scheme in dynamic mode for long integrations, but in static mode the comparison between model and observations was not qualitatively sensitive to whether the diagnostic model or the reference prognostic cloud scheme is used. Both explain a significant portion of the observation variance, except (for both) in the tropics at 11 µm. The diagnostic scheme explains a slightly lower portion of the variance, but has smaller biases due to different tuning. Further, the two schemes should have similar sensitivities, since one of them is a simplified version of the other one (this aspect will be documented elsewhere). As a consequence, it seems that the diagnostic model could be used in the 4D-Var observation operator in lieu of the prognostic model, which would avoid the introduction of cloud variables in the 4D-Var control vector.

Secondly, it is clear that the channels that are the most impacted by clouds are very nonlinear for perturbations of the order of the current background errors in global NWP. They cannot be processed directly within the current ECMWF 4D-Var system, where only weak nonlinearities are tolerated. In addition, large errors were shown for these channels, both in terms of bias (which could be removed) and of standard deviations. Nonlinearities actually reveal the deficiency of the NWP background in resolving the ambiguity of the cloud-affected radiances in terms of temperature and humidity information. These observations can be pre-processed by a local nonlinear retrieval method, as is done for the geostationary atmospheric-motion vectors (e.g. Schmetz *et al.* 1993).

In contrast, the observation operator showed much more linear and accurate behaviour for some of the upper tropospheric channels, at 4.5, 6.3 and 14.3 µm. It is worth emphasizing two features of the approach. First, accuracy is achieved in these channels despite a lower spatial resolution compared to the observations, which seems to indicate that the representativeness error is not a significant issue here. Second, the focus is on
temperature and humidity fields and not on cloud variables, since the latter are diagnosed from the former. An obvious advantage is that temperature and humidity analysis increments are likely to improve the forecast far away from the analysis (e.g. Marécal and Mahfouf 2002). On the other hand, the accuracy of the present observation operator is still limited and only part of the information of the observations can be extracted.

The conclusion of our assessment is that there is a potential benefit in assimilating cloud-affected satellite radiances at 4.5, 6.3 and 14.3 µm (for instance, the 6.3 µm channel on board all the geostationary satellites) *directly in 4D-Var*. This would avoid blending 4D-Var and local retrieval methods to exploit these channels in the presence of cirriform clouds. Scientific developments to the current 4D-Var systems may still be needed, for instance to improve the estimation of background-error statistics, or to harmonize the resolution of the observations and the variable model resolutions within the incremental formulation. This is a concern for all types of assimilated observations.

A similar study will be performed for a selection of microwave channels in the presence of clouds and rain (Moreau *et al.* 2004). Channels in strong water-vapour absorption bands, for instance at 22.235 GHz or 183.31 GHz, are well modelled (Chevallier and Bauer 2003) and may be sufficiently linear for direct 3D- or 4D-Var assimilation.

## ACKNOWLEDGEMENTS

This work was done in the context of the Satellite Application Facility on Numerical Weather Prediction which is co-sponsored by EUMETSAT. M. Doutriaux-Boucher at EUMETSAT coordinated the set-up of the Radiance Data from Clouds product. The kind help from P. Watts and A. McNally (ECMWF) in dealing with the AIRS data was very much appreciated. Comments from E. Andersson, A. Hollingsworth, A. Simmons and J.-N. Thépaut at ECMWF, from R. W. Saunders at the Met Office and from one anonymous reviewer helped to improve the discussion.

## REFERENCES

<table><tr><td>Chevallier, F. and Bauer, P.</td><td>2003</td><td>Model rain and clouds over oceans: Comparison with SSM/I observations. <i>Mon. Weather Rev.</i>, <b>131</b>, 1240–1255</td></tr><tr><td>Chevallier, F. and Kelly, G.</td><td>2002</td><td>Model clouds as seen from space: Comparison with geostationary imagery in the 11 µm channel. <i>Mon. Weather Rev.</i>, <b>130</b>, 712–722</td></tr><tr><td>Chevallier, F., Bauer, P., Kelly, G., Jakob, C. and McNally, T.</td><td>2001</td><td>Model clouds over oceans as seen from space: Comparison with HIRS/2 and MSU radiances. <i>J. Climate</i>, <b>14</b>, 4216–4229</td></tr><tr><td>Chevallier, F., Bauer, P., Mahfouf, J.-F. and Morcrette, J.-J.</td><td>2002</td><td>Variational retrieval of cloud profile from ATOVS observations. <i>Q. J. R. Meteorol. Soc.</i>, <b>128</b>, 2511–2526</td></tr><tr><td>Courtier, P., Thépaut, J.-N. and Hollingsworth, A.</td><td>1994</td><td>A strategy for operational implementation of 4D-Var, using an incremental approach. <i>Q. J. R. Meteorol. Soc.</i>, <b>120</b>, 1367–1388</td></tr><tr><td>Dee, D. P. and da Silva, A. M.</td><td>2003</td><td>The choice of variable for atmospheric moisture analysis. <i>Mon. Weather Rev.</i>, <b>131</b>, 155–171</td></tr><tr><td>Derber, J. and Bouttier, F.</td><td>1999</td><td>A reformulation of the background error covariance in the ECMWF global data assimilation system. <i>Tellus</i>, <b>40A</b>, 1–25</td></tr><tr><td>Eyre, J. R.</td><td>1991</td><td>‘A fast radiative transfer model for satellite sounding systems’. ECMWF Technical Memorandum 176, European Centre for Medium-Range Weather Forecasts, Shinfield Park, Reading, Berkshire RG2 9AX, UK</td></tr><tr><td>Harris, B. A. and Kelly, G.</td><td>2001</td><td>A satellite radiance bias correction scheme for radiance assimilation. <i>Q. J. R. Meteorol. Soc.</i>, <b>127</b>, 1453–1468</td></tr><tr><td>Ide, K., Courtier, P., Ghil, M. and Lorenc, A.</td><td>1997</td><td>Unified notation for data assimilation: Operational, sequential and variational. <i>J. Meteorol. Soc. Japan</i>, <b>75</b>, 181–189</td></tr><tr><td>Janisková, M., Mahfouf, J.-F., Morcrette, J.-J. and Chevallier, F.</td><td>2002</td><td>Linearized radiation and cloud schemes in the ECMWF model: Development and evaluation. <i>Q. J. R. Meteorol. Soc.</i>, <b>128</b>, 1505–1528</td></tr></table>
<table>
  <tr>
    <td>Köpken, C., Kelly, G. and Thépaut, J.-N.</td>
    <td>2004</td>
    <td>Assimilation of Meteosat radiance data within the 4D-Var system at ECMWF: Assimilation experiments and forecast impact. Q. J. R. Meteorol. Soc., 130, (in press)</td>
  </tr>
  <tr>
    <td>Lorenc, A. C., Ballard, S. P., Bell, R.S., Ingleby, N. B., Andrews, P. L. F., Barker, D. M., Bray, J. R., Clayton, A. M., Dalby, T., Li, D., Payne, T. J. and Saunders, F. W.</td>
    <td>2000</td>
    <td>The Met Office global three-dimensional variational data assimilation scheme. Q. J. R. Meteorol. Soc., 126, 2991-3012</td>
  </tr>
  <tr>
    <td>Marécal, V. and Mahfouf, J.-F.</td>
    <td>2002</td>
    <td>Four-dimensional variational assimilation of total column water vapor in rainy areas. Mon. Weather Rev., 130, 43-58</td>
  </tr>
  <tr>
    <td>McNally, A. P.</td>
    <td>2002</td>
    <td>A note on the occurrence of cloud in meteorologically sensitive areas and the implications for advanced infrared sounders. Q. J. R. Meteorol. Soc., 585, 2551-2556</td>
  </tr>
  <tr>
    <td>McNally, A. P. and Watts, P. D.</td>
    <td>2003</td>
    <td>A cloud detection algorithm for high spectral resolution infrared sounders. Q. J. R. Meteorol. Soc., 129, 3411-3423</td>
  </tr>
  <tr>
    <td>Moreau, E., Bauer, P. and Chevallier, F.</td>
    <td>2003</td>
    <td>Variational retrieval of rain profiles from spaceborne passive microwave radiance observations. J. Geophys. Res., 108(D16), 4521, doi:10.1029/2002JD003315</td>
  </tr>
  <tr>
    <td>Moreau, E., Lopez, P., Bauer, P., Tompkins, A. M., Janisková, M. and Chevallier, F.</td>
    <td>2004</td>
    <td>Variational retrieval of temperature and humidity profiles using rain rates versus microwave brightness temperatures. Q. J. R. Meteorol. Soc., 130, (in press)</td>
  </tr>
  <tr>
    <td>Nash, J.</td>
    <td>2003</td>
    <td>'Review of test results on the accuracy of radiosonde relative humidity sensors'. Pp. 117-123 in Proceedings of the ECMWF/GEWEX workshop on humidity analysis, 8-11 July 2002, European Centre for Medium-Range Weather Forecasts, Shinfield Park, Reading, Berkshire RG2 9AX, UK</td>
  </tr>
  <tr>
    <td>Rabier, F., McNally, A. P., Andersson, E., Courtier, P., Undén, P., Eyre, J., Hollingsworth, A. and Bouttier, F.</td>
    <td>1998</td>
    <td>The ECMWF implementation of three-dimensional variational assimilation (3D-Var). II: Structure functions. Q. J. R. Meteorol. Soc., 124, 1809-1829</td>
  </tr>
  <tr>
    <td>Rabier, F., Järvinen, H., Klinker, E., Mahfouf, J.-F. and Simmons, A.</td>
    <td>2000</td>
    <td>The ECMWF operational implementation of four-dimensional variational assimilation. I: Experimental results with simplified physics. Q. J. R. Meteorol. Soc., 126, 1143-1170</td>
  </tr>
  <tr>
    <td>Räisänen, P.</td>
    <td>1998</td>
    <td>Effective longwave cloud fraction and maximum-random-overlap clouds—a problem and a solution. Mon. Weather Rev., 126, 3336-3340</td>
  </tr>
  <tr>
    <td>Saunders, R. W., Brunel, P., Chevallier, F., Deblonde, G., English, S. J., Matricardi, M. and Rayer, P. J.</td>
    <td>2002</td>
    <td>'RTTOV-7 - Science and validation report'. Met Office Forecasting Research Technical Report, 387, available from URL http://www.metoffice.com/research/nwp/publications/papers/technical_reports/sa.html</td>
  </tr>
  <tr>
    <td>Schmetz, J. and Turpenheim, O.</td>
    <td>1988</td>
    <td>Estimation of the upper-tropospheric relative humidity field from METEOSAT water vapor image data. J. Appl. Meteorol., 27, 889-899</td>
  </tr>
  <tr>
    <td>Schmetz, J., Holmlund, K., Hoffman, J., Strauss, B., Mason, B., Gaertner, V., Koch, A. and van De Berg, L.</td>
    <td>1993</td>
    <td>Operational cloud-motion winds from Meteosat infrared images. J. Appl. Meteorol., 32, 1206-1225</td>
  </tr>
  <tr>
    <td>Simpson, J. and Wiggert, V.</td>
    <td>1969</td>
    <td>Models of precipitating cumulus towers. Mon. Weather Rev., 97, 471-489</td>
  </tr>
  <tr>
    <td>Stokes, G. M. and Schwartz, S. E.</td>
    <td>1994</td>
    <td>The Atmospheric Radiation Measurement (ARM) Program: Programmatic background and design of the cloud and radiation testbed. Bull. Amer. Meteorol. Soc., 75, 1201-1221</td>
  </tr>
  <tr>
    <td>Tiedtke, M.</td>
    <td>1989</td>
    <td>A comprehensive massflux scheme for cumulus parametrization in large-scale models. Mon. Weather Rev., 117, 1779-1800</td>
  </tr>
  <tr>
    <td>Tiedtke, M.</td>
    <td>1993</td>
    <td>Representation of clouds in large-scale models. Mon. Weather Rev., 121, 3040-3061</td>
  </tr>
  <tr>
    <td>Tompkins, A. M. and Janisková, M.</td>
    <td>2003</td>
    <td>'A cloud scheme for data assimilation: Description and initial tests'. ECMWF Technical Memorandum 410, European Centre for Medium-Range Weather Forecasts, Shinfield Park, Reading, Berkshire RG2 9AX, UK</td>
  </tr>
</table>