## Abstract

Visible satellite images contain high-resolution information about clouds that would be well suited for convective-scale data assimilation. This application requires a forward operator to generate synthetic images from the output of numerical weather prediction models. Only recently have 1D radiative transfer (RT) solvers become sufficiently fast for this purpose. Here computationally efficient methods are proposed to increase the accuracy and consistency of an operator based on the Method for Fast Satellite Image Synthesis (MFASIS) 1D RT. Two important problems are addressed: the 3D RT effects related to inclined cloud tops and the overlap of subgrid clouds. It is demonstrated that in a rotated frame of reference, an approximate solution for the 3D RT problem can be obtained by solving a computationally much cheaper 1D RT problem. Several deterministic and stochastic schemes that take the overlap of subgrid clouds into account are discussed. The impact of the inclination correction and the overlap schemes is evaluated for synthetic 0.6-*μ*m SEVIRI images computed from operational forecasts of the German-focused COSMO (COSMO-DE) Model for a test period in May–June 2016. The cloud-top inclination correction increases the information content of the synthetic images considerably and reduces systematic errors, in particular for larger solar zenith angles. Taking subgrid cloud overlap into account is essential to avoid large systematic errors. The results obtained using several different 2D cloud overlap schemes are very similar, whereas small but significant differences are found for the most consistent 3D method, which accounts for the fact that the RT problem is solved for columns tilted toward the satellite.

## 1. Introduction

Clear sky satellite radiances have become the most important observations for numerical weather prediction (NWP) systems, and there is a rising interest to use also cloud-affected radiances. Satellite observations are not only utilized for data assimilation (see Bauer et al. 2011a,b, and references therein) but also for model evaluation (see e.g., Bikos et al. 2012; Heinze et al. 2017) and retrieval algorithm validation (see e.g., Bugliaro et al. 2011). All these applications require a forward operator, which computes synthetic satellite images from the NWP model state. Forward operators have to be sufficiently accurate and consistent with the assumptions made in the NWP model (see the discussion in Senf and Deneke 2017). Moreover, for operational purposes, forward operators also have to be very fast and are strongly constrained in terms of memory consumption. Suitable forward operators for thermal infrared and microwave satellite instruments have been available for some time and have been continuously improved. In most of the weather prediction centers, either the RTTOV radiative transfer package (Saunders et al. 1999) or the Community Radiative Transfer Model (Han et al. 2006) is used. In contrast, sufficiently fast forward operators for solar channels are only in an early development state.

Thermal satellite channels provide mainly temperature and humidity information with high spatial and temporal resolution and have proven to provide valuable observations for data assimilation (DA). While clouds also affect infrared radiances and this cloud effect could be exploited for data assimilation (see, e.g., Harnisch et al. 2016), more detailed and complementary information about clouds is contained in visible satellite images. However, these images have not been exploited for data assimilation so far, mainly because generating their model equivalents was computationally too expensive. This is related to two fundamental differences between visible and infrared images. First, radiances in the solar part are dominated by multiple scattering, while scattering can usually be neglected in the thermal part. Therefore, computationally cheap methods developed for thermal infrared channels cannot be applied for visible channels. Second, infrared radiances are determined by thermal emission of the atmosphere and the surface, whereas visible radiances are determined by the reflection of directional radiation from the sun by the atmosphere and the surface. The position of the sun is not required for the computation of synthetic infrared images and assuming a plane-parallel atmosphere is a sufficiently good approximation in this case. In contrast, the sun–satellite geometry is decisive and 3D radiative transfer (RT) effects are much more important for the visible channels.

Accurate RT solvers taking 3D effects and multiple scattering into account are available for the visible spectral range [e.g., the Monte Carlo code for the physically correct tracing of photons in cloudy atmospheres (MYSTIC) 3D Monte Carlo solver; see Mayer B. 2009], but they typically take several CPU days to generate one satellite image. Forward operators based on standard 1D RT methods like the one developed by Kostka et al. (2014) in the framework of the Hans-Ertel Centre for Weather Research (HErZ) project (Weissmann et al. 2014) are less accurate and faster, but they are still too slow for operational purposes. Only recently sufficiently fast 1D RT methods for visible wavelengths based on lookup tables have become available (Jonkheid et al. 2012; Scheck et al. 2016a). While these methods allow for efficient treatment of multiple scattering, they do not yet provide a full solution for the problem of generating accurate and consistent synthetic visible satellite images. The two presumably most important aspects still missing are 3D RT effects and the consistent treatment of subgrid-scale clouds.

In this work we discuss and evaluate methods to improve the accuracy and consistency of a forward operator based on the fast 1D RT method of Scheck et al. (2016a) without strongly degrading its performance. The first topic we addresses is a 3D RT effect related to the fact that the radiance observed for a cloud-top surface changes when it is tilted toward or away from the sun. The second topic is the overlap of subgrid clouds. While subgrid cloud water content assumptions consistent with the NWP model have already been used in Kostka et al. (2014), the overlap of the subgrid clouds was not taken into account. Several articles have discussed how to include cloud overlap in large-scale models with the aim to improve heating rates and surface irradiances. However, it is not clear what is the most consistent and efficient way to implement cloud overlap in the context of synthetic visible satellite images generated from convection-permitting models with higher resolution. Here we compare several methods and discuss the differences between stochastic (STO) or deterministic (DET) implementations and the impact of 3D effects related to the slant satellite viewing angle.

The rest of the paper is structured as follows: In section 2 we discuss the NWP model, the RT method, the satellite data, and the test periods used in this work. Section 3 contains a description of the newly implemented operator features that take cloud-top inclination and subgrid cloud overlap into account. In section 4 we compare images generated by the operator to Monte Carlo 3D RT results, and in section 5 we compare synthetic images based on operational forecasts to 0.6-*μ*m observations of the SEVIRI instrument on Meteosat. Last, section 6 contains a summary and conclusions.

## 2. Data and methods

### a. COSMO forecasts

The NWP model states for which we compute synthetic satellite images are 3-hourly operational forecasts of the German Weather Service [Deutscher Wetterdienst (DWD)]. These forecasts have been generated using the nonhydrostatic limited-area COSMO community model in its German-focused convection-permitting model configuration COSMO-DE. This model configuration will also be employed for future DA experiments we intend to perform with the Kilometre-Scale Ensemble Data Assimilation (KENDA) system (Schraff et al. 2016) using an implementation of the local ensemble transform Kalman filter (LETKF; see Hunt et al. 2007). The model domain consists of 421 × 461 grid points with a horizontal grid spacing of 2.8 km and 50 layers in the vertical. The model domain covers Germany, Switzerland, and Austria, and parts of their neighboring countries. The model allows for deep convection on the grid scale and has a parameterization for shallow convection (Baldauf et al. 2011). The heating and cooling rates due to radiation are computed according to Ritter and Geleyn (1992) using a two-stream solver. In this process subgrid clouds are assumed to exist in grid cells where the relative humidity exceeds a height-dependent critical value (see Senf and Deneke 2017).

