# The surface-boundary layer connection across spatial scales of irrigation-driven thermal heterogeneity: An integrated data and modeling study of the LIAISE field campaign

Mary Rose Mangan, Oscar Hartogensis, Aaron Boone, Oliver Branch, Guylaine Canut, Joan Cuxart, Hugo J de Boer, Michel Le Page, Daniel Martínez-Villagrasa, Josep Ramon Miró, et al.

► To cite this version:

Mary Rose Mangan, Oscar Hartogensis, Aaron Boone, Oliver Branch, Guylaine Canut, et al. The surface-boundary layer connection across spatial scales of irrigation-driven thermal heterogeneity: An integrated data and modeling study of the LIAISE field campaign. Agricultural and Forest Meteorology, 2023, 335, pp.109452. 10.1016/j.agrformet.2023.109452. hal-04271460v2

HAL Id: hal-04271460

https://hal.science/hal-04271460v2

Submitted on 6 Nov 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /></figure>

# The surface-boundary layer connection across spatial scales of irrigation-driven thermal heterogeneity: An integrated data and modeling study of the LIAISE field campaign

Mary Rose Mangan<sup>a,*</sup>, Oscar Hartogensis<sup>a</sup>, Aaron Boone<sup>b</sup>, Oliver Branch<sup>c</sup>, Guylaine Canut<sup>d</sup>, Joan Cuxart<sup>e</sup>, Hugo J. de Boer<sup>f</sup>, Michel Le Page<sup>g</sup>, Daniel Martínez-Villagrasa<sup>e</sup>, Josep Ramon Miró<sup>h</sup>, Jeremy Price<sup>i</sup>, Jordi Vilà-Guerau de Arellano<sup>a</sup>

## ARTICLE INFO

**Keywords:**

Surface heterogeneity

Scaling surface fluxes

Boundary-layer dynamics

Heat & moisture advection

## ABSTRACT

Irrigation in semi-arid regions induces thermal heterogeneity across a range of spatial scales that impacts the partitioning of energy at the surface, the development of the atmospheric boundary layer, and the bidirectional interactions between the atmosphere and the surface. In this analysis, we use data from the Land Surface Interactions with the Atmosphere in the Iberian Semi-Arid Environment (LIAISE) experiment combined with a coupled land-atmosphere model to understand the role of the scales of irrigation-induced, thermal heterogeneity on the surface fluxes and consequently, the development of the diurnal convective boundary layer. The surface heterogeneity is characterized by Bowen ratios that range from ~0.01 in the irrigated areas to ~30 in the non-irrigated areas; however, the observed boundary-layers dynamics in both locations are similar. In this analysis, we address the questions of how the surface fluxes impact the development of the boundary-layer dynamics and how the boundary layer influences the diurnal cycle of surface fluxes. To interpret the observations, we introduce a heterogeneity scaling scheme where length scales range from local scale (~100 m) to regional scale (~10 km) to investigate the role of scale on surface representation in numerical models and to address the discrepancy between surface observations and their representation in weather and climate models.

We find that at the surface, both the available energy and its partitioning depend on spatial scale. The observed boundary-layer properties can be explained through the composite of surface fluxes at the regional scale. Surface fluxes at the local scales are unable to replicate the observed boundary layer — even when including large-scale contributions. We find that non-local boundary layer processes like advection are important for partitioning energy at the local scale. We explore the connection between surface fluxes and the development of the boundary layer and the potential non-local effects on boundary-layer development.

## 1. Introduction

The Earth's surface is highly heterogeneous over a variety of spatial scales, which presents a challenge for interpreting the physical processes that govern both the partitioning of energy at the surface and the development of the atmospheric boundary layer. Historically,

land-atmosphere interactions have been studied by observing the one-directional impact of the surface on the atmosphere, or vice versa. Recently, there has been an effort to advance observations of land-atmosphere interactions by co-locating surface flux stations with continuous boundary-layer observations (Beamesderfer et al., 2022; Helbig

<footnote><sup>a</sup> Meteorology and Air Quality Group, Wageningen University and Research, Wageningen, The Netherlands</footnote>

<footnote><sup>b</sup> Météo-France/CNRS, Toulouse, France</footnote>

<footnote><sup>c</sup> Institute of Physics and Meteorology, University of Hohenheim, Stuttgart, Germany</footnote>

<footnote><sup>d</sup> Météo-France, Toulouse, France</footnote>

<footnote><sup>e</sup> Department of Physics, University of the Balearic Islands, Palma, Spain</footnote>

<footnote><sup>f</sup> Copernicus Institute of Sustainable Development, Environmental Sciences, Universiteit Utrecht, Utrecht, The Netherlands</footnote>

<footnote><sup>g</sup> CESBIO, Centre d'Etudes Spatiales de la Biosphère, Univ. de Toulouse, CNRS, CNES, IRD, UPS, INRAE, Toulouse, France</footnote>

<footnote><sup>h</sup> Meteorological Service of Catalonia, Barcelona, Spain</footnote>

<footnote><sup>i</sup> Met Office, Exeter, United Kingdom</footnote>

<footnote>* Corresponding author.
E-mail address: maryrose.mangan@wur.nl (M.R. Mangan).</footnote>
et al., 2021). By coupling surface flux observations with boundary-layer observations, we can disentangle processes that govern land-atmosphere dynamics. Although Beamesderfer et al. (2022) and Helbig et al. (2021) call for long-term measurement sites for quantifying land-atmosphere interactions, intensive and comprehensive experimental campaigns provide insight on the bidirectional land-atmosphere feedbacks across spatial scales. Therefore, we use data from the Land Surface Interactions with the Atmosphere over the Iberian Semi-Arid Environment (LIAISE) campaign (Boone et al., 2021) to investigate both the bidirectional land-atmosphere interactions in a thermally heterogeneous environment.

In this analysis, we address how the impacts of heterogeneity are felt across spatial scales. We are motivated both by process understanding of the bidirectional feedbacks between the surface and the boundary-layer and examining the usefulness of applying an idealized mixed-layer model to a heterogeneous area. Therefore, we aim to address the following questions about process understanding for the LIAISE campaign: (1) How does the thermally heterogeneous surface control the boundary-layer dynamics across spatial scales, and (2) What is the impact of the boundary-layer dynamics on the diurnal variability of the observed surface fluxes? Although the study is motivated by physical process understanding of land-atmosphere interactions, the methodology is motivated by the need to address the mismatch between how observations and models represent the physical world. We introduce a site-specific framework for modeling the scales of surface heterogeneity, which could be applied systematically to other regions to investigate the role of the scale of surface heterogeneities on the ABL. In idealized boundary-layer models low resolution limits the representation of the land surface and the atmospheric boundary layer. By comparing the impacts of the resolution of the heterogeneities between measurements and models, we can evaluate more comprehensively the parametrizations of land surface models in heterogeneous areas.

Because land-surface models consider subgrid scale surface heterogeneity by making a composite of either the land surface properties within the grid cell (parameter aggregation) or the fluxes above the surface from separate non-interacting tiles (flux aggregation), there is a discrepancy between the land surface modeling community's needs and the measurements of surface fluxes. Land-atmosphere exchange has been measured by eddy-covariance systems for decades (Swinbank, 1951; Baldocchi et al., 1988; Helbig et al., 2021), and it remains one of the most common observational sources for verification of land surface models. Although through networks like FluxNet (Baldocchi et al., 2001) and ICOS (Kadygrov et al., 2015), there is high global coverage of surface flux measurements, each eddy-covariance system measures fluxes with a relatively small footprint. This means that each tower measures the fluxes from a small area within the grid cell of the global scale model. Moreover, eddy-covariance towers typically measure below the first model level of a weather or climate model. For these reasons, measured surface fluxes do not necessarily represent the same flux that land surface models aim to reproduce; however, these eddy-covariance systems are commonly used as model verification for the land surface models. To address the spatial disparity between the measurements and models, we created maps of surface fluxes and properties (Section 4.1) to replicate the surface and its fluxes at a regional spatial scale.

There have been multiple studies with large-eddy simulation and with experimental campaigns to quantify the impacts of surface heterogeneities on the atmospheric boundary layer. Large-eddy simulations have been used to study heterogeneity scaling in the convective boundary layer (Patton et al., 2005; Shen and Leclerc, 1994, 1995; van Heerwaarden et al., 2014), quantifying regional scale impacts of surface heterogeneity including secondary circulations (Raasch and Harbusch, 2001), and the development of an internal boundary layer (Bou-Zeid et al., 2004). In addition to large-eddy simulation, field campaigns which include observations that span scales of heterogeneity have

been used to study heterogeneous land surfaces, including the recent CHEESEHEAD experiment (Butterworth et al., 2021), the BLLAST experiment (Lothon et al., 2014), and the CloudRoots Campaign (Vilà-Guerau de Arellano et al., 2020). Nonetheless, it remains difficult for field measurements to span the same scales as the surface heterogeneity; therefore, modeling experiments have been used to connect the local scale measurements in heterogeneous regions to their mesoscale impacts on the boundary layer.

To support the interpretation of multi-scale observations from the LIAISE experiment, we used numerical experiments performed by a conceptual mixed-layer model (CLASS) (Vilà-Guerau de Arellano et al., 2015). In so doing, we could study the dominant processes at each scale and study whether new and emerging processes become dominant at different scales. We have deliberately chosen a conceptual, coupled land-atmosphere model for this study to reduce the complexities of the LIAISE domain so that we investigate the essential processes that govern land-atmosphere interactions in the thermally heterogeneous environment. Furthermore, it allowed us to replicate systematically how a land surface model handles subgrid scale surface heterogeneity at different spatial resolutions.

The paper is structured as follows. In Section 2, we will introduce the LIAISE campaign as it is suited for our multi-scale analysis due to the strong contrast of surface properties. We subsequently propose a spatial scaling scheme to represent local, landscape and regional processes (Section 3). In this way, we can replicate how a land surface model handles heterogeneous land surfaces depending on scales, and how these scales translate to a model grid cell. Next, we will describe our methods to upscale our measurements to represent each model scale, and we will describe how we replicate spatial heterogeneity using a zeroth-order model (Section 4). We will show the results from the conceptual model with the measurements and the ERA5 reanalysis model (Section 5). Finally, we will bring the scales together to discuss the bidirectional impacts between boundary-layer dynamics and surface fluxes in a heterogeneous environment (Section 6).

## 2. Experimental and modeling approach

In this study, we used data from the LIAISE campaign and a mixed-layer coupled land-atmosphere model to investigate surface fluxes at different scales of heterogeneity. Results from observations and the mixed-layer model were compared to the 0.25° resolution European Centre for Medium-Range Weather Forecasts' (ECMWF) reanalysis product, ERA5 (Hersbach et al., 2020). We use ERA5 as a prototypical ABL for the modeling efforts because its land surface model partitions energy in a similar way as the CLASS model. In this way, we can make a direct comparison between our modeling study and how global models handle surface heterogeneity, especially at the regional scale where the horizontal domain matches that of ERA5.

