Surrogate cloud fields generated with the iterative amplitude adapted Fourier transform algorithm

A new method of generating two-dimensional and three-dimensional cloud fields is presented, which share several important statistical properties with real measured cloud fields.Well-known algorithms such as the Fourier method and the Bounded Cascade method generate fields with a specified Fourier spectrum. The new iterative method allows for the specification of both the power spectrum and the amplitude distribution of the parameter of interest, e.g. the liquid water content or liquid water path. As such, the method is well suited to generate cloud fields based on measured data, and it is able to generate broken cloud fields. Important applications of such cloud fields are e.g. closure studies. The algorithm can be supplied with additional spatial constraints which can reduce the number of measured cases needed for such studies. In this study the suitability of the algorithm for radiative questions is evaluated by comparing the radiative properties of cloud fields from cloud resolving models of cumulus and stratocumulus with their surrogate fields at nadir, and for a solar zenith angle of 0◦ and 60◦. The cumulus surrogate clouds ended up to be identical to the large eddy simulation (LES) clouds on which they are based, except for translations and reflections. The root mean square differences of the stratocumulus transmittance and reflectance fields are less than 0.03% of the radiative budget. The radiances and mean actinic fluxes fit better than 2%. These results demonstrate that these LES clouds are well described from a radiative point of view, using only a power spectrum together with an amplitude distribution.


Introduction
In order to investigate radiative transfer (RT) in the cloudy atmosphere, not only a RT model is needed, but also realistic twodimensional (2-D) or three-dimensional (3-D) cloud fields. In this paper, a new method is presented to generate cloud fields -based on temporally and spatially restricted measurementswhich can be used as input for RT models.
No current measuring technique exists which can provide instantaneous 3-D distributions of cloud properties such as the liquid water content (LWC) and number concentration. The best measurements currently available are 2-D, e.g. a horizontal field of radiances, or a time-height cross-cut of LWC. Therefore, the needed 3-D cloud fields are often simulated by cloud-resolving models (e.g. Scheirer and Macke, 2001). However, it would be dangerous to rely only on results based on model clouds alone; cloud-resolving models have to parameterize the cloud physics and chemistry, have a limited resolution and are often forced with very simple boundary conditions. Thus, a careful validation of the structure of model clouds is needed, and an independent method to model cloud fields directly from measurements is mandatory.
One such alternative are fractal cloud generators. The arguably simplest example would be cloud fields constructed with the Fourier method; these cloud have a power-law power spectrum and random phases. Kew (2003) developed a Fourier-based algorithm dedicated to ice clouds which includes fall streaks. Benassi et al. (2004) designed the tdMAP method to make fractal cumulus clouds. This method generates a tree, in which every height level represents a certain scale; the tree is pruned to generate cloud-free sections. Schertzer and Lovejoy developed methods to make multifractal clouds. Multifractals allow for a more accurate description of cloud structure than the monofractal Fourier spectrum; see Schertzer et al. (2002) for a clear introduction. A popular method is the Bounded Cascade algorithm that creates fractal cloud fields with a discontinuous structure (Cahalan, 1994). Di Giuseppe and Tompkins (2003) present an advanced algorithm using the Fourier transform, which has many features of fractal cloud generators, even if the chosen power spectrum is not a power law, but a gamma function.
The above fractal cloud generators (except for Di Giuseppe and Tompkins) have a structure that follows a perfect power law. Even though fractal mathematics describes several aspects of the observed cloud structure well, cloud fields also have features that are not fractal and are not scale-free. Common examples are gravity and lee waves in cloud fields, cloud streets and sometimes very regular hexagonal Bénard-type convective cells. Other, more rare examples are orographic clouds, contrails and ship tracks. Even Von Kármán vortices have been observed in cloud fields behind islands. Additionally, one can expect random deviations from scaling at large scales, as the fractal structure is a statistical property. At large scales, i.e. at scales close to the length of the measurement, the estimate of, e.g., the mean variance is inaccurate. Thus, deviations from the power-law scaling are to be expected. Most fractal cloud generators, however, do not take this into account and generate cloud fields that scale perfectly up to the largest scales.
Apart from fractal algorithms, many methods have been developed to generate cloud fields that are based directly on measurements. Examples are Scheirer and Schmidt (2004) and Los and Duynkerke (2000), who have generated stratocumulus clouds based on in situ LWC measurements and Liou et al. (2002) who have generated 3-D ice clouds based on ground-based lidar and radar measurements and on satellite observations. The structure of these cloud fields, which are based directly on measurements, has clear nonrealistic features. The stratocumulus clouds of Los and Duynkerke, for example, have a constant LWC perpendicular to the flight direction and the ice clouds of Liou et al. have a fixed cloud top and base height perpendicular to the wind. However, the advantage of this group of methods is that their cloud fields correspond closely to the measurements. This is important in closure experiments, where the clouds macro-, microphysical and optical properties are measured and modeled in an attempt to see whether they fit together, or there is still some physics missing. Moreover, such fields can be used to bridge the scale gap between various types of measurements and between measurements and models. This paper introduces surrogate cloud fields by combining ideas from fractal cloud fields with realistic structure and fields based directly on measurements. The method was designed to generate cloud fields for closure and comparison studies. Our method can also be applied as a user-friendly algorithm to generate fractal fields, but the cloud fields do not have to be selfsimilar. Thus, the method has the flexibility to generate cloud fields with a measured power spectrum. Moreover, these surrogate cloud fields also have the measured cloud water (height) distribution. A version of the algorithm interpolates scanning liquid water path (LWP) measurements without smoothing the field. In addition, a cloud mask can be included in the algorithm. Evans and Wiscombe (2004) developed an algorithm to make surrogate cumulus field, with a similar concept in mind, i.e. they also used measured power spectra and LWC height distributions. However, they describe the structure in terms of the power spectrum of a binary cloud mask; thus, the LWC structure inside the cloud is not taken into account. The Evans and Wiscombe algorithm can be seen as the specialist for the creation of surrogate cumulus fields from profiling measurements. By contrast, this paper will show that our iterative method is best used for (broken) stratiform clouds and is able to include non-statistical measurements, such as the measured values themselves (interpolating the rest of the field), and a cloud mask.
A more or less complete description of a cloud field requires at least one-and two-point statistics. One-point statistics relate to the values themselves: mean, standard deviation, etc.; the most complete description being the probability density function (PDF) or its empirical counterpart, the amplitude distribution. In this paper, the amplitudes denote either the LWC or the LWP values. The importance of one-point statistics in RT is exemplified by the success of the Independent Pixel Approximation (IPA; Cahalan, 1994). Thus, it is an important feature of our cloud generator that it reproduces this important cloud water distribution as precisely as possible.
Clouds made with the standard Fourier method, which uses complex Fourier coefficients with random phases, have a PDF that is on average close to Gaussian. Bounded Cascade clouds have a lognormal-like PDF. Such methods are, consequently, not able to replicate a measured distribution. Three-dimensional cloud fields will generally have a distribution with a significant number of zeros, as a result of the variability of the cloud boundaries, i.e. the distribution will neither be Gaussian nor lognormal like. Cumulus fields will have distributions with a significant number of zeros even in the 2-D case. Consequently, methods with a smooth PDF will only be able to generate 2-D stratiform cloud fields of LWP, or optical depth. As a workaround one can set all LWP values below a certain threshold to zero, but this will affect the power spectrum. Such a change does not necessarily represent a problem in theoretical work (Di Giuseppe and Tompkins, 2003), but is not desired in studies utilizing measured or fractal power spectra. The tdMAP algorithm and the multifractal generators are able to indirectly influence the PDF of their output cloud fields.
Two-point statistics describe how two points of a time series -or field -relate to each other, e.g. how strongly they are correlated. The structure of a stationary linear time series, i.e. a time series made by a linear dynamical process, is fully described by its autocorrelation function, or equivalently by its power spectrum. Cloud dynamics are known to be nonstationary and nonlinear, so more elaborate statistics may be required. This paper will show that the power spectrum and the amplitude distribution are not able to describe some nonlinear structures found in cloud fields. However, from the viewpoint of the radiative cloud properties, they provide a very accurate description of the cloud liquid water structure. The importance of combining one-and two-point statistics is also stressed in the paper of Evans and Wiscombe (2004).
The nonlinear dynamics community has developed methods to generate time series that have both a predefined amplitude distribution and a power spectrum. The most used recent method is the iterative amplitude adapted Fourier transform (IAAFT) algorithm by Schmitz (1996, 2000), which will be adapted in this study. Section 2 explains this algorithm in detail. The IAAFT method for time series will be generalized to 2-D and 3-D cloud fields in Section 3. The main objective of this paper in relation to the IAAFT method is to make 2-D LWP fields from one-dimensional (1-D) LWP measurements and 3-D LWC fields from 2-D LWC fields (i.e. LWC profiles as a function of time). In Section 4, 3-D LWC surrogate fields are made from 3-D cumulus and stratocumulus LWC fields generated by a Large Eddy Simulation (LES) model. Using RT calculations, we will show that the surrogate fields and the originals have very similar radiative properties.

