Operational numerical weather prediction systems currently only assimilate infrared and microwave satellite observations, whereas visible and near-infrared reflectances that comprise information on atmospheric clouds are not exploited. One of the reasons for that is the absence of computationally efficient observation operators. To remedy this issue in anticipation of the future regional Kilometer-Scale Ensemble Data Assimilation (KENDA) system of Deutscher Wetterdienst, we have developed a version that is fast enough for investigating the assimilation of cloudy reflectances in a case study approach. The operator solves the radiative transfer equation to simulate visible and near-infrared channels of satellite instruments based on the one-dimensional (1D) discrete ordinate method. As input, model output of the operational limited-area Consortium for Small-Scale Modeling (COSMO) model of Deutscher Wetterdienst is used. Assumptions concerning subgrid-scale processes, calculation of in-cloud values of liquid water content, ice water content, and cloud microphysics are summarized, and the accuracy of the 1D simulation is estimated through comparison with three-dimensional (3D) Monte Carlo solver results. In addition, the effects of a parallax correction and horizontal smoothing are quantified. The relative difference between the 1D simulation in “independent column approximation” and the 3D calculation is typically less than 9% between 0600 and 1500 UTC, computed from four scenes during one day (with local noon at 1115 UTC). The parallax-corrected version reduces the deviation to less than 6% for reflectance observations with a central wavelength of 810 nm. Horizontal averaging can further reduce the error of the 1D simulation. In all cases, the bias is less than 1% for the model domain.
Extending the use of satellite radiances for numerical weather prediction (NWP) is a high priority at many forecast centers. While the assimilation of satellite radiances has led to some of the greatest increases in forecast skill that have been achieved during the last decade, the current use of satellite radiances is still restrictive with only a small fraction of the available observations being included in a data assimilation (DA) process. Particularly, a better exploitation of cloud- or precipitation-affected satellite measurements could bear great potential for further improvements of weather forecasting (Bauer et al. 2011a). These data specifically provide information from overcast regions that are typically sensitive regions with great importance for NWP (McNally 2002). In particular, information linked to cloud variables and precipitation could help to improve the forecast of convective precipitation, which is one of the key targets for regional high-resolution limited-area models.
The assimilation of radiances that are affected by clouds or precipitation is, however, much more difficult than in clear air (Errico et al. 2007). Crucial reasons for this are the complexity and nonlinearity of the relevant forward operators that increase substantially in the presence of water in the condensed or frozen phase (see, e.g., Bennartz and Greenwald 2011). Such forward operators (also called observation operators), which compute the model equivalent for the respective observation types, are vital parts of modern DA systems. For variational DA systems, their linearized and adjoint versions are also required, while for ensemble DA systems the forward operator itself is sufficient.
For satellite radiances, the forward operator includes a radiative transfer (RT) model that computes the radiances that would be measured by the satellite instrument for a given atmospheric state. In the presence of clouds, RT computations can become very demanding (Liou 1992), especially in the solar spectral range. However, a crucial requirement for developing a DA system that can deal with cloudy radiances is a sufficiently fast and reliable RT model for the respective wavelengths.
So far, most of the radiance assimilation efforts (including those concerning cloud affected measurements) were made for global models (i.e., synoptic scale) and were focused on radiation in the microwave (MW) or infrared (IR) spectral bands (Bauer et al. 2011b). In some respect, the situation is easiest for the MW spectrum, where clouds are rather transparent and only very thick water clouds and rain significantly impair the ability to undertake quantitative retrievals. As a consequence, the corresponding RT operator is more linear than for IR radiances and an all-sky approach has been successfully adopted at the European Centre for Medium-Range Weather Forecasts (Bauer et al. 2010).
For IR radiances, RT computations are more nonlinear and very sensitive to the input cloud variables. For this reason, assimilation methods have been developed that intend to “subtract” the influence of clouds on the RT computations in order to assimilate the same fields as for clear-air assimilation despite the presence of clouds (McNally 2009; Pavelin et al. 2008; Pangaud et al. 2009), rather than exploiting the cloud information contained in the cloudy radiances. The temperature and humidity fields constrain the occurrence of clouds to a certain extent, but the full observed information on clouds is not directly assimilated.
A central task for limited-area models is to produce a more accurate short-term forecast of clouds and precipitation. For the initialization of such models, the explicit exploitation of cloud information therefore has higher priority than for global models. One of the most fundamental problems in this context is to improve location errors, that is, situations where observed clouds are displaced or completely missing in the model (or where model clouds have no counterpart in the observations). Some recent work has shown that variational DA methods (while showing skill in improving properties of correctly located model clouds) have strong limitations in such situations and often a cloud mask is employed for explicitly limiting the assimilation to cases where model clouds and observed clouds are sufficiently close (Polkinghorne and Vukićević 2011; Seaman et al. 2010; Okamoto 2013; Chevallier et al. 2004; Stengel et al. 2013, 2010). An interesting method for tackling such limitations was developed by Renshaw and Francis (2011). Another approach involves ensemble DA methods, which seem to be less severely affected by this problem (Otkin 2010, 2012a,b; Zupanski et al. 2011).
While most of the radiance assimilation experiments so far have focused on the IR and MW radiances, forward operators that also include the visible (VIS, 390–700 nm) and near-infrared (NIR, 0.7–5 μm) spectral range have also been developed (e.g., Greenwald et al. 2002, 2004, Evans 2007).
In this paper we present another forward operator that is also suitable for radiances in this spectral range and that can be used in the experimental regional Kilometer-Scale Ensemble Data Assimilation (KENDA) system of Deutscher Wetterdienst (DWD), which is based on a local ensemble transform Kalman filter (LETKF; Hunt et al. 2007). More precisely, the operator is designed to enable the KENDA system to assimilate data from the geostationary platform Meteosat, which are available with a high temporal resolution.
If the aim is to exploit cloud information, then it seems natural to draw the attention to the VIS and NIR spectra even though the corresponding RT computations are comparably complex. VIS and NIR observations provide a wealth of cloud information and by this a much earlier detection of convective activity than, for example, radar observations, which are sensitive to larger droplets only. Given the major focus of convective-scale models to forecast convective precipitation for comparably short-lead times (typically a few hours to one day), these are seen as a promising data source to represent convective activity correctly already at early stages. VIS and NIR channels also saturate less quickly than IR for water clouds and by this they contain more information on the optical thickness and the related cloud water content, where the IR would provide only a yes/no information and the cloud-top temperature. For this reason, remote sensing of optical thickness and effective radius is only done during daytime using the solar channels.
Another advantage of VIS channels is that low cumulus clouds are better distinguishable from the surface signal, since they are usually much brighter than the surface, whereas in the IR low clouds are hardly distinguishable from the surface due to their similar brightness temperatures. Finally, compared to IR channels, VIS and NIR are less sensitive to thin cirrus clouds and may therefore also provide information about clouds below thin cirrus that would be hidden in the IR. The resolution of VIS and NIR satellite observations of typically a few kilometers also matches well with the grid spacing of current regional models. The Spinning Enhanced Visible and Infrared Imager (SEVIRI) aboard the satellites of Meteosat Second Generation (MSG), for example, has a resolution of 3 km at the subsatellite point. MW and most IR observations that are currently assimilated at NWP centers [such as the Atmospheric Infrared Sounder (AIRS), the Infrared Atmospheric Sounding Interferometer (IASI), the Special Sensor Microwave Imager (SSM/I), etc.] in contrast are well matched with the grid spacing of global models. The goal for convective-scale data assimilation systems should therefore be to include VIS and NIR in addition to MW and IR channels, as the different observation types are in many ways complementary.
In the past, many decisions with respect to wavelength selection and assimilation strategy were made with regard to variational DA systems that are extremely demanding concerning the possible linearization of the forward operator, as nonlinearities can prevent the convergence of the minimization of the cost function. Lately, many operational centers started to develop DA systems based on ensemble Kalman filter (EnKF) methods for their limited-area models. While these also make assumptions about the linearity of the assimilation problem, they are expected to be more robust with respect to the occurrence of nonlinear effects (Kalnay et al. 2008). Since the assimilation of cloud information is a high priority for these models, we believe that the direct assimilation of VIS and NIR radiances yields a great potential. However, no operational global or regional NWP model assimilates such observations, and also assimilation experiments exploring the impact of these wavelengths seem to be extremely rare and, to our knowledge, all in the context of variational DA systems (where no significant positive impact could be demonstrated; this, however, could be linked to the inability of such systems to correct for location errors; see Polkinghorne and Vukićević 2011).
The paper is structured as follows. Section 2 introduces the configuration of the operational limited-area Consortium for Small-Scale Modeling (COSMO) model used at DWD (COSMO-DE) and its relevant output for the RT calculations. Furthermore, the concept of RT and the particular solvers applied in this article are described. In section 3, important parameterizations used in the forward operator are summarized. These include the total liquid and ice water content calculated from both grid-scale model variables and assumptions about the subgrid-scale cloud water mixing ratios (liquid and frozen). In addition, the parameterizations of effective scattering radii of water droplets and ice crystals in clouds are given. Section 4 describes the preprocessing parallax correction that is applied to simulate 1D RT in columns tilted toward the satellite to account for the slant viewing angle. The accuracy assessment based on the comparison of 1D and 3D results is presented in section 5 and a summary is given in section 6.
This section provides a description of the configuration of the operational limited-area model COSMO-DE used at DWD, the processing of its output to synthetic satellite images using forward operators, and the main properties of the employed 1D and 3D RT solvers used in this study.
a. Meteorological model and data
The forecast fields used to simulate synthetic satellite images are produced by the COSMO community model (Baldauf et al. 2011). The COSMO model has been used for operational numerical weather prediction at DWD since 1999. The convection-permitting model configuration COSMO-DE has been operational since April 2007. The model domain has a horizontal grid spacing of 2.8 km and consists of 421 × 461 grid points. The area covers Germany, as well as Switzerland, Austria, and parts of the other neighboring countries of Germany. In the vertical, it consists of 50 terrain-following model layers. The lowest level lies 10 m aboveground, and the model top lies 22 km above mean sea level. The model explicitly resolves deep convection, while shallow convection is parameterized (Baldauf et al. 2011).
The VIS and NIR operators use the model output of temperature, pressure, mixing ratios of humidity, cloud liquid water, cloud ice, and snow, as well as cloud fraction in each layer and the base and top heights of shallow convective clouds. In addition, the temporally constant parameters orography, geometrical height of model layer boundaries, latitude, and longitude are input for the operator. As a case study, 22 June 2011 has been chosen and output fields from 3-h forecasts at 0600, 0900, 1200, 1500, and 1800 UTC have been used for the simulations. This is a particularly interesting day from the meteorological point of view, since on 22 June 2011 a well-developed cold front at the leading edge of an upper-level trough passed Germany. A strong jet streak at 500 hPa overlapped with low-level instability, providing favorable conditions for deep convection. Heavy rain, hail, strong winds, and a tornado were observed in central Germany. Satellite imagery of this event is provided in section 5. On such a day, the assimilation of VIS and NIR channels could be particularly beneficial by identifying the convective activity better and at an early stage. Though visible channels can be used to identify convective activity early in its development, their utility for data assimilation and as a monitoring tool is limited only to daytime hours. Infrared channels, however, can provide information about convective development both day and night. Their simultaneous use would be optimal.
b. Radiative transfer models
As a tool to simulate RT for solar radiation, the software package libRadtran by Mayer and Kylling (2005) is applied. It contains the uvspec model, a command-line-based executable to solve the RT equations. Input files are used to concisely define an atmospheric scene in terms of profiles of water and ice clouds represented by their liquid water content (LWC), ice water content (IWC), surface albedo, trace gases, aerosol, pressure, and temperature. In combination with information about microphysical cloud properties, such as the effective radii of scattering particles, the corresponding optical properties are searched for in lookup tables. The parameterizations used to calculate LWC, IWC, and the corresponding effective radii are described in section 3. Subsequently, the optical properties given in terms of the extinction coefficient, the single-scattering albedo, and the scattering phase function are passed on to the RT solver, which calculates reflectances. Finally, a postprocessing step takes into account the extraterrestrial solar spectrum, including earth–sun distance variations, to determine the final output (as chosen by the user; in our case, reflectance).
The libRadtran software includes several RT solvers of varying complexity and degree of approximation. In the context of this study, two solvers are applied. The first one is the 1D solver based on the Discrete Ordinate Radiative Transfer model (DISORT) by Stamnes et al. (1988), modified and translated into C code by Buras et al. (2011), which is used in our proposed forward operator. The second one is the Monte Carlo code for the physically correct tracing of photons in cloudy atmospheres (MYSTIC) 3D solver (Emde and Mayer 2007; Mayer 2009; Buras and Mayer 2011), which is used as “model truth.”
Each solver provides a numerical solution to the radiative transfer equation (Chandrasekhar 1960):
where I denotes the radiance for a certain location and direction, β is the volume extinction coefficient, ω is the single-scattering albedo, B(T) is the Planck function, and P(Ω, Ω′) is the scattering phase function determining the probability of scattering from a beam direction Ω′ to Ω. For the case at hand, where the focus lies on RT in the solar channels, the emission given by the last term involving B(T) is negligible for VIS and comparably small for the NIR channel used in this study. At longer wavelengths, however, thermal emission becomes more important.
The 1D solver DISORT solves Eq. (1) in a horizontally homogeneous plane-parallel atmosphere1 by discretizing into a finite amount of angular streams s on which the scattering integral is evaluated in terms of Gaussian quadrature. For this purpose, the scattering phase function is expanded into a finite series of Legendre polynomials.2 The RT equation is solved in each of the nz atmospheric layers with constant optical properties. Thus, a total number of 2 × (s × nz) equations has to be evaluated, where continuity requirements for the radiance field need to be satisfied at the level interfaces. In the presented examples, nz is set to 50 and the number of angular streams s is set to 16. The 1D solver is sufficiently fast for case study purposes in an offline DA calculation. Nevertheless, having a computation time of approximately 5–10 min per scene over the whole model domain (run on 37 processors), it is still beyond the limitations of an operational ensemble DA system.
The Monte Carlo solver MYSTIC is a probabilistic approach to the solution of Eq. (1). It traces model photons on their way through the atmosphere. Scattering and absorption in the atmosphere and reflection and absorption at the ground are accounted for. At each interaction point, the properties (type of extinction process, scattering angles in the case of scattering, etc.) are drawn randomly using the respective cumulative probability density and the Mersenne Twister (MT 19937) random number generator (Matsumoto and Nishimura 1998). The length of a path in between interacting grid boxes can be calculated by integrating the extinction coefficient along the path until the optical depth drawn randomly from the inverse Lambert–Beer probability density is reached. For each scattering process, the same scattering phase function as for the DISORT solver is used for randomly choosing the scattering angle. These steps are repeated for a large number of model photons. MYSTIC has been validated in an extended model intercomparison project [Intercomparison of Three-Dimensional Radiation Codes (I3RC)] in Cahalan et al. (2005), where the agreement between the individual models was typically on the 1% level.
For our application, we are interested in satellite radiances (or equivalently reflectances), which are difficult to obtain from standard Monte Carlo simulations, because the photons rarely hit the detector, let alone come from the direction of viewing. Therefore so-called variance reduction techniques are used that increase the efficiency by several orders of magnitude. We use the backward Monte Carlo approach, where photons are generated in the final outgoing direction at top of the atmosphere and travel backward. At each interaction with the atmosphere or surface, a local estimate is performed; that is, the probability that the photon scatters/reflects toward the sun and is not extinct on its subsequent way through the atmosphere is calculated. The sum of all local estimates yields the correct result for the radiance measured by the satellite, as can be proven with the von Neumann rule (Marchuk et al. 1980). For a detailed description of the local estimate technique, see Mayer (2009). Because of convergence problems arising when using the local estimate technique in the presence of clouds, we also use the set of variance reduction techniques [Variance Reduction Optimal Options Method (VROOM)] described in Buras and Mayer (2011).
The main uncertainty of MYSTIC is the statistical photon noise (roughly proportional to ), which is small provided that the number of photons N is large enough. For the purpose of this study, the 3D RT simulations will be considered as model truth against which the results of the 1D operator are verified. The big disadvantage of the Monte Carlo method is certainly the excessively large amount of computer time required to obtain a result with a small statistical error (t ~ N ~ σ−2). Therefore, it remains a good research tool for producing very realistic simulations; however, its capability for operational applications—for example, observation operators for cloudy satellite radiances—is very limited with current computer systems. An example of the computational time in the cases at hand is about 12 h per scene, run on 37 processors.
For the parameterization of molecular absorption, the Low-Resolution Atmospheric Transmission (LOWTRAN) band model by Pierluissi and Peng (1985) has been applied as adopted from the Santa Barbara DISORT Atmospheric Radiative Transfer (SBDART) code by Ricchiazzi et al. (1998). Thus, a three-term exponential fit is used for the transmission, which implies that one simulation corresponds to three solutions of the RT equation for one spectral increment. Standard precalculated Mie lookup tables are used for scattering by water droplets. The scattering tables are based on the algorithm described in Wiscombe (1996). For the scattering of radiation by nonspherical ice crystals, the parameterizations by Baum et al. (2005a,b, 2007) are used. Since the main concern of the present work is the effect of clouds on solar radiation, aerosols have been neglected at the current stage.
Within this article, the calculated radiance is converted to reflectance, defined by
where E0 denotes the extraterrestrial flux and θ0 is the solar zenith angle (SZA). For the sake of clarity, we have explicitly included the dependencies on viewing angles (zenith angle θ and azimuth angle ϕ).
Because of unresolved processes in the model, assumptions about subgrid-scale contributions to liquid and frozen cloud water have to be implemented as parameterizations in the forward operator besides approximations about the sizes of scattering particles.
a. Liquid and ice water content
The input parameters to the forward operator are the grid-scale fields of pressure P, temperature T, and the mixing ratios of humidity QV, liquid cloud water QC, cloud ice QI, and snow QS. Model fields of cloud fraction CLC as well as the base height and top height of shallow convective clouds are also input for the forward operator. Since the COSMO model resolves deep convection, the corresponding mixing ratios are contained in the grid-scale fields in contrast to the treatment of shallow convection, which is parameterized as a subgrid-scale process. The cloud-related input variables (QC, QI, and QS) are all grid-scale quantities. To include the impact of subgrid processes in the calculations of radiation, the COSMO model uses a subgrid parameterization that derives the respective cloud variables and used in the model’s radiation scheme. To derive the input quantities for the RT solver, the VIS and NIR forward operator largely follows this subgrid scheme. The only difference is that the forward operator replaces the input variable QI by a mixed variable . This slightly revises the separation between ice and snow carried out by the COSMO model, whose radiative interaction has been tuned with respect to thermal radiation only. In the following, we have chosen κ = 0.1 (which should be well within the uncertainty related to the partitioning between ice and snow). Although we are aware of the fact that this particular choice of κ is rather heuristic, a sensitivity study determining an optimal choice of this parameter goes beyond the scope of this work. For convenience, Table 1 summarizes relevant variables and their meanings, which are used in the following description of the COSMO model’s subgrid scheme.
In the latter, the grid-scale input variables QC and QI only serve to specify lower bounds for the subgrid variables and of in-cloud water mixing ratios (liquid and frozen) from which the radiatively active quantities are, respectively, derived. Apart from these lower bounds, and are determined by
the assumption that the subgrid in-cloud water Qsgs is half a percent of the saturation value, that is, Qsgs = 0.005Qsat, and
the partitioning of Qsgs, which is done through a simple temperature-dependent coefficient fice, that is, and .
As seen from Eq. (A6), the coefficient fice decreases linearly from a value of 1 for temperatures below −25°C to 0 at −5°C (and above). This coefficient is also used in the definition of the effective saturation value Qsat, which is a linear combination of the saturation values over liquid water and ice ; see Eqs. (A5) and (A3) for definitions.
It has to be noted that the Qsgs variable described above represents only one part of the subgrid variations, which are parameterized in the COSMO model. A second type of subgrid variability that the subgrid scheme accounts for stems from shallow convective clouds (which are also parameterized in the COSMO model). For this cloud type, Qcon = 0.2 g kg−1 has been chosen for the in-cloud cloud water mixing ratio (liquid and frozen) except for very large values of Qsat (with Qsat > 20 g kg−1) for which 1% of Qsat is assumed for Qcon. As above, the partitioning of Qcon into liquid and ice clouds ( and , respectively) is also determined by the coefficient fice.
Relating the in-cloud variables to the effective, radiatively active variables and requires a partitioning of the total cloud fraction = CLC/100 into a shallow convective part con (which is related to Qcon) and the remaining part ( − con), which is related to Qsgs. Following the COSMO model’s subgrid scheme, con is diagnosed from the total height of shallow convective clouds, as given in Eq. (A7) of the appendix. One can write the radiatively active total mixing ratios as
from which the corresponding values (g m−3) of LWC and IWC are given by
where ρ is the density of humid air and ρd is the density of dry air (g m−3). The densities are determined using the ideal gas equation of state [Eq. (A1)]. In the last step on the right of Eq. (4), the fact that ρ can be approximated by ρd at sufficiently low temperatures was used (which holds for the temperature range where ice processes are active in this scheme). For the RT simulations, a plane-parallel assumption is made, which implies that the cloud condensate determined by Eqs. (3) and (4) is constant within a grid box.
b. Microphysical parameterizations
Once the total LWC and IWC from both grid-scale and subgrid-scale quantities have been calculated, further assumptions concerning the associated cloud microphysics have to be made. In particular, the effective radii of the scattering particles of solar radiation need to be estimated.
Following the assumptions in Bugliaro et al. (2011), the effective radii of water droplets in clouds are parameterized depending on LWC (g m−3), droplet number concentration N (m−3), and water density ρ ≈ 106 g m−3 at 4°C. The parameterization for the effective radius reads
where is the ratio between the volumetric radius of droplets and the effective radius. For all examples given, N = 1.5 × 108 m−3 is chosen according to Bugliaro et al. (2011) and the value of k = 0.67 is chosen for mainly continental clouds according to Martin et al. (1994). Lower and upper limits on the effective radii of water droplets are taken to be 1 and 25 μm, respectively, since we are primarily concerned about cloud droplets. Larger droplets, such as rain drops, are neglected.
For ice crystals, a parameterization of randomly oriented hexagonal columns described in Bugliaro et al. (2011) is used, who adopted from Wyser (1998) and McFarquhar et al. (2003). Similar as for water droplets, the effective radii of ice crystals in cirrus clouds depend on IWC (g m−3) and temperature T (K) as given by
Effective radii of the scattering ice particles calculated by Eqs. (6) are determined (μm). They are restricted to values between 20 and 90 μm.
4. Parallax correction
In this section, a grid transformation on the input variables LWC, IWC, , and used by the RT solver is described that corrects the error due to the slant satellite viewing angle through the atmosphere. The correction is referred to as parallax correction.
Each grid box, defined by the indices (i, j, k) representing longitude, latitude, and altitude, respectively, is shifted horizontally by (Δi, Δj) pixels. The Δi, Δj need to be chosen such that they correct the parallax. For this purpose, the shift should be
for the latitudinal direction, where ϕ is satellite azimuth angle, θ is the satellite zenith angle, and Δz is the altitude of the upper boundary of the grid box (see Fig. 1).
For the longitudinal direction, the shift should be
We discretize the shifts by dividing (Δx, Δy) by the grid resolution of 2.8 km and finally compute the rounded integers (Δi, Δj). To give an example of typical shifts, we have calculated the average over the model domain in each layer. According to Eqs. (7) and (8), the shifts are proportional to the height Δz. Hence, they increase linearly from at the ground to at 10 km in the y direction. In the x direction, the shifts are much less and only increase from at the ground to at 20 km. The smaller adjustments of Δi are because the longitude of the satellite position (in our case, at 9.5°E) lies within the model domain.
The transformation mapping the input variables from the old to the new grid is thus carried out according to
run over all grid boxes (i, j, k), where X refers to the three-dimensional arrays containing the variables LWC, IWC, and , and to their values on the new grid. Using the transformed grid to simulate RT in independent column approximation (ICA) takes the effect of the satellite viewing angles into account, however, with the advantage of using the faster 1D RT solver instead of the computationally expensive 3D RT solver. In section 5, the results, including the parallax correction, are compared to the uncorrected 1D operator results.
5. Accuracy assessment
a. Experimental setup
As mentioned above, 22 June 2011 has been chosen for the case study to assess the 1D operator accuracy. The 3-h forecast fields of COSMO-DE are used to simulate synthetic satellite images in 3D and 1D at 0600, 0900, 1200, 1500, and 1800 UTC. For this case study, observations are simulated for the SEVIRI instrument aboard the Meteosat-8 satellite of MSG. Nonetheless, the forward operator introduced here is not limited to this particular instrument.
The satellite viewing angles on each individual pixel of the COSMO-DE domain are accounted for. To have a direct comparison between 3D and 1D RT, additional simplifications are made to ensure that no error is introduced due to different treatments in the calculations. The simplifications made are that the model levels, as well as the solar angles, are kept constant over the scene at a particular time. Therefore, a constant SZA is assumed throughout the whole domain (corresponding to the pixel in the middle of the scene with latitude 50.8° and longitude 10.4°). Given that we are only interested in the accuracy of the 1D operator as compared to a “perfect” 3D simulation, this slightly unrealistic model representation is acceptable. In both 1D and 3D calculations, aerosols have been ignored. For our purpose (i.e., improving the location and structure of clouds in a weather forecasting model) this seems acceptable, since in the large majority of cases aerosols have a subdominant effect on VIS and NIR radiation compared to cloud water and ice. In addition, the operational COSMO-DE forecasts do not contain aerosols. Any usage of aerosols would thus be a crude estimation from which we do not expect significant benefit.
To avoid errors due to boundary effects, a smaller grid of 390 × 420 pixels is used for the evaluation of the accuracy. The first reason for this is that the MYSTIC simulations use periodic boundary conditions that would introduce an error in our model truth at the boundaries. Second, COSMO-DE forecasts are integrated with lower-resolution 7-km COSMO boundary conditions [COSMO Europe (COSMO-EU) model–domain]. These introduce a kind of “driving” error at the edges of the model domain due to possible inconsistencies between COSMO-EU and COSMO-DE fields, which also requires that the edges are neglected in future assimilation experiments. Removing 26 pixels in the north, 15 in the south, 15 in the west, and 16 in the east of the original COSMO-DE domain ensures that at least 42 km are cut off of each boundary.
The 3D MYSTIC simulations have been carried out with N = 3 × 104 photons per pixel. In the cases at hand, the MYSTIC simulations have an uncertainty of about 1%–1.5%. This estimated range for the standard deviation includes clear and cloudy scenes and the considered wavelengths.
To quantify the relative difference between 3D and 1D simulations, we use the following formula:
where the sums are calculated over all pixels of the relevant domain and Rij is the reflectance in pixel (i, j). Unless stated otherwise, the term relative difference refers to the quantity defined in Eq. (10). Similarly, the relative bias is given by
Another measure commonly used is the root-mean-square error (RMSE), which we normalize with the mean 3D reflectance , yielding the quantity
where nx and ny denote the number of pixels in the i and j directions of the relevant model domain.
By looking at different times of the day, the dependence of the relative difference on the SZA is determined in Table 2. The table shows the results of the relative difference defined in Eq. (10) obtained using different settings for the simulation. For completeness, the corresponding solar azimuth angle (SAA) at each time is also given in the table (0° corresponds to the southern direction and the angle increases clockwise). The “ICA” stands for the plain independent column approximation on 2.8-km resolution, “parallax” denotes the 1D solver applied to the parallax-corrected fields on 2.8-km resolution, “3 × 3 mean” is a moving average of the parallax-corrected version where the reflectance in each pixel is calculated by taking the moving average over 3 × 3 pixels (centered in the respective pixel), and “5 × 5 mean” denotes a moving average over 5 × 5 pixels.
An example of the 3D and 1D operator output of a full COSMO-DE scene is depicted in Fig. 2. Comparing the two simulations, one can easily distinguish the main differences. Cloud shadows become apparent in the 3D simulation in this afternoon scene at 1500 UTC with a SZA of 50°. These can obviously not be captured by the 1D operator.
Overall, the parallax correction improves the plain ICA result by about 2%. Taking the moving average over 3 × 3 pixels smoothes the field and therefore eliminates errors due to small horizontal displacements which results in a further improvement by 1%–2%. Going to smoothing over 5 × 5 pixels results in yet another small improvement. Between 0600–1500 UTC, the relative difference is smaller than 9% in all cases, while at 1800 UTC it increases significantly to over 20% in the nonaveraged cases. This strong increase in the differences is a result of the large SZA of 78°, which leads to larger cloud shadows than in the earlier scenes. A sensitivity study, in which we artificially changed the SZA for the 1800 UTC case to 50° (the value at 1500 UTC), revealed that the difference is not very sensitive to the type of clouds involved. We conclude that for the assimilation of cloudy VIS and NIR reflectances, one might want to discard observations with an SZA larger than 70° or adjust the errors in the assimilation system unless further corrections are applied. The absolute value of the relative bias is very small (less than 0.6%) for all simulated cases (Table 3). For the readers more familiar with RMSE statistics, the same results in terms of a normalized RMSE [see Eq. (12)] are provided in Table 4.
To provide an example of the corresponding results for the SEVIRI channels VIS006 in the visible with a central wavelength of 635 nm and NIR016 in the near-infrared with the central wavelength at 1.64 μm, 3D and parallax-corrected 1D simulations have been carried out at 1500 UTC. Figure 3 shows the corresponding 3D operator output reflectance fields. For channel VIS006, the relative difference is 6.1% with a bias of −0.4% and for channel NIR016 it is 7.0% with a bias of −1.2%. We conclude that the accuracies are of similar magnitude for the two VIS channels, while the NIR channel is slightly less accurate. The model cloud fraction at 1500 UTC is depicted in Fig. 4. When comparing it to the RT simulations in Figs. 2 and 3, it can be seen that the VIS channels mostly represent the lower- and medium-height (400–800 hPa) water clouds. The NIR channel is a good discriminator between ice clouds (<400 hPa), which appear dark because ice absorbs radiation stronger than liquid water at 1.6 μm and the water clouds, which appear bright. In particular, the thunderstorm cells can well be detected in the NIR. This may be a desirable feature, since it provides additional information on the location of clouds while making a clear distinction between high ice clouds and low/medium water clouds.
Figure 5 depicts the relative differences in reflectance between the 3D and the 1D calculations of channel VIS008 at 1200 UTC in each pixel (i, j) of the evaluated domain as an example of the effect of the parallax correction. Without the correction, large differences are present near the edges of cloud structures. These differences are substantially reduced by applying the parallax correction in the 1D calculation. As a comparison, Fig. 5 also contains the 3D and parallax-corrected 1D reflectance fields. It seems that the most severe relative differences occur at higher latitudes, in particular at sharp northern cloud edges where ice clouds are present. A reasonable explanation for this is the fact that the southern position of the sun at noon produces the largest shadows north of the high clouds.
In addition, we separately analyzed areas where the differences are largest, that is., at cloud edges. For this investigation, we have applied a threshold considering only those pixels in which the difference between the 3D and 1D reflectances |ΔR| > 0.1. For these pixels with a large difference, the effect of the parallax correction is stronger and the mean relative difference between 3D and ICA reduces from 31% to 23% with the parallax correction at 1200 UTC.
A histogram of the relative differences between 3D and parallax-corrected 1D reflectances at 1200 UTC for channel VIS008 is depicted in Fig. 6. The Gaussian fits included show that the differences deviate somewhat from a Gaussian distribution. One obvious reason is the higher peak around zero, which arises from clear-sky and homogeneously cloudy regions. In such regions, the 1D simulation is nearly a perfect method and as good as the 3D simulation. Hence, the height of the peak depends on the cloud cover of the simulated scene. Also, there are some events at large multiples of the standard deviation that broaden the Gaussian fit. In particular, about 2% of the differences are outside of the 3σ range. The observed deviation from Gaussian error statistics is, however, expected.
To demonstrate how the synthetic scenes simulated from model output look compared to real observations from MSG-SEVIRI, we provide a time sequence of observations and simulations on 22 June 2011 in Fig. 7. The SEVIRI observations of channels VIS008 over the diurnal cycle are depicted in the top row, the middle row displays the 3D simulations from 3-h forecast fields, and the bottom row shows the parallax-corrected 1D simulations from 3-h forecast fields. Overall, both 1D and 3D synthetic satellite images look realistic with the exception of the 1800 UTC scene, where missing shadow effects in the 1D operator lead to unrealistic structures. These missing shadow effects also led to large mean deviations of 1D and 3D results (Table 2).
On this particular day, the model forecasts contain substantially more clouds than the observations, particularly in the morning scenes. These discrepancies are clearly higher than the observation error and the estimated operator error, and thus reflect errors in the representation of clouds in the model forecasts. The developed forward operator can therefore also be used as a tool to identify potential model weaknesses. To evaluate this in more detail, however, requires the systematic comparison of a longer time period, as the interpretation of individual scenes may be misleading. Such an evaluation of systematic and stochastic differences for a longer period with the goal to identify model deficiencies in the representation of clouds is ongoing and will be subject of a follow-on publication.
Furthermore, Fig. 7 illustrates the differences between synthetic images from the 1D and 3D operators, which strongly depend on the SZA of the respective scene. For a smaller SZA (around 0900 or 1200 UTC; local noon is around 1115 UTC), it is hard to tell the difference between the two. With increasing angles at, for example, 0600 and 1500 UTC, shadow effects become more obvious in the output of the 3D operator and at 1800 UTC they lead to comparably large differences, as described before.
The largest deviations of observed and simulated imagery clearly result from the different location (or existence) of clouds in the model forecast, and reality and correcting these errors is therefore the main intention for assimilating such observations. Pixels that are cloud free in both the model forecast and reality lead to comparably similar results, reflecting that other operator error sources, for example, albedo or aerosol assumptions, are second-order effects. Different cloud types in the forecast and reality, for example, a semitransparent cirrus cloud instead of an opaque water cloud with very high reflectance values, can obviously lead to differences, but nevertheless these are still much smaller than the signal of cloudy-sky versus clear-sky values.
6. Summary and outlook
This article introduces an observation operator for visible and near-infrared satellite reflectances. The operator is intended as a fast enough tool to study the impact of directly assimilating visible and near-infrared observations within the experimental KENDA–COSMO system of DWD (or other DA systems that do not require a linearized and adjoint operator). Since particularly water clouds have a clearer signal at these wavelengths, it seems to be a natural extension to include such observations as a valuable source of cloud information. In addition to introducing the technical aspects of the forward operator, we have evaluated its accuracy with respect to a computationally expensive Monte Carlo radiative transfer model.
Moreover, a parallax correction is introduced that corrects 1D simulations for the slant path of radiation through the atmosphere toward the observing satellite. The accuracies of the independent-column calculation and its parallax-corrected version are evaluated by comparison to 3D Monte Carlo simulations. The latter are considered “perfect” model simulations due to their ability to account for arbitrarily complex cloud structures and corresponding shadow effects. Furthermore, the effect of horizontal averaging of the 3D and 1D reflectance fields over both 3 × 3 pixels and 5 × 5 pixels is evaluated to investigate the sensitivity of operator accuracy to resolution. The input fields are 3-h forecasts of the limited-area COSMO model at 0600, 0900, 1200, 1500, and 1800 UTC 22 June 2011.
In summary, relative differences between 0600 and 1500 UTC are about 6%–8% without parallax correction for the visible channel VIS008 of MSG-SEVIRI with a central wavelength of 810 nm. Including the parallax correction in the 1D calculations improves these results to about 4%–6%. The horizontal averaging over 3 × 3 and 5 × 5 pixels gives a further improvement to a difference of less than about 5% and less than about 4.5%, respectively. This is because the averaging cancels out some of the horizontal variations on small scales. Since the effective model resolution is lower than the grid spacing, similar smoothing routines might be relevant for future assimilation experiments to reduce representativeness errors. In addition, given the deficiency of current models to capture every individual convective system, assimilating such observations at a reduced resolution may be a desirable approach. As examples, the differences in the two VIS and NIR channels of the SEVIRI instrument, VIS006 and NIR016, have also been evaluated at 1500 UTC of the same day. The results for VIS006 are similar to those for VIS008, while NIR016 is about 1% less accurate.
At 1800 UTC, the differences turn out to be substantially larger than between 0600 and 1500 UTC due to the larger SZA leading to an increase in cloud shadows. In the absence of further corrections that can account for these 3D effects in the faster 1D simulations, one can draw the conclusion that for the assimilation of VIS and NIR satellite reflectances, it is only sensible to assimilate when solar zenith angles are smaller than about 70°. Otherwise, they would need to be assimilated with a suitable adaption of the errors in the assimilation system.
Another error source that is currently neglected are aerosols. During certain conditions—for example, volcanic eruptions, large forest fires, or desert storms—the radiative impact of aerosols may be of similar magnitude (in the VIS and NIR spectral range) as that of clouds. While in central Europe such events are very rare and/or of very small horizontal extent for the operational practice, it could be useful to develop methods (using, e.g., a combination of different channels) by which the data assimilation system can differentiate such signals from those of clouds. Quality control methods that prevent the assimilation of such data if the probability of a contamination is particularly high could also be a solution. Similar strategies may have to be employed for the treatment of snow surfaces, whose radiative signal can be similar to that of low-level clouds in the considered frequency range.
A more general limitation of the forward operator accuracy is the simplified one-moment microphysics scheme, which computes particle size and density from a single cloud water variable (for liquid and frozen clouds, respectively). In reality there is more variability in these parameters depending on cloud age and cloud type. For the key issue of correcting location error (i.e., mismatches between the locations of observed and modeled cloud), these errors are probably not very decisive. For improving, for example, the ice water content in clouds, the adequacy of the employed microphysics scheme might need to be revisited.
In future studies, the 1D forward operator presented here shall be applied in the KENDA–COSMO system of DWD to study the impact of directly assimilating reflectance observations of MSG-SEVIRI solar channels. The presented 1D operator is sufficiently fast for such case study purposes in an offline calculation as opposed to the 3D operator (which runs on 37 processors with a computation time of about 12 h per scene). Nevertheless, a computation time of approximately 5–10 min per scene over the whole model domain (run on 37 processors) is beyond the limitations of an operational ensemble DA system. Thus, a second objective for future research is to test methods to accelerate RT in the VIS and NIR spectral range and to assess the respective loss in accuracy. We are currently working on radiation schemes that are more than two orders of magnitude faster than 16-stream DISORT, using alternatively a strongly modified 2-stream approach or a lookup table. The implementation and test of such solvers is ongoing research. In addition to assimilation experiments, the observation operator can also be used as a tool to identify model weaknesses, in particular, concerning the representation of clouds.
This study was carried out in the Hans-Ertel-Centre for Weather Research. This research network of universities, research institutes, and Deutscher Wetterdienst is funded by the BMVBS (Federal Ministry of Transport, Building and Urban Development). P. M. K. would like to thank Annika Schomburg, Christian Keil, and Axel Seifert for the discussions and help concerning COSMO-related issues. P. M. K. is indebted to Luca Bugliaro for providing relevant satellite data.
In this appendix, relevant formulas and physical constants used in the operator calculations are summarized. Note that the definitions used in the parameterizations of subgrid-scale quantities are adopted from the subgrid scheme of the COSMO model code and are stated here for completeness only.
With pressure P and temperature T given, the densities are determined through the equation of state for ideal gases:
For the gas constant of dry air, one can plug in the value Rd = 287.05 m3 Pa kg−1 K−1and for water vapor, it is given by Rυ = 461.51 m3 Pa kg−1 K−1.
The saturation vapor pressure over water and ice is given by the Magnus formulas (Sonntag 1990)
respectively, where the particular constants have been adopted from the COSMO model code. The approximate temperature ranges of validity of the Magnus formula lie in between −45° and 60°C over water and in between −65° and 0.01°C over ice. Furthermore, the saturation mixing ratios can be calculated as
from which one can derive the relative humidity φ = Qtot/Qsat using the total humidity mixing ratio Qtot = QV + QC + QI. For x one can plug in either water or ice.
In the case of a mixed state the gas constant is, strictly speaking, not a constant but rather depends on pressure and temperature. It is given by
and lies between Rd and Rυ. Another equivalent way to treat a mixed state is to use the virtual temperature TV = TR/Rd, where the gas constant Rd is in fact kept constant.
In the following, some definitions are introduced that are used in the parameterizations summarized in section 3. The total saturation mixing ratio is defined as a sum of water and ice contributions by
where the ice fraction is defined as
In addition to the mixing ratios, the COSMO model uses cloud fractions. The shallow convective cloud fraction in the subgrid scheme of the model is defined by
where the magnitude depends on the heights of the shallow convective clouds, is the top height, and is the base height. The latter fields are model output (m). Terms and are nonzero where the convection scheme produces shallow convective clouds. If the height of the considered layer lies between and , then Eq. (A7) is applied; otherwise, we set = 0.
This refers to a horizontally infinitely extended model atmosphere with parallel layers in which optical properties only vary vertically.