## Abstract

Results are presented of two large-eddy simulation (LES) runs of the entire year 2012 centered at the Cabauw observational supersite in the Netherlands. The LES is coupled to a regional weather model that provides the large-scale information. The simulations provide three-dimensional continuous time series of LES-generated turbulence and clouds, which can be compared in detail to the extensive observational dataset of Cabauw. The LES dataset is available from the authors on request.

This type of LES setup has a number of advantages. First, it can provide a more statistical approach to the study of turbulent atmospheric flow than the more common case studies, since a diverse but representative set of conditions is covered, including numerous transitions. This has advantages in the design and evaluation of parameterizations. Second, the setup can provide valuable information on the quality of the LES model when applied to such a wide range of conditions. Last, it also provides the possibility to emulate observation techniques. This might help detect limitations and potential problems of a variety of measurement techniques.

The LES runs are validated through a comparison with observations from the observational supersite and with results from the “parent” large-scale model. The long time series that are generated, in combination with information on the spatial structure, provide a novel opportunity to study time scales ranging from seconds to seasons. This facilitates a study of the power spectrum of horizontal and vertical wind speed variance to identify the dominant variance-containing time scales.

## 1. Introduction

The quality of the representation of boundary layer turbulence in weather and climate models, as well as our understanding of the phenomenon, has greatly improved over the past decades. A number of observational supersites have played a key role in this regard. One such supersite is the Cabauw Experimental Site for Atmospheric Research (CESAR) located at Cabauw, the Netherlands. It features a 213-m high meteorological tower equipped with a wide range of measurement facilities. The tower was erected in 1972, and its capabilities have been continuously extended with new observational equipment, while great care was also taken to ensure the continuity of existing observations (Monna and Van der Vliet 1987; Beljaars and Bosveld 1997). Currently, the site is no longer characterized by the meteorological tower alone, but also features, for instance, a wide range of advanced ground-based remote-sensing facilities (e.g., Russchenberg et al. 2005). As such, the site also functions as a test bed for new measurement techniques.

CESAR observations were used in the development of a number of parameterizations (see, e.g, van Ulden and Wieringa 1996). The site was instrumental in a number of large intercomparison projects aimed at investigating existing parameterizations and developing new ones (e.g., Chen et al. 1997; Crewell et al. 2004; Bosveld et al. 2014).

In attempting to fully understand the turbulent phenomena of interest, boundary layer observations are often complemented by numerical models, ranging from regional weather models to high-resolution LES. A popular methodology is to set up an idealized numerical case study on the basis of observations from a measurement campaign. Such a numerical study has a number of advantages: 1) It represents an assessment of the capability of the current state-of-the-art models to simulate the measured boundary layer state. 2) The model describes an “ideal” experiment, since all input and output are known (in contrast to the actual measurement campaign). 3) The simulations can be repeated at will, allowing researchers to test a certain hypothesis, for instance, by slightly altering the simulation setup and studying the differences.

Examples of such studies are the large-eddy simulations of the Barbados Oceanographic and Meteorological Experiment (BOMEX; Siebesma et al. 2003) based on measurements by Holland and Rasmusson (1973), simulations of the Beaufort and Arctic Seas Experiment (BASE; Curry et al. 1997; Kosović and Curry 2000), or the GEWEX Atmospheric Boundary Layer Studies (GABLS; Holtslag et al. 2013) that has focused on the stable boundary layer. Accordingly, many theories and parameterizations have been tested against such datasets (e.g., Siebesma and Cuijpers 1995; Cheng et al. 2002; Heus and Jonker 2008). However, downsides to such an approach have also been identified (e.g., Neggers et al. 2012). First, these specific cases might not form a representative dataset to base models on. Second, the cases might not include those situations that are actually most troublesome for the models.

For these reasons, the KNMI Parameterization Testbed (KPT) was designed: a setting where data are gathered from observations, global/regional models, and single-column models (SCM) on a daily basis (Neggers et al. 2012). In the KPT, output from large-scale models (i.e., regional or global) is used to initiate SCM, which can then be compared to the available observational datasets. The key difference here is the statistical line of approach: by generating a dataset that covers long periods of time (i.e., from months to years), a statistically representative set of weather situations is formed.

Until recently, such long time datasets were unattainable in the realm of high-resolution large-eddy simulation (LES) modeling. The excessive costs related to LES runs precluded these models from covering large periods of time. However, recent computational developments now allow the first attempts to bridge that gap. Although computational costs still impose some severe limitations, we present two LES datasets that cover an entire year (2012) in a single, continuous simulation situated around the Cabauw supersite, based on the input of the regional weather model available in the KNMI parameterization test bed. One LES run was performed for its relatively large domain (25 km × 25 km), and the other was performed for its relatively high resolution (25 m). Such LES studies possess the advantages of a numerical case study, but mitigate the disadvantages relating to the representativeness of the simulations.

The rest of this paper proceeds as follows: Section 2 describes the setup of the LES runs. The resulting dataset is characterized and compared to observations in section 3. We continue the validation of the runs in section 4, in which we compare cloud properties with observations. Then, in section 5, we investigate the turbulent time scales present in the LES runs by recreating the classic wind speed spectrum of van der Hoven (1957), based on LES data, and compare this to multiple observational sources and the large-scale “parent” model. Prompted by the latter results, we finally use section 6 to illustrate a possible avenue of future study that might benefit from a dataset as introduced here. In this section, we revisit Taylor’s hypothesis of frozen turbulence based on the temporal and spatial data.

## 2. Simulation