The IAAFT algorithm
This section will describe the original IAAFT algorithm Schmitz, 1996, 2000) and a small alteration of this algorithm to make it more accurate. In a first reading of this paper, this section can be omitted. It is sufficient that the reader knows that given a power spectrum and a distribution, the iterative algorithm is able to generate a time series which has these statistical properties. In Section 3, the algorithm will be extended to multidimensional power spectra and fields.
The IAAFT algorithm was developed for nonlinearity testing (Theiler and Prichard, 1996;Schreiber and Schmitz, 2000). This is a test which is used in the nonlinear sciences to distinguish between time series originating from linear and nonlinear dynamical systems. As one is only interested in dynamical nonlinearities and not in static nonlinearities, the null-hypothesis of this statistical test is that the time series originates from a linear dynamical process which is modified by a nonlinear static filter.
The IAAFT algorithm is explained mathematically below and illustrated with a LWP time series in Fig. 1. On our webpage (see Section 6), the algorithm is clarified using pseudo Matlab code.

Spectral adaptation
The measured time series (Fig. 1, panel 1) is denoted by m n . The algorithm starts with a random shuffle of the data points (panel 2). In each iteration, i, the Fourier spectrum is adjusted first (panel 3) and then the amplitudes (panel 4). Based on the original time series the power spectrum is calculated as To produce a time series with the desired power spectrum {x 1, i }, the Fourier transform of the initial or the iterated time series {x 2, i −1 } is calculated (S k ) and the magnitudes of its coefficients are replaced by those of the original time series (M 2 k ). The phases, ϕ k = S k /|S k |, remain unaltered. Thus, the complex Fourier coefficients of {x 1 } are given by Consequently, the time series in panels 1 and 3 have the same power spectrum; the difference between their structures is due to differences in their distributions.

Amplitude adaptation
In the second step, the amplitudes are adjusted based on their ranking.
As this amplitude adaptation will alter the power spectrum, both steps are repeated until a convergence threshold is reached and the surrogate time series is returned (panel 5).