### 2.1. Site description: LIAISE domain

The LIAISE field campaign took place in the Ebro River Valley in Catalonia, Spain during May through October 2021. The LIAISE experiment was designed to improve understanding of land-atmosphere interactions in a thermally heterogeneous environment (Boone et al., 2021, 2019). Atmospheric flow in the LIAISE domain is complex for three reasons: (1) it is located in a nearly closed valley in the Ebro River Basin, (2) the Mediterranean Sea is located approximately 70 km to the southeast which induces a land-sea circulation, and (3) within the LIAISE domain, there is a strong thermal heterogeneity (Fig. 1). Fig. 1 indicates the land use classes in the LIAISE study domain. The western part of the domain contains annual, irrigated crops like maize and alfalfa and irrigated fruit trees. The eastern part of the domain contains rainfed vegetation consisting mainly of winter cereal crops, vineyards and rainfed orchards. The box inside of right plot in Fig. 1 represents the extent of the ERA5 grid cell in the LIAISE domain.
<figure><img src="image_4.png" /><figcaption>**Fig. 1.** The LIAISE campaign occurred in the Ebro River Valley in Northeast Spain in the summer of 2021. The left panel shows the aerial view of the LIAISE domain (image from GoogleEarth). The inset in the left picture is the location of the experimental site in the Iberian Peninsula. The right panel is the 100 m land cover map from CORINE Land Cover Product (Buttner, 2014). The box in the land use map is the extent of the ERA5 grid cell in the LIAISE domain.</figcaption></figure>

For the purposes of our study, we were most interested in the ther-
mal heterogeneity within the LIAISE domain that arises from irrigation
applied in one area with a length scale on the order of 10–100 km.
There were measurements in a number of locations near the sharp
boundary between the irrigated and the semi-arid areas. In the LIAISE
campaign, there were three “supersite” locations: Mollerussa (mixed
orchards), La Cendrosa (alfalfa field) and Els Plans (fallow field with
rainfed, natural vegetation) (Fig. 1). In addition to the three supersites,
there was a network of nine surface energy budget (SEB) stations
in each of the predominant crop types in the LIAISE domain. All
measurement locations were located within 10 km of the wet-dry
boundary.

## 2.2. Experimental data: LIAISE campaign

The LIAISE field campaign included a short-term observation period from 15 July through 30 July 2021. At Els Plans (rainfed fallow) and La Cendrosa (irrigated alfalfa), measurements spanned scales from the leaf level to the boundary-layer level. Data collected at each supersite include ecophysiology measurements, SEB stations and boundary-layer measurements including 50 m towers, tethered balloons and hourly radiosondes during the daytime. Across the entire LIAISE domain, there were aircraft measurements of the boundary layer (e.g., turbulent fluxes) and the surface (e.g., solar induced fluorescence and soil moisture). Additionally, there was a network of nine SEB stations in the predominant vegetation covers in the LIAISE domain. For more information about the extent of the measurements from the LIAISE experiment, see Boone et al. (2021).

In Table 1, we show the overview of the data used in this study. At the surface, we used data from the network of the nine SEB stations and the lowest level of flux measurements at Els Plans (2 m) and La Cendrosa (3 m). At the surface stations, there were eddy-covariance systems, four-stream radiometers and ground heat flux measurements. The sampling times for the systems varied from 10–20 Hz for the eddy-covariance to 60 s for the radiative fluxes, however, all SEB data was block averaged to 30 min. The fetch during convective conditions for the surface stations varied from 50 to 150 m downwind, depending on the measurement and crop heights, which fell within the field boundaries in the predominant wind direction for all fields. For boundary-layer measurements, we used hourly radiosondes where we derived the boundary layer height with the parcel method (Kaimal

and Finnigan, 1994), the mixed-layer mean potential temperature and specific humidity and the entrainment jumps in the potential temperature and specific humidity at the top of the boundary layer (Fig. 1 Conzemius and Fedorovich, 2007; Driedonks and Tennekes, 1984). The mixed-layer scalar means were calculated with the average of the radiosonde values of potential temperature and specific humidity below the mixed-layer height. The entrainment jumps in scalars were found by first determining the lapse rate of the scalar in the free atmosphere. Using that lapse rate, we extrapolated it to find the expected value at the top of the mixed-layer. The jump is defined as the difference between the value at the mixed-layer top and the mean mixed layer value. Finally, we used data from the automated weather stations, which will be explained in more detail in Section 4.1.

We averaged a three day period from 20 July through 22 July to represent a “composite day” of the LIAISE experiment. Both the synoptic situation and the surface fluxes were similar over this period. There was anticyclonic flow at the surface with a thermal low building to the west of the study site. Over the course of this three day period, local conditions in the LIAISE domain were slowly getting hotter and drier. For these reasons, we created a composite LIAISE composite day by taking the average of the surface fluxes and boundary-layer properties over 20–22 July for use in the mixed-layer model. By using a composite day, we were able to better capture a typical situation for the LIAISE domain instead of modeling a situation that was heavily infused by random extremes due to measurement limitations or non-local events. Furthermore, creating a composite day allowed us to have a robust procedure for handling missing data.

At both La Cendrosa and Els Plans, there were SEB stations which directly measure the components of the energy budget: net radiation (R<sub>n</sub>), sensible heat flux (H), latent heat flux (LE), ground heat flux (G). The average energy budget non-closure from 7 to 17 UTC at these sites varied from 33% in La Cendrosa to 15% in Els Plans. We forced energy budget closure using the method from Twine et al. (2000), which preserves the observed Bowen ratio, for comparison with the CLASS model. Additionally, there were measurements of three-dimensional wind, temperature, humidity, air pressure, soil temperature and soil moisture. Boundary-layer measurements at both sites include a 50 m tower with three-dimensional wind, temperature and moisture measurements at 10, 25 and 50 m above ground level. Moreover, at both sites, there were hourly radiosondes launched starting at 06:00 LT
<table><caption>**Table 1**
An overview of the data from the LIAISE field campaign used in this study including the nominal sampling frequency and the averaging time.</caption><thead><tr><th></th><th>Instrument</th><th>Derived variable</th><th>Locations</th><th>Sampling frequency</th><th>Averaging time</th></tr></thead><tbody><tr><td><strong>Surface - Fluxmaps</strong></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Sonic Anemometers</td><td>Campbell Sci. CSAT3<br>Gill WindMaster<br>R.M. Young 81000</td><td>H</td><td>9</td><td>10-20 Hz</td><td>30 min</td></tr><tr><td>Gas Analyzers</td><td>Campbell Sci. Irgason<br>EC150<br>LICOR-7500<br>Krypton KH2O</td><td>LE</td><td>9</td><td>10-20 Hz</td><td>30 min</td></tr><tr><td>Radiation</td><td>Hukseflux NR-01<br>Kipp & Zonnen CNR4</td><td>R<sub>n</sub></td><td>9</td><td>1-60 s</td><td>30 min</td></tr><tr><td>Surface Heat Flux</td><td>Hukseflux HFP01</td><td>G</td><td>9</td><td>1-60 s</td><td>30 min</td></tr><tr><td><strong>Boundary layer</strong></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Radiosondes</td><td>Vaisala RS92-SGP</td><td>h, θ, q</td><td>2</td><td>1 h</td><td></td></tr><tr><td><strong>Meteorological stations</strong></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>AWS</td><td>Vaisala HMP155<br>R.M. Young 05103</td><td>T<sub>2m</sub>, q<sub>2m</sub> U<sub>10m</sub></td><td>15</td><td>1 min</td><td>30 min</td></tr></tbody></table>

(UTC+2) through 19:00 LT during all three of the composite days. Ecophysiological measurements, including stomatal conductance, CO₂ and light-response curves, and leaf-area index were taken at La Cendrosa on 17 and 19 July.

To demonstrate the contrast between the irrigated and rainfed areas, we show Fig. 2 which includes the diurnal cycle of the surface energy balances for both Els Plans and La Cendrosa during the composite day and the radiosondes launched in both sites at 15:00 LT (UTC+2). Hereafter, all of the time series will be presented in UTC instead of local time (UTC+2) because in the study domain, solar noon is approximately 12 UTC. In Fig. 2a and b, the surface energy balances for La Cendrosa and Els Plans are displayed. The observed albedo between the two sites were comparable ($\alpha$ ≈ 0.23), but the surface temperature, and therefore, outgoing longwave radiation was higher for Els Plans than for La Cendrosa, so there was more available energy at La Cendrosa than Els Plans. Much of the available energy at La Cendrosa was partitioned into latent heat flux, and the observed sensible heat flux became negative after 13:00 UTC which indicates stable conditions in the surface layer. Conversely, at Els Plans, much of the available energy was partitioned into sensible heat flux, while latent heat flux remained near zero the entire day. In Fig. 2c and d, the potential temperature and specific humidity measured via radiosondes were averaged together from 13:00 UTC during the three days. At La Cendrosa (green), the surface was cooler and wetter than at Els Plans (yellow), but Els Plans had a better defined mixed-layer. The boundary-layer heights were approximately the same regardless of land cover; however, we observe a thin stable boundary-layer in the lowest ~200 m at La Cendrosa. Furthermore, we observe that above 500 m, there is dry air intrusion in the La Cendrosa radiosonde above the locally wet surface layer, while the Els Plans radiosonde shows a well-mixed profile. This supports the idea that there is a blending height above which the influence of the landscape scale diminishes, as suggested in Fig. 3.

## 2.3. Coupled land-atmosphere model