The period from 28 May to 30 June 2016 is used to evaluate the operator. For some additional tests we consider also a shorter period in 2012 (10–28 June).

### b. SEVIRI observations

The Spinning Enhanced Visible and Infrared Imager (SEVIRI) aboard the geostationary Meteosat Second Generation (MSG) satellites operated by EUMETSAT provides full-disc observations in 11 spectral channels in the visible and infrared spectral range every 15 min. Here we focus on the 0.6-*μ*m satellite images from *Meteosat-10*, which was located at 0° longitude during the 2016 test period. The spatial resolution of the visible and near-infrared channels is 3 km × 3 km at the subsatellite point and about 6 km × 3 km in the COSMO-DE domain.

To derive albedo values required for the RT calculation, we use also the SEVIRI Clear Sky Reflectance Map (CRM) product (EUMETSAT 2015, 266–267), which is computed as the average reflectance value for all cloud-free cases found in the past seven days. If no cloud-free cases occurred in this period, then data from the next earlier 7-day period or a climatological value is used. The clear sky product is generated for every day for 1200 UTC and on Wednesdays also for every even hour between 0600 and 2000 UTC. We derive a time-dependent albedo from the CRM by starting from a value of zero and iteratively adjusting the albedo value until the resulting reflectance agrees with the clear sky reflectance. Linear interpolation is used to compute albedo values for times where no CRM is available. The algorithm used to compute the CRM product often causes patches with significantly too-high reflectances. Using the minimum CRM value within the test period (at the same time of the day) for each pixel strongly reduces these artifacts.

### c. Radiative transfer and optical properties

The operator version of Kostka et al. (2014) made use of the discrete ordinate method (DISORT; Stamnes et al. 2000) to compute reflectances. Here we use the newly developed Method for Fast Satellite Image Synthesis (MFASIS; Scheck et al. 2016a), which is based on compressed reflectance lookup tables (LUTs) computed with DISORT. The 21-MB LUT discussed in Scheck et al. (2016a) is used. MFASIS is about four orders of magnitude faster than DISORT and only slightly less accurate. As in Kostka et al. (2014), a parallax correction is performed to take the satellite viewing angle into account. NWP model data from columns tilted toward the satellite are used as input for the 1D RT solver [the tilted independent column approximation (TICA) approach; see Wapler and Mayer 2008].

For each atmospheric column, MFASIS requires the vertically integrated optical depths and mean effective particle sizes of water and ice clouds as input variables. These variables are computed in the same way as in Kostka et al. (2014). The effective droplet radius is determined following Martin et al. (1994) and for the effective ice particle size we use the parameterization of Wyser (1998). Water and ice contents are converted to optical depths using the parameterizations of Hu and Stamnes (1993) and Fu (1996), respectively. We used subgrid cloud water and ice contents and cloud fractions that are consistent with the assumptions used in the internal RT code in the COSMO Model. In the current version, aerosols are not yet taken into account.

## 3. Operator improvements

### a. Cloud-top structure

Synthetic satellite images in the visible and near-infrared obtained by 1D RT methods show much less structure within the clouds than images generated using 3D RT methods or real satellite images (see Fig. 6 in Kostka et al. 2014), in particular for large solar zenith angles (SZAs). The missing structure can thus be considered as a systematic error that is mainly related to 3D RT effects. Three-dimensional RT methods like Monte Carlo ray tracing are orders of magnitude too slow for operational purposes. However, the 3D RT effect dominating the cloud structure is rather simple and can be approximately taken into account with a computationally cheap method. The clouds are darker where the cloud-top surface is tilted away from the sun, and they are brighter where the surface is tilted toward the sun. This effect is simply related to the fact that the same number of photons from the sun is distributed over a large or a smaller cloud surface area. The orientation of the cloud surface with respect to the satellite is also important, as it determines whether the surface is visible from the satellite and which solid angle it occupies.

This variation in brightness is a local effect in the sense that to first order it depends on only the direct radiation from the sun and the local morphology of the cloud. Only when the cloud surface is tilted so far away from the sun that the diffuse illumination of the cloud is of similar importance as the direct radiation, nonlocal effects must not be neglected anymore and the RT problem becomes much more complicated. However, for not-too-extreme SZAs, this is a rather rare case in geostationary satellite images.

To quantify the change in reflectance caused by a tilting of the cloud surface by a zenith angle into azimuthal direction , it is useful to consider a rotated frame of reference, in which the normal vector of the cloud surface points toward the zenith. For the moment we ignore all clear sky effects (e.g., Rayleigh scattering) and assume that only sunlight reflected from the ground and the clouds contributes to the radiance observed by the satellite. Under these assumptions and as illustrated in Fig. 1, the radiance emanating from the tilted cloud measured for a sun zenith angle and a satellite zenith angle *θ* should be the same as the radiance emanating from a ground-parallel cloud for modified sun and satellite zenith angles and , respectively. The sun and satellite positions relative to the cloud surface normal vector are the same in Figs. 1a and 1b, so the scattering angle does not change. Also, the length *L* of that part of the viewing path that lies within the cloud is the same. The optical thickness of the uniform cloud with extinction coefficient *β* in the nonrotated frame is given by

As *L* and *β* have the same value in the rotated frame, the optical thickness there is given by , where .

The problem in Fig. 1b is nearly a 1D RT problem, except for the fact that the ground is tilted. We assume that this difference can be taken into account by using a modified surface albedo . The radiance in the case of a cloud tilted by an angle into direction can thus be written as the radiance for a ground-parallel cloud with modified parameters,

where the modified viewing angles and depend on the original angles and the inclination angles, respectively.

Using the definition of reflectance,

where is the extraterrestrial solar flux, and Eq. (2) can also be written in terms of reflectances,

where . Without clouds () and atmosphere, the reflectance is equal to the surface albedo, , which requires that .