Convergence
As a measure of the accuracy ( ) of the IAAFT algorithm, we use the average of the changes that must be made to the iterated time series at each step, relative to the total standard deviation (σ ) of the time series: The final accuracy of the surrogates depends on the number of values, i.e. for longer time series, a better fit can be achieved (Schreiber and Schmitz, 1996). The algorithm typically converges to a positive non-zero accuracy. Repeating the algorithm with the same statistical parameters, it will converge to another accuracy factor. This illustrates that the algorithm finds local minima. Consequently, the accuracy can be improved by repeating the calculation and choosing the most accurate surrogate. Typical calculation times on a normal PC for the IAAFT method range from a few seconds for a time series of 4096 samples to about an hour for large 3-D fields with 5 × 10 6 data points. For the final surrogate, one can select the time series either after the last Fourier adaptation (x 1 , with an identical power spectrum), or after the last amplitude adaptation (x 2 , with the same values as the original). We have selected the latter, in light of the fact that small errors in the LWC values are probably more important than small errors in its structure, especially for LWC fields with values close to zero.
In the engineering community the IAAFT algorithm was recently discovered independently and compared to their older methods (Masters and Gurley, 2003). The IAAFT algorithm was found to be most accurate for the cases studied. Lewis and Austin (2002) used an algorithm similar to one of those older algorithms to create fractal clouds with a lognormal distribution.

Stochastic algorithm
Using the normal IAAFT algorithm the accuracy is not sufficient for cumulus fields. This is apparent in cumulus fields as they are very sensitive to small deviations; here one observes single pixel clouds with low water content in the cloud-free portions, which have a relatively strong impact on the radiative properties of the cloud field.
We developed a new, more accurate IAAFT algorithm by slightly altering the normal scheme, namely, by starting the algorithm with a stochastic stage. In this first stage, only one-fifth of the values in the amplitude adaptation were substituted; a sen-sitivity study showed that any value between a few percent and 40% performs well. As this is done in every iteration, in the end, the amplitude distribution of the surrogate is still very close to the measured time series. However, the surrogate does not have exactly the same amplitude distribution. Consequently, a second stage is applied using the standard IAAFT algorithm. This second stage starts with the surrogate from the first stage as first guess instead of white noise. We call this algorithm the stochastic iterative amplitude adapted Fourier transform (SIAAFT) algorithm and it is used in Section 4.2 to generate cumulus surrogates from LES cumulus fields.

Surrogate cloud fields of measurements
In this section, we will use the IAAFT method to make cloud fields with power spectra and amplitude distributions from measured clouds. First, we will make surrogates of the same dimension as their original measurements in Sections 3.1 and 3.2, to illustrate the strength and the limitations of the approach. Then, in Sections 3.3 and 3.4 we will demonstrate how to generate a 2-D surrogate field from a 1-D measurement, or a 3-D field from a 2-D measurement. In this case, the higher dimensional power spectrum is calculated from its measured lower-dimensional power spectrum by assuming that the field is horizontally isotropic. Section 3.5 shows that for scanning measurements one can estimate the anisotropic power spectrum and that the algorithm can also be used for interpolation. The typical problems of spectral analysis also apply to the IAAFT method. Hence, use of a sufficient number of data points and appropriate windows before application of the Fourier transforms is important. Windows are needed in Fourier analysis as the time series are assumed to be periodic, i.e. it is assumed that after the last value of the measured time series, the time series resumes with the first measured value. Consequently, if the first and the last value of a measured time series are very different, spurious high frequencies are introduced that correspond to this discontinuity. However, since windows also have their disadvantages, we have opted not to use windows, but to carefully select our time series and fields to have similar values at the beginning and end.

1-D LWP surrogate
As a first example, we utilized a 1-Hz LWP time series measured with the MICCY 22-channel microwave radiometer (Crewell et al., 2001) on September 5, 2001 at Cabauw, The Netherlands, during the BALTEX Bridge Campaign (BBC). The observed stratocumulus cloud had a cloud base at 2 km, with some wisps of cloud at 1 km, and the high LWP periods in the time series are associated with virga from the cloud. The wind speed at 2 km at the beginning of the measured time series was 5.5 m s -1 , as estimated from radiosonde data. We used the wind speed to convert the temporal scale to a spatial scale. There were some gaps in the data due to automatic calibration of the microwave radiometer; these were bridged by linear interpolation.
From this LWP time series (Fig. 2a) we created a 1-D surrogate ( Fig. 2b) that exhibits a similar structure. In the LWP measurement there are peaks, which are associated with the fall streaks, in the high LWP part of the time series, but the low LWC section has much less variability. The surrogate has much variability in the low LWC section of the time series. Visually comparing the surrogate with the measured time series is important, as the human eye is good at recognizing patterns and spotting differences. Figure 2c and 2d show that the amplitude distributions and power spectra of the measurement and the surrogate are (almost) identical. In this case, we achieved an accuracy of 0.3% with 4096 data points after 3 s of calculation on a desktop Pentium IV PC.