The exponential increase in computational resources is now taken for granted. However, this growth is increasingly dependent on the addition of extra processor cores and less on improved processor speed. The consequence is that parallel computing (i.e., distribution of the work among multiple processor cores) has become the standard. Consequently, the number of computational cells used in simulations has increased in order to fully utilize the extra available cores. Whereas this has naturally led to steady improvements in domain and resolution, long time simulations have remained difficult as the time span cannot be computed in parallel.

We have recently ported our large-eddy simulation model to run completely on the graphical processing unit (GPU or video card), resulting in GPU-resident Atmospheric Large-Eddy Simulation (GALES). GPU residency, in this context, means that the computational data (i.e., the three-dimensional fields of momentum, temperature, etc.) remain on the GPU card. This approach can be contrasted with the “accelerator” approach by Michalakes and Vachharajani (2008) in WRF, where only selected routines are performed on the GPU. Instead, GALES performs *all* of its three-dimensional computations on the GPU. This avoids the communication bottleneck between GPU and CPU and allows for a more optimal use of the GPU’s computing power [see Schalkwijk et al. (2012) for details].

GALES is able to employ multiple GPUs to perform large-domain runs (Schalkwijk et al. 2015), but can also significantly improve time-to-solution performance (i.e., the time one has to wait for a given simulation to finish) for smaller simulations that reside on one GPU only. As a result, GALES can perform much longer simulations than previously feasible. This allowed us to perform a continuous simulation of 1 yr of boundary layer evolution around the Cabauw observational supersite, the Netherlands. We will refer to this run as YOGA-2012 (year of GALES, 2012) or simply YOGA.

### a. LES setup

GALES is based on the Dutch Atmospheric Large-Eddy Simulation, described extensively in Heus et al. (2010). For the setup of YOGA-2012, several modifications were implemented. Following Böing et al. (2012), the Boussinesq approximation was replaced by an anelastic approximation to account for changes in density in the vertical, allowing the extension of the domain to a 13-km altitude. Additionally, a simple, single-moment ice microphysics scheme is employed based on Grabowsky (1998) and Tomita (2001). Since this scheme is keyed toward tropical deep convective events, the autoconversion method proposed by Sundqvist (1978) was implemented, identical to the implementation in the European Centre for Medium-Range Weather Forecasts (ECMWF, cycle 31r1) model. This method also employs a simple modification to account for the Bergeron–Findeisen process and is expected to improve the representation of autoconversion in the midlatitudes in comparison to the original scheme of Grabowsky (1998).

To illustrate the setup of YOGA-2012, it is useful to separate phenomena by their typical length scale. Subfilter-scale phenomena occur on scales smaller than the LES filter scale. They are not explicitly calculated in the LES; the subgrid model is responsible for their representation. Resolved-scale phenomena are explicitly resolved by the LES. They are larger than the LES grid size, but smaller than the LES domain size. Large-scale phenomena are phenomena that occur on larger scales than the LES domain size. A formal decomposition of the interaction between these scales is presented in Soong and Ogura (1980) and Grabowski et al. (1996).

In GALES, subfilter-scale motions are treated through eddy viscosity/diffusivity fluxes, modeled as a function of the subfilter-scale turbulent kinetic energy. Their treatment in our model follows the initial ideas of Deardorff (1972) and is described in more detail in Heus et al. (2010).

In the simulations performed here, the large-scale motions are accounted for by coupling the simulation to the Dutch operational weather and climate forecasting model Regional Atmospheric Climate Model, version 2.1 (RACMO). The atmospheric dynamics in RACMO are based on the High Resolution Limited Area Model (HIRLAM); the physical processes are adopted from the ECMWF model. RACMO is extensively described in van Meijgaard et al. (2008).

In YOGA-2012, the LES is performed with periodic boundary conditions in the horizontal directions and thus is statistically homogeneous. This facilitates the conservation of resolved turbulent structures on the relatively small domain, but hinders the representation of domain-scale gradients. The following subsections describe how the large-scale phenomena are treated in YOGA-2012.

#### 1) Temperature and humidity

The atmospheric scalar fields that the LES integrates in time are *θ*_{l}, *q*_{t}, *q*_{r}, and : liquid water potential temperature, total (nonprecipitative) specific humidity, precipitative water–specific humidity, and the root of the subgrid turbulent turbulent energy, respectively. The total amount of water is given by *q*_{t} + *q*_{r}. The distinction between the two is as follows: *q*_{t} is assumed to follow the atmospheric flow, whereas *q*_{r} develops motions relative to the flow (i.e., precipitates). The microphysics scheme in the LES is responsible for the conversion between *q*_{t} and *q*_{r}.

Consider a conserved variable *φ* ∈ {*θ*_{l}, *q*_{t}}, the conservation equation of which can in general be written as follows:

where **u** = (*u*, *υ*, *w*) is the wind vector, and *S*_{φ} represents source terms that may apply. Molecular diffusion is neglected. The equation solved by the LES follows from filtering this equation with a filter length *λ* [extensively treated in, e.g., Wyngaard (2004)]. The resulting equation for the filtered variable is as follows:

where the subscript *h* is used to indicate the horizontal components. The first two terms represent turbulent transport on the resolved and subfilter scale, respectively, and are written in flux form using the anelastic approximation:

where *ρ*_{0}(*z*) is a time-independent base state density that varies in the vertical only. The third and fourth terms in Eq. (2) represent the vertical and horizontal components of transport on large scales. We employ the values of RACMO, denoted with superscript *R*, to represent the large-scale effects.