Up to now we have ignored the clear sky contributions to the reflectance. The latter can become important for large zenith angles, which means that the optical path through the clear sky atmosphere is large. For geostationary satellite images, *θ* is usually not very large and even for the clear sky contribution to reflectances in the 0.6-*μ*m channel we are focusing on is typically at most a few . However, also for small zenith angles *θ* and , the modified zenith angles and associated with tilted clouds can be large. If we used the standard MFASIS LUT, which includes clear sky effects, to determine in such a case, then an erroneously strong clear sky signal would contribute to the reflectance. Therefore, it is necessary to correct Eq. (4) such that it includes the clear sky signal for the original zenith angles and a cloud signal computed using the modified zenith angles. To separate cloud- and clear sky contributions, an additional MFASIS LUT was computed that does not contain any clear sky contributions. In the following, we will denote reflectances computed with this cloud-only LUT as and reflectances including clear-sky contributions are denoted as . The additional table allows us to estimate the reflectance for a tilted cloud in a Rayleigh atmosphere as

where the reflectance correction term

accounts for the inclination of the cloud and *f* is a function of the optical depth that reduces the correction to zero for thin clouds. This is necessary, as only thick clouds behave similar to solid surfaces, in the sense that they follow Eq. (5) for . For thin clouds there is no well-defined surface close to which most of the incoming solar photons are scattered and thus the inclination correction is not valid any more. A comparison with 3D Monte Carlo simulations (see section 4) suggest that should be scaled with the logarithm of the optical depth. Therefore, we adopt

and empirically choose and .

In principle another correction term accounting for deviation of the actual cloud-top height from the one used in the calculation of the MFASIS LUT could be added to Eq. (5). This would require an additional LUT dimension and for the 0.6-*μ*m channel we focus on in this work, the reflectance correction would be very small. Therefore, we omit this term here. For channels more strongly affected by water vapor absorption or Rayleigh scattering, it could be necessary to include such a term.

Before Eq. (5) can be evaluated, one needs to define a cloud-top surface, to compute its tilt angles and to calculate the modified zenith angles and . For each tilted column of the model grid, we determine the height where the optical depth is 1. The local cloud-top surface is assumed to be a plane that is determined by a least squares fit to the points in the current column and the eight neighboring columns. For the COSMO-DE grid considered here, the eight neighboring columns are located ±2.8 km away from the central column in the two grid directions. Neighbor columns with a total optical depth less than one are excluded from the plane fit. If fewer than five points are available, we consider the orientation of the plane not to be sufficiently well defined and assume that the plane is parallel to the ground. Using the normal vector of the plane obtained from the least squares fit, , one can compute the tilt zenith angle . To determine the modified angles and , we first have to define the unit vectors pointing at the satellite and the sun as

where *ϕ* and are the satellite and sun azimuth angles, respectively. The unit vectors , , and define a Cartesian basis for a frame of reference in which one of the dimensions is aligned with the cloud normal vector. The second vector, , could be chosen as any vector perpendicular to , because the 1D RT problem is invariant with respect to rotation around the vertical axis. The rotation matrix transforms to in the nonrotated frame (Fig. 1a). As is visible in Fig. 1b, a rotation in the opposite direction is required to transform the sun and satellite vectors into the modified sun and satellite vectors in the rotated frame, that is, and Here, and are defined analogously to Eq. (8), but with the modified (primed) angles. The modified zenith angles can thus be computed as

Finally, in the presence of water and ice clouds with surface inclinations that are in general different, we need a way to combine two correction terms of the form Eq. (6). As MFASIS cannot by design treat ice and water clouds completely independently, an approximation is required here. We simply assume that the water cloud correction should be scaled down where the ice cloud correction is more important and adopt

where and are the column-integrated water and ice cloud optical depths, respectively. The correction terms and are computed using the azimuthal angles [Eq. (9)] relative to the water and ice cloud surfaces, respectively.

### b. Point spread function

The radiance attributed to a pixel of a SEVIRI image also depends weakly on the radiance emanating from the neighboring pixels, as described by the point spread function (PSF) of the instrument. At the nadir the SEVIRI point spread function has been determined to be

where *d* is the distance in pixels and (Deneke and Roebeling 2010). By convolving the synthetic reflectance field with this function, the smoothing influence of the point spread function can be taken into account.

The albedo based on the SEVIRI clear sky reflectance product contains already the effects of the point spread function. Therefore, cloud-free regions should be excluded from the smoothing. This is achieved by masking out pixels where

in the pixel and the eight neighboring pixels.

Applying the point spread function is particularly important when the cloud-top inclination is taken into account. Otherwise, sharp features with a too-high contrast related to locally tilted cloud tops appear. Applying the PSF blurs these sharp features and leads to minimum and maximum reflectance values that are in better agreement with the observations.

### c. Cloud overlap schemes

NWP and climate models cannot resolve clouds on scales smaller than their grid length. As the influence of these clouds is not negligible, the models rely on parameterizations to account for these unresolved clouds on subgrid scales. Often these parameterizations generate cloud water and cloud ice content in cells where the grid-scale cloud water variable is zero but the relative humidity is higher than a height-dependent critical value (for an overview, see Quaas 2012). For global and convection-permitting models these partially cloudy grid cells, in which the subgrid cloud fraction variable *γ* is between zero and one, constitute an important contribution to the total cloud cover.

When several partially cloudy grid cells are present in the same column of the model grid, the radiative transfer calculation requires assumptions on how subgrid clouds overlap. The operator version discussed in Kostka et al. (2014) made use of the subgrid cloud water and ice content parameterizations implemented in COSMO. However, in contrast to the RT implementation in COSMO, it was assumed that the subgrid clouds fill the whole grid cell, that is, . This simplification meant that the issue of cloud overlap could be ignored. The new operator version discussed here is more consistent and it takes the cloud cover parameterization used for the radiative transfer code in COSMO into account.

The simplest overlap assumptions are that the cloud overlap is either random or maximum. For maximum cloud overlap, the total cloud cover of *n* layers with cloud cover is assumed to be

For two randomly overlapping clouds, the probability to encounter at least one of them at a randomly chosen location is . The probability to miss both of them is thus . Generalizing this result for *n* layers yields

where is the mean total cloud cover of many different realizations of the randomly overlapping clouds. A slightly more complex and often used overlap assumption that is also implemented in the internal RT code in COSMO (Ritter and Geleyn 1992) is the “maximum random” cloud overlap. Clouds in adjacent layers are assumed to follow the maximum overlap rule, whereas clouds separated by a cloud-free layer are assumed to overlap randomly. A formula approximating the mean total cloud overlap for this case is given by

