The Benefits and Challenges of Downscaling a Global Reanalysis With Doubly-Periodic Large-Eddy Simulations

Global reanalyzes like ERA5 accurately capture atmospheric processes at spatial scales of 𝐴𝐴  (10) km or larger. By downscaling ERA5 with large-eddy simulation (LES), LES can provide details about processes at spatio-temporal scales down to meters and seconds. Here, we present an open-source

large domain and/or grid nesting to allow turbulence to fully develop at the inflow boundaries, making the simulations expensive (Heinze et al., 2017a).The periodic boundary LES setup-which we will employ-circumvents this problem by using doubly-periodic lateral boundary conditions, only coupling LES to the large-scale weather by applying horizontally homogeneous but time and height varying large-scale forcings.These forcings typically contain the advective tendencies of heat, moisture, and momentum, the large-scale vertical (subsidence) velocity, and the geostrophic wind components (e.g., Heinze et al., 2017b;Neggers et al., 2012;Schalkwijk et al., 2015;van Laar et al., 2019).Although such a setup clearly has shortcomings in complex large-or mesoscale conditions, it allows for smaller domains.This makes the simulations computationally cheaper, or allows the use of a finer computational mesh and hence resolve turbulence on smaller scales.As we will demonstrate, such high resolution LESs are needed to capture high frequency interactions between solar radiation and clouds, which is a key process required to advance our understanding of land-atmosphere interactions (e.g., Gentine et al., 2019;Vilà-Guerau de Arellano et al., 2023).
Previous work on realistic periodic boundary LESs was often done using limited area host models (e.g., Heinze et al., 2017b;Neggers et al., 2012;Schalkwijk et al., 2015).One advantage of these models over a global model like ERA5 is the increased horizontal resolution.Nevertheless, we have chosen to use ERA5 for a number of reasons.First, ERA5 has a global coverage, and therefore does not limit the LES simulations to a certain geographical region.Second, ERA5 covers a long time period from 1950 to within 5 days from real time.And third, and perhaps most important, the ERA5 data can easily be accessed through the Copernicus Data Store (CDS [https:// cds.climate.copernicus.eu]),bypassing the need to request data from various national weather services in case it is not openly available.Another commonality in most previous work is that the underlying code used to preprocess the host model data was closed-source.In order to improve scientific transparency and reproducibility, we released our code as an open-source Python package named the "Large-eddy simulation and Single column model-Large-Scale Dynamics," or (LS) 2 D (pronounced as "el es too di") in short.This Python package contains a number of routines meant to simplify and automate all the steps required to generate the input for doubly-periodic LESs.
The goal of this paper is threefold.First, we document (LS) 2 D by describing the coupling method between ERA5 and LES, and the procedures in the (LS) 2 D Python package (Section 2).Although (LS) 2 D has already proven its usefulness and skill in a number of studies (Kreuwel et al., 2020;Mol, van Stratum, et al., 2023;Ražnjević et al., 2022;Tijhuis et al., 2023;Veerman et al., 2022), this previous work mostly focused on clear or shallow convective conditions.Therefore, the second goal of this paper is to quantify the added value of doubly-periodic LES with respect to ERA5 for (a) dry and shallow convective conditions, and (b) more challenging conditions.We do this by simulating the full month of August 2016 over Cabauw (the Netherlands) with (LS) 2 D and MicroHH (van Heerwaarden et al., 2017), and validating the results with the Cabauw surface, tower, and remote sensing observations (Sections 4.1 and 4.2).The third goal is to discuss the pros and cons of using a doubly-periodic LES setup by comparing a key process-the surface solar irradiance variability-with detailed observations from the Baseline Surface Radiation Network (BSRN, Knap, 2022;Mol, Knap, & van Heerwaarden, 2023).With a sensitivity experiment consisting of six different 1-month LESs, we address the impact of domain size and resolution on the ability of LES to capture different surface solar irradiance time scales (Section 4.3).Finally, we discuss our findings, and conclude the paper in Section 6.

Methods
To account for the influence of the large scale weather acting on the local LES domain, (LS) 2 D creates a one-way coupling between ERA5 and LES.The coupling method that we apply is similar to the methods described by Heinze et al. (2017b), Neggers et al. (2012), andSchalkwijk et al. (2015).In this approach, the initial conditions (atmosphere and soil) and a number of time and height varying large-scale processes acting on the LES domain, are derived from routine output of (in our case) ERA5.These processes, partly in the form of a tendency of a state variable like temperature, humidity, or wind, are then added to the prognostic LES equations, as shown in Equations 1 and 2: 10.1029/2023MS003750 3 of 14 (2) Here, ψ is a generic scalar, in MicroHH either the liquid water potential temperature (θ l ), the total specific humidity (q t ), or other scalars, and u i is a vector with the horizontal wind components (u, v, 0).As MicroHH uses a saturation adjustment scheme for cloud water and ice, these quantities are implicitly included in q t .(LS) 2 D also provides the individual specific humidity phases for models which have for example, prognostic cloud water and/ or ice.Variables with a superscript LS are variables from ERA, and variables with a tilde denote the filtered LES fields.The angle brackets indicate a horizontal averaging operation, either over a number of ∼30 × 30 km 2 grid points in ERA5, or the entire LES domain.Depending on the local and meteorological conditions, the spatial averaging of the host model may be required to prevent processes from being represented twice, when they are both included in the large-scale forcings from the host model, but also resolved by LES itself.
The coupling between the resolved ERA5 variables and the turbulent LES variables consists of a number of terms: • The advective tendency contains the resolved advective tendency from the host model.Extracting this tendency directly from the dynamical core of IFS would potentially be more accurate, but as this is highly impractical, this term is approximated offline from the three-dimensional fields.• The subsidence term contains the interaction between the large-scale vertical velocity from ERA5, and the turbulent LES fields.• The source term can contain any external forcing, for example, radiative heating rates for LES simulations without interactive radiation.• The Coriolis term contains the influence of the large-scale pressure gradient on the horizontal wind components in LES, through the geostrophic wind.• Finally, the relaxation term is a safety measure for long experiments, which prevents the LES model from deviating too far from the host model, by nudging the horizontal mean state of LES to the mean state of ERA5, on a timescale τ n .
In these equations, f c is the Coriolis frequency f c = 2 Ω sin(ϕ), where Ω = 7.2921 • 10 −5 rad s −1 and ϕ is the latitude.The two horizontal geostrophic wind components are denoted by   LS g;j , and are defined as: where ∂Z LS /∂x i are the horizontal gradients of the geopotential height on constant pressure levels, and g is the gravitational acceleration.
The momentum tendency in the IFS model contains an additional term to account for the drag from subgrid-scale orography, which has been omitted from Equation 2.
All terms from Equations 1 and 2 are calculated on the native (137) ERA5 model levels, as this provides the most detail in the vertical.The geostrophic wind components are an exception, which are calculated on pressure levels, and next interpolated to the model levels.All terms additionally vary in time.They are, however, applied as horizontally homogeneous quantities in the LES domain.Although it is in principle possible to calculate and apply spatially varying tendencies, this approach would be difficult to unite with the choice for periodic boundary conditions in LES.
The initial conditions for LES (θ l , q t , u, v) are directly interpolated from the ERA5 model levels onto the vertical LES grid.In line with the large-scale forcings, these initial conditions are optionally spatially averaged over a number of ERA5 columns, and next applied as horizontally homogeneous quantities in the LES domain.
The gradients ∂ϕ LS /∂x i are approximated from the ERA5 data with either second or fourth-order accurate centered finite differences: 10.1029/2023MS003750 4 of 14 This is a user choice, after which the chosen method is applied consistently in both space and across different variables.

Python Package
The main task of (LS) 2 D is to automate all the steps required to generate the input for doubly periodic LES.To explain the typical workflow with (LS) 2 D, we will step-by-step walk through a simplified (LS) 2 D Python script, shown in Figure 1.
After installing the Python package, (LS) 2 D is available with import ls2d.To simplify passing the case settings to (LS) 2 D, most settings are gathered in a Python dictionary.Next, the first step is to download the ERA5 data.(LS) 2 D always first checks whether the required ERA5 data has already been download before, to prevent unnecessary duplicate downloads.If the ERA5 data needs to be downloaded, (LS) 2 D supports two methods, either using the Copernicus Data Store (CDS, openly available), or the Meteorological Archival and Retrieval System (MARS, requires ECMWF supercomputer access).In both cases, and especially when using CDS, the queuing time can be substantial, as the model level data is only stored in the tape archive.Therefore, when using the CDS option, (LS) 2 D will store the unique CDS request IDs once the downloads are submitted, and stop the Python script.On a subsequent launch of the Python script, download_era5() will check if there are active request, if so check if they are finished, and if this is the case, download the NetCDF files from the CDS server.When using MARS to download the ERA5 data, the MARS retrievals are submitted using the SLURM workload manager.
Once the required ERA5 NetCDF files are available, (LS) 2 D has the Read_era5() routine which reads the required files, and calculates some derived properties like the model level pressure and height, and state variables in other units.Next, calculate_forcings() calculates the required terms from Equations 1 to 2. The n_av option allows the user to average the initial conditions and large-scale forcings over ±n_av ERA5 grid points, and the method argument switches between the second and fourth-order accurate finite differences.
After specifying a vertical LES grid (in this case a stretched grid, where the kth level has a grid spacing Δz = Δz 0 (1 + α) k ), the get_les_input() routine interpolates all vertical profiles to the specified LES grid.The resulting les_input object is an Xarray data set (Hoyer & Hamman, 2017) containing all the required LES input, and additional information related to surface variables like for example, the roughness lengths, leaf area index (LAI), and vegetation and soil types from ERA5.For LES simulations without an interactive land-surface model, (LS) 2 D also provides the time varying surface sensible and latent heat fluxes.Each variable contains attributes describing the variable and its units.This data set can easily be saved in NetCDF format using xarray's to_netcdf() method, or further processed in the (LS) 2 D script, and saved into the LES model specific input format.

LES Model
We used MicroHH for all LES simulations in this study.Note however that (LS) 2 D is not specifically designed for MicroHH, and can (and has) been used with other similar LES models.The core of the model is described in detail by van Heerwaarden et al. (2017).Over the last years, the model has been extended with various new physics parameterizations required for simulations of realistic weather.This includes the RTE-RRTMGP radiative transfer model (Pincus et al., 2019), an interactive land-surface model based on HTESSEL (Balsamo et al., 2009), and a single moment ice microphysics scheme (Tomita, 2008).Describing all model settings is practically unfeasible for such complex LES setups.Therefore, we provide a limited description of the primary settings required to understand the nature of the experiments.The LES code 10.1029/2023MS003750 5 of 14 and input are available as supplementary material, and all options can be found in our online documentation (https://microhh.readthedocs.io).
The LES domain size and horizontal resolution varied between the different experiments, as specified in Table 1.In the vertical direction, we used a stretched grid with a grid spacing of 20 m near the surface, stretched over 192 levels up to ∼18 km height.Although the largest domain approaches the size of 2 × 2 ERA5 grid points, all simulations used (LS) 2 D forcings from only a single ERA5 grid point.For advection, we used the fifth order scheme from Wicker and Skamarock (2002), and a non-dynamic Smagorinsky subfilter-scale model.The radiative fluxes, calculated with RTE-RRTMGP over the full 3D fields, were updated every 60 simulation seconds.The experiments used the interactive HTESSEL land-surface model, which predicts the surface skin temperature and humidity, and calculates the surface momentum and scalar fluxes using Monin-Obukhov similarity theory.We used a spatially homogeneous surface, with short grass with a vegetation fraction of 95% and a LAI of 2.6, and the medium-fine soil type from IFS (ECMWF, 2016).Finally, the nudging timescale τ n (Equations 1 and 2) was set to 3 hr for all model levels.
Our experiments cover the full month of August 2016.Each month is simulated as 31 individual days, and each simulated day is started at 22:00 UTC (00:00 Central European Summer Time) the previous day, with an integration time of 26 hr.In the post-processing of the statistics, the first 2 hr of each simulation are discarded as spin up, resulting in continuous 1 month long time series.The spinup duration of 2 hr is relatively short, but typically sufficient to create relatively smooth time series without discontinuities, without spending excessive amounts of computing time on the spinup period.
Included in Table 1 are estimates of the computational costs, expressed in Simulated Days Per Day.These estimates are provided for both the central processing unit (CPU) and graphics processing unit (GPU) versions of MicroHH, using either single (four byte) or double (eight byte) precision floating point numbers.All benchmarks were performed on the Snellius supercomputer, using either a single compute node with 128 AMD EPYC 7H12 CPU cores, or a single NVIDIA A100 GPU.As a reference: the most expensive setup-the XL experiment in double precision on the CPU-costs a total of ∼130,000 core hr for a single month.

Post-Processing
We used model output from virtual observation sites (individual LES grid points, without any spatial or time averaging), sampled at a 5-s frequency, to calculate most LES statistics.For these individual columns, radiative transfer was also diagnosed at a 5-s frequency.For the comparison with observations, we averaged the LES output over a 10-min window, in line with the time averaging of the Cabauw observations.Non-averaged individual column statistics were used for the results in Section 4.3, to preserve the highest amount of variability.

Results
The results section starts by providing a first impression of the typical weather variability over Cabauw, and the ability of both the (LS) 2 D-MicroHH combination (hereafter simply referred to as LES) and ERA5 itself, to capture this variability (Section 4.1).Next, we examine the skill of all models by statistically comparing LES and ERA5 for a full 1 month period with the Cabauw observations (Section 4.2).Finally, in Section 4.3 the surface solar irradiance variability is studied.

Characterization of Diurnal Variability
The weather over Cabauw is often highly variable, with different weather types within the time span of a few days.To give a first impression of the ability of ERA5 and LES to capture this variability, we present an overview of a short time period with varying weather conditions controlled by a nearby low pressure area (11-16 August) in Figure 2. In general, both LES and ERA5 perform visually similar (the actual statistics are provided in the next section).
Both models capture most of the variation in the surface solar (SW ↓) and longwave (LW ↓) irradiance (Figures 2a  and 2b), and the surface upward longwave radiation (LW ↑, Figure 2c).Only on clear nights (12 → 13 and 14 → 15 August) both models overestimate LW ↑, indicating that the modeled surface temperatures are too high.
The surface sensible (H) and latent (L v E) heat fluxes (Figures 2d and 2e) are in line with the observations as well.
The Bowen ratio β ≡ L v E/H ≈ 1/3 is typical for Cabauw, and correctly reproduced by ERA5 and the land-surface model in LES.As a result, the 10-m potential temperature, specific humidity, and wind speed (Figures 2f, 2g, and 2h) closely follow the observations.However, both models fail to capture some of the fast observed fluctuations in wind speed.This is perhaps expected in ERA5, with its coarse spatial and temporal resolution, but LES should in theory capture these fluctuations.The surface rain rate (Figure 2i) during the frontal passage was weak at around 0-2 mm per hour, which is reproduced by the microphysics schemes in both IFS and LES, although in both models precipitation stops a bit early.Most of the variability in cloud cover (Figure 2j) is captured by both models.On the night from the 12-13 August, and interesting shift between the cloud cover in ERA5 and LES is visible.The transition from overcast to clear conditions occurs earlier in LES, while the opposite transition from clear to overcast is delayed.This behavior is likely caused by coupling LES to ERA5 only through the domain mean state.Our hypothesis for this behavior is discussed in Section 5. On the last day with shallow cumulus convection, LES is much closer to the observed cloud fraction compared to ERA5.The skill of both LES and ERA5 in predicting the cloud cover is studied in more detail in Section 4.2.2.
Overall, the (LS) 2 D coupling between ERA5 and LES successfully introduces most of the day-to-day variability observed in reality into LES, both for the 5 day period shown in this section, and for the entire month of August 2016 (not shown).In the next section, the skill of LES and ERA5 is further analyzed by statistically comparing both models to the Cabauw observations for a full 1-month period.

Surface and Tower Observations
The statistical analysis will focus on the mean biases of LES and ERA5.To distinguish between biases in different parts of the diurnal cycle, the statistics have been calculated over three hourly periods.For the LES simulations, we focus on the highest resolution (S-HR) and largest domain (XL) experiments.
Figures 3a and 3b show the statistics for the surface shortwave radiation.The mean bias in the incoming radiation in LES is significant at a maximum of 60-75 W m −2 .In contrast, ERA5 has a mean bias of only −1.5 W m −2 .The bias in the outgoing shortwave radiation in LES is almost entirely caused by the bias in the incoming radiation, indicating that the albedo used in LES (0.24%) is accurate.The ERA5 grid point for Cabauw, however, has an albedo of around 17%, resulting in a maximum negative bias in the outgoing shortwave radiation of −40 W m −2 .The large positive bias in the solar irradiance in LES is most likely caused by an underestimation of clouds in certain weather regimes, which is discussed in more detail in the following sections.As a result of the underestimation of clouds, the surface longwave irradiance (Figure 3b) is underestimated in LES.However, compared to the solar irradiance, the biases are much smaller at ∼−5 to −10 W m −2 .The biases in the surface longwave outgoing radiation (Figure 3c) are similar in LES and ERA5.During the night, both models overestimate the outgoing longwave radiation, indicating that the surface temperatures are too high.
The biases in the solar radiation in LES (incoming) and ERA5 (outgoing) are not compensated by the biases in longwave radiation, and therefore both models have a net excess of energy at the surface (Figure 3e).As a result, both the sensible (Figure 3f) and latent (Figure 3g) heat fluxes are overestimated.The timing of the overestimation differs between the two variables: in LES, the sensible heat flux almost perfectly follows the net radiation bias, but the bias in the latent heat flux is delayed.During the night, the sensible and latent heat fluxes in LES are too negative, but these small biases are likely within the measurement uncertainty of eddy-correlation flux measurements at night (e.g., Bosveld et al., 2020;de Roode et al., 2010).ERA5 does a slightly better job at predicting the sensible heat flux, but has a positive evaporation bias throughout the diurnal cycle.
As expected, the LES biases in the 10-m temperature (Figure 3h) and specific humidity (Figure 3i) follow the pattern of the surface sensible and latent heat fluxes.For the 10 m temperature, this results in an overestimation of the diurnal amplitude, as the model is too cold at night and too warm during the day.This is in contrast with ERA5, which underestimates the amplitude of T10 m , which is a well known problem in IFS (Sandu et al., 2013). 10.1029/2023MS003750 9 of 14 The positive evaporation biases in both LES and ERA5 directly translate to relatively small (+0.1-0.2 g kg −1 ) biases in the 10 m specific humidity.Finally, the biases in the 10 m wind speed (Figure 3j) show a similar pattern in both LES and ERA5, with a positive bias during the night, and a negative bias during the day.ERA5 has the largest (negative) bias, especially during the day, and in the LES experiments, the high-resolution experiment reduces the biases at night.
In summary, LES typically performs either similarly or worse than ERA5.For the short-and longwave irradiances, the biases are related to the representation of clouds in LES, which will be studied in more detail in the next section.

Clouds
The Cabauw site has multiple instruments measuring cloud properties, including an LD-40 ceilometer (Bosveld et al., 2020).The latter provides time series of both cloud presence and cloud base at a 15-s frequency, which allows for a direct (statistical) comparison with the high frequency (5 s) single column output from MicroHH.From the observations and LES, a cloud fraction is derived by first calculating a cloud mask (criteria: cloud present in the ceilometer output, positive liquid water path in LES), and second calculating the cloud fraction as the fraction of time that clouds are present in a 1-hr period.From ERA5, we use the total cloud cover variable, which is based on the prognostic sub-grid cloud cover, corrected with a generalized cloud overlap assumption (ECMWF, 2016).
Figure 4 shows the statistics for the LES experiments and ERA5, binned in 10% intervals as a function of the observed cloud fraction.From these results, we can conclude that the influence of the domain size (S vs. XL) and resolution (S vs. S-HR) in LES is small, as all LES experiments show similar biases.For low cloud fractions of ≤30% (0-2 okta), the bias in LES is close to zero, whereas ERA5 has a positive bias of ∼10%-25%.This improvement in LES can be expected as LES explicitly resolves these small clouds.For intermediate cloud fractions of 30%-50% (3-4 okta), LES underestimates the cloud fraction with ∼15%, where ERA5 has a positive bias of ∼10%.For the highest cloud fractions in the range 50%-100% (5-8 okta) both LES and ERA5 underestimate the cloud fraction, although the bias in LES is much larger at ∼30%-40%.
From these results, it is evident that in doubly periodic LESs, with the domain sizes that we used, cloudiness is underestimated for observed cloud fractions larger than 30%.The potential reasons behind these biases, and possible solutions to overcome these biases, are further discussed in Section 6.

Surface Solar Irradiance Variability
The previous sections statistically compared the time averaged model output to observations, in order to reveal potential biases in LES or ERA5.This section will extend that analysis, but with a focus on quantities where LES can potentially add value over ERA5, as it explicitly resolves the interaction between individual clouds and radiation.We should note that our simulations use the two-stream approximation, and therefore exclude 3D-radiative effects (e.g., Kassianov et al., 2011;Villefranque & Hogan, 2021).The impact of these 3D radiative effects has previously been studied with MicroHH and (LS) 2 D by Tijhuis et al. (2023), Veerman et al. (2022), and will not be further examined here.
One variable which is of particular interest is the surface solar irradiance, which at Cabauw (BSRN station) is measured at a one-second frequency.As explained by Mol, van Stratum, et al. (2023), solar irradiance variability and cloud properties are closely related.Therefore, we use this quantity to analyze both the high frequency interactions between clouds and radiation, but also the low frequency interactions, which are likely influenced by the biases in cloudiness (Section 4.2.2).In order to capture irradiance variability at a time scale down to seconds, this section uses the diagnosed radiative fluxes at a 5 s frequency.

Shadow Duration and Size Distributions
The first part of the analysis will focus on the shadow duration and size.In both LES and the observations, we define a shadow as a continuous period where the difference between the surface solar irradiance (global horizontal irradiance [GHI]) and its clear sky (CS) value is more than 50 W m −2 , that is, GHI cs -GHI > 50 W m −2 .In LES the CS irradiance is output by the RTE-RRTMGP model, for the observations we use the CS product from McClear (Lefèvre et al., 2013).The shadow duration is next translated to an approximate shadow length by multiplying it with the 200 m wind speed.This wind speed is likely less than the wind speed in the cloud layer(s), but is the highest observed wind speed available in the Cabauw observations.Figures 5a and 5b shows the probability density functions (PDFs) of the shadow duration and size.All LES experiments are included to examine the influence of domain size and horizontal resolution.Both the observations and LESs exhibit a power-law relationship with a slope around −5/3 across a wide range of scales.The power-law slope is in line with the findings of Mol, van Stratum, et al. (2023), who studied similar shadow properties using 10 years of Cabauw BSRN data.The similarity in the PDF slope indicates that our chosen period is not anomalous in terms of the shadow (or cloud) properties.shadow .At the smallest time and length scales, the influence of the horizontal resolution is obvious.At the highest resolution (S-HR) the variability is close to the observations, but as the grid spacing increases, less variability is captured.This is most clearly visible in the shadow lengths, where the LES simulations at Δ = 50 m and Δ = 100 m seem to drop off at scales below ∼2Δ.At intermediate time and length scales-between 30 s > t > 10 min, equivalent to 200 m < l < 5,000 m-the highest resolution experiment almost perfectly follows the observations, while the lower resolution experiments slightly overestimate the shadow occurrences.At time scales beyond ∼10 min and ∼5 km, all LES experiment show a sharp drop-off.This is to some extent related to the domain size, which limits the maximum cloud size and therefore the shadow durations and sizes.However, in all experiments except S, shadow lengths are missing which could fit into the LES domain.Only at the largest time and length scales do the LES experiment again converge with the observations.As these length scales are much larger than the domain sizes, this is simply the result of an overcast cloud deck which is advected multiple times through the periodic boundary conditions.

Solar Irradiance Power Spectrum
The PDFs from the previous section analyzed the surface solar irradiance as a discrete (on/off) quantity.To get more insight in the full extent of the solar variability, this section studies power spectra of the BSRN observations, ERA5, and all LES experiments.
The power spectral density (PSD) is obtained from time series using fast Fourier transforms.As these PSDs tend to be noisy on the original frequencies, we averaged the PSDs over logarithmically increasing bin sizes (e.g., Schalkwijk et al., 2015).For the LES experiments, we additionally averaged the binned PSDs over the five column locations.
Figure 6 shows the power spectral densities.The energy at time scales (τ) corresponding to the diurnal cycle and daylight period is clearly visible, and captured correctly by the LES experiments.ERA5 is close to the observations, but underestimates the variability slightly at τ > 12 hr.At τ < 12 hr, ERA5 clearly underestimates the variability at all time scales that can potentially be resolved with the hourly output.The LES spectra also drop off at around τ = 7 hr, although less than ERA5, and only converge with the observations around τ = 10-15 min.The larger LES domains (M, L, XL) resolve more energy at time scales between 10 min and 2 hr, but still do not capture the full variability seen in the observations.At time scales below 10 min, LES can capture the full variability down to the 10-s scale, but only in the highest resolution experiment.When the horizontal resolution is increased from 25 to 50 m or 100 m, the resolved variability at time scales less than a minute quickly drops off.
The results from Figures 5 and 6 are in line with the findings from the previous sections.At time scales larger than 10-15 min-corresponding to a spatial scale of ∼10 km-variability in the surface solar irradiance is underestimated as the result of a misrepresentation of large clouds or conditions with large cloud fractions in LES.However, these results also demonstrate the added value of downscaling ERA5 with LES, as such an LES setup is capable of capturing solar irradiance variability from time scales of 10 min down to 10 s.

Discussion
Our validation with observations showed that (LS) 2 D, combined with MicroHH, captures most of the day-to-day variability in the weather for the 1 month period that we selected.However, the results also indicate that such a doubly-periodic LES is not a general purpose tool, suitable for all weather conditions.The validation of the cloud properties and solar irradiance variability against a comprehensive set of surface and remote sensing observations, revealed that these LESs underestimate the cloud fraction for observed cloud fractions >30% (Section 4.2.2), and the size of cloud structures >5 km (Section 4.3.1).This is most likely caused by our approach, where LES is only coupled to the large-scale model through a set of horizontally mean large-scale forcings.Consider for example, a cloud layer that is correctly predicted by the host model, at a height where LES has not developed any spatial variability in temperature or humidity.If this cloud layer is advected into the LES model only through the horizontal mean state, then all grid points at that height in LES will be either saturated or not, resulting in a cloud deck in LES that essentially behaves as an on/off switch.And if the cloud deck in the host model consists of broken clouds, the relative humidity in LES will likely stay below 100%, and LES will not capture any clouds.Similar issues have recently been reported by Jansson et al. (2022) in a superparameterization setup.A setup where LES is coupled to the host model at the lateral boundaries might solve these issues, as this can result in the advection of spatial variability in for example, temperature or humidity into the LES domain.In addition, more detailed observations like for example, vertical profiles of cloud fraction-which were not available for the period that we studied-might be beneficial to better understand these issues in LES.
A second issue that we did not discuss in detail is the formal validity of our experiments in the nocturnal boundary layer (NBL).Resolving sufficient turbulence in a weak to moderately stable boundary layer requires a grid spacing of   (1) m (e.g., Beare et al., 2006).In addition, LES might simply not be the correct tool when studying more strongly stable conditions with little or no turbulence.These resolution requirements and other limitation make it difficult to unite a valid LES setup for nighttime conditions, with a domain large enough for (deep) convection.An LES setup with two-way grid nesting could potentially solve some of these issues, by letting the high-resolution domain feed back the correct mean thermodynamic state to its larger parent domain.Alternatively, for studies where the NBL is of secondary importance, the early morning biases could be ignored, as they are unlikely to influence processes the following day (van Stratum & Stevens, 2015, 2018).

Summary and Conclusion
This paper presented (LS) 2 D: an open-source Python package designed to simplify downscaling ERA5 with doubly-periodic LES.With a number of 1 month long simulations over Cabauw (the Netherlands), consisting of various sensitivity experiments on domain size and resolution with the MicroHH LES model, we demonstrated both the benefits and challenges of using such a realistic but doubly-periodic LES setup.
Overall, the combination of (LS) 2 D and MicroHH manages to capture most of the day-to-day variability in the weather.However-as discussed in the previous section-the model setup has difficulties with capturing conditions with large cloud fractions and/or structures.This results in mean biases in the short-and longwave irradiances, which propagate into other quantities such as the surface heat fluxes.Although this does not seem like a ringing endorsement for our methodology, we believe that our downscaling approach is still beneficial for studying processes that are internally resolved by LES, but poorly resolved or parameterized by the host model.This way, the host model sets the large-scale environment, in which LES resolves the small scale processes.A key example, important for both weather and climate, is shallow convection, which is fully parameterized by most large-scale weather and climate models, but explicitly resolved by LES.For these conditions, our comparison with the Cabauw observations showed that LES improves on its host model when predicting for example, cloud cover.
Being able to resolve shallow convection accurately has an important implication, as it enables the use of this model setup to study the complex interplay between radiative transfer, land-surface processes, turbulent transport, and moist convection.As we have demonstrated, a sufficiently high resolution LES setup manages to capture solar irradiance variability across a wide range of temporal scales, all the way down to a timescale of seconds.This is important for correctly modeling land-atmosphere interactions, but also provides the potential to use LES for forecasting for example, irradiance variability for solar energy applications, for which information on high frequency variability is crucial (Kreuwel et al., 2020).
In summary, the downscaling of ERA5 with (LS) 2 D and LES has proven to be a useful tool to advance our understanding on the interplay between several key atmospheric processes, when applied to meteorological appropriate conditions.

Figure 4 .
Figure 4. Statistics cloud cover: cloud cover statistics for ERA5 and three of the large-eddy simulation experiments, averaged in 10% bins of the observed cloud fraction.

Figure 5 .
Figure 5. Cloud shadow duration and size: Normalized probability density functions (PDFs) of the cloud shadow duration (a) and size (b) from the Baseline Surface Radiation Network (BSRN) observation and all large-eddy simulation (LES) experiments.In the bottom row, the mean slope-typically consisting of a x −5/3 slope(Mol, van Stratum, et al., 2023)-has been divided out to improve readability.

Figure 6 .
Figure 6.Solar irradiance spectra: Power spectra of the surface solar irradiance, comparing ERA5 and large-eddy simulation with the 1 s Baseline Surface Radiation Network (BSRN) observations.The gray shading is the original (not averaged) power spectral density (PSD) of the BSRN data.
Name Horizontal size (km) Resolution (m) SDPD CPU SDPD GPU Note.The last two columns provide an estimate of the average simulated days per day (SDPD) for the setup used in this paper, on 128 AMD EPYC CPU cores, or a single NVIDIA A100 GPU.The two numbers provided for each setup indicate the SDPD when running MicroHH using double (8 byte)/ single (4 byte) precision floating point numbers.