The RACMO runs that drive YOGA-2012 are constrained by observations. Every 24 h, RACMO restarts on the basis of a state constructed through data assimilation, but GALES continues uninterrupted (YOGA-2012 contains no intermediate cold restarts). To constrain GALES from drifting away from the observed state, the LES state is slowly relaxed to the state of RACMO on a time scale of *τ* = 6 h: this is the fifth term in Eq. (2). The time scale *τ* is chosen to be small enough for the LES to adapt to RACMO on the time scale at which RACMO performs data assimilation, yet long enough to allow the LES to develop some measure of independent turbulence (Neggers et al. 2012). Note that relaxation occurs with respect to the slab-averaged value 〈*φ*〉, to avoid affecting the magnitude of turbulent fluctuations (Heus et al. 2010).

Last, no GPU-resident full radiation scheme is available yet, and the coupling to a CPU-based radiation scheme would slow the computation down up to a point where a year-long simulation becomes infeasible. Therefore, we apply the RACMO radiative forcings [RACMO employs the GCM version of the Rapid Radiative Transfer Model (RRTMG) radiation module] to GALES. An additional advantage of this method is that this improves the comparability between the regional model and GALES; the disadvantage is that the radiative tendencies applied to GALES are detached from its actual thermodynamic state. Whereas the latter effect is anticipated to be small for clear air situations, it might cause a bias in cloudy cases. For instance, the cloud-top radiative cooling tendency that RACMO has may be applied to cloud-free regions in GALES if RACMO has a higher cloud top. Likewise, if GALES forms clouds that are absent in RACMO, no consistent radiative cooling is applied to the cloud field. Since both effects are anticipated to lead to cloud breakup in GALES, they might result in a net bias of reduced cloudiness.

#### 2) Momentum

The momentum equations can be written as follows:

where *p* is the pressure, *ρ* is the density, and Ω_{i} is the angular velocity vector representing Earth’s rotation; **F**_{u} represents other forces on the momentum equation including gravity and hydrometeor drag from precipitation.

The LES-filtered momentum equations that we have employed in YOGA-2012 are similar to Eq. (2), but include a number of additional effects:

where one can immediately recognize the first six terms from Eq. (2).

Large-scale (LS) pressure gradients are accounted for through the term . The geostrophic wind is acquired from RACMO. By expanding the pressure gradient term with the help of the environmental hydrostatic state, buoyancy forces are extracted from the pressure term, resulting in a modified pressure term with , where is the pressure difference with the environmental state, and *e* is the subgrid kinetic energy. Buoyancy forces are represented through fluctuations in the virtual potential temperature *θ*_{υ} and act only in the vertical direction; represents the vertical unit vector.

#### 3) Surface processes

The turbulent drag and exchange of scalars between the surface and atmosphere are parameterized in the surface model. For the simulations in this manuscript, we employed a land surface model to parameterize the sensible and latent evaporative heat fluxes at the surface.

The land surface model of DALES, described in Heus et al. (2010), was modified in order to mimic the Tiled ECMWF Scheme for Surface Exchanges over Land (TESSEL) (Viterbo and Beljaars 1995; van den Hurk et al. 2000) that is used in RACMO. The surface fluxes are calculated using horizontally averaged properties in order to create a horizontally homogeneous surface flux. This is consistent with the statistically horizontally homogeneous setup and increases comparability to large-scale models.

The surface model’s energy input is provided by a net radiative input term *Q*_{net}. The surface model divides this energy over the turbulent heat fluxes and a ground flux *G* such that the surface energy balance

is satisfied, where is the sensible heat flux and is the evaporative heat flux, with *ρ* as the density, *c*_{p} as the specific heat, and *L*_{υ} as the latent heat of vaporization. Since GALES applies the radiation of RACMO, *Q*_{net} represents an externally applied forcing to the LES. Surface model parameters (e.g., surface roughness) are used in accordance with RACMO. The prognostic variables in the surface model are the soil temperature *T*_{soil} and water content *ϕ*_{soil}. A relaxation term is added to these variables to prevent drift, equivalent to that in Eqs. (2) and (5).

### b. Processing of RACMO output

The LES runs were driven by RACMO forecasts that restart every 24 h on the basis of an observed state. To facilitate the study of turbulent time scales, however, we have performed the YOGA simulations as single, continuous runs, without interruptions. For this reason, some processing of the RACMO data was required.

First, all RACMO data used for driving the LES were averaged to an hourly series to ensure smooth forcings. Then, the resulting series were interpolated onto the LES time step, which adapts to the flow (Heus et al. 2010) and typically ranges between 0.5 and 5 s.

Second, the RACMO dataset was converted to a continuous time series. Each RACMO forecast that we employed starts at 1200 UTC, but there is significant overlap between forecasts. We interpolate the RACMO forecasts between consecutive days, smoothly switching from forecast to forecast between 1200 and 0000 UTC. This procedure is illustrated in Fig. 1. Between 1200 and 0000 UTC, the fraction *f*_{i} that combines the RACMO data of 2 days linearly increases between 0 and 1:

where indicates the *i*th day (*i* = 0, 1, 2, …) of RACMO input. The fraction *f*_{i} is given by

where *t*_{i} the time into day *i* since midnight. We have initialized the runs at *i* = 0 and *t*_{0} = 12 h, that is, 1200 UTC 31 December 2011, and used the first 12 h as initialization (*f*_{0} = 1 throughout initialization).

In the comparisons between RACMO and the YOGA simulations throughout this paper, the RACMO results have been concatenated to a continuous series using the same procedure. For this reason, and the reason that the driving RACMO run is set up differently (e.g., lower resolution) than it would be in an operational setting, the comparison should not be interpreted as an assessment of the quality of RACMO. Rather, the RACMO results provide a reference state to the LES.

### c. Computational choices

YOGA-2012 was set up to simulate a domain of 25 × 25 km^{2} at 100-m grid spacing in the horizontal directions. Vertically, the grid spacing increases as