2-D LWC profiles
In this section, we create 2-D vertical LWC surrogates from 2-D LWC profile retrievals. This illustrates the quality of the surrogates and is a first step towards making 3-D LWC fields from such 2-D LWC profiles, which is discussed in Section 3.4. To generate these surrogates, we used a LWC profile time series, i.e. a 2-D time-height field (Fig. 3a). These LWC profiles were derived using an optimal estimation technique ) that combines microwave radiometer brightness temperatures, radar data, radiosonde profiles, lidar ceilometer cloud-base heights, and 2 m temperature and humidity measurements, with a priori information from a microphysical cloud model. Because of gaps in the data, we had to linearly interpolate the field to a regular grid, resulting in some loss of variance. The winds were strong with speeds reaching 9.7 m s -1 at cloud height as measured by the wind profiler at 8.5 h UTC. In this paper, all references to times are expressed in decimal hours. Since the temporal resolution of this retrieval was 10 s, the horizontal spatial resolution of the LWC profiles is 97 m. Apparently this poor spatial resolution makes the cloud LWC field appear somewhat noisy. However, its power spectrum shows good scaling behavior at all scales, indicating that measurement and algorithmic noise, which tend to be white, are not observed.
The LWC field is clearly anisotropic and we would like the surrogate to have the same height-dependent LWC amplitude distribution. In this case, we did not use one sorted vector with all of the LWC values, as with the basic algorithm, but rather, we performed this operation separately for every height level; that is, we utilize a 2-D sorted LWC field. With this approach, the cloud top and base height distributions are automatically preserved.
For the LWC profiles, one typically observes trends, e.g. an adiabatically increasing LWC, that lead to spurious high frequencies in the power spectra. To prevent this problem, we subtract the mean LWC profile from the 2-D field itself before any calculations are made and then add this mean profile to the surrogate at the end. The resulting structure is thus defined globally, whereas the amplitudes are defined per height level.
The structure of the 2-D surrogate (Fig. 3b) looks similar to the original LWC profiles (Fig. 3a), but there are some differences. Especially strong vertical structures, like the fall streaks in and below the cloud, are not as striking in the surrogate as they are in the measurement.