While more sophisticated cloud overlap assumptions exist (see, e.g., Di Giuseppe and Tompkins 2015, and references therein), we restrict ourselves here to the maximum-random assumptions that are consistent with the internal two-stream RT code in COSMO. The response of a global model to changes between these assumptions that were used in the computation of heating rates has been investigated by Morcrette and Jakob (2000). Here we are interested in convective scales and the impact of taking the slant satellite viewing angle into account, which can be ignored in global models. We implemented four different methods to take maximum-random cloud overlap into account. These methods differ in the assumed subgrid cloud size distribution, the required computational effort, and whether the stochastic nature of the subgrid clouds and the slant satellite viewing angle are considered. In addition, an implementation of the random overlap assumption will serve as a reference. In each of these methods, the column containing partially cloudy layers is subdivided into *N* subcolumns in which the cloud fraction is either 0 or 1 in each layer. The reflectance for the column is computed as the weighted average of the reflectances for the subcolumns. Examples of subcolumn cloud configurations generated by the different methods are given in Fig. 2 for a column containing four partially cloudy layers.

The first method is deterministic and aims to estimate the mean reflectance resulting from all subgrid cloud realizations that are allowed according to the cloud overlap rules. Method DET is the “streams” approach (Matricardi 2005) used in the RTTOV package. Here the total cloud cover is computed for each layer according to Eq. (15). Adopting a dimensionless coordinate that spans the grid cell, the cloud in layer *k* is assumed to extend from to . Each cloud is thus placed at the rightmost location that does not cause the total cloud cover to exceed (see Fig. 2a). The grid cell is then divided into *N* streams or subcolumns, such that every layer in every subcolumn is either completely cloud covered or cloud free. This process can lead to a rather high number of subcolumns. In the worst case, subcolumns have to be considered, where is the number of model levels. The increased memory requirements associated with a high number of subcolumns is a disadvantage of the DET method. A fixed subcolumn number version that avoids this problem is discussed in Scheck et al. (2016b). Not all physically reasonable cloud configurations are considered in the DET approach. For instance, cases where the cloud in km km] overlaps with the cloud in km, 6 km] for the example in Fig. 2 would not be prohibited by the maximum-random overlap rules, but they are not accounted for. DET thus underpredicts the degree of overlap of nonadjacent clouds. In contrast, the internal two-stream code in COSMO accounts for this case. Therefore, reflectances computed with the DET method are not fully consistent with the two-stream code, even though they exactly reproduce the total cloud cover given by Eq. (15).

In addition to this deterministic cloud overlap scheme, several stochastic schemes were implemented. These schemes generate only one random subgrid cloud realization compatible with the overlap rules. It should be noted that by design these stochastic schemes are not differentiable, which makes them unsuitable for variational DA methods. However, there is a general trend toward DA methods that do not require an adjoint model and exact differentiability and use information derived from ensembles instead (see, e.g., Gustafsson et al. 2018). These methods are based on the assumption that changes in the model state and the operator output can be described approximately by a linear relation. This may not be the case for the reflectance of a single pixel that depends on random numbers used in a stochastic scheme. However, one can use the standard approach of averaging over several pixels to obtain a “superobservation” that is nearly free of stochastic noise and depends on a sufficiently linear way on the model state.

Methods STO-N and STO-C (refer to the appendix) use a uniform decomposition of the column into *N* subcolumns, where, in contrast to the DET scheme, *N* can be chosen a priori. Clouds are placed into the subgrid cells according to nondeterministic rules. Details on the algorithms used in these methods are given in the appendix. Here we discuss only the most important differences. Method STO-N computes the probability that layer *k* is filled with a cloud independently for each subcolumn *i*. Averaged over many realizations (see Fig. 2c for an example), STO-N yields the correct mean cloud cover in each layer, and the total cloud cover given by Eq. (15) is reproduced. STO-N takes all cloud overlap configurations allowed by the maximum-random overlap rules into account and is in this sense more consistent than the DET method. The fact that all subcolumns are treated independently has two consequences. First, subgrid clouds spanning several adjacent subcolumns form only by chance and for the average subgrid cloud size goes to zero for . Second, computing the average reflectance for *n* realizations with *N* subcolumns requires exactly the same operations as computing the reflectance for one realization with subcolumns. Therefore, also single realizations converge toward the ensemble mean for . The fact that all realizations converge toward the same reflectance values implies also that the spread vanishes for . STO-N thus generates a spread that is related to only discretization errors, not to the uncertainty in the subgrid cloud configuration. In contrast to the noncontinuous clouds in STO-N, method STO-C generates in each layer *k* one continuous subgrid cloud (see Fig. 2c) with a cloud fraction that agrees with (except for a discretization error ). For , the cloud fraction and the subgrid cloud size thus converge toward . However, STO-C does not exactly reproduce the total cloud cover as defined by Eq. (15). In contrast to STO-N, the ensemble spread converges toward a finite value for that quantifies the uncertainty in the reflectance related to the unknown subgrid cloud configuration. Method STO-R is implemented very similar to STO-C, but it assumes the random overlap rules. For this purpose the cloud position in each layer is chosen randomly (see Fig. 2d).

Finally, the last cloud overlap scheme, STO-C3D, is also stochastic and takes into account that the maximum-random overlap assumption, as used in the internal RT code in COSMO, applies for vertical columns, whereas the operator makes use of TICA (see Wapler and Mayer 2008) to account for the satellite viewing angle *θ*. A group of maximum overlapping cloud layers should be stacked on top of each other in the vertical direction, not along the tilted column. Descending in the tilted column by one layer of height *z* causes a horizontal displacement away from the satellite. To compensate for this displacement, the subcolumn clouds in each layer have to be shifted toward the satellite by the same distance. The fact that the compensating displacement must occur in a certain direction implies that this problem is intrinsically three-dimensional. Therefore, two subcolumn dimensions are required. We consider a bundle of subcolumns that are tilted by *θ* toward the satellite, into a direction that is given by the projection of the satellite vector onto the surface (see Fig. 3). The correction for the tilt has to be applied in this *x* direction but not in the perpendicular direction . Clouds are placed into the layers of these subcolumns in way that is very similar to method STO-C (for details see the appendix), except for two differences. First, within a group of maximally overlapping layers, the clouds are additionally shifted by in the *x* direction. Second, for two-dimensional cloud layers an assumption has to be made about the horizontal cloud aspect ratio. Here we use the simplest choice, an aspect ratio of one.

It should be noted that according to the reciprocity principle not only is the slant viewing path important but also the slant path of the photons arriving from the sun, which means in particular that cloud shadows must be taken into account. With the current version of MFASIS this is not yet possible. Thus, STO-C3D is more consistent than the other schemes, but it does not yet provide a full solution for the problem. However, STO-C3D would provide a full solution in this sense for thermal channels, for which the direct radiation from the sun is irrelevant.

## 4. Comparison with 3D Monte Carlo results

### a. Cloud-top inclination

To verify that the cloud-top inclination correction discussed in section 3a improves the results, synthetic MFASIS images generated with and without the correction are compared to reference results obtained with the 3D Monte Carlo code MYSTIC (Mayer B. 2009).