starting at Δ*z*_{0} = 30 m and increasing with *α* = 0.9% per model level *k* to roughly 120 m at 13-km altitude. Although relatively coarse for current LES standards, some resolution must be sacrificed for long time series to remain feasible. A relatively large domain was chosen in order to pick up larger-scale motions and to keep large-scale convective phenomena from experiencing the domain limits. This domain is sufficient to resolve relatively intense convective events (e.g., cumulus congestus, scattered showers) since the tropopause typically lies at roughly 10-km height, but organized deep convective phenomena and fronts are naturally out of scope. The consequence of the larger domain, however, is that all small-scale turbulence is poorly resembled in YOGA-2012, most notably in the nocturnal boundary layer, where YOGA is heavily reliant on the subgrid scheme.

To better resolve these motions and thereby check the resolution dependence, YOGA-2012 is complemented by a second high-resolution simulation to which we will refer as YOGA-HR-2012 or YOGA-HR. This run is equivalent to YOGA-2012 in setup and forcings but is performed on a much finer grid, using a limited domain size. YOGA-HR employs a 4.8 × 4.8 km^{2} domain, using 25-m horizontal grid spacing. The vertical grid spacing increases from Δ*z*_{0} = 8 m at the surface, with *α* = 0.7%, to 40 m at the top of the domain, which now lies at 3.6 km. The domain of this simulation is clearly insufficient to represent large-scale cloud structures and even a significant portion of low clouds. Whereas 25 m is still too coarse to sufficiently resolve the stable boundary layer, turbulence is better represented in this simulation, providing more insight in, for example, turbulence close to the surface and in the morning and afternoon transitions.

Contrary to frequent practice, the LES domain is not translated with the mean wind (Galilean transformation). The continually changing wind in YOGA would require frequent adjustments to the transformation, which would probably cancel its beneficial effects. Moreover, the absence of translation eases the emulation of local ground-based measurements.

### d. Data output

The storage of the six, prognostic, three-dimensional fields at each time step throughout the year would require the storage of 2.5 petabytes per simulation. As this amount of data is hard to store and even harder to retrieve and wield, much of the data were processed during the simulation, writing only the resulting statistical properties to disk.

Additionally, we have stored both the full three-dimensional prognostic fields at low time-resolution (once a day) and high-resolution time series of these variables at low spatial resolution (four selected locations throughout the domain). The three-dimensional cloud field (liquid water–specific humidity, along with in-cloud values for *w*, *θ*_{l}, and *q*_{t}) was stored every 900 s, as the sparsity of this field allows for a very efficient compression. Altogether, the total output was reduced to roughly 0.5 terabytes per simulation.

## 3. Characterization

To provide a visual cue of the YOGA datasets, images of the Cabauw webcam are combined with cloud volume renderings of the YOGA-HR and YOGA dataset in Fig. 2. The images were created for several sample days at 1200 UTC to show a hint of the diversity of weather situations that were encountered during the runs. The cloud visualizations were made using GALES’ native volume renderer.

The difference between YOGA and YOGA-HR is immediately visually apparent in Fig. 2. The small domain of YOGA-HR seems to represent a “zoom” image of the YOGA cloud field. As a result, the limited vertical extent of YOGA-HR sometimes causes problems when tall clouds reach the upper-domain limit (second panel). Alternatively, the middle panel is an example of a case in which YOGA-HR completely lacks the high-altitude stratiform cloud deck, causing the simulation to look like a different state altogether, although actually the webcam, YOGA-HR, and YOGA all share a low-level cumulus cloud layer (remember that the domain top lies at 3.6 km for YOGA-HR and 13 km for YOGA).

Overall, both runs were able to represent a diverse range of weather situations. Even potentially problematic situations like the passage of frontal systems are represented by the LES, although their evolution is represented in time rather than in space (i.e., they arrive in and depart from the entire domain instantaneously). Of course, this is not to say that the LES represents such systems truthfully (this can hardly be expected), but it is encouraging that it remains stable and finds some form of representation for these conditions.

To provide a more quantitative evaluation of YOGA, we compare the average wind speed, temperature, and humidity between the simulation results and observations from the meteorological tower in Fig. 3. The figure shows the year-averaged diurnal cycle of wind speed , temperature *T*, and total specific humidity *q*_{t} 140 m above the surface. Simulation values are domain and time averaged; observation values are time averaged. For reference, the RACMO series that drives YOGA is also shown.

The difference in resolution between YOGA and YOGA-HR is most visible in the top panel of Fig. 3, where YOGA has problems following the diurnal cycle of absolute wind speed. In fact, YOGA even represents a worse match with observations than RACMO does. YOGA-HR performs much better in this regard and is the only simulation that shows a diurnal cycle that resembles the one found in observations (although slightly exaggerated). This diurnal cycle persists even though YOGA-HR is continually being nudged toward the RACMO values, which indicates that turbulent time scales this close to the surface dominate over the nudging time scale. Possibly, the diurnal cycle would be even more pronounced without nudging.

The resolution of YOGA effectively implies that the simulation is largely dependent on the subgrid formulation at 140-m height. Figure 3 indicates problems in the subgrid scheme in representing weather conditions that, as YOGA-HR shows, can be represented through explicit resolution of turbulence.

The diurnal cycle of temperature (middle panel) is closely captured by both simulations, but is admittedly straightforward. The diurnal cycle of humidity is more interesting, as humidity peaks at sunset and dawn and dips around noon. These features are presumably caused by the interplay between surface fluxes and boundary layer growth. In the morning, increasing latent heat fluxes cause a buildup of moisture in a shallow boundary layer. Then, strong boundary layer growth spreads the accumulated humidity over the increased boundary layer height, while entrainment causes further drying. As these processes diminish toward the evening, humidity increases to a second peak. The double-peak structure is reproduced by all simulations. The simulations are close also in absolute terms, differing at most 0.2 g kg^{−1} with observations. The YOGA simulations are slightly moister than RACMO, and YOGA-HR again has the smallest mean error.