2-D LWP surrogate fields
In the next two sections, we will illustrate an important application of the algorithm, the ability to create cloud fields that have one dimension more than the measurement. For isotropic fields the power spectrum can be converted from any dimension to another, see e.g. Christakos (1992, chap. 6). For fractal spectra Austin et al. (1994) gives this relation for 1-D to 2-D and Kew (2003) for 1-D to 3-D. Evans and Wiscombe (2004) use an optimization method to find a 3-D horizontally isotropic power spectrum that has the same 2-D power spectrum as the original.
For our purposes we will restrict the discussion to the easy case that the field is isotropic and the distance equal in both directions. Then we can create a 2-D power spectrum from a 1-D one by rotating, as well as rescaling the 1-D Fourier coefficients around the origin in 2-D k-space (Fig. 4a). The rotation makes the coefficients independent of direction. The rescaling is necessary so that the power integrated over a circle for a certain frequency is equal to the 1-D power at that frequency. For a n-dimensional hyperphere the equation reads (Christakos, 1992)  The temporal scale of the time series was converted into this spatial scale using a wind speed of 5.5 m s -1 , taken from a nearby radiosonde. The high peaks in the LWP values are due to the fall streaks in the virga or very light rain; thus, they would hardly be visible in a photo or satellite measurement of this cloud field. The accuracy was at 5‰. The computation required about 1.5 h of CPU time for this 2-D field with about 2 × 10 7 data points.
with A n being the integration area, θ a unit vector, and the inner product k·θ the distance from the origin. The subscript n denotes the n-dimensional power spectrum that is calculated from the 1-D power spectrum S 2 1 . As the size of a circle is proportional to the frequency, this means one has to divide all coefficients by their wave number. As a result, the squared Fourier coefficients of the isotropic 2-D LWP field are given by S 2 (k x , k y ) ∝ S 2 (k)/k with k = k 2 x + k 2 y being the magnitude of the 2-D wavenumber. In the algorithm, k is rounded to the nearest wavenumber of the 1-D measurement. After this step the total variance of the power spectrum is scaled to the total variance of the amplitude distribution.
If k x and k y both have the value of the maximum wavenumber of the original time series, the maximum isotropic wavenumber is √ 2 times the maximum 1-D wavenumber. This means the rotation will not produce spectral values to high k x and k y , i.e. the grey upper right area in Fig. 4a. We chose to fill this area with the power at the Nyquist frequency, i.e. we assumed a flat spectrum at these frequencies. A more elegant method would be to extrapolate the power spectrum to higher frequencies using its self-similarity. However, we expect the results of either approach to be similar.
In Fig. 2d, one can see that the average 1-D power spectrum calculated from the rows and columns of the 2-D surrogate field is almost the same as the power spectrum of the 1-D measurement, only the typical Fourier noise is missing.
If the measured LWP vector has N values, the 2-D LWP field will have N × N values. Hence, the vector with the sorted LWP values, used in the amplitude adaptation, needs N × N entries; every measured value has N copies in this vector. After these extra preparations, the basic algorithm from Section 2 was followed again.
One of the resulting surrogates of this calculation is shown in Fig. 4b. This 2-D surrogate was calculated using the 1-D time series presented in Section 3.1. To make the calculation, we had to average the data to 2-s resolution, due to memory limitations. The LWP field in Fig. 4b may look noisier than the clouds we see in everyday experience. However, the eye does not perceive a LWP field, but rather, a radiance field. This radiance field is smoother than the LWP field at small scales. Because of light scattering, the radiance is representative for a region, not just a point. This is called radiative smoothing and reduces much of the variability in the radiance field at scales smaller than the cloud depth. Furthermore, fall streaks contribute significantly to the variance of the LWP field, but are not important for visible radiation.
We constructed the 2-D surrogate cloud field using a 1-D time series, which had a much smaller number of data points. From standard stationary statistics, we have learned that this is acceptable. The 1-D time series represents a noisy estimator of the larger 2-D field. However, cloud time series are not stationary (Davis et al., 1996(Davis et al., , 1999. The mean, width and all other moments and statistical properties of a time series may not be defined for a nonstationary time series or inhomogeneous fields. The fact that most power, i.e. variability in LWP is at the large scales means that the amplitude distribution will never converge. A larger sample from this cloud field would most likely have larger LWP variance and a broader LWP distribution, while a small sample from a non-stationary data set will generally have a smaller width of the amplitude distribution. As an illustration, consider the standard deviation of the surrogate LWP cloud field illustrated in Fig. 4b. The standard Tellus 58A (2006), 1 deviation of the entire 2-D LWP field is 109 g m −2 , while the average standard deviation of 1-D line measurements of this field is only 90 g m −2 . Thus, if we had collected a zenith pointed measurement from this spatially correlated cloud field, we would have found a 17% smaller width of the LWP distribution. This can be understood by imaging a 1-D measurement through the maximum in the field. Due to the strong spatial correlations, the other values will also be higher than average and the standard deviation consequently smaller for such a 1-D measurement.
The broadening of the amplitude distribution to larger values of cloud liquid water is not that important for RT as the relations between the radiative fluxes and LWP levels saturated at high LWP. Thus the radiance is insensitive to errors in LWP in this region. However, the broadening of the amplitude distribution to lower values has the potential to be important and should be carefully studied to estimate the size of possible biases.

3-D LWC surrogate fields
To generate a 3-D surrogate field based on the 2-D LWC profiles, we apply the methods developed in Sections 3.2 and 3.3. Again we remove the mean LWC profile from the measurement and use a 2-D height-dependent sorted LWC field in the amplitude adaptation. The 2-D time-height-based Fourier coefficients are subsequently cylindrically rotated around the vertical wavenumber axis and scaled like in Section 3.2 to make a 3-D horizontally isotropic power spectrum.
A 3-D surrogate field made from the 2-D LWC retrieval presented in Section 3.3 is illustrated in Fig. 5. The three 2-D fields actually shown in Fig. 5 are the average values of the 3-D field in one of the three spatial dimensions. Note that the horizontally integrated fields look much smoother than the measured 2-D LWC field, as an average field should. The cross-sections do have a structure similar to the measurement (Fig. 3c).  Fig. 3c. The calculation of this 3-D surrogate with 5 × 10 6 data points took about an hour on a desktop PC and its accuracy was 5 × 10 −4 .

2-D LWC field with local forcing
In the future, we would like to apply the IAAFT method to create surrogate cloud fields based on ground-based scanning microwave radiometer measurements. A typical scanning measurement would be to rotate the instrument on its azimuth-axis, at a fixed elevation angle, while the clouds drift by on the wind. In this way one would measure LWP values in a spiral pattern at cloud height, assuming that the cloud thickness is less than the cloud height.
For the example in this subsection, we used an enhanced version of the algorithm that also takes the LWP values measured on the spirals as its third input. In an extra iterative step, between the spectral adaptation and the amplitude adaptation, the algorithm forces the LWP values at the measured locations, to correspond with the measurement.
One advantage of a scanning measurement is that the 2-D LWP power spectrum can be estimated without assuming isotropy. From the LWP values on the spiral, the 2-D periodic autocorrelation function is estimated. With the Fourier transform the autocorrelation function is converted into the 2-D anisotropic power spectrum. The LWP distribution is estimated based on the LWP values on the spirals, i.e. if the full field has x times more values, the LWP distribution contains every measured LWP value x times.
In this subsection, we will simulate an idealized scanning measurement on a LES maritime stratocumulus field, which was simulated by Duynkerke et al. (2004). These LES clouds were modeled for the conditions during the FIRE I campaign with a resolution of 50 m in the horizontal (52 grid boxes) and 10 m in the vertical (112 grid boxes). Based on this LWC field, the LWP field was calculated. The scanning pattern is plotted on top of the LWP field (see Fig. 6a). The measured values, which correspond to 16.5% of all values, were used to calculate the surrogate shown in Fig. 6b. Fig. 6. A simulation of a LWP measurement on a LES stratocumulus field. The statistics and the local LWP values of the LES LWP field (Fig. 6a) are used to make the surrogate LWP field (Fig. 6b).
This version of the algorithm can be seen as an interpolation algorithm that does not smooth the field, but maintains the structure of the measurement. In future work, we intend to investigate the influence of the fraction of points and the scanning pattern on the accuracy of the cloud field. It might be possible to find a more optimal method to estimate the PDF where the local density of the number of data points is taken into account. Furthermore, we would like to apply it to measured data (LWP and LWC profiles) and compare the results to measured radiative fluxes.

Evaluation
In this section, we assessed the appropriateness of using power spectra (linear spatial correlations) and LWC height distributions to create surrogate cloud fields for closure or comparison studies. Some examples in Section 3 illustrated that these statistics do not capture all aspects of the structure of clouds. We do not know a priori whether the missing nonlinear structures are importnt for RT or not. This question can only be answered for a specific application. In this paper, we aim to study it for the distribution of the radiances, irradiances and actinic fluxes. The irradiances are important for the radiative budget (i.e. climate and weather models), while the radiances are important for remote sensing purposes. Actinic fluxes, the radiance falling on a point integrated over all spatial angles, are important for photochemistry.
To assess the quality of the surrogate clouds, we utilized two sets of LES clouds, comprised 31 stratocumulus LWC fields and 52 cumulus LWC fields and created surrogate clouds from each LES cloud. The difference in the radiative properties between the LES original and its surrogate is then used as a measure of the quality of the surrogates. Next to the normal IAAFT surrogates, we made Fourier surrogates and PDF surrogates. The Fourier surrogates are calculated by using the power spectrum of the original together with random Fourier phases. As always the mean height profile is added afterwards. These Fourier surrogates have, on average, normally distributed LWC values. The PDF surrogates are created by randomizing the LWC values at each height level. These PDF surrogates will thus have exactly the same LWC height distribution, but no spatial correlations.
One of the maritime stratocumulus clouds (Duynkerke et al., 2004) used is shown in Figure 6. Drop sizes were calculated for each cloud grid box by assuming a monodisperse distribution (i.e. 1 drop size) with 300 drops per cm 3 . The average reflectance is 0.66 and the average optical depth is 47. The clouds are without gaps and their depth varies from 30 to 460 m and is on average 300 m. These LES clouds were modeled for the conditions during the FIRE I campaign with a resolution of 50 m in the horizontal (52 grid boxes) and 10 m in the vertical (112 grid boxes).
The cumulus cases were originally modeled for the ARM program and represent the diurnal cycle of cumulus over land (Brown et al., 2002). Until the middle of the morning there are no clouds. The cloud cover grows to a 31% (16% cloud reflectance) in the afternoon. After this the number and cloud cover decreases, but the clouds become deeper (almost 5 km) in the middle of the evening. The run stops at midnight. The clouds have a resolution of 100 m in the horizontal and 112 m in the vertical. The number of grid boxes is 66 by 66 horizontally with 122 height levels. The average reflectance is 8% (assuming 1000 drops per cm 3 ).
The radiation calculations for all sets of clouds have been made assuming the wavelength of the incoming monochromatic solar radiation to be 550 nm with a solar zenith angle (SZA) of 60 • . Additionally the radiative properties of the IAAFT surrogates where also tested for a SZA of 0 • . The surface albedo was set to zero, periodic boundary conditions were applied and the number of photons amounted to 10 7 . The cloud field was assumed to be in a vacuum.
The upward and downward irradiances were calculated using the Leipzig Monte Carlo Model (LMCM; García and Trautmann, 2003). The actinic fluxes were calculated with SHDOM (Evans, 1998) using N µ = 16 and N ϕ = 32, a splitting accuracy of 0.03 and a solar flux of 1.83 W m −2 nm −1 . For easier comparison with the other 2-D fields, the 3-D actinic flux field has been vertically integrated to a 2-D integrated actinic flux (IAF) field.

Stratocumulus
For each of the LES cloud fields, we created two surrogate fields. An example is shown in Figs. 7a and 7b. The surrogate stratocumulus fields showed poor convergence with an average accuracy of 12%, with the worst value at 18% (see Section 2). Two surrogates from the same original can differ greatly in accuracy; the algorithm finds a local minimum.
Stratocumulus and altocumulus often display beautiful cell structures, similar to Bénard convection. The input LES clouds show such features. Their 3-D IAAFT surrogates show these much less. For example the LES cloud in Fig. 7a has clearer line patterns than its IAAFT surrogate in Fig. 7b. The same phenomenon can be seen in surrogates made from 2-D LWP fields and 2-D radiance fields with Bénard cells. It was even clearer in a calculation made with one high-resolution LES stratocumulus which had about 2 × 10 7 data points (Schroeter and Raasch, 2002).
To validate the radiative properties of the surrogate clouds, we compared the radiative fields of all pairs of LES clouds and their corresponding surrogates; see Fig. 8 for one example of an LES cloud and its three different surrogates. The main visible difference between the Fourier and the IAAFT surrogate is that the latter captures the extreme values better. The PDF surrogates are very similar to homogeneous clouds and show little variance.
With an analysis of variance (ANOVA) we compared the differences between two sets with RT results from the Monte Carlo models at 60 • SZA. The first set is the result for the LES cloud and its surrogate, the second set of results is from same LES cloud which is calculated twice. The deviations in the mean optical properties of the IAAFT surrogates are not statistically significant. The Fourier surrogates and the PDF surrogates did show significant deviations in their mean optical fields. Thus, in this case the remaining differences can be attributed to the differing structure of these surrogates. For the actinic fluxes, calculated with SHDOM, a root mean square (RMS) error cannot be determined.
the relative RMS error of all radiative properties of the IAAFT surrogate stratocumuli are less than 1.4% compared to the LES clouds they are based on. In case of the mean optical depth only differences for the Fourier surrogate are possible. For the irradiance and the radiance the Fourier surrogates deviate most, as here the right PDF is most important. The too low transmittance for the PDF surrogates is analogous to the bias found in plane parallel homogeneous clouds. The actinic flux is too high for the Fourier surrogates and too low for the PDF surrogates.
To compare the full distributions of the optical properties, the values of the original fields and their surrogates were sorted and plotted against each other, see Fig. 10. This results in one line for every set of originals and surrogates, which should ideally lie on the dotted line that indicates y = x. If the lines run flatter this means that the surrogates have less variance in this range of values. Random shades of grey have been used to make it easier to follow individual lines. To further compare the distributions we have also made scatterplots of the standard deviation, and the minimum of the fields. These are not shown here, but are listed in Table 1.
The distribution of the optical depth is best captured by the IAAFT surrogates (Fig. 10a). The standard deviation of the PDF surrogates has a bias of 71% and consequently also the extreme values are off. The bias of the standard deviation of the Fourier and the IAAFT surrogates is only 2 and 1%, respectively. The distribution of the PDF and the Fourier surrogates is much more symmetric as those of the originals or the IAAFT surrogates.
The deviations of the distribution of the irradiances are very similar to those of the optical depth. However, remarkable is that the standard deviations of the reflectance are too small, showing biases for the Fourier and the PDF surrogates of 13 and 37%, respectively.
The standard deviations of the 2-D actinic flux fields of the Fourier and the PDF surrogates are too low; their biases are 25 and 62%. In the actinic flux fields of the original LES cloud (and the IAAFT surrogates) the maximum flux is always at the top. However, in about half of the PDF and Fourier surrogates the maximum actinic flux is not at the highest level, but one level (10 m) lower. This is likely due to the difference in the cloud top structure. Here it should be noted that SHDOM linearly interpolates between the grid points, and thus reduces the number of clear boxes in the PDF surrogates.
For all optical properties the power spectra of the PDF surrogates where flat. The spectra of the IAAFT and the Fourier surrogates followed that of the originals, except for a scaling factor in case of differences in the standard deviation of the fields.

Cumulus
For the cumulus clouds, the accuracy turned out to be very important, but the average accuracy of the standard IAAFT algorithm Tellus 58A (2006), 1 was only 7%. These fields still have too many small clouds that have not merged to form one or a few larger ones. As a result, a bias towards a larger reflectance and albedo occurs in the surrogates. The radiative properties of the best converged surrogates of each pair of cumulus surrogates correspond closer to those of the originals than those of the worst converged surrogates. By repeating the calculation and choosing the best converged surrogates, we were able to reach an accuracy of 4‰. However, even these surrogate cumulus clouds still showed a bias.
Therefore, we applied the Stochastic Amplitude Adapted Fourier Transform algorithm (see Section 2). This algorithm provided a much better convergence. In the first run, 39 of the 52 surrogates converged fully, by repeating the calculation all surrogates could be made to converge fully. Surprisingly, the surrogate cumulus clouds are in the end identical to their underlying LES clouds, except for translations or reflections, see Fig. 7 for an example. Thus, the radiative properties of the LES cumulus clouds and their surrogates will be the same as well. In this case, the amplitude distribution and the power spectrum are able to completely describe the cloud structure.
That the surrogates are identical, may be understood by recognizing that the cumulus fields are almost binary fields; that is, two thirds of the fields' standard deviation is determined by the cloud 'mask'. Specifically, if we set the LWC of the cloudy boxes (i.e. LWC > 0.1 g m −3 ) to the average LWC of all the cloudy boxes, the resulting binary field has, on average, 68% of the variability of the original field. By calculating the autocorrelation function for all permutations of a few short sparse correlated binary time series, we found that the autocorrelation function is able to constrain the structure to a great degree. We expect that the cumulus surrogates will not be identical to their originals for cumulus clouds with a higher cloud cover or for cumulus fields with a higher resolution.
We did not create Fourier or PDF surrogates from these sparse cumulus fields, as they would be completely inappropriate in this case. Both types of surrogates would have been much more homogeneous and stratiform as the LES cumulus fields.

Discussion
This section will discuss the limitations and strengths of the IAAFT surrogate clouds based on the examples from Section 3 and the quantitative evaluation in Section 4. This then leads to a discussion on the applications that could benefit from these fields and further work that would be needed.
The IAAFT method does have limitations. The main limitation is that the nonlinearities of its surrogates are static; they do not originate from nonlinear dynamical equations. This is probably why the algorithm is not fully able to handle fall streaks, Bénard cells and other line-like structures. On the other hand the inclusion of the PDF makes the description with a Fourier spectrum much more powerful, as can be seen by the surrogate fields of the sparse cumulus fields.
A theoretical limitation we did not show is that the IAAFT algorithm cannot generate time series with asymmetric dynamics, e.g. time series with slow rises and fast drops. Formulated more precisely, on average the dynamics of the IAAFT-generated surrogate time series will be symmetric. This is a natural consequence of the fact that the overlying structure of the time series is described by an autocorrelation function which contains only information on the strength of the correlation between two points as a function of their distance.
These limitations appear to pose no serious problem for issues of RT. The analysis in Section 4 showed that the radiative properties of the IAAFT surrogate clouds are very close to those of their original LES clouds. This result, if found to be generally applicable, greatly simplifies the measurement and analysis of cloud structure. It is expected that the improvement of the IAAFT surrogates over the Fourier surrogates and the PDF surrogates is even larger for thin or broken stratocumulus clouds. A remarkable result is that the surrogate stratocumulus cloud fields that lose their nonlinear Bénard cells still match the original clouds closely with respect to their radiative properties.
In order to execute a closure study, the uncertainties due to limited data and its measurement errors should be closely examined. In our evaluation we had error-free and complete information on the 3-D structure of the original cloud fields. Even the most comprehensive cloud measurement cannot provide this quantity and quality. For closure studies, comparing physical and radiative properties of a cloud field, the limiting factor is thus expected to be the amount and the quality of the data, not the algorithm or the statistical description of the cloud field.
In the evaluation of the sparse cumulus fields, we have seen that the convergence needed to be very high to get good results for RT. This also means that the initial estimate of the power spectrum should be very accurate. This accuracy cannot be reached in a real measurement and as a result a cumulus surrogate based on PDF and spectrum only, will show biases in its albedo. For closure studies this problem can be solved by including a cloud mask from an imager. Preliminary studies with a ver-sion of the algorithm where every iterative step a cloud mask was applied, showed good results and smooth transitions at the cloud boundaries.
As the Fourier phases are kept intact during the spectral adaptation, the algorithm changes the structure of the fields only little. This allows one to impose the PDF on the field, but also for spatial constraints, like a cloud mask or the local scanning values (see Section 3.5). Other spatial constraints we are working on are coarse grained mean LWP or coarse grained cloud fraction which could be derived from satellite measurements. Likely, fields of cloud top or base height could also be included.
To test the general applicability of the conjecture that our statistical description is sufficient for RT, a more thorough analysis would be necessary. More cloud types should be studied. It should be noted that the structure of the LES stratocumulus fields used in the analysis did not have much variance in LWP and were very well behaved compared to real cloud fields. For instance, there were no multiple cloud layers and the cloud boundaries were relatively smooth. This well-behaved nature is due to the idealized atmospheric profiles and the periodic boundary conditions used in LES modeling. Accordingly, using cloud fields from other types of cloud-resolving models and multifractal or wavelet-based cloud generators would provide an additional test. Still, this study using surrogates based on LES clouds gives confidence that these further tests will give positive results.
One could argue that for developing retrievals and parameterizations, it is undesirable that all surrogate cloud fields from one set of statistics have the same radiative properties. Ideally, such retrievals and parameterizations should be based on as many realizations as possible. However, for closure studies, this property is a clear vice, and for the development of retrievals and parameterizations one can each time use different measurements to base the surrogates upon. We have, for example, routinely made surrogate cloud fields from the 2 months of measurements of the BBC campaign .
Compared to fractal cloud generators, the main advantage of the IAAFT algorithm is that its clouds are based on measured statistics and not on fractal idealizations of the measured structure. The multifractal (Schertzer et al., 2002) and wavelet-based (tdMAP) generators have the advantage that they have a more comprehensive description of the structure, but the evaluation in Section 4 of the IAAFT surrogate clouds suggested that for many applications our statistical description is sufficiently accurate. A study that would show that or when this holds would help our fundamental understanding of 3-D RT.
However, we would rather not compare our method to fractal cloud generators, but rather to algorithms that create clouds based on measurements. Compared to these methods the advantage of the IAAFT surrogate clouds is that its structure is more realistic. Only the algorithm of Evans and Wiscombe (2004) has a similar degree of realism. This algorithm describes the cloud structure in terms of a binary cloud mask in stead of LWC. Thus, the LWC at the IAAFT cloud boundaries will be smoother.
Furthermore, due to the iterative nature and the phase preserving spectral adaptation, the IAAFT algorithm is more flexible and can more easily integrate spatial constraints by additional measurements.

Conclusions and outlook
The IAAFT (Iterative Amplitude Adapted Fourier Transform) algorithm is able to generate realistic cloud fields for use in RT calculations. These cloud fields are described by a power spectrum and by a height-dependent LWC distribution. As these synthetic or 'surrogate' clouds can be based on the full statistics from observations, they are ideally suited for empirical studies. An important application of the surrogate clouds will be closure studies. Such studies will especially benefit from future versions of the algorithm that can include various spatial constraints.
Although the surrogate clouds do not capture the entire nonlinear structure of clouds, the evaluation using LES clouds showed that their radiative properties are (almost) equal to the LES clouds. For this evaluation gentle homogeneous stratocumulus and sparse cumulus fields were used.
These sparse cumulus fields were identical to their LES original, i.e. they were completely described by the two statistical parameters of the IAAFT algorithm. Thus, the cumulus surrogates had exactly the same radiative properties as the LES cumulus clouds. We expect that this result is not valid for cumulus clouds with a higher cloud cover or made from a LES cloud with a higher resolution. In practice the power spectrum can not be measured with sufficient accuracy, real surrogates will thus not be identical, but rather biases in cumulus surrogates are to be expected. As these biases are mainly due to wisps of cloud in the cloud free sections, these biases may be reduced by introducing a measured cloud mask.
The radiative budget of stratocumulus clouds is also matched well by the surrogates. The RMS differences of the transmittance and reflectance are less than 0.03% of the budget. The radiances and actinic fluxes fit better than 2%. These numbers should be seen as an order-of-magnitude estimation, as they will be different for other cloud fields. The results for the radiances and the irradiances were calculated with a Monte Carlo (MC) model and the differences can be attributed to the stochastic error of the RT models. Simpler PDF surrogate clouds, which had no spatial correlations, had much larger differences and these differences could no longer be attributed to the MC models, showing that structure is important. Fourier surrogate clouds, with only the same spatial correlations, also had much larger differences as the IAAFT surrogates, showing that an accurate LWC distribution is important. These results are only valid for the relatively homogeneous LES stratocumulus fields. Suggestions are given in the discussion to study the generality of this conjecture.
In the future, we aim to focus on utilizing spatial constraints, expanding the work in Section 3.5 on iterating local LWP values. For example, one could combine ground-based cloud structure measurements with horizontal 2-D cloud-top height or cloud mask field as measured by an airplane or satellite.
We have performed 10-Hz scanning measurements with the microwave radiometer MICCY and the cloud radar MIRACLE during the BBC and BBC2 campaigns. The high temporal resolution and the spatial sampling of these scanning measurements will allow better estimates of the LWP distributions of inhomogeneous cloud fields than standard zenith measurements and are ideal for closure studies using surrogate clouds. Schreiber and Schmitz (2000) have shown that it is possible to iterate two time series simultaneously (i.e. breath rate and heart rate), while also taking into account their cross power spectrum. So one could, for example, generate surrogate cloud fields using LWC and effective radius statistics simultaneously. Alternatively, one could utilize cloud LWP and cloud-top height statistics simultaneously. However, cloud-boundary structure is much more difficult to describe than cloud-liquid water. Even the relatively homogeneous altocumulus cloud in this study, has two cloud layers in some parts and many clouds have gaps, i.e. no cloud-top height. Many clouds will therefore not have a useable time series of cloud-top height, and in many cases we will have to use fractal geometry to describe their structure. Unfortunately, algorithms to calculate fractal measures do not have inverse functions, like the Fourier transform; hence, iterative algorithms to generate surrogate clouds cannot be constructed.
This problem can be solved by using a global search algorithm, which also promises to converge better. Schreiber and Schmitz (2000) use simulated annealing to search for time series with almost arbitrary statistical properties. We are working on an evolutionary search algorithm to find cloud fields that have specified properties with respect to cloud-liquid water, cloud base and top height and number of layers. This method is slow, but flexible and may be able to handle (multi-)fractal, non-linear and non-stationary cloud structure measures.
The code of the IAAFT algorithms and many more examples are available on our website: http://www.meteo.unibonn.de/victor/themes/surrogates/.

Acknowledgments
We would like to express our thanks to Klaus Lehnertz from Department of Epileptology of the University of Bonn for introducing us to surrogate time series, Michael Schröter and Siegfried Raasch from the Institute of Meteorology and Climatology of the University of Hannover for the high-resolution LES stratocumulus field. This research was carried out in the framework of the 4-D-clouds project, which is sponsored by the German Ministry of Research, BMBF, in the AFO2000 Program on Atmospheric Research.