For this comparison output from the Icosahedral Nonhydrostatic (ICON) model runs with a 150-m grid resolution performed in the framework of the High Definition Clouds and Precipitation for Advancing Climate Prediction [HD(CP)^{2}] project was mapped onto a regular latitude–longitude grid with a grid spacing of about 1 km and used to compute input data for the RT. Details on the model run settings and a model evaluation can be found in Heinze et al. (2017). Several simplifications were adopted for the RT computations. Cartesian geometry and constant viewing and sun angles are assumed. As MYSTIC does not support subgrid clouds, the latter were removed. For this purpose, in grid cells with a cloud fraction higher than 0.5 the cloud fraction was set to one and in all other cells to zero. Because of the computational effort required for 3D Monte Carlo simulations, we considered only 7 days (listed in Fig. 7a). The cloud situation varies strongly from day to day and also within the domain and includes both water and ice clouds. For each day only the 1200 UTC model output was used, but we considered three different sun angle combinations that roughly correspond to the situation at 0600, 0900, and 1200 UTC in June. In the following 0600 UTC means , 0900 UTC means , and 1200 UTC means , where is counted clockwise from south. Per day and angle combination synthetic images for two rectangular regions with a size of about 250 km × 250 km were computed, resulting in cases.

From the synthetic images for a 0600 UTC example case shown in Fig. 4, it is obvious that much more structure is visible in the clouds when the inclination correction is applied. Moreover, the pixels that become brighter as a result of the correction seem to be strongly correlated with the bright pixels of the Monte Carlo solution. It is a common feature of all cases that the mid- to high-reflectance regions become more similar to the reference solution. This effect is also visible in the reflectance histograms computed for all the 0600, 0900, and 1200 UTC cases shown in Fig. 5. For all solar angle settings, the peak at a reflectance of about 0.75 (a typical value for opaque clouds) that is present in the 1D results becomes smaller and broader when the inclination correction is used. In particular, for reflectances larger than 0.75, the histogram curve is much closer to the reference solution when the correction is applied. Both the maximum reflectance and the slope of the histogram for reflectances higher than about 0.7 fit much better to the MYSTIC results when the correction is used. To quantify the differences between the histograms, we define the *histogram error * as the area between the histogram curve for a reflectance field and the one for a reference solution divided by the area under the reference solution histogram,