The vertical structure of the boundary layer is shown in Fig. 4, year averaged for daytime and nighttime conditions. Daytime, in this manuscript, corresponds to the time between 1 h after sunrise and 1 h before sunset. Similarly, nighttime is the time between 1 h after sunset and 1 h before sunrise. Since this procedure places more weight on the summer months for daytime data and on the winter months for nighttime data, the profiles in Fig. 4 are not directly comparable to Fig. 3.

The wind speed profile is rather sensitive to the resolution and turbulence model, especially in stable conditions (e.g., Beare et al. 2006). Generally, the 30-m vertical resolution of YOGA is insufficient to resolve the boundary layer structure, especially at night. In fact, wind profiles of YOGA do not satisfy Monin–Obukhov similarity theory (MOST) near the surface (Businger et al. 1971; Dyer 1974). Although GALES imposes MOST at the surface, the presently used subfilter-scale model seems incapable of realizing a transition that satisfies MOST between the surface and the well-resolved (and indeed MOST satisfying) regime. This important issue is currently under investigation. RACMO, featuring parameterizations better suited to coarse resolution, performs better. At 8-m vertical resolution, YOGA-HR can sustain steeper gradients and performs slightly better than RACMO.

For *θ*_{l} and *q*_{t}, the agreement between simulations and observations is very good over the tower heights. There is remarkably little difference between YOGA, YOGA-HR, and RACMO (the difference is at most 0.5 K and 0.2 g kg^{−1} throughout the boundary layer), which indicates there are few consistent differences between the models. Close to the surface, YOGA lacks resolution to resolve the steep surface gradients in *θ*_{l} and *q*_{t}, as it exhibits weaker gradients than found in observations. Therefore, its better numerical match for daytime humidity is probably a case of compensating errors. In general, YOGA-HR shows a better vertical structure and is slightly warmer during daytime, which seems to correspond to measurements. The gradients in the warm surface layer are still underestimated though, even at 8-m vertical resolution.