We used the atmospheric mixed layer, slab model Chemistry Land-surface Atmosphere Soil Slab model (CLASS) (https://classmodel.github.io/; Vilà-Guerau de Arellano et al. (2015)) to help us interpret the LIAISE composite day. By using a coupled land-atmosphere mixed-layer model, the complexities of the LIAISE campaign can be simplified such that topography and advection were prescribed. To further disentangle the situation, we were interested mainly in the daytime before the sea breeze arrives in the LIAISE domain. The sea breeze arrives between 14:30 and 15:45 UTC during the composite day. Although we show the results of the model until 18:00 UTC, our discussion is focused on processes that occur before the arrival of the sea breeze. By using this model, we could delineate the bidirectional impact of the

land surface on the atmosphere and the atmosphere on the land-surface fluxes.

### Land Surface Representation

The surface layer model was based on Monin-Obukhov similarity theory. At the land surface, the vegetation was represented by a big-leaf model, where the partitioning of energy was done with the Penman-Monteith equation. We used the Jarvis-Stewart land-surface model to parameterize the vegetation surface resistance (van Heerwaarden et al., 2009). The evolution of both soil moisture and temperature used the force-restore method (Noilhan and Mahfouf, 1996; Noilhan and Planton, 1989).

The land surface was represented in the CLASS model with static variables — like leaf area index (LAI), vegetative fraction ($c_{veg}$), and soil moisture. Most of the variables to describe the land surface are consistent between the scales. For example, we assumed that the soil properties including wilting point and field capacity are relatively constant throughout the domain. Soil properties were estimated using soil maps from Institut Cartogràfic i Geològic de Catalunya (2009-2020). Soil thermal conductivity was used as a tuning parameter for capturing the ground heat flux for each scale. Initial surface temperature were taken from observations from the early morning. The variables to describe the land surface were determined by creating surface maps which are described in more detail in Section 4.

The CLASS model had two soil layers: the top one which responded to the atmosphere dynamically, and the bottom one which represented the root zone. In the bottom, reservoir soil layer, the temperature and moisture were constant throughout the day. There was slow diffusion of heat and moisture from the reservoir soil layer to the top soil layer. In this analysis, the reservoir soil layer had soil moisture that is set to field capacity for all scales, and the top soil layer is set to observed values at each scale. This means that the top soil layer controls the surface evaporation and the partitioning of energy at the surface, and the reservoir soil layer controls plant transpiration. Thus, the plants modeled with CLASS are never water stressed. In the dry scales, where plants were likely water stressed, this effect was taken into account by altering the vegetative cover (Section 4.3).

### Mixed Layer

In the CLASS model, the mixed layer prognostic model was described by Tennekes (1973). Initial conditions for mixed-layer characteristics (e.g. potential temperature, specific humidity, mixed-layer height) were prescribed based on the mean radiosondes from the composite day (see Section 4.2). The model calculated the time-evolving mixed-layer properties every five seconds. For example, the potential temperature ($\theta$) was calculated with

$$ \frac{\delta\theta}{\delta t} = \frac{\overline{w'\theta' - w'\theta'^{e}}}{h} + ad\upsilon_{\theta} \quad (1) $$
<figure><img src="image_5.png" /><figcaption>**Fig. 2.** LIAISE Composite day: 20–22 July 2021. (a) The surface energy balance for La Cendrosa averaged across composite day, and (b) The surface energy balance for Els Plans averaged across the composite day. (c) Potential temperature from a radiosonde launched at 13:00. (d) Specific humidity measured from the radiosondes launched at 13:00 during the composite day. In all subsequent figures, the solid line is the mean for all composite day and the shaded area represents one standard deviation from the mean.</figcaption></figure>

where $\overline{w'\theta'}$ is the surface heat flux, $\overline{w'\theta'_e}$ is the entrainment heat flux and $adv_\theta$ is the temperature advection, and $h$ is the mixed-layer depth. The surface fluxes impacted the heating (and moistening) of the mixed-layer, as well as the growth of the mixed-layer. In turn, the mixed-layer properties controlled the gradient between the surface and the atmosphere which impacts the surface fluxes.

The LIAISE domain was characterized by its different scales of heterogeneity, so a modeling scheme was developed to reflect these scales (see Section 3). The CLASS model is single column slab model with added advection, so each spatial scale is represented by one vertical column where the surface conditions are changed to reflect the composite land cover at that scale. The CLASS model incorporated larger scale forcing by including advection terms for momentum, temperature and moisture (Vilà-Guerau de Arellano et al., 2015, 2020). In order to better capture the observed boundary layer in the CLASS model runs, large scale advection of temperature and moisture were included based on a network of Automated Weather Stations (AWS) operated by the Servei Meteorològic de Catalunya. More information on this procedure can be found in Section 4.2.

### 3. Spatial scaling scheme

In the LIAISE domain, surface heterogeneities occur across a range of spatial scales. We have defined three scales of heterogeneities in order to compare the results of the mixed-layer model with the surface fluxes that were constructed from local measurements during the LIAISE campaign. In this way, we can both quantify the impact of scale on modeled and measured fluxes and evaluate how the surface fluxes impact the development of the boundary layer across each scale. We have defined three scales: regional (~10 km), landscape (~1 km) and local (~100 m). The regional scale consists of wet and dry landscape scales, and within the landscape scale, there are alfalfa and fallow local scales to represent individual fields. Existing methods to characterize heterogeneities focus on size of the heterogeneity (van Heerwaarden et al., 2014) or the structure of the heterogeneity (Bou-Zeid et al., 2020); however, much of the scaling research has been focused on one

<figure><img src="image_6.png" /><figcaption>**Fig. 3.** The schematic representation of the scales in the LIAISE and how they interact with each other with height in the atmospheric boundary layer. The local scale is impacts the surface layer. Above the surface layer, there is a blending zone which is most impacted by the landscape scale. Above the blending height, the regional scale controls the boundary layer. The LIAISE campaign has instrumentation to measure with height in the boundary layer: flux towers measure in the surface layer, tethered balloons measure in the blending zone, and radiosondes and aircraft measure above the mixed layer.</figcaption></figure>

scale or type of heterogeneity — not how the type of heterogeneity depends on the scale.

The largest scale in the LIAISE domain is the regional scale, which has a length scale on the order of 10s km. On this scale, the heterogeneity is from a large, single wet patch surrounded by dry land. This scale represents the extent of the ERA5 grid cell that is shown in Fig. 1. This is classified as a Type III — large individual patch class of heterogeneity by Bou-Zeid et al. (2020). Within the regional scale, there are landscape scales within both the wet and the dry patches each with a length scale on the order of 1 km. Within this scale, there are heterogeneities
between fields — in both moisture and surface roughness. This type of heterogeneity is considered Type IV unstructured heterogeneity as defined by Bou-Zeid et al. (2020). The smallest scale is the local scale, which is on the order of 100 m. It represents individual fields: La Cendrosa alfalfa for the irrigated (alfalfa) local scale and Els Plans natural vegetation for the dry (fallow) local scale. We assume that the local scale is a statistically homogeneous area. We propose our scaling scheme to be complementary to the one proposed by van Heerwaarden et al. (2014) in which they defined heterogeneity scaling as macroscale, mesoscale or microscale depending on the size of the heterogeneous patch compared to the domain. The regional scale is a macroscale heterogeneity, the landscape scales are mesoscale heterogeneities and the local scales are microscale heterogeneities.

Fig. 3 is an abstract representation of the horizontal scales and how they interact with each other in the boundary layer. The impact of the local scales are felt near the surface. The landscape scales, which capture the heterogeneity between fields in both the irrigated and rainfed areas, are felt above the surface layer where impacts from local fields are blended together. Above the blending height near the top of the boundary layer, the impacts of the regional scale are felt in the atmosphere. Depending on the type of measurements, we observe fluxes that are representative of different scales. These are the first-order effects: each horizontal scale feeds a vertical scale in the atmosphere. The second-order effect is that the vertical scales interact with each other in the boundary layer and communicate downwards to impact the surface fluxes.

In order to replicate the spatial scales using the single pixel mixed-layer model, we defined a composite land surface for each scale based on a spatial average of surface characteristics. This is analogous to the parameter aggregation scheme used in land surface models. We imposed the calculated advection term in the local scale cases to replicate the boundary layer. We will discuss the observational data preparation that was necessary to replicate these scales both in the CLASS model and with observational data in the next section.

## 4. Data integration and upscaling

In order to evaluate the behavior of the CLASS model representative of different spatial scales, we needed to upscale the surface fluxes and surface properties, measured at local scales, to match the regional and landscape scales. The purpose of this data integration is to combine networks of measurements to estimate a spatial distribution of surface properties and fluxes. By upscaling the data, we could directly compare measurements to both the CLASS and ERA5 models. In this section, we introduce the mixed layer properties and surface parameters that constrain the CLASS model. We also introduce the data products that are used to verify the model.

In Section 4.1, we describe the data used to constrain and verify the land surface representation in CLASS. We describe the upscaling of surface flux data to represent surface fluxes across the regional and landscape scales. By representing the land surface as a composite of realistic surface, we are using parameter aggregation to describe the heterogeneity, but we are verifying the product with flux maps, which represents the flux aggregate approach to describe the heterogeneity. In Section 4.2, we describe the data approach for constraining the mixed layer. We introduce an approach to calculate boundary-layer advection from a network of automatic weather stations so that we can replicate the boundary layer at the local scales. In Section 4.3, we summarize the data approach and outline the experimental design based on spatial scale.

### 4.1. Land surface data: Surface maps

To create the land surface differences between the scales, we adjusted vegetation properties and the top layer of soil moisture based on observations. The spatially aggregated land surface constrains the

CLASS model. We used the diurnal cycle of observed surface energy budget components to verify the model. In order to prepare the input land surface parameters and the surface energy budget data for model verification for the regional and landscape scales, we created surface flux and representation maps. Using these maps, we created a dataset that represents the mean surface properties at the landscape and regional scales.

The LIAISE campaign included a network of SEB stations across the predominant crop types in the area during the summertime. In total, there were nine different SEB stations used for this analysis, including stations located in alfalfa, maize and fallow fields, fruit, and nut orchards, and a vineyard. All SEB stations were processed uniformly using EddyPro (LI-COR Biosciences; Fratini and Mauder (2014)). We gap-filled the data using available data from the other composite days. In addition to uniform processing of the network of eddy-covariance stations, a 100 m resolution land cover map from the Sistema de Información Geográfica de Parcelas Agrícolas (SIGPAC) was used to provide the spatial extent of the crops. The SIGPAC crop cover map was supplemented by the 100 m CORINE Land Cover product for the urban areas and the water bodies (Buttner, 2014) (see Fig. 1). Crop types for SIGPAC were reclassified to match those of the predominant vegetation types from the SEB station network. Measured fluxes were applied to the corresponding crop type within the LIAISE region to create flux maps. See Appendix A for more information on the reclassification method for the SIGPAC crop cover maps. Although the 2020 crop cover map was used, we assume that the total distribution of crop types is similar between 2020 and 2021.

At the regional scale, 21.3% of the land area was urban. Because there were no SEB stations in an urban area in the LIAISE campaign, we modeled the expected surface fluxes using the method described in Román-Cascón et al. (2021). We assumed that the urban land had a Bowen ratio ($\beta$) of 5, an emissivity of 0.92 and an albedo of 0.15 (Grimmond and Oke, 1999; Lemonsu et al., 2004). G was assumed to be 10% of $R_n$ during daytime. Fluxes of sensible and latent heat flux were solved by iteratively updating surface temperature using measured incoming shortwave radiation, air temperature, humidity observed at the grass SEB site. The same procedure was used for the fluxes over a water surface (which accounts for 0.3% of the regional surface area). In that case, the assumed $\beta$ was 0.1, an emissivity of 0.98 and an albedo of 0.08. G was assumed to be 30% of $R_n$ during daytime.

Fig. 4 displays an example of the flux maps for latent heat flux for the LIAISE domain on 20 July 2021 at 14:00 UTC. The entire map represents the LIAISE regional scale. The dashed line is the separation between the wet and dry landscape scales. The regional fluxes were calculated as a spatial average of the fluxes in the entire domain, and the wet and dry landscape fluxes were calculated as a spatial average from the area inside and outside of the dashed line in Fig. 4 respectively. Consequently, we were able to derive time series of regional, dry and wet landscape scale energy budget components.

In the wet landscape, the latent heat flux was as high as 400 W m⁻², but this depends on the crop cover at the local scale. Within the wet landscape scale, there were urban areas near Mollerussa which provide relatively high sensible heat flux compared to the rest of the landscape. In the dry landscape, the fields that were fallowed in the summer, like Els Plans, had high Bowen ratios. The orchards in the dry landscape provided more latent heat flux compared to the rest of the landscape. There were maize fields in the north of the LIAISE region, which moistened the dry landscape scale. Within the landscape scales, it is evident that there can be strong differences between fields.

In addition to using the land cover map to create flux maps for the components of the surface energy balance, we have created maps based on the ecopyschological measurements to estimate vegetative properties at the landscape and regional levels, including LAI, $c_{veg}$, and stomatal conductance, which are used to prescribe the surface at the regional and landscape scales. See Appendix A for more details about the flux and surface map products.
<figure><img src="image_7.png" /><figcaption>Fig. 4. An example flux map for latent heat flux for 20 July 2021 at 14:00 UTC. The dashed line represents the boundary between the wet and dry landscape scales. The total grid cell represents the LIAISE regional scale as well as the ERA5 grid (0.25° × 0.25°) cell for the LIAISE domain. The dashed line represents the separation between the wet (west) and dry (east) landscape scales. The local scales are represented by the points at La Cendrosa (alfalfa local) and Els Plans (fallow local) respectively. The box in the top right corner has the spatial standard deviation for each spatial scale.</figcaption></figure>

## 4.2. Mixed layer data: Model initialization & advection

Unlike the input parameters for the land surface representation of the CLASS model, the mixed-layer properties change during the day. We inputted mixed-layer characteristics (e.g. mean mixed-layer potential temperature, mean mixed-layer specific humidity and mixed-layer height) at the start of the model run (6:00 UTC). During the LIAISE composite day, there were hourly radiosondes between 04:00 UTC and 17:00 UTC that were launched from both La Cendrosa and Els Plans. The radiosondes launched at 06:00 UTC at both sites were averaged together to calculate the initial mixed-layer properties. After providing the initial conditions, the CLASS model calculated the mixed-layer characteristics, including mean potential temperature, mean specific humidity and mixed-layer height. The dynamically changing mixed-layer characteristics were verified using hourly radiosondes from both the wet and the dry areas. We calculated the mixed-layer height with the parcel method approach from the radiosondes (Kaimal and Finnigan, 1994).

For all scales, the initial boundary-layer profile from the soundings at 06:00 UTC was used. The CLASS model is insensitive to the initial conditions in the range the measurement uncertainties after the first hour. The only difference in the mixed-layer input between scales was the diurnal advection terms. At the regional and landscape scales, there was weak advection that corresponds to synoptically driven hot, dry westerlies to the region during the day. In the late afternoon (after 15:00 UTC), we prescribed slightly cool and moist advection to represent the sea breeze. At the regional and landscape scales, we assumed that the boundary-layer development was primarily formed within the region; however, at the local scale, this weak synoptic advection was insufficient to describe the observed boundary layer. Without advection, the surface fluxes from the alfalfa local scale yielded

a boundary-layer height of ~800 m, while the surface fluxes from the fallow local scale yielded a boundary-layer height ~1700 m. This implies that the boundary layer is not formed locally over either scale. Therefore, we have calculated advection of moisture and heat for the alfalfa and fallow local scale cases using an AWS network operated by the Servei Meteorològic de Catalunya.

In order to calculate advection, we have selected fifteen AWS locations in the LIAISE domain (Fig. 5a). During the day, the wind direction was predominantly from the west. In the late afternoon when the sea breeze arrives, the wind direction shifted so that it is predominantly from the east. This means that during the day, air is advected from the hot and dry semi-arid steppes at the center of the Ebro basin through the irrigated area. As the air is advected through the wet landscape, the air mass moistens and cools. After the air mass is modified by the wet landscape, then it advects over the hot and dry natural vegetation area in the LIAISE domain. Therefore, we expect relatively warm and dry advection across the wet landscape and relatively cool and wet advection into the dry landscape during the day. When the sea breeze arrived in the late afternoon, it introduced cool, moist air from the sea. At that time, the air mass was modified due to a relatively hot and dry surface. Based on these assumptions, we have selected the stations in blue and white in Fig. 5a to calculate the local advection at the alfalfa local scale, and the red and white stations in Fig. 5a to calculate the local advection at the fallow local scale.

We calculated the advection using

$$adv_x = \bar{U} \frac{dX}{dr} \qquad (2)$$

where $\bar{U}$ is the mean wind speed, $\frac{dX}{dr}$ is the gradient of the scalar (e.g. potential temperature or specific humidity) between stations that align in the mean wind direction where $dr$ is the distance between stations. We calculated the advection between each station that falls within the mean wind direction for each 10-min interval and averaged them to calculate the mean advection for both Els PLans (fallow) and La Cendrosa (alfalfa) (Fig. 5). We input the advection of heat and moisture terms hourly in the CLASS model runs, and between the updated values, the advection terms are linearly interpolated. See Appendix B for more details about the advection calculation using the AWS network.

A typical diurnal cycle of advection was calculated from an average of data for all of July 2021. We chose to use the monthly mean diurnal cycle of advection in the CLASS model because it smoother than the advection for the composite day, but shows the same pattern and magnitude. Fig. 5b is the mean diurnal cycles of temperature and Fig. 5c is the mean diurnal cycle of moisture advection. The diurnal cycle of advection terms are reasonable based on the large-scale forcing observed during the LIAISE campaign. In the wet area, there was warm and slightly dry advection in the mid-day, as the synoptic forcing was from westerlies from inland in the Iberian peninsula. The mid-day temperature advection across the wet-dry boundary was not appreciable, but there was strong moist air advection. In the mid-afternoon – after 15:00 UTC – when the sea breeze arrived to the LIAISE domain, the advection terms for both the wet and dry areas look similar. The advected air mass was cool and moist, which corresponds to what is expected from the sea breeze.

## 4.3. Data constraints and verification for CLASS model

The CLASS model cases for the scales were heavily constrained by data observed locally and upscaled using the technique described in Section 4.1. The mixed-layer profile and the state of the atmosphere were initialized once at 06:00 UTC. We defined the surface at each scale using measurements from the LIAISE campaign. The mixed layer was initialized at the start of the model run, but advection was a dynamically changing parameter. At the local scale, we defined advection with the AWS network described in Section 4.2. At the regional and landscape scales, advection was weak and used only as a tuning
<figure><img src="image_8.png" /><figcaption>Fig. 5. (a) The Automated Weather Stations from SMC used to calculate the advection in the local alfalfa and fallow fields. The blue points are stations used only for the wet advection, the red points are used only for dry advection, and the white points are stations used in both the wet and dry fields. The base map shows the elevation of the LIAISE domain above sea-level (European Digital Elevation Model, version 1.1). (b) Diurnal Cycle of temperature advection for the wet field (green) and the dry field (yellow). (c) Diurnal cycle of moisture advection for the wet field (green) and the dry field (yellow). The black dashed line indicates noon, and the gray dot-dashed line indicates the approximate time that the sea breeze starts.</figcaption></figure>

parameter to capture the sea breeze. We assumed that there is no large-scale subsistence, because we have selected days with weak synoptic forcing. According to ERA5, wind divergence is on the order of $10^{-6}$ s⁻¹. In the CLASS model, the subsistence velocity is defined as the wind diverge times the boundary layer depth (Vilà-Guerau de Arellano et al., 2015). If we assume that the boundary layer height is on the order of 1000 m over the day, we find the subsistence velocity is on the order of $10^{-3}$ m s⁻¹, which is an order of magnitude lower than the entrainment velocity. For all cases, the CLASS model was verified using either local surface fluxes or the aggregated fluxes to the landscape and regional scales for the surface performance. Hourly radiosondes launched from Els Plans and La Cendrosa were used to verify the mixed-layer model performance.

Table 2 summarizes the differences between the spatial scales introduced in Section 3 and how they are replicated using the CLASS model. The columns indicating scale view and model view show the spatial extent of the heterogeneity in each scale and how that is represented using LAI and $c_{veg}$ in the CLASS model. The columns indicating the model initialization and verification focus on the differences between the scales. For the verification, the variables are the same, but the representative areas differ.

## 5. Results

After developing the scaling scheme to integrate the observations in a conceptual modeling framework, and combining measurements to create a dataset to represent fluxes and advection across these scales, we can evaluate the results of the CLASS model runs across these scales. The results are presented with a bottom-up approach. In Section 5.1, we display the results for the surface energy components at all scales. In Section 5.2, we display the mixed-layer development at the regional and landscape scales.

### 5.1. Surface energy balance

#### Observations

To verify the results of the CLASS model, we focused on the surface energy budget components for each of the scales. We have observed surface fluxes at all scales via direct eddy-covariance at the local scales and the composite surface fluxes from the flux maps at the regional

and landscape scales. At the local scale, we measured a $\beta$ of ~30 for the dry location and a $\beta$ of ~0.01 for the wet location. At larger scales, these extremes are tempered: we observed $\beta$ of ~0.6 and ~2.7 for the wet and dry landscape scales respectively. At the regional scale, the observed Bowen ratios converged on a $\beta$ of ~1.5. There were different surface fluxes depending on the spatial scale. We mimicked this in the CLASS model using parameter aggregation of the land surface based on the surface maps (Section 4.1).

#### Regional and Landscape Scales

At the regional scale, the net radiation was similar across the scales, however, the dry landscape had higher midday ground heat flux compared to the regional and the wet landscape scales (Fig. 6). This indicates that there is more available energy in the irrigated areas than the non-irrigated areas. ERA5 overestimated the ground heat flux at the regional scale, so it provided less available energy than there was in reality.

The partitioning of the available energy changed across the scales. In the wet landscape, approximately $2/3$ of the available energy was contributing to latent heat flux compared to sensible heat flux. Conversely, in the dry landscape, approximately $3/4$ of the available energy was partitioned into the sensible heat flux. At the regional scale, there was slightly more energy given to sensible heat flux than latent heat flux, but it fell somewhat in the middle of the wet and the dry landscapes. The model results for the latent heat flux match the observations reasonably for the wet landscape and regional scales; there was a slight dry bias in the regional scale model. Conversely, at the dry landscape scale, there was a dry bias in the latent heat flux, which corresponded to an overestimation of sensible heat flux at the same scale. This was likely due to the top layer of soil moisture that is prescribed to be drier than reality at this scale. The CLASS model overestimated the sensible heat flux relative to observations at all scales.

Compared to the CLASS model at the regional and landscape scales, ERA5 performs worst in all components of the surface energy balance. Although, the results of net radiation were reasonable, it overestimated ground heat flux, so ERA5 provided less available energy than observed. Moreover, ERA5 partitioned this energy almost entirely into the sensible heat flux. The sensible heat flux from ERA5 seems to be a reasonable match for the dry landscape scale, however, ERA5 captured next to no latent heat flux in the LIAISE regional scale.
<caption>**Table 2**</caption>

The scale definition and numerical scheme for the conceptual model. The scale-view shows the horizontal extent of each scale. The model-view demonstrates how the model observes the land surface based on the LAI and $c_{veg}$. The input parameters are the ones derived from data that change between scales, while the verification parameters are dynamically changing from observations described in the “Data Source” column.

<table><thead><tr><th rowspan="3">Scale</th><th rowspan="3">Scale-View</th><th rowspan="3">Model-View</th><th colspan="4">Input Parameters</th><th colspan="2">Verification Parameters</th><th rowspan="3">Data Source</th></tr><tr><th colspan="3">Surface</th><th rowspan="2">Mixed Layer</th><th rowspan="2">Surface</th><th rowspan="2">Mixed Layer</th></tr><tr><th>LAI</th><th>c<sub>veg</sub></th><th>SM</th><th>Adv</th></tr></thead><tbody><tr><td>Regional (~10 km)</td><td>LIAISE</td><td></td><td>1.25</td><td>0.5</td><td>0.2</td><td>Synop.</td><td>R<sub>n</sub>, G<br>H, LE</td><td>h<br>θ, δθ<br>q, δq</td><td>Radiosondes<br>Surface Maps: Regional</td></tr><tr><td></td><td>Wet</td><td></td><td>3</td><td>0.75</td><td>0.2</td><td>Synop.</td><td>R<sub>n</sub>, G<br>H, LE</td><td>h<br>θ, δθ<br>q, δq</td><td>Radiosondes<br>Surface Maps: Landscape</td></tr><tr><td>Landscape (~1 km)</td><td>Dry</td><td></td><td>0.75</td><td>0.7</td><td>0.2</td><td>Synop.</td><td>R<sub>n</sub>, G<br>H, LE</td><td>h<br>θ, δθ<br>q, δq</td><td>Radiosondes<br>Surface Maps: Landscape</td></tr><tr><td></td><td>Alfalfa</td><td></td><td>3</td><td>1.0</td><td>0.21</td><td>AWS</td><td>R<sub>n</sub>, G<br>H, LE</td><td>h<br>θ, δθ<br>q, δq</td><td>Surface Fluxes</td></tr><tr><td>Local (~10 m)</td><td>Fallow</td><td></td><td>0.01</td><td>0.1</td><td>0.1</td><td>AWS</td><td>R<sub>n</sub>, G<br>H, LE</td><td>h<br>θ, δθ<br>q, δq</td><td>Surface Fluxes</td></tr></tbody></table>

<figure><img src="image_9.png" /><figcaption>**Fig. 6.** The components of the energy budget for the LIAISE composite day for both data (dots), CLASS model runs (black lines), and ERA5 (gray line). The data for the landscape (yellow and dark green) and regional scales (green-yellow) were calculated via the flux maps and the SEB network.</figcaption></figure>

## Local Scales

The local scales of the CLASS model were compared to the field level measurements in both Els Plans and La Cendrosa (Fig. 7). The advection term described in Section 4.2 was prescribed here for both the wet and fallow local scales respectively (black lines in Fig. 7). The brown lines represent the local scale model case without advection applied for reference. At the local scale, there was higher net radiation at the wet site than the dry site, which the CLASS model was unable to capture. Between the higher outgoing longwave radiation and ground heat flux in the dry site, overall, there was on average 271 W m⁻² h⁻¹ of additional energy at the irrigated site compared to the dry site.

At the local scale, the surface fluxes were more extreme than at the landscape and regional scales. The alfalfa local case had almost all of the energy partitioned into the latent heat flux. At the alfalfa local scale, the sensible heat flux even became negative in the afternoon, which the CLASS model was able to capture. Conversely, at the fallow local scale, all of the energy was partitioned into the sensible heat flux, and there was little measured or modeled latent heat flux. The model case for the fallow local scale overestimated the sensible heat flux because the extra available energy that the model prescribes is added to this term. Like the regional and landscape cases, ERA5 overestimated the ground heat flux, however, it matched both the sensible and latent heat flux for the fallow local scale better than any other scale.
<figure><img src="image_10.png" /><figcaption>Fig. 7. The components of the energy budget for the LIAISE composite day for both data (dots), CLASS model runs including advection (black lines), the CLASS model runs without advection (brown lines) and ERA5 (gray line). The observations and the CLASS model runs are only for the local scale.</figcaption></figure>

## 5.2. Mixed-layer development

### Observations

During the LIAISE composite day, the boundary layer grew to a maximum value of approximately 1500 m in both the wet and the dry areas (Fig. 8), but there is a clear difference in the morning growth of the boundary layer. In the dry area, the boundary layer grew faster in the early morning than in the wet area. This can be explained by the magnitude of the sensible heat flux as shown in Figs. 7 and 8. It is evident from Fig. 8 that the maximum observed mixed-layer height in the dry area is 150 m higher than the maximum in the wet area, but the timing of the maximum boundary-layer height is different between landscapes. The mixed layer above the dry area tended to be warmer and drier than the wet area. On average, the mixed layer in the dry area was 1.3 K warmer and 0.87 $kg^{-1}$ drier than the mixed layer in the wet area.

The maximum boundary-layer height occurred around 13:00 and 15:00 UTC for the wet and dry sites respectively before the sea breeze arrived and stunted the boundary layer growth. The temperature followed a similar pattern: the mixed layer heated up throughout the morning and reached its maximum just before the sea breeze arrived in the late afternoon. In the morning, the mixed layer was humid. As the boundary layer grew, it dried out for both the wet and the dry areas. When the sea breeze arrived, the mixed layer humidified.

### Regional and Landscape Scales

In Fig. 8 we show the results of the mixed-layer height and mean potential temperature and specific humidity for the regional, landscape

and local scales, ERA5 and the observations for the mixed layer over the LIAISE composite day. The modeled regional scale fell between the observed mixed layers in the wet and the dry areas. Although we do not have observations of the mixed layer for the regional scale, we assumed that the regional scale mixed-layer characteristics would fall somewhere between those measured in both the wet and dry areas as it is driven by the sensible and latent heat fluxes in the regional scale (Fig. 6). Therefore, the regional scale CLASS model captures the observed regional scale mixed layer. The landscape scale CLASS model cases demonstrated the differences observed in the wet and dry landscape; the wet landscape had a cooler and wetter mixed layer than the dry landscape.

ERA5 performed poorly for the boundary-layer height as it predicted a maximum boundary-layer height of 2870 m. This relates to its partitioning of the surface energy budget: nearly all of the available energy was prescribed to sensible heat flux, which was used to grow the boundary layer. This result can be replicated with the CLASS model when we switched off the interactive land surface model if the same $\beta$ was prescribed. ERA5 captured the mixed-layer potential temperature well compared to the dry landscape observations, but it tended to dry out the mixed layer more than observed in either landscape.

## 6. Discussion: Integration of spatial scales

As shown in the previous section, the surface fluxes in a heterogeneous domain depend on the scale of the heterogeneity; however, we need to integrate processes across spatial scales to understand the processes that govern land-atmosphere interactions. Therefore, we
<figure><img src="image_11.png" /><figcaption>Fig. 8. The components of the mixed layer for the LIAISE composite day for both data (dots), CLASS model runs for the landscape and regional scales (black lines) and the local scale (brown lines), and ERA5 (gray line) for the regional and landscape scales.</figcaption></figure>

will first discuss the limitations of our observational approach before exploring the feedbacks between surface fluxes and boundary-layer dynamics from our numerical experiments. Finally, we will discuss briefly the differences between ERA5 and our regional scale case.

We used observational data observed at the local scale to define the spatial scales in the CLASS model and to verify the dynamic behavior of the CLASS model. Nonetheless, some of the assumptions we made to create the flux maps warrant discussion. We assumed that surface characteristics, including LAI and the soil moisture can be linearly averaged over the domain. Some parameters – like soil composition and moisture – may be non-linear in space, so our procedure may cause errors in the observations that are used as model input. Another potential limitation in the data approach is that there is a mismatch between what experimentalists measure with a single SEB station and how modelers use the data for model verification. Measurements are typically in the surface layer – well below the lowest grid cell of a model – and there could be problems with the constant flux layer assumption, as we observed in the alfalfa local scale. Furthermore, single point measurements have relatively small footprints compared to that of a model, so measurements are heavily impacted by surface conditions. It is unlikely that the single SEB measurement is representative of all of that crop type in the LIAISE domain. We addressed this mismatch in the data upscaling section (Section 4), but our method is limited through the assumption that all fields of a given crop type behave the same.

At the regional scale, we replicated surface fluxes using a composite of surface characteristics including soil moisture, vegetation cover and leaf area index. The observed regional scale $\beta$ was ~1.5, which was the necessary $\beta$ to model the mixed-layer height using the conceptual model. At the local scales, microscale advection was required to correctly capture the high latent heat fluxes in the alfalfa local scale. At this scale, the sensible heat flux became negative in the afternoon, which we were able to capture with the CLASS model. We also found that this negative heat flux only occurs locally. It has been observed only in the roughness sublayer of the La Cendrosa alfalfa field. This suggests that non-local processes are impacting the local scale, including advection. Because this is such a local phenomenon, it implies there is an internal boundary layer that forms in the heaviest vegetative areas within the wet landscape scale. This could mean that energy from the regional or landscape scales impact the local energy budget.

At the landscape scales, the surface fluxes were less extreme than local scales; however, they retained more of the characteristics of the different surfaces than the regional scale. This is because the landscape scales are not entirely wet or dry: there was a mix of crop types across both scales. For example, fluxes from fruit and nut trees – both with and without irrigation – have higher $\beta$ than the annual crops like maize

or alfalfa. For that reason, the crop mixture moderates the fluxes in both the wet and the dry landscapes. In the wet landscape scale, there are drier fields and even fallow fields which drove the boundary layer development.

Because the sensible heat flux is one of the most important factors in determining boundary-layer height (Ball, 1960), the mixed layer results for each scale changed based on the surface fluxes. We found that a modeled $\beta$ of ~1.8 replicates the correct mixed-layer characteristics defined by the boundary layer growth and the diurnal variability of state variables. This is just slightly higher than what we measured at the regional scale, which implies that the boundary layer on the regional scale is formed via a composite of surface fluxes from the LIAISE region. This is supported by the fact that the mean of the surface fluxes for wet and dry landscape cases is approximately the surface fluxes from the regional scale. The regional scale had a boundary-layer height that is 10% higher and 8% lower than the wet and dry landscape scales respectively. Therefore, even at the landscape scale, a composite of surface characteristics is able to reproduce the growth of the boundary layer.

We did not capture the sea breeze well at the end of the model period, which is due to the assumption of no synoptic scale subsidence. In the CLASS model, the boundary layer can only shrink if (1) the sensible heat flux is negative or (2) there is a contribution of the large-scale subsidence. When the sea breeze arrived, the observed sensible heat flux at the regional scale was positive, so even including advection, we could not capture the collapse of the mixed-layer height without imposing subsidence. One should be aware of the sea breeze when interpreting the mixed-layer results at the end of the afternoon.

At the local scale, surface fluxes were more extreme; however, we found that the fallow local case also well represents the boundary-layer height and temperature. The wet landscape case showed a marked difference in boundary-layer height (30% lower), temperature (2 K cooler), and specific humidity (0.3 g kg$^{-1}$ wetter) than the regional case. This illustrates that the regional scale boundary layer is more influenced by the presence of the extremely dry fields than those of the extremely wet fields. Furthermore, because we added advection in the local cases and the mixed layer for the alfalfa local scale does not match that of the regional scale, we can infer that the mixed layer – even at the local scales – is formed by a composite of surface fluxes in the region instead of the boundary layer being advected from downwind. This could mean that if the dry region were to be irrigated in the future, there would be a pronounced difference in the regional boundary layer compared to now.

The idea that there is an internal boundary layer over the wet region and that the regional boundary layer is formed by a composite
of surface fluxes was supported by the radiosonde observations in Fig. 2. The radiosonde over La Cendrosa, which was characterized by stable thermal stratification, gradually dries and warms with height as it began to observe fluxes that originate from the regional scale. Moreover, the radiosonde from Els Plans observed a mixed layer that is both cooler and wetter than its surface. With this analysis, it is not apparent if there is a traditional internal boundary layer that is formed on the local and landscape scales, or if it better matches the idea of a blending height where the regional measurements converge within the boundary layer.

At the local scale when we added advection (Fig. 7), we could quantify how much the non-local boundary-layer processes impact the surface fluxes in both the alfalfa and fallow fields. With advection in the alfalfa local scale, the latent heat flux was 10% higher and the sensible heat flux was 13% lower than the case without advection. The mean differences do not tell the entire story: the inclusion of advection allows the slight temporal shift in the latent heat flux compared to net radiation observed in Figs. 2 and 7. It also ensures that the sensible heat flux becomes negative in the afternoon locally. At the fallow local scale, the inclusion of advection did not appreciably change the latent heat flux; however, it changed how quickly the surface cools down during the sea breeze, which impacted the partitioning between sensible and ground heat flux in the late afternoon.

Mixed-layer theory has no applicability during times of stable stratification, which has been observed in the afternoon at the alfalfa local scale. However, because the stably stratified air near the surface is topped by a convective boundary layer, we apply the CLASS model to capture the convective layer above it. For that reason, we cannot capture the correct sensible heat flux nor the correct boundary layer stratification at this scale without the presence of advection. At the alfalfa local scale, the main limitation of applying a mixed-layer model is that one is unable to capture the internal boundary layer that was observed.

In order to understand why ERA5 performed badly in the LIAISE domain, we should consider that by using the CLASS model, we replicated a land surface model that prescribes heterogeneous surfaces using the parameter aggregation method, while ERA5 uses the flux aggregation approach. We found that ERA5 best matches the local surface fluxes from Els Plans (fallow local scale), although its spatial extent was that of the regional scale. We hypothesize that the reason ERA5 performs badly in the domain is that it fails to capture the subgrid scale heterogeneity due to irrigation. The average soil moisture across all tiles is 0.1 m⁻³, which is similar to that measured in the top 5 cm of Els Plans. The lack of soil moisture insures that even if the vegetation were to be parameterized correctly, there is not enough water in the model to correctly partition the surface fluxes. Capturing irrigation (or the lack thereof) is a well-documented weakness in weather and climate models in arid and semi-arid agricultural areas (Alexander et al., 2022; Lawston et al., 2015, 2020; Qian et al., 2020). However, we found that the lack of irrigation is a small part of the differences between ERA5 and our regional scale case. When we ran the regional case using the soil moisture from ERA5, we found that the β increases from ~1.8 to ~2.2 (compared with β<sub>ERA5</sub> ~17.9). The difference in the mixed-layer characteristics is negligible between the regional case and the regional case run with ERA5 soil moisture. Instead, the difference in the boundary-layer development must be explained through a combination of the vegetation characteristics and the available soil moisture.

## 7. Summary and conclusions

In this study, we introduced a site-specific framework to investigate how observations of the boundary-layer dynamics connect to the surface energy budget across spatial scales. We used comprehensive observations of surface and atmospheric observations from the two-week LIAISE campaign in July 2021. The LIAISE domain was characterized by an extreme surface heterogeneity; there was a sharp

contrast between an irrigation area (~10 km) and the semi-arid area. The LIAISE experiment offered a unique possibility to study both how the surface heterogeneity controls on the boundary-layer dynamics and how the boundary-layer dynamics controls the diurnal variability of surface fluxes across spatial scales. We interpreted the scaling scheme by using a conceptual mixed-layer land-atmosphere model coupled to various land surfaces characterized by extreme heterogeneity. This study offered a unique opportunity to determine the reliability of land surface models in a heterogeneous environment. The introduced framework for interpolating spatial scales of heterogeneity is a promising method for verifying the performance of land surface models in heterogeneous areas.

By combining observations and systematic numerical results, we were able to quantify the relationship between surface properties and boundary-layer dynamics at local (~100 m), landscape (~1 km) and regional (~10 km) scales. Our aim was to connect the “large-patch” type of heterogeneity from the regional scale with the “unstructured” type of heterogeneity at the landscape scales and the “statistically homogeneous microscale” heterogeneity at the local scales. We raised two research questions to disentangle the bidirectional impacts between the land surface and the boundary-layer dynamics specifically for the LIAISE domain.

### 1. How does a heterogeneous surface control the boundary-layer dynamics across spatial scales?

The observed β in the alfalfa (wet) field was ~0.01 while the observed β in the fallow (dry) field was ~30, but the observed evolution of boundary-layer characteristics were similar despite the extreme surfaces. Using a land-atmosphere coupled model, we found that a β of ~1.8 was required to form the observed boundary-layer characteristics. This was slightly higher than the observed β at the regional scale (β~1.5) at the LIAISE regional scale. Using the observed surface fluxes, where the local surface fluxes are more extreme, the modeled mixed layer differed from the observed boundary-layer growth. The local fallow scale mean boundary-layer height was 10% higher than the regional scale and the local alfalfa scale mean boundary-layer height was 34% lower than the regional scale. This indicates that the observed boundary layer is formed from the aggregated regional landscape where the dry landscape has slightly more weight than the wet landscape, and there was an internal boundary layer that forms in the wet landscape and local scales. Therefore, advection was necessary to describe the boundary layer at local scales.

### 2. What is the impact of the boundary-layer dynamics on the diurnal variability of the surface fluxes?

In analyzing the influence of boundary-layer dynamics on surface fluxes, we found that another non-local effect plays a key role in reproducing the observations with the model: the advection of heat and moisture. We determined advection by calculating gradients of temperature and moisture from an AWS network in the LIAISE domain. During the day, there was warm, dry air advection within the wet landscape and cool, moist air advection across the wet-dry boundary. In the late afternoon after the sea breeze arrives, both the irrigated and dry landscapes experienced cool and moist advection. At the local scale, the inclusion of advection of heat and moisture was important for capturing the both the magnitude and timing diurnal cycle of surface fluxes. In the alfalfa (wet) local scale, advection allowed for more latent heat flux than the surface alone would allow, especially in the afternoon. At the landscape and regional scales, it appears that the boundary-layer dynamics played a less important role than the surface in controlling the partitioning of energy at the surface.

By using a coupled land-atmosphere model constrained by observations, we could disentangle the controls that the boundary layer exerts on the surface and those that the surface exerts on the boundary-layer dynamics. In thermally heterogeneous environments, land-atmosphere
interactions are complex: there are non-local drivers of both the atmospheric boundary layer and the partitioning of energy. We introduced a site-specific scaling framework to address the role that spatial scale plays in a subgrid scale heterogeneity, which can offer a template for future studies. We find that local scale surface observations are insufficient for explaining the boundary-layer dynamics at any scale. Instead, the observed boundary layer is formed via composite fluxes of sensible and latent heat over a horizontal extent of 10 km. In the LIAISE region, which is characterized by strong thermal heterogeneity, we found that the regional surface properties developed the boundary-layer dynamics, but the boundary layer feedback on the heat and moisture surface fluxes was a much less clear connection. In the future, this framework for surface heterogeneity can be used to evaluate the drivers of latent heat flux to further investigate the impacts of the boundary-layer dynamics on the surface fluxes.

## Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Data availability

Data will be made available on request.

## Acknowledgments

We would like to thank all of organizers, hosts and participants of the LIAISE campaign. In particular, we would like to acknowledge Martin Best, Joaquim Bellvert, Jennifer Brooke, Jan Polcher and Pere Quintana for their work on the LIAISE steering committee, and Henk Snellen, Getachew Adnew, Marc Castellnou, Jerry Chen, Kevin van Diepen, Kim Faassen, Raquel Gonzalez Armas, Wouter Mol, Robbert Moonen, Ruben Schulte and Gijs Vis for their work on the LIAISE-NL team during the experiment. Pere Quintana also provided the crop cover maps in the Ebro River Valley.

We would like to thank the two anonymous reviewers for their comments which greatly improved the quality of this paper.

This PhD project partly was supported by the appointment of Jordi Vila as Chair of the Meteorology and Air Quality Group of Wageningen University, The Netherlands.

## Appendix A

The SIGPAC land use map was reclassified to match the measured land cover types in the LIAISE domain (Fig. 9 and Table 3). Table 3 shows the original crop types and how they were reclassified to match the LIAISE experiment, while Fig. 9 shows the reclassified land use map that has been used to make the flux maps. There were a number of assumptions made to reclassify the actual crop types into the LIAISE land use types:

* **Cereal** In the LIAISE region, there are multiple types of cereal crops grown during the winter in the rainfed area including wheat, barley and oats. We assumed that by July, all of these cereal crops had been harvested, so that a fallow field or one with dry stubble remained. Therefore, we use the measurements from Els Plans to represent all cereal crops.

* **Maize** In the irrigated region, corn accounted for 20% of the landscape (Table 3). There were two stations measuring in corn fields during the LIAISE campaign, so these sites were averaged together to provide a more robust measure of the variability of maize fields in the region during the LIAISE campaign. We assumed with maize, most fields were in approximately the same growth stage and that they were all pre-senescent.

<figure><img src="image_12.png" /><figcaption>Fig. 9. The reclassified 100 m SIGPAC crop cover map for the LIAISE regional scale.</figcaption></figure>

* **Alfalfa** In the irrigated area, alfalfa accounted for approximately 15% of the landscape (Table 3). The LIAISE campaign fell within one growing cycle of alfalfa: the campaign began about a week after harvest and five days after the first irrigation. Therefore, we measured one alfalfa growing cycle from nearly bare soil to full crop cover. In order to properly replicate the temporal variability of growing stage in the regional alfalfa, we apply a growing stage to each alfalfa pixel and use the observed surface fluxes from the La Cendrosa alfalfa field at that stage. Because the surface plays an important role in determining both the available energy and its partitioning, we take all surface components instead of recalculating fluxes based on a Bowen ratio. This means that we have to assume that daily variability in synoptic and boundary layer dynamics are constant throughout the LIAISE domain, which is not necessarily true.

* **Fruit Orchards** We assume that the energy partitioning between different types of fruit trees (e.g. apple, pear, olives) are similar. There are also fruit trees located in both the irrigated and rainfed areas, but our reference apple orchard is partially irrigated; however, there were weighing lysimeter in both irrigated and non-irrigated apple trees in the orchard which were used to correct the Bowen ratio of the eddy-covariance system for the non-irrigated area.

* **Nut Orchards** We measured in one non-irrigated almond orchard during the LIAISE campaign. Because most of the almonds were all located in the non-irrigated area, it is fair to assume that all almond trees behave similarly.

* **Vineyards** Like almond orchards, most of the vineyards were located in the non-irrigated area where the vineyard surface energy budget (SEB) station was located. Like the almond trees, the one station is likely representative of all of the vineyards.

* **Water** There is a lake in the irrigated region. In this study, we modeled the energy budget for the water. In the future, energy budget components measured directly from the lake will be available. The flux maps will be updated to reflect this.

In addition to creating surface flux maps, we have also created products using the leaf-level ecological measurements including leaf area index, vegetative fraction, stomatal conductance, surface resistance,
<caption>**Table 3**</caption>

Reclassification scheme from SIGPAC crop cover to match the measurements of LIAISE SEB stations. The re-classified crop type is the selected land cover that matches a surface energy budget (SEB) station. The left three columns display the proportion of the total land cover that each crop type makes up across the regional and landscape scales.

<table><thead><tr><th>Re-classified crop</th><th>SEB station</th><th>SIGPAC crop</th><th>%Regional</th><th>%Wet landscape</th><th>%Dry landscape</th></tr></thead><tbody><tr><td>Cereals</td><td>Els plans - Natural vegetation</td><td>Soft Wheat, Colza, Barley, Oat, Triticale, Vetch, Peas</td><td>36.1</td><td>22.7</td><td>50.0</td></tr><tr><td>Maize</td><td>Boldu maize and IRTA maize</td><td>Corn</td><td>12.5</td><td>20.0</td><td>4.4</td></tr><tr><td>Alfalfa</td><td>La Cendrosa - Alfalfa</td><td>Alfalfa</td><td>8.5</td><td>14.6</td><td>2.0</td></tr><tr><td>Grass</td><td>IRTA - Grass</td><td>Festuca grass, Ray-grass</td><td>2.7</td><td>4.4</td><td>0.7</td></tr><tr><td>Vineyards</td><td>Verdu - Vineyard</td><td>Vineyards</td><td>1.0</td><td>0.1</td><td>2.0</td></tr><tr><td>Fruit trees</td><td>IRTA - Apple</td><td>Olive, Pear, Peach, Nectarine, Apricot, Date, Apple, Other fruit trees</td><td>13.3</td><td>16.0</td><td>10.3</td></tr><tr><td>Nut trees</td><td>Prexiana - Almonds</td><td>Almond trees</td><td>1.7</td><td>0.7</td><td>3.0</td></tr><tr><td>Urban</td><td>Modeled</td><td>N/A</td><td>21.2</td><td>19.6</td><td>23.3</td></tr><tr><td>Water</td><td>Modeled</td><td>N/A</td><td>0.27</td><td>0.5</td><td>0.0</td></tr></tbody></table>

and soil respiration. Although these products are not shown, they were used to create the composite land surface at different scales used in the CLASS model. In the future, we hope to maps for net ecosystem exchange, soil moisture and temperature from the data collected in the LIAISE network of SEB stations. The data from the LIAISE Unified Eddy-Covariance processing as well as the flux maps will be available to the LIAISE community and the public through the LIAISE database (https://liaise.aeris-data.fr/).

## Appendix B

We must make a number of assumptions in order to use the network of AWS locations to calculate boundary-layer advection:

1. We calculate advection with 2 m temperature and humidity and 10 m wind velocity to represent of the boundary-layer advection. This assumption will induce errors in our advection estimates because 2 m observations are heavily influenced by the surface conditions.

2. For the alfalfa local scale, we assume that the advection is the average of the advection across the wet landscape. Therefore, we have selected AWS stations in the irrigated part of the LIAISE domain to calculate advection (blue and white locations Fig. 5a).

3. For the fallow local scale, we assume that the advection is that which crosses from the wet landscape into the dry landscape. Therefore, we have selected AWS near the boundary of the wet-dry boundary (white and red locations in Fig. 5a).

4. We assume that the 10 m wind speed and direction is constant on the landscape scale.

5. The stations for each scale are aligned by mean wind direction for a 10-min interval. The advection term is calculated using

$$adv_x = \bar{U} \frac{dX}{dr} \qquad (3)$$

where $\bar{U}$ is the mean wind speed, $\frac{dX}{dr}$ is the gradient of the scalar (e.g. potential temperature or specific humidity) between stations that align in the mean wind direction where $dr$ is the distance between stations.

6. Stations that fall within a 30° window from the mean wind direction of another station are used to calculate the gradient of the scalar. All combinations of stations that fall within the 30° window of the mean wind direction are averaged together to find the advection for a 10 min interval.

7. The advection term is calculated from a mean of all scalar gradients from all appropriate station combinations for a given wind direction.

## References

Alexander, G.A., Holmes, H.A., Sun, X., Caputi, D., Faloona, I.C., Oldroyd, H.J., 2022. Simulating land-atmosphere coupling in the Central Valley, California: Investigating soil moisture impacts on boundary layer properties. Agricult. Forest Meteorol. 317, 108898. http://dx.doi.org/10.1016/j.agrformet.2022.108898.

Baldocchi, D., Falge, E., Gu, L., Olson, R., Hollinger, D., Running, S., Anthoni, P., Bernhofer, C., Davis, K., Evans, R., Fuentes, J., Goldstein, A., Katul, G., Law, B., Lee, X., Malhi, Y., Meyers, T., Munger, W., Oechel, W., U, K.T.P., Pilegaard, K., Schmid, H.P., Valentini, R., Verma, S., Vesala, T., Wilson, K., Wofsy, S., 2001. FLUXNET: A new tool to study the temporal and spatial variability of ecosystem-scale carbon dioxide, water vapor, and energy flux densities. Bull. Am. Meteorol. Soc. 82 (11), 2415-2434. http://dx.doi.org/10.1175/1520-0477(2001)082<2415: FANTTS>2.3.CO;2.

Baldocchi, D.D., Hincks, B.B., Meyers, T.P., 1988. Measuring biosphere-atmosphere exchanges of biologically related gases with micrometeorological methods. Ecology 69 (5), 1331-1340. http://dx.doi.org/10.2307/1941631.

Ball, F.K., 1960. Control of inversion height by surface heating. Q. J. R. Meteorol. Soc. 86 (370), 483-494. http://dx.doi.org/10.1002/qj.49708637005.

Beamesderfer, E.R., Buechner, C., Faiola, C., Helbig, M., Sanchez-Mejia, Z.M., Yáñez-Serrano, A.M., Zhang, Y., Richardson, A.D., 2022. Advancing cross-disciplinary understanding of land-atmosphere interactions. J. Geophys. Res.: Biogeosci. 127 (2), e2021JG006707. http://dx.doi.org/10.1029/2021JG006707.

Boone, A., Bellvert, J., Best, M., Brooke, J., Canut-Rocafort, G., Cuxart, J., Hartogensis, O., e Moigne, P., Miró, J.R., Polcher, J., Price, J., Quintana Seguí, P., Wooster, M., 2021. Updates on the international land surface interactions with the atmosphere over the iberian semi-arid environment (LIAISE) field campaign. GEWEX News 31 (4), 17-21.

Boone, A., Best, M., Cuxart, J., Polcher, J., Quintana, P., Bellvert, J., Brooke, J., Canut-Rocafort, G., Price, J., 2019. Land surface interactions with the atmosphere over the iberian semi-arid environment (LIAISE). GEWEX News 29 (1), 8-10.

Bou-Zeid, E., Anderson, W., Katul, G.G., Mahrt, L., 2020. The persistent challenge of surface heterogeneity in boundary-layer meteorology: A review. Bound.-Layer Meteorol. 177 (2), 227-245. http://dx.doi.org/10.1007/s10546-020-00551-8.

Bou-Zeid, E., Meneveau, C., Parlange, M.B., 2004. Large-eddy simulation of neutral atmospheric boundary layer flow over heterogeneous surfaces: Blending height and effective surface roughness. Water Resour. Res. 40 (2), http://dx.doi.org/10.1029/2003WR002475.

Butterworth, B.J., Desai, A.R., Townsend, P.A., Petty, G.W., Andresen, C.G., Bertram, T.H., Kruger, E.L., Mineau, J.K., Olson, E.R., Paleri, S., Pertzborn, R.A., Pettersen, C., Stoy, P.C., Thom, J.E., Vermeuel, M.P., Wagner, T.J., Wright, D.B., Zheng, T., Metzger, S., Schwartz, M.D., Iglinski, T.J., Mauder, M., Speidel, J., Vogelmann, H., Wanner, L., Augustine, T.J., Brown, W.O.J., Oncley, S.P., Buban, M., Lee, T.R., Cleary, P., Durden, D.J., Florian, C.R., Lantz, K., Riihimaki, L.D., Sedlar, J., Meyers, T.P., Plummer, D.M., Guzman, E.R., Smith, E.N., Sühring, M., Turner, D.D., Wang, Z., White, L.D., Wilczak, J.M., 2021. Connecting land-atmosphere interactions to surface heterogeneity in CHEESEHEAD19. Bull. Am. Meteorol. Soc. 102 (2), E421-E445. http://dx.doi.org/10.1175/BAMS-D-19-0346.1.

Buttner, G., 2014. CORINE land cover and land cover change products. In: Manakos, I., Braun, M. (Eds.), Land Use and Land Cover Mapping in Europe: Practices & Trends. Springer, Dordrecht, the Netherlands, pp. 55-74.
Conzemius, R., Fedorovich, E., 2007. Bulk models of the sheared convective boundary layer: Evaluation through large eddy simulations. J. Atmos. Sci. 64 (3), 786–807.
http://dx.doi.org/10.1175/JAS3870.1

Driedonks, A.G.M., Tennekes, H., 1984. Entrainment effects in the well-mixed atmospheric boundary layer. Bound.-Layer Meteorol. 30 (1), 75–105. http://dx.doi.org/10.1007/BF00121950.

Fratini, G., Mauder, M., 2014. Towards a consistent eddy-covariance processing: an intercomparison of EddyPro and TK3. Atmos. Meas. Tech. 7 (7), 2273–2281.
http://dx.doi.org/10.5194/amt-7-2273-2014.

Grimmond, C.S.B., Oke, T.R., 1999. Aerodynamic properties of urban areas derived from analysis of surface form. J. Appl. Meteorol. Climatol. 38 (9), 1262–1292.
http://dx.doi.org/10.1175/1520-0450(1999)038<1262:APOUAD>2.0.CO;2

van Heerwaarden, C.C., Vilà-Guerau de Arellano, J., Moene, A.F., Holtslag, A.A.M., 2009. Interactions between dry-air entrainment, surface evaporation and convective boundary-layer development. Q. J. R. Meteorol. Soc. 135 (642), 1277–1291. http://dx.doi.org/10.1002/qj.431.

van Heerwaarden, C.C., Mellado, J.P., Lozar, A.D., 2014. Scaling laws for the heterogeneously heated free convective boundary layer. J. Atmos. Sci. 71 (11), 3975–4000.
http://dx.doi.org/10.1175/JAS-D-13-0383.1.

Helbig, M., Gerken, T., Beamesderfer, E.R., Baldocchi, D.D., Banerjee, T., Biraud, S.C., Brown, W.O.J., Brunsell, N.A., Burakowski, E.A., Burns, S.P., Butterworth, B.J., Chan, W.S., Davis, K.J., Desai, A.R., Fuentes, J.D., Hollinger, D.Y., Kljun, N., Mauder, M., Novick, K.A., Perkins, J.M., Rahn, D.A., Rey-Sanchez, C., Santanello, J.A., Scott, R.L., Seyednasrollah, B., Stoy, P.C., Sullivan, R.C., de Arellano, J.V.-G., Wharton, S., Yi, C., Richardson, A.D., 2021. Integrating continuous atmospheric boundary layer and tower-based flux measurements to advance understanding of land-atmosphere interactions. Agricult. Forest Meteorol. 307, 108509.
http://dx.doi.org/10.1016/j.agrformet.2021.108509.

Hersbach, H., Bell, B., Berrisford, P., Hirahara, S., Horányi, A., Muñoz-Sabater, J., Nicolas, J., Peubey, C., Radu, R., Schepers, D., Simmons, A., Soci, C., Abdalla, S., Abellan, X., Balsamo, G., Bechtold, P., Biavati, G., Bidlot, J., Bonavita, M., De Chiara, G., Dahlgren, P., Dee, D., Diamantakis, M., Dragani, R., Flemming, J., Forbes, R., Fuentes, M., Geer, A., Haimberger, L., Healy, S., Hogan, R.J., Hólm, E., Janisková, M., Keeley, S., Laloyaux, P., Lopez, P., Lupu, C., Radnoti, G., de Rosnay, P., Rozum, I., Vamborg, F., Villaume, S., Thépaut, J.-N., 2020. The ERA5 global reanalysis. Q. J. R. Meteorol. Soc. 146 (730), 1999–2049. http://dx.doi.org/10.1002/qj.3803.

Institut Cartogràfic i Geològic de Catalunya, 2009-2020. GT IV. Soil map 1:25,000.
URL: https://www.icgc.cat/en/Public-Administration-and-Enterprises/Downloads/Geological-and-geothematic-cartography/Soil-cartography/GT-IV.-Soil-map-1-25-000.

Kadygrov, N., Broquet, G., Chevallier, F., Rivier, L., Gerbig, C., Ciais, P., 2015. On the potential of the ICOS atmospheric CO₂ measurement network for estimating the biogenic CO₂ budget of Europe. Atmos. Chem. Phys. 15 (22), 12765–12787.
http://dx.doi.org/10.5194/acp-15-12765-2015.

Kaimal, J.C., Finnigan, J.J., 1994. Atmospheric Boundary Layer Flows: Their Structure and Measurement. Oxford University Press, New York.

Lawston, P.M., Santanello, J.A., Hanson, B., Arsensault, K., 2020. Impacts of irrigation on summertime temperatures in the Pacific Northwest. Earth Interact. 24 (1), 1–26.
http://dx.doi.org/10.1175/EI-D-19-0015.1.

Lawston, P.M., Santanello, J.A., Zaitchik, B.F., Rodell, M., 2015. Impact of irrigation methods on land surface model spinup and initialization of WRF forecasts. J. Hydrometeorol. 16 (3), 1135–1154. http://dx.doi.org/10.1175/JHM-D-14-0203.1.

Lemonsu, A., Grimmond, C.S.B., Masson, V., 2004. Modeling the surface energy balance of the core of an old mediterranean city: Marseille. J. Appl. Meteorol. (1988-2005) 43 (2), 312–327.

Lothon, M., Lohou, F., Pino, D., Couvreux, F., Pardyjak, E.R., Reuder, J., Vilà-Guerau de Arellano, J., Durand, P., Hartogensis, O., Legain, D., Augustin, P., Gioli, B., Lenschow, D.H., Faloona, I., Yagüe, C., Alexander, D.C., Angevine, W.M., Bargain, E., Barrié, J., Bazile, E., Bezombes, Y., Blay-Carreras, E., van de Boer, A., Boichard, J.L., Bourdon, A., Butet, A., Campistron, B., de Coster, O., Cuxart, J., Dabas, A., Darbieu, C., Deboudt, K., Delbarre, H., Derrien, S., Flament, P., Fourmentin, M., Garai, A., Gibert, F., Graf, A., Groebner, J., Guichard, F., Jiménez, M.A., Jonassen, M., van den Kroonenberg, A., Magliulo, V., Martin, S., Martinez, D., Mastrorillo, L., Moene, A.F., Molinos, F., Moulin, E., Pietersen, H.P., Piguet, B., Pique, E., Román-Cascón, C., Rufin-Soler, C., Saíd, F., Sastre-Marugán, M., Seity, Y., Steeneveld, G.J., Toscano, P., Traullé, O., Tzanos, D., Wacker, S., Wildmann, N., Zaldei, A., 2014. The BLLAST field experiment: Boundary-Layer Late Afternoon and Sunset Turbulence. Atmos. Chem. Phys. 14 (20), 10931–10960. http://dx.doi.org/10.5194/acp-14-10931-2014.

Noilhan, J., Mahfouf, J.F., 1996. The ISBA land surface parameterisation scheme. Glob. Planet. Change 13 (1), 145–159. http://dx.doi.org/10.1016/0921-8181(95)00043-7.

Noilhan, J., Planton, S., 1989. A simple parameterization of land surface processes for meteorological models. Mon. Weather Rev. 117 (3), 536–549. http://dx.doi.org/10.1175/1520-0493(1989)117<0536:ASPOLS>2.0.CO;2

Patton, E.G., Sullivan, P.P., Moeng, C.-H., 2005. The influence of idealized heterogeneity on wet and dry planetary boundary layers coupled to the land surface. J. Atmos. Sci. 62 (7), 2078–2097. http://dx.doi.org/10.1175/JAS3465.1.

Qian, Y., Yang, Z., Feng, Z., Liu, Y., Gustafson, W.I., Berg, L.K., Huang, M., Yang, B., Ma, H.-Y., 2020. Neglecting irrigation contributes to the simulated summertime warm-and-dry bias in the central United States. Npj Clim. Atmos. Sci. 3 (1), 1–10.
http://dx.doi.org/10.1038/s41612-020-00135-w.

Raasch, S., Harbusch, G., 2001. An analysis of secondary circulations and their effects caused by small-scale surface inhomogeneities using large-eddy simulation. Bound.-Lay. Meteorol. 101 (1), 31–59. http://dx.doi.org/10.1023/A:1019297504109.

Román-Cascón, C., Lothon, M., Lohou, F., Hartogensis, O., Vilà-Guerau de Arellano, J., Pino, D., Yagüe, C., Pardyjak, E.R., 2021. Surface Representation Impacts on Turbulent Heat Fluxes in WRF(v.4.1.3). preprint, Atmospheric sciences, http://dx.doi.org/10.5194/gmd-2020-371.

Shen, S., Leclerc, M.Y., 1994. Large-eddy simulation of small-scale surface effects on the convective boundary-layer structure. Atmos.-Ocean 32 (4), 717–731. http://dx.doi.org/10.1080/07055900.1994.9649519.

Shen, S., Leclerc, M.Y., 1995. How large must surface inhomogeneities be before they influence the convective boundary layer structure? A case study. Q. J. R. Meteorol. Soc. 121 (526), 1209–1228. http://dx.doi.org/10.1002/qj.49712152603.

Swinbank, W.C., 1951. The measurement of vertical transfer of heat and water vapor by eddies in the lower atmosphere. J. Atmos. Sci. 8 (3), 135–145. http://dx.doi.org/10.1175/1520-0469(1951)008<0135:TMOVTO>2.0.CO;2.

Tennekes, H., 1973. A model for the dynamics of the inversion above a convective boundary layer. J. Atmos. Sci. 30 (4), 558–567. http://dx.doi.org/10.1175/1520-0469(1973)030<0558:AMFTDO>2.0.CO;2.

Twine, T.E., Kustas, W.P., Norman, J.M., Cook, D.R., Houser, P.R., Meyers, T.P., Prueger, J.H., Starks, P.J., Wesely, M.L., 2000. Correcting eddy-covariance flux underestimates over a grassland. Agricult. Forest Meteorol. 103 (3), 279–300.
http://dx.doi.org/10.1016/S0168-1923(00)00123-4.

Vilà-Guerau de Arellano, J., van Heerwaarden, C.C., van Stratum, B.J.H., van den Dries, K., 2015. Atmospheric Boundary Layer: Integrating Chemistry and Land Interactions. Cambridge University Press, New York, New York.

Vilà-Guerau de Arellano, J., Ney, P., Hartogensis, O., de Boer, H., van Diepen, K., Emin, D., de Groot, G., Klosterhalfen, A., Langensiepen, M., Matveeva, M., Miranda-García, G., Moene, A.F., Rascher, U., Röckmann, T., Adnew, G., Brüggemann, N., Rothfuss, Y., Graf, A., 2020. CloudRoots: integration of advanced instrumental techniques and process modelling of sub-hourly and sub-kilometre land-atmosphere interactions. Biogeosciences 17 (17), 4375–4404. http://dx.doi.org/10.5194/bg-17-4375-2020.