Here, *i* is the reflectance bin index and we consider 140 reflectance bins evenly spaced between reflectance values of 0 and 1.4. For all cases the histogram errors of the MFASIS solution with respect to MYSTIC is strongly reduced when the inclination correction is used (Fig. 6). Also, the reflectance RMSE with respect to the MYSTIC results should be reduced by the correction. As Fig. 7 shows, this is the case for nearly all of the RT experiments. However, the reduction of the reflectance RMSE is much less significant than the one for the histogram error. This is related to the fact that 3D RT effects not related to the cloud-top inclination contribute to and, at low to midreflectances, dominate the total reflectance error. These effects cause also significant deviations of the MFASIS curves from the MYSTIC curves in Fig. 5 at low to midreflectances. For larger SZA one would expect 3D effects to be more important and indeed the RMSE and the reflectance bias (Fig. 8 are significantly larger for the 0600 UTC case than for the 0900 and 1200 UTC cases. It should be noted that if these effects are approximated by future developments, then probably the parameters and from Eq. (7) could be chosen differently and the histogram shape could be improved further without increasing the RMSE.

Larger differences between MFASIS and MYSTIC often result from the violation of the assumptions made for independent column 1D RT or the inclination correction. In 1D RT it is assumed that the clouds extend to infinity in horizontal direction. No photons can enter or leave the cloud through the cloud’s sides. In reality and in 3D Monte Carlo simulations, there is a photon flux through the cloud’s sides, which can lead to a positive or negative reflectance bias, compared to 1D RT. In particular, clouds with a high aspect ratio (height divided by horizontal size) and intermediate optical thickness are problematic in this respect, as the cloud’s sides contribute a larger fraction to the total cloud surface area and photons entering through the sides have a higher probability of reaching the cloud top than in very thick or very thin clouds. Moreover, in the independent column approach, clouds in one column cannot cast shadows on other columns. In particular for large SZAs and clouds with a high aspect ratio, not accounting for the extinction of direct radiation by clouds can lead to a positive bias, which can also be seen here in Fig. 8. Finally, the inclination correction assumes there is a sufficiently well-defined cloud surface and that the inclination angle is smaller than 90°. These assumptions can be invalid, for example, in pixels where mainly cloud sides are visible.

While these problems exist and should be addressed by future developments, our results indicate that for most optically thick clouds and clouds with a not-too-extreme aspect ratio, 1D RT together with the inclination correction is a reasonable and computationally much cheaper approximation to Monte Carlo 3D results and thus presents a significant improvement compared to pure 1D RT. It should also be noted that the seven test days used here were chosen to be simulated in the framework of the HD(CP)^{2} project because they are characterized by strong convection or other extreme weather events. Therefore, cases that should be less sensitive to the missing 3D RT effects discussed above—for example, mainly stratiform clouds—are underrepresented in this sample. The mid- to low-reflectance bias discussed above should thus on average be less severe for samples randomly drawn from the climatology (see section 5).

### b. Cloud overlap

It would be interesting to test also the overlap schemes discussed in section 3c by means of a comparison with Monte Carlo results. Such a comparison could indicate how important it is to take not only the slant viewing path into account but also the slant path to the sun for partially cloudy columns. This would require constructing an ensemble of high-resolution model states with known cloud overlap statistics below a certain spatial scale, to compute MYSTIC images for this ensemble and to compare them to MFASIS images computed for a version of the model state that has been coarsened to that scale. However, because of the bias related to missing 3D effects discussed above, in particular the missing cloud shadows, such a comparison would favor the overlap implementation that best compensates for the existing bias and not the one providing the best representation of the cloud overlap statistics for the high-resolution data. Therefore, we refrain from a comparison with MYSTIC here and will just estimate the potential impact of taking the slant viewing path into account in the next section.

## 5. Comparison with SEVIRI observations

### a. Setup and metrics

The impact of the operator improvements described in section 3 is investigated using 0.6-*μ*m SEVIRI observations and model equivalents computed from COSMO-DE forecasts for a 34-day test period in May–June 2016. To reduce potential lateral boundary effects, we consider only the part of the COSMO-DE domain that is at least 12 grid cells away from the boundaries for the evaluation of the results. Moreover, to avoid problems with snow and steep orography, we exclude the region with longitudes between and and latitudes less than , which contains the Alps.

Several metrics are considered to evaluate how well the synthetic satellite images agree with the observed ones. We define the cloudiness *C* as the fraction of pixels in which the reflectance is higher than a critical value *R*_{c} of 0.2, which is an upper limit for the clear-sky reflectance in the considered domain for the 0.6-*μ*m SEVIRI channel. This simple definition is not perfect, as it may miss some thin clouds over darker ground. However, the error should be rather small and this definition has the advantage that cloudiness can be computed in a consistent way for both observed and synthetic images. The cloudiness bias is the difference between cloudiness values of the observation and the model equivalents. The reflectance bias is defined as the difference between the mean reflectances for the observations and the model equivalents. The RMS reflectance error would be an obvious choice to measure how strongly the synthetic images deviate from the observed ones. However, this quantity is likely to be dominated by displacement errors, that is, the fact that clouds in the model are not at the same positions as in reality. As this study is not focused on displacement errors, which are related to the NWP model and not the operator, we consider a quantity that is not influenced by these errors. The reflectance histogram for a given satellite image does not depend on the exact cloud positions, but it still contains information about how much of the domain is covered by thin and dense clouds and is also sensitive to 3D RT effects. The histogram error defined in Eq. (16) is used to measure the differences in the histograms for observed () and synthetic () reflectance fields. For one given time, the histogram error can still be strongly dominated by errors in the model state, for example, when the cloudiness is wrong because a front enters the domain a few hours later in the model run than in reality. Averaged over many instants of time these error should be small, if the model does not have a significant bias. Therefore, reflectance histograms and computed for the full 34-day test period are used in Eq. (16). However, as the systematic error of the operator is likely to depend on the SZA, we still distinguish between different times of the day here and compute using data from all days but only one time of the day. In reality models have biases, which can also affect (see section 5b). As discussed in section 3, the inclination correction is applied only for columns containing sufficiently dense clouds, where the reflectance should not strongly depend on the cloud overlap assumption. Therefore, we discuss the impact of the inclination correction and the choice of the cloud overlap scheme independently in the following sections.

### b. Cloud-top inclination

An example of images computed with and without taking the cloud-top inclination (CTI) into account for a 0600 UTC case can be found in Fig. 9. To make it easier to distinguish clouds from the ground, not only the 0.6-*μ*m but also the 0.8*-μ*m channels have been used to generate these images (see figure caption). Without CTI, dense clouds are visible as uniformly white areas without features (Fig. 9c). With the inclination, much more structure is visible in the synthetic image (Fig. 9b). Although some variability on the finest scales is still missing, the synthetic image is qualitatively much more similar to the observation than the pure 1D result. For instance, a band of high clouds starting at the middle of the left domain boundary is barely visible without the inclination correction and much easier to see with the correction. A similar feature is present in the SEVIRI image. Moreover, it is now possible to distinguish convective from stratiform clouds and to see boundaries between low and high clouds. For smaller SZAs, the impact of the CTI correction is weaker but still present.

Taking the CTI into account also leads to improved reflectance distributions, in particular for larger SZA. As is visible from the reflectance histograms for 0600 and 1800 UTC in Fig. 10, the slope of the histogram curve at high reflectances () is too steep for 1D RT but becomes less steep and shows better agreement with the observations at high reflectances when the CTI correction is applied. However, the improvements in the histograms are less pronounced than in comparison between MFASIS and MYSTIC (see Fig. 5). When the reflectance correction terms in Eq. (10) are multiplied by a factor of 3 (the curves marked CTI in Fig. 5), the deviations from the histogram for the observations become significantly smaller. This difference is probably related to the fact that for the MFASIS–MYSTIC comparison, both images were based on the same dataset with a resolution of about 1 km, whereas now we are comparing SEVIRI images based on nature to synthetic images based on model results with an effective resolution of about 15 km (Bierdel et al. 2012). The gradients of all model fields and also of the cloud-top height will be less steep in the COSMO state than in nature, which directly affects the inclination correction. The gradients related to structures larger than the effective model resolution should be present in the model state but are too weak. Sharpening the gradients, perhaps with a more sophisticated method than the multiplication with a constant factor used here, could be a way to compensate for this model deficiency and lead to synthetic images with increased information content that is more useful for DA. The optimal sharpening parameters will be model and resolution dependent and must be determined by means of tuning experiments. There are also other potential reasons for why the inclination correction without the factor of 3 seems to be too weak. The constants and in Eq. (7) may not be determined optimally and could lead to an amplitude of the correction term that is too small. Including further 3D effects, as discussed in section 4, could facilitate a better determination of these parameters or an improved replacement for Eq. (7). Furthermore, the least squares plane fit used to determine the angles relative to the cloud surface (see section 3) leads to relatively smooth results and may underestimate gradients.

While the slope of the histograms in Fig. 10 is improved by the inclination correction at high reflectances, in particular for 0600 and 1800 UTC, other problems are still present. The dense cloud peak is still located at too-high reflectances for the model results. Moreover, for 0600 and 1800 UTC there are too many clouds in the model, which is indicated by ground peaks at that are less pronounced than in the observations for these times. These systematic differences are independent of whether the CTI correction is applied, and they are also independent of the type of cloud overlap assumed (see next section). The ground peak problem could be related to a deficiency of the COSMO Model or the ICON model that was used to compute the boundary conditions (e.g., a microphysics problem) or to problems with the data assimilation, but it could also be simply an artifact related to a not sufficiently long test period. Further investigations are required to determine the cause of the problem and synthetic SEVIRI images may prove useful for this purpose.

Figure 11a shows more quantitatively that in particular for large SZA is reduced significantly by the CTI correction. The CTI correction thus has the potential to extend the daily time span in which the systematic error is sufficiently low to allow for an assimilation of visible satellite images. The cloudiness bias is not sensitive to the CTI (Fig. 11b) and the reflectance bias is in most cases reduced (Fig. 11c), except for 1800 UTC. In this case the positive bias could possibly be reduced by taking cloud shadows into account. It should be noted that looking just at the metrics presented in Fig. 11 is a bit misleading, as it shows primarily advantages of the CTI correction for large SZAs. However, even at noon often additional information about the cloud morphology is visible in the synthetic images. This increased information content could be valuable for data assimilation.

### c. Cloud overlap

The deterministic and stochastic cloud overlap methods discussed in section 3c differ in important aspects. Only the methods DET and STO-N adhere to the total cloud cover formula [Eq. (15)]; only STO-C, STO-R, and STO-C3D reproduce the desired cloud fraction profile in each realization; and only STO-C3D takes the slant satellite viewing angle into account. For certain profiles of the cloud fraction and optical depths, the methods can in some cases yield strongly differing results, as the idealized examples in section 3c demonstrate. Here, we investigate how much influence the choice of the overlap method has on synthetic satellite images of the test period, that is, for profiles that can be regarded to be more realistic. As discussed in section 3c, the 3D scheme is more consistent than the 2D schemes, but it is still restricted in the sense that it does not account for cloud shadows (like all the other schemes). To minimize this effect, we will mainly focus on only the 1200 UTC cases in the following paragraphs.

We will first consider only results obtained for subcolumns (except for DET, where the number of subcolumns is determined by the cloud fraction profile) for which discretization errors are negligible, as we will show later. Only as references, we consider operator results for which only grid-scale clouds are taken into account (GRD) and results where, as in Kostka et al. (2014), a subgrid cloud fraction of was adopted for all grid cells containing subgrid cloud water or ice (CF1). As is visible from the time evolution of cloudiness (Fig. 12a), grid-scale clouds alone lead to cloudiness values that are much too low compared to the observations. When subgrid clouds are added but their overlap is not taken into account (the CF1 case), significantly too-high cloudiness values result. If the random or maximum-random cloud overlap method is used, then the model cloudiness is much closer to the observed one. It is thus essential to take the overlap of the subgrid clouds into account. These conclusions could also be drawn from the reflectance histograms shown in Fig. 13. Taking cloud overlap into account causes the histogram for the synthetic images to agree much better with the one for the observations.

Despite all of the differences between the different overlap methods, the results for the cloudiness are rather similar. The average cloudiness bias in the test period depends more strongly on the time of the day than on the overlap method (Fig. 12b). The 2D stochastic maximum-random overlap methods and the deterministic method lead to nearly identical cloudiness biases (Fig. 12b) and reflectance histograms (not shown). The 2D random and the 3D maximum-random images are somewhat brighter than the 2D maximum-random images. However, these relatively small differences in global metrics are somewhat misleading. To get a clearer view on the differences between the overlap methods and separate systematic from random deviations, we performed additional experiments. For each of the stochastic methods, an ensemble consisting of 100 different realizations was computed for a shorter test period (10–28 June 2012) and only at 1200 UTC. An example of the ensemble mean reflectance computed with STO-C is displayed in Fig. 14a. For this case the difference between the STO-R and the STO-C ensemble means is shown in Fig. 14b. It is obvious that for some parts of the scene, the difference of the ensemble mean reflectances can exceed 0.1. However, only in about of the domain, mainly at the boundaries of dense clouds and in broken clouds fields, do the ensemble means differ significantly so that the total impact on the reflectance histograms seems weak. A better way to quantify such local systematic differences between two reflectance distributions and is to consider, for example, the 95th percentile of the deviation . In Table 1 we list the RMS difference, the bias, the 95th percentile, and the maximum reflectance deviation of the ensemble means obtained with the stochastic methods with respect to results obtained with the deterministic DET method. While the RMS differences and the biases are rather small, the 95th percentile of the deviation indicates that locally the reflectance distributions obtained for STO-C3D and STO-R differ by several from the DET values. These local differences are larger than the SEVIRI calibration uncertainty (which is a few percent; see Meirink et al. 2013) and thus could be significant. The results for the maximum-random STO-C3D method are in fact closer to the random overlap STO-R results than to the maximum-random DET results. The 95th percentile deviation of STO-C3D from STO-R is only , whereas the deviation from DET is . Taking the slant viewing path into account thus has a similar impact as changing the fundamental assumption about the overlap rules. While the exact number certainly would be different if cloud shadows were also taken into account in STO-C3D (see discussion in section 3c), it is unlikely that this would cancel out the impact of the slant viewing path. In contrast to the differences between DET and STO-C3D, and in terms of the deviation, the differences between DET and STO-C or STO-N are small. Neither the fact that not all possible cloud configurations are taken into account in DET nor the incompatibility of STO-C with Eq. (15) thus seems to have an important impact.

For data assimilation the ensemble mean value and the reflectance spread are of interest. The spread generated by a stochastic method quantifies the uncertainty related to the unknown subgrid cloud configuration. However, like the mean reflectance, the reflectance spread also can be affected by discretization errors related to the finite number *N* of subcolumns used in the stochastic schemes. As discussed in section 3c, the STO-N method implicitly assumes that the subgrid cloud sizes go to zero for and that in this limit the spread vanishes. This behavior is visible in Fig. 15a, which shows that not only the error of the mean reflectance but also the mean and maximum spreads decline continuously with an increasing number of subcolumns. The spread generated for finite *N* is thus just a discretization artifact. In contrast, the subgrid clouds in the STO-C method are as large as possible and the spread converges toward a nonzero value for . Figure 15b shows that the mean and maximum spreads decline faster than for STO-N at small *N* and do not change significantly for . Neither of the two assumptions about the subgrid cloud size is well justified. In broken cumulus cloud fields, the cloud size distribution often follows power laws [see, e.g., Heinze et al. (2017), and references therein] and thus lies between these two extreme cases. This suggests that the spread obtained using STO-C and a high number of subcolumns is an upper limit for the spread one would obtain for a more realistic subgrid cloud size distribution. A typical example for the spread distribution for an ensemble of 100 members computed with the STO-C method and is shown in Fig. 14c. In only about of the pixels does the spread exceed and the maximum values are limited to a few . While a bias of such a magnitude is relevant, as we have argued in the discussion about mean reflectance differences above, a random noise of this amplitude probably is not. In the assimilation of visible reflectances probably “superobbing,”—that is, averaging the reflectance over a certain radius, like the effective model resolution—will be applied and reduce the spread down to levels far below the assumed observation error. Therefore, the spread caused by the unknown subgrid cloud distribution seems not to be very relevant for DA. In principle, a deterministic method would therefore be sufficient. However, it still seems advantageous to use the stochastic STO-C3D method, as it is the most consistent approach and can result in reflectances that are significantly different from results obtained with the DET method. For both the mean spread and the RMS reflectance error obtained with STO-C3D are about , which should be sufficiently small to be negligible, compared to the assumed SEVIRI observation errors. Finally, it should be noted that the reflectance lookup process becomes a factor of *N* slower when *N* subcolumns are used. However, this process is still not the overall bottleneck for , because, for example, computation of the 3D fields of optical properties and retrieving data from the tilted columns still require more time.

## 6. Summary and conclusions

We have proposed efficient methods to increase the accuracy and consistency of a forward operator for visible satellite images based on a fast 1D RT method. Two problems were addressed: how to model the variation of reflectance related to the inclination of cloud-top surfaces and how to account for the overlap of subgrid clouds.

Using a transformation into a frame of reference aligned with the tilted cloud top, we reduce the 3D RT problem to a nearly one-dimensional problem that can be approximately solved with a fast 1D RT method. A comparison with 3D Monte Carlo RT calculations based on high-resolution model runs confirms that the inclination correction reduces errors, in particular for dense clouds. A comparison of synthetic images based on operational COSMO Model forecasts with 0.6-*μ*m SEVIRI observations show that the information content of synthetic images is increased significantly and that the reflectance histograms are improved when the cloud-top inclination is taken into account. In this case a tuning factor had to be introduced to compensate for the fact that the lower model resolution leads to gradients that are too weak. The additional information visible in synthetic images computed with the cloud-top inclination 3D correction could be a significant advantage for DA. Further studies are required to investigate how this information can be exploited.

The comparison with 3D Monte Carlo RT calculations also demonstrates that several other 3D effects still not accounted for can cause significant errors. Photon transport through the cloud sides can lead to a positive of negative reflectance bias and is especially problematic for clouds with intermediate optical depths. Moreover, cloud shadows were not taken into account, which leads to a positive reflectance bias. In particular for high clouds or larger SZA, a way to simulate cloud shadows could reduce errors and provide additional features that could be useful for DA. A preliminary version of a cloud shadow approximation was already used to generate synthetic satellite images for Heinze et al. (2017). For more accurate approaches, methods to estimate the diffuse radiation similar to the one of Wissmeier et al. (2013), who considered surface irradiances, could also be useful for top-of-atmosphere reflectances. The development of methods accounting for further 3D RT effects in an efficient way could be facilitated once fast 3D RT methods (see, e.g., Jakub and Mayer 2015) are applied within the NWP model and provide information about the diffuse radiation field.

Ignoring the problem of subgrid cloud overlap causes large errors in the cloudiness. Several different implementations of maximum-random cloud overlap yielded similar results and a much better agreement with the observed cloudiness. A 3D maximum-random approach taking the slant satellite viewing path into account is the most consistent method and leads to results that are similar to those of 2D schemes for the random overlap assumptions. The 3D scheme is still limited in the sense that it does not take cloud shadows into account. However, the results still support the conclusion that taking the slant viewing path into account can be of similar importance as changing the fundamental overlap assumption. The reflectance spread generated by stochastic methods is probably too small to be relevant for DA. As our goal was to maximize consistency with the RT code in COSMO, we did not consider newer suggestions for cloud overlap assumptions that are often based on a decorrelation length (see Hogan and Illingworth 2000) and vary with the seasons or parameters, like vertical wind shear (Di Giuseppe and Tompkins 2015) or the SZA (Tompkins and Di Giuseppe 2007). However, it should be relatively easy to adapt the 3D stochastic overlap method to more complicated assumptions. A second subgrid effect that could also be included in a future version of the stochastic scheme is the subgrid variation of water content, which is already taken into account in broadband radiation schemes (see, e.g., Pincus et al. 2003) and was found to be important by Shonk et al. (2010).

## Acknowledgments

The study was carried out in the Hans-Ertel Centre for Weather Research (Simmer et al. 2016). This German research network of universities, research institutes, and DWD is funded by the Federal Ministry of Transport and Digital Infrastructure (BMVI Grant DWD2014P8). Data generated in the research program High Definition Clouds and Precipitation for Advancing Climate Prediction [HD(CP)^{2}], founded by the Federal Ministry of Education and Research in Germany (BMBF), was used in this study.

### APPENDIX

#### Stochastic Cloud Overlap Algorithms

The stochastic methods STO-C, STO-N, and STO-R use for each model grid column (tilted toward the satellite) the profiles of subgrid cloud fraction , optical depths , , and effective particle radii , for water and ice clouds as input, where is the layer index with denoting the uppermost layer. From these input profiles, the methods generate *N* subcolumns with profiles , where *i* is the subcolumn index and *X* denotes any of the optical depths and particle radii. When layer *k* in subcolumn *i* is filled with a cloud, this is indicated by , otherwise . The column reflectance is then obtained by computing the reflectances for all subcolumns and averaging the results. Often several identical subcolumns exist, for which the reflectance is computed only once.

The stochastic methods must determine according to the maximum-random overlap rules such that averaged over a high number of realizations, the number of subcolumns in which layer *k* is filled, , converges toward so that on average the desired cloud fraction profile is reproduced.

In method STO-N [used already in Marquart and Mayer (2002) and corresponding to the stochastic maximum-random cloud generator described in Räisänen et al. (2004)] each subcolumn is filled independently starting from . For each layer a random number is generated and *f* is set according to

Here the fill probability depends on , , and , and it is defined by the following rules:

.

The first rule applies for empty layers, the second for the topmost layer of a group of cloudy layers, and the last two rules for cloudy layers directly below other cloudy layers. As can be shown easily, the expectation value of is indeed for each of these cases.

In methods STO-C and STO-R continuous clouds are generated, which can span several subcolumns, that is, for each layer the subcolumn indices and exist with

where the modulo operator is included to account for periodic boundary conditions. The indices are chosen such that the actual cloud fraction in each realization does not deviate by more than a discretization error from the desired value , and the cloud width , which is averaged over many realizations, has an expectation value of . For STO-R, which is based on the random overlap assumptions, is chosen randomly for each layer. For STO-C is chosen randomly only for the uppermost layer of each group of nonempty (i.e., ) layers. For further layers in these groups, is chosen randomly but under the constraint that the overlap with the layer directly above is maximized.

Finally, method STO-C3D is a generalization of STO-C, which takes the slant satellite viewing angle into account. For this purpose a bundle of subcolumns tilted by an angle *θ* toward the satellite (see Fig. 3) is filled with clouds in a similar way as in the method STO-C. For each layer *k*, the smallest integer number with is determined. An block of subgrid cells is filled with clouds. Starting from a randomly chosen corner of this block, or cells at the block edge are removed such that averaged over many realizations the desired cloud cover is reproduced. The cloud in the topmost layer and the clouds below empty layers are placed at a random position. Clouds below nonempty layers are in a first step placed randomly, but with the constraint that the overlap with the cloud in the layer above is maximized. In a second step their position is corrected by shifting them by an integer number of indices or in the *x* direction (i.e., toward the satellite) such that the expectation value of the shift is .

## REFERENCES

*Quart. J. Roy. Meteor. Soc.*, https://doi.org/10.1002/qj.3179, in press.

*ERCA 2008—From the Human Dimensions of Global Environmental Change to the Observation of the Earth from Space*, C. Boutron, Ed., Vol. 1, EDP Sciences, 75–99.

## Footnotes

© 2018 American Meteorological Society. For information regarding reuse of this content and general copyright information, consult the AMS Copyright Policy (www.ametsoc.org/PUBSReuseLicenses).