The evolution of the boundary layer is ultimately dominated by a simple balance between surface fluxes and sources (e.g., Lilly 1968; van Driel and Jonker 2011; Schalkwijk et al. 2013). Therefore, differences in the profile means might be directly related to the magnitude of the surface fluxes. Figure 5 compares the climatology of the LES surface fluxes with the year 2012 of the long-term Cabauw observational surface dataset [see Beljaars and Bosveld (1997) for details]. The observational data have been corrected to ensure the closure of the surface energy balance by increasing the sensible and latent heat fluxes, keeping *H*/*E* constant, until the surface energy balance is satisfied. The shaded regions around the observational data (solid lines) indicate uncertainty estimates.

The net energy input at the surface is provided by *Q*_{net}. Since GALES applies the radiation of RACMO, this term is equal for all simulations. The land surface model of the LES ascertains the distribution of that energy over the sensible heat flux *H*, the latent (evaporative) heat flux *E*, and the ground flux *G*. The modeled *Q*_{net} overestimates the observed values, most notably in the summer months. Latent heat flux coincides well with the observations, whereas sensible heat flux is overestimated.

This implies that the land surface scheme favors a higher Bowen ratio in Cabauw than observed. This effect may be because of irrigation measures that are unaccounted for, which wet the soil in dry periods. Therefore, were *Q*_{net} to better conform to measurements, the sensible and latent heat fluxes would not necessarily improve.

This is true for all models, which agree well on the surface flux distribution. This was expected, as both the LES runs and RACMO utilize the same land surface parameterization. The agreement on surface fluxes is also consistent with the close intermodel agreement on the average values provided in Figs. 3 and 4.

## 4. Cloud properties

In conclusion to the previous section, we find that the boundary layer processes are better resolved, and generally closer to observations, in YOGA-HR than they are in YOGA. This should come as no surprise, given the fourfold difference in resolution (100- vs 25-m horizontal and 30- vs 8-m vertical resolution at surface level). The added value of YOGA-2012 lies in its much larger domain (25 vs 5 km), which is mostly relevant for the representation of moist convection.

The most robust method of observing cloud cover at Cabauw is arguably the ceilometer. This device (Vaisala LD-40) determines, at high temporal resolution, at which height (if any) a cloud is found directly above the device. It thus provides a sequence of cloud heights *z*_{c,i} with *i* = 1, 2, 3, …, *N* within a time frame *T*, where the number of observations *N* = *fT* is dependent on the measurement frequency *f* and the time window.

The ceilometer output can be converted to an average cloud cover *C*_{T}(*z*_{max}) for clouds found below *z*_{max} by counting the number of hits satisfying 0 < *z*_{c,i} < *z*_{max}:

where *H*(*z*) is the Heaviside function.

In the YOGA runs, we have implemented several “virtual” ceilometers, which determine *z*_{c,i} at four locations spread evenly over the domain. In the LES, *z*_{c,i} is defined as the height of the lowest grid cell that has nonzero liquid water content in the “ceilometer” column. This allows a direct comparison to Cabauw observations. This is impossible in RACMO, however, as the cloud cover is parameterized, such that we must use RACMO’s statistical, parameterized values for comparison.

Therefore, Fig. 6 compares the yearly averaged diurnal cycle of total cloud cover *C*_{T}(∞) and of low cloud cover *C*_{T}(3 km) with the parameterized total and low cloud cover in RACMO, respectively. The cloud cover was calculated for *T* = 900 s, but is relatively insensitive to the value of *T*. Figure 6a shows that the total cloud cover compares very favorably with observations for both YOGA and RACMO. The domain height of YOGA-HR is insufficient for a fair comparison and is left out.

But even the domain of YOGA might be too small to resolve high cloud dynamics. Therefore, Fig. 6b focuses on the *low* cloud cover, that is, *C*_{T}(3 km). Both RACMO and YOGA seem to underestimate the low cloud cover. Consequently, comparison with Fig. 6a suggests that these models thus overestimate the *high* cloud cover. This conclusion must be met with caution, however, since the ceilometer cloud cover estimates become less trustworthy with height. The significant overestimation by YOGA-HR of the daytime cloud cover is probably due to numerical artifacts as the clouds approach the top domain edge.

The performance of low cloud predictions is further assessed in a simple “quality matrix” in Table 1. The table separates “cloud states” in three classes: cloudless, broken clouds, and an uninterrupted cloud deck. The low cloud cover is abbreviated *c* = *C*_{T}(3 km) for *T* = 900 s. If all simulations were to reproduce the observations in terms of cloud state, we would find 100% at the diagonal positions and 0% in all other positions.

Table 1 suggests that RACMO creates broken clouds too often, in cases where observations indicate cloud-free or stratiform clouds. YOGA and YOGA-HR are more successful in representing clear-sky situations; a clear sky is reproduced in roughly three-quarters of the times in both YOGA and YOGA-HR. Failure to reproduce a clear-sky situation is often related to delayed cloud breakup or accelerated cloud formation. However, the success of representing clear sky is partially related to an overall shortage of cloudiness in both YOGA runs.

A broken cloud deck is reproduced more often in YOGA than in YOGA-HR (53% against 44%), whereas the full cloud deck is reproduced more often in YOGA-HR. This might have been expected, as the broken cloud deck is associated with cumulus convection, the representation of which requires a relatively large and high domain, whereas the full cloud deck is associated with stratiform clouds, requiring high resolution.

Radiative cooling is also an important effect in stratiform clouds, so the noninteractive radiation employed in both runs is a weakness that might preclude the LES from sustaining a stratiform cloud deck (for instance, misaligned, radiative, cloud-top cooling may cause cloud breakup, as discussed in section 2), potentially causing an underestimated cloud cover. Another potential reason for the underestimation of cloud cover by the LES runs is the limited variability due to the domain size. Since cloud formation often requires sufficient variance, the failure to resolve variance at scales larger than the domain size may artificially limit cloud formation.

### Sensitivity

For low clouds, Table 1 indicates that the performance of the YOGA runs is worst in the case of broken clouds (cumulus clouds). This is somewhat unexpected because LES models are generally believed to perform well for these types of clouds (e.g., Stevens et al. 2001; Siebesma et al. 2003). To test the sensitivity with regard to the setup, we first selected a manageable subset of days that are characteristic for surface-driven cumulus clouds. By requiring that broken clouds are present [0 < *C*_{T}(3 km) < 0.8] during more than 50% of the daytime, and that the sensible heat flux is above average (*H* > 40 W m^{−2}), we selected 10 convectively driven cumulus days. These days were rerun individually using a setup identical to YOGA, but with the following differences:

The cloud cover for low clouds, for the selected 10 days, is shown in Fig. 7. Figure 7a shows the cloud cover for the selected days for the original cases and compares it with RACMO. In general, the results for this subset are in line with Fig. 6b; YOGA creates a higher cloud cover than RACMO, but still less than observed. It is therefore unsurprising that Fig. 7b shows an improved (increased) cloud cover as the relaxation is removed and decreased cloud cover as the relaxation strength increases. The sensitivity to the relaxation time scale *τ* is not very strong; however, as *τ* ranges between 1 h and ∞, the cloud cover changes by roughly 10%. Note that the reference line (dashed blue) also slightly differs between Figs. 7a and 7b, since Fig. 7b shows reruns of individual days (i.e., not made continuous as in section 2b).

The full radiation scheme of RRTMG has been coupled to the LES for reruns 3 and 4. Since this radiation scheme is not available for the GPU, the radiation calculations were performed on the CPU once every 10 simulated minutes. This setup was yet infeasible for the full YOGA run because of the computational constraints. When coupled to RACMO in setup 3, the use of an interactive radiation scheme has little effect in terms of cloud cover. Only when the LES is run independently (setup 4), the radiation has a significant effect, as it can now significantly alter the mean thermodynamical state. That the cloud cover in setup 4 is closest to observations is encouraging, since it indicates that the more independent the LES becomes, the better it matches observations.

## 5. Spectral comparison

The year-long LES runs of this study provide a novel opportunity to study time scales ranging from seconds to seasons on the basis of a single, uninterrupted large-eddy simulation run. This allows us to revisit the power spectrum of horizontal wind speed, which was first constructed by van der Hoven (1957). In his classic paper, he created the power spectrum by piecing together various portions of the spectrum, retrieved from different sources at different periods of time. We can now recreate these spectra on the basis of LES by a straightforward Fourier transformation of the LES data.

Figure 8a shows the temporal power spectrum over the full year of 2012, averaged over the four virtual towers, both for YOGA and YOGA-HR. The light blue shaded region shows the fast Fourier transform of the YOGA-HR dataset; the thick lines represent an average over exponentially increasing bin sizes. For reference, the RACMO data that serves as input to the LES model are also shown (solid red line).

Since no single observational device at Cabauw provides the full year of wind velocities from second to year scale, the spectra are compared to a couple of sources. The cup anemometer (solid black line) provides reliable long-term wind speed measurements, but lacks accuracy in the high-frequency domain. Therefore, an estimate of the high-frequency variance is provided by sonic anemometer data. The sonic anemometers were not operational throughout the entire year of 2012, and part of the data had to be discarded because of interference from the tower wake, resulting in a fractured dataset. The sonic anemometer spectrum was therefore constructed by averaging the spectra of a number of continuous time series of at least 12 h (~80 of such time series were found, scattered throughout the year). The dataset is probably not fully representative, as indicated by the discontinuity between the cup and sonic anemometer data, but should give an idea of the general trend.

The YOGA runs coincide relatively well with observations in both the high-frequency and low-frequency regimes. At the low-frequency end, the spectrum is dominated by daily and seasonal time-scale influences. At the high-frequency end, the spectrum is dominated by boundary layer turbulence. In the intermediate regime, roughly between 12 and 0.5 h, the observations suggest a *ω*^{−5/3} falloff from the low-frequency regime, associated with the turbulent energy cascade.

This −5/3 falloff is not witnessed in the spectra of YOGA and YOGA-HR, which show a “gap” in the intermediate regime. This gap is found where the resolved turbulence scales of the YOGA simulations transition into the time scales of the input forcings. At time scales larger than roughly 12 h, the simulation input provides sufficient information. The LES itself can create turbulent time scales up to roughly *L*/〈*U*〉, with *L* as the domain size, amounting to 1.5 h for YOGA and 900 s for YOGA-HR, for 〈*U*〉 ≈ 5 m s^{−1}. Since the RACMO input too features little variance at scales smaller than 12 h, the LES has no variance source at these scales and thus shows too little variance. For an LES run to be able to sufficiently resolve the variance at these scales, we require *L* ≃ 200 km. Although LES runs of that scope are now becoming possible for a period of a day or so (e.g., Khairoutdinov et al. 2009; Schalkwijk et al. 2015), they are presently infeasible for longer periods.

Finally we plot in Fig. 8b the spectrum of vertical velocity. For large time scales, one observes a power law (~*ω*^{+1}) covering more than four decades. This implies that the integral time scale of the vertical velocity at 180-m height is well defined and smaller or of the order of 1 h. This corresponds with the notion that at low levels in the atmosphere vertical movement is mainly produced by boundary layer processes, whereas these motions are limited by the presence of the earth surface.

A discussion on the time and length scales at which turbulent transport occurs in the YOGA runs is provided in a companion paper (J. Schalkwijk et al. 2014, unpublished manuscript) in which we investigate the cospectra and .

## 6. Taylor’s hypothesis of frozen turbulence

In the previous section, we concentrated on spectral characteristics in the time domain. Because of the nature of the forcing of the model we could span a very large range of scales. In the spatial domain however, the horizontal spectral scales are cutoff at the domain size because of the periodic boundary conditions. Here, we look into the relation between temporal and spatial fluctuations. Taylor’s hypothesis of frozen turbulence is often invoked (either implicitly or explicitly) to relate temporal characteristics of turbulent quantities to spatial characteristics. The hypothesis is that if the advective velocity *U* is much larger than the turbulence velocity *u*′, the changes in time of a quantity *ϕ* at a fixed point “are simply due to the passage of an unchanging pattern of turbulent motion” (Taylor 1938, p. 478), that is,

allowing one to relate the temporal to the spatial characteristics.

The YOGA simulations provide a dataset that describes a unique realization of turbulence that is (at least statistically) representative of the observed turbulence in Cabauw. Since both data with high spatial resolution and data with high temporal resolution are available, we might explore the difference of time-based and spatial-based determination of turbulent properties. Recent studies have explored the relation between the spatial and temporal view on turbulence by applying Taylor’s hypothesis in the context of shear-driven turbulence in the channel flow (Moin 2009; Del Álamo and Jiminéz 2009) and in the atmospheric boundary layer (Higgins et al. 2012).

The spatial field of the vertical velocity *ϕ* = *w* was compared with the time series in four virtual “towers” in the LES, daily at 1200 UTC. As the full fields are not stored for hardware reasons, averaging is performed over the four towers and over the year, as denoted by the overbar .

The most direct test of Taylor’s hypothesis is the direct application of Eq. (11); we might directly calculate the relative RMS error between the left- and right-hand side of this equation

where *σ*_{ϕ} is the standard deviation of *ϕ*: . Note that *δ*_{ϕ}(*τ*) is closely related to the correlation between *ϕ*(*x*, *t* + *τ*) and *ϕ*(*x* + *Uτ*, *t*); *δ*_{ϕ}(*τ*) = 0 corresponds to a correlation of 1, and *δ*_{ϕ}(*τ*) = 1 corresponds to a correlation of 0.

The advective velocity *U* is estimated as the slab-averaged velocity and has an angle *α* = tan^{−1}(〈*υ*〉/〈*u*〉) with an east–west direction (*α* = 0 for westerly wind). Figure 9a shows the resulting RMS error *δ*_{w}. The error rapidly increases with *τ*, especially at low heights (<100 m) where the error saturates within 1–2 min. At larger heights, some correlation remains even after 5 min, since here larger-scale motions exist that have longer lifetimes and thus are better suited to Taylor’s hypothesis.

In a recent paper, Higgins et al. (2012) investigated Taylor’s hypothesis in terms of the autocorrelation function *R*_{ϕ}:

which is a generalized formulation of *δ*_{ϕ}(*τ*). It represents the correlation between *ϕ*(*x*, *t*) and *ϕ*(*x* + *ζ*, *t* + *τ*) for any combination *ζ*, *τ*. Note that for *ζ* = *Uτ*, *R*_{ϕ} can be directly related to *δ*_{ϕ}(*τ*):

The autocorrelation *R*_{w}(*τ*, *ζ*) is shown in Figs. 9b and 9c at 33 and 101 m above the surface, respectively. By plotting *R*_{w} as a function of *Uτ* and *ζ*, the correlation is shown in terms of length scales instead of time scales, using the mean velocity at the given height to convert the time scales in length scales.

The correlation is strongest on the diagonal, *ζ* = *Uτ*, as expected. In comparing the top panel with the lower panels of Fig. 9, note that the difference between heights is amplified in the lower panels. Indeed, by plotting the correlation in terms of length scales, we amplify the decreased correlation at lower heights with the fact that wind velocities rapidly diminish toward the surface. This results in a much smaller decorrelation length scale close to the surface.

The small decorrelation length might be misleading if one is interested more in the statistical properties of the spatial structure than in the representation itself. Therefore, we check whether the reconstructed (using temporal measurements) spatial representation of *ϕ* is statistically equivalent to the actual spatial representation of *ϕ*. This is investigated by comparing the spatial energy spectrum *E*_{ϕ}(*k*) with the reconstructed spatial energy spectrum . Note that this is in fact the context in which Taylor (1938) posed the hypothesis that now bears his name. The spatial energy spectrum decomposes the variance of a signal *ϕ* into wavenumbers *k*, that is,

Figure 10 compares *E*_{ϕ}(*k*) (solid lines) with (dashed lines) for *ϕ* ∈ {*w*, *θ*_{l}, *q*_{t}} at heights between 8 and 100 m. The left panels show a direct comparison, and the right panels show the ratio . The spatial and reconstructed spectra are roughly equal for all variables at length scales between 300 m and 1 km. For smaller length scales (larger wavenumbers), the reconstructed spectra drop significantly faster than the spatial spectra. This is presumably because of the numerical issues since the spectra deteriorate close to the numerical resolution. The difference between *E*_{ϕ} and at small scales increases with height, which is related to the increase in advective velocity *U* with height, which results in a poorer high-frequency sampling.

There is also a significant difference between reconstructed and spatial spectra at large scales. This is because of the slow transients in time that are not represented spatially. The statistically horizontal homogeneous fields cannot represent the large scales related to the diurnal cycle that are found in the time (and thus also in the reconstructed) spectra. This is especially evident in the spectra of *θ*_{l}, but *w* and *q*_{t} also undergo a diurnal cycle, although smaller in magnitude.

Results from section 5 suggest that the decrease of turbulence intensity at larger scales may simply be caused by the limited domain size, in which case the temporal fluctuations may even be more accurate. On the other hand, effects related to the diurnal cycle are governed by a steady time scale that cannot be converted to a length scale using any advective velocity *U*. Therefore, we might wonder whether a similar simulation setup, but with much larger domain size, would (and should) show much larger *spatial* scales.

All in all, given the temporal shortcomings of the YOGA setup observed in Fig. 8 and the obvious spatial shortcomings related to the limited domain size and resolution, our analysis turns out to be not so much a test of Taylor’s hypothesis but rather an illustration of the spectral issues one encounters in the type of simulations ventured in this study. As computational resources increase, future large-scale simulations might shed further light on the subject.

## 7. Conclusions

Computational developments now allow LES runs that cover a time span of a year. The YOGA and YOGA-HR runs presented in this manuscript serve as illustrations of this principle. As computational resources continue to increase, long time runs will become relatively cheaper and easier to perform and therefore also more common.

The runs presented here can be much improved. For one, the resolution of the runs, even for YOGA-HR, is insufficient to resolve the turbulence near the surface or in the stable boundary layer. Second, the lack of interactive radiation precludes the investigation of the cloud feedback on radiation. Also the microphysics model, responsible for the formation of rain, might require further attention.

Nevertheless, several preliminary conclusions can be drawn from the YOGA runs. First and foremost, it is encouraging that the LES philosophy of explicit simulation of turbulence is stable enough to simulate an entire year of conditions that vary from stable boundary layers to deep convective events. Specifically, we found that the part of the boundary layer that is explicitly resolved coincides best with observations. This might act as “proof of concept” for future studies that might further improve on the representation of modeled processes.

Furthermore, the simulation setup utilized in this study can be generalized to a predictive approach, since simulations complete faster than real time and can be forced by predictive models. Therefore, the overall quality and stability are important factors that encourage future research in the application of LES in forecasting, for instance in the form of a “super-parameterization” in large-scale models (as in, e.g., Grabowski 2001).

A spectral decomposition of horizontal wind speed variance indicates that the larger scales (>10 h) are sufficiently represented, in comparison to observations, through the coupling with a weather model. The small (<1 h) scales are created by the LES itself and also seem reasonable. The observations show that substantial energy is present in the intermediate regime, which is insufficiently resolved by the YOGA runs. This regime is of importance for the short-term weather forecasting that is pursued with present-day mesoscale atmospheric models.

The lack of resolved variation in these scales is related to the limited domain size of the YOGA runs; to resolve all these intermediate scales one would require a domain of roughly 200 km, which is presently out of reach computationally. The vertical velocity variance, however, shows no intermediate-scale gap and is remarkably constant in intensity over all frequencies.

## Acknowledgments

This work would not have been possible without the KNMI Parameterization Testbed (KPT). We are thankful for the continuous efforts of the KNMI team responsible for maintaining this project. Of particular value has been Dr. Roel Neggers for his initial work on the test bed, Dr. Erik van Meijgaard for generating and providing the advective forcings from RACMO, and Dr. Bert van Ulft for continuously maintaining the test bed. We are also indebted to the CESAR team in general and Ir. Henk Klein Baltink in particular, who continue to ensure the availability of indispensable observational datasets.

## REFERENCES

*IEICE Trans. Commun.,*

**E88-B,**2252–2258.

*Bull. Amer. Meteor. Soc.,*in press.