Connecting Land-Atmosphere Interactions to Surface Heterogeneity in CHEESEHEAD 2019

The Chequamegon Heterogeneous Ecosystem Energy-balance Study Enabled by a High-density Extensive Array of Detectors 2019 (CHEESEHEAD19) is an ongoing National Science Foundation project based on an intensive field campaign that occurred from June-October 2019. The purpose of the study is to examine how the atmospheric boundary layer responds to spatial heterogeneity in surface energy fluxes. One of the main objectives is to test whether lack of energy balance closure measured by eddy covariance (EC) towers is related to mesoscale atmospheric processes. Finally, the project evaluates data-driven methods for scaling surface energy fluxes, with the aim to improve model-data comparison and integration. To address these questions, an extensive suite of ground, tower, profiling, and airborne instrumentation was deployed over a 10×10 km domain of a heterogeneous forest ecosystem in the Chequamegon-Nicolet National Forest in northern Wisconsin USA, centered on the existing Park Falls 447-m tower that anchors an Ameriflux/NOAA supersite (US-PFa / WLEF). The project deployed one of the world’s highest-density networks of above-canopy EC measurements of surface energy fluxes. This tower EC network was coupled with spatial measurements of EC fluxes from aircraft, maps of leaf and canopy properties derived from airborne spectroscopy, ground-based measurements of plant productivity, phenology, and physiology, and atmospheric profiles of wind, water vapor, and temperature using radar, sodar, lidar, microwave radiometers, infrared interferometers, and radiosondes. These observations are being used with large eddy simulation and scaling experiments to better understand sub-mesoscale processes and improve formulations of sub-grid scale processes in numerical weather and climate models.

Land-atmosphere exchanges of energy, water, and carbon influence weather and climate. The 86 biological processes that mediate these exchanges with the atmosphere occur at multiple spatial 87 and temporal scales, necessitating a variety of cross-scale observational platforms. Accurate 88 accounting of land-atmosphere interactions is critical for improving the predictive performance 89 of numerical weather and climate models. Unfortunately, there is a persistent mismatch between 90 the scales of observations and models. This scale mismatch is problematic because natural 91 environments exhibit substantial heterogeneity in their surface characteristics, which means that 92 observations are not always accurate reflections of the entire model grid cell. Furthermore, the 93 atmosphere is strongly influenced by nonlinear two-way interactions with radiation, land cover, 94 and soil, so that the spatial and temporal scaling of surface fluxes is fundamental to assessing the 95 parameterizations used in atmospheric models to represent land-atmospheric interactions. 96 97 The notion that land surface heterogeneity influences the surface energy balance, and the 98 resulting atmospheric responses, emerged from early model simulations showing the importance 99 of soil moisture, vegetation, albedo, roughness, and heating on the atmosphere (Garratt 1993; 100 Mahrt 2000; Betts et al. 1996;Charney 1975;Avissar 1995;Pielke et al. 1998 responses scale linearly or non-linearly and whether they differ for dry versus moist dynamics 104 (Raupach and Finnigan 1995). Modeling studies on this topic have been developed from limited 105 sets of observations of prior field experiments and from specialized modeling domains using 106 simplified boundary conditions (e.g., Kang et al., 2007;Hill et al., 2008Hill et al., , 2011Zhu et al., 2016). 107 From these previous studies, scaling laws have been derived based on numerical simulations 108 (van Heerwaarden et al. 2014; Rihani et al. 2015), but a systematic regional-scale observational 109 experiment that quantifies the multi-scale nature of sub-grid scaling and patterning has never 110 been fully realized (Steinfeld et al. 2007). 111 112 An issue related to how heterogeneity influences transport processes in the ABL is the energy 113 balance closure problem. This refers to an observed tendency in eddy covariance ( instrumentation reliability to test hypotheses on spatial heterogeneity and atmospheric feedbacks. 136 The two main research objectives for the CHEESEHEAD19 experiment were to 1) investigate 137 causes of energy balance non-closure over heterogeneous ecosystems and 2) to address the 138 problem of scaling surface energy fluxes. 139 140 There is currently no definitive answer as to what is responsible for energy balance non-closure. 141 The project was designed specifically to test the hypothesis that heterogeneity is responsible for 142 generating organized (sub-)mesoscale structures that are not resolved by traditional EC methods. 143 144 Various theories suggest that "spatial" EC, where multiple towers are combined to estimate the 145 mesoscale contribution to the total flux, could be used to analyze this contribution and "close" 146 the energy balance (Steinfeld et al. 2007; Mauder et al. 2008b The EC tower network consisted of 17 towers from the NSF Lower Atmosphere Observing  213  Facility (LAOF) ISFS, two additional towers, and the central tall Ameriflux tower. Ground-214  based measurement of vegetation occurred at 41 plots in the domain, plus an additional 10 plots  215 for measuring phenology. Airborne spectroscopy imaging was used to map leaf chemistry and 216 canopy properties. 217 218 The suite of atmospheric profiling instruments included the LAOF Integrated Sounding System 219 (ISS; Fig. 2c) and the UW SPARC system (Fig. 2a) Towers sampled three-dimensional wind velocity, temperature, moisture, and CO2 at 20 Hz to 276 measure land-atmosphere fluxes (τ, HS, HL, FCO2). Each tower also measured net radiation, soil 277 heat flux at 5 cm depth (and soil temperature profile, heat capacity, and moisture to determine 278 soil heat storage), and a 3-level air temperature and humidity profile to estimate canopy heat 279 storage. A majority of the sites were forested and had flux instruments mounted 33 m AGL (Fig.  280 3; Table S1). Instruments for wetland, grass, and lake sites were mounted between 1 -3 m AGL 281 to maintain consistent vegetation within the flux footprint. Tower placement within the 10 × 10 282 km study domain followed a stratified random grid pattern, taking into account practical 283 considerations including distance to road, suitable gap in trees for a tower, USFS-owned land, 284 etc. Individual towers were an average of 1.4 km from their nearest neighboring tower and an 285 average of 3.5 km from the tall tower. This meant that under certain conditions (e.g., high wind 286 speeds, stable stratification) several of the towers shared overlapping flux footprints; a favorable 287 condition for applying some of the data-driven scaling methods used in the project. Additionally, 288 the semi-random placement meant that the towers were not chosen by distributing the towers in 289 the centers of the most homogeneous areas of the various land cover types. Thus, within the 290 individual footprint of each tower there was often spatial variability in vegetation height and type 291 (deciduous vs. evergreen). While this can complicate analyses of flux measurements, it generates 292 more representative data from these types of mixed forests. Furthermore, we expect it will 293 enhance the ability of the data-driven methods for estimating domain-wide fluxes. 294 295

Fig. 3. (a) EC tower SW1 -an example of the 33 m AGL telescoping towers 296 deployed by NCAR ISFS and (b) EC instruments mounted at the top of tower 297
SW2.

299
A suite of high-quality radiation sensors was deployed in the ISS field as a complement to the 300 net radiometers installed on each flux tower. The full suite included high-quality upwelling and 301 downwelling broadband surface radiation measurements to determine the surface radiation 302 budget, as well as ancillary measurements of meteorological parameters, photosynthetically 303 active radiation (PAR), and clouds as described in Table S2. Radiation measurements are 304 manually screened and then processed through an automated data quality procedure (Long and  305 Shi calculated. Measurements of cloud properties will allow us to quantify their impacts on the 310 radiative and turbulent heat fluxes to better understand the two-way coupling between cloud-311 radiative interactions and boundary layer evolution, and to investigate the effect on EC non-312 closure. 313 314 A smaller suite of radiation, cloud, and surface meteorological measurements were deployed at 315 the Prentice and Lakeland Airports, approximately 45 km south and east from the ISS field, 316 respectively ( Fig. 1), to characterize the larger spatial scale inhomogeneities. These 317 measurements include downwelling shortwave and longwave irradiance as well as diffuse and 318 direct components of shortwave irradiance (Table S2); sufficient information to derive cloud 319 radiative effects and fractional sky cover using the Radiative Flux Analysis method described 320 above. Ceilometers deployed at the two airport sites provided additional cloud and boundary 321 layer information. 322 323 radar wind profilers by providing higher temporal and vertical resolution measurements than the 434 radars, but the radars were able to profile winds several km higher than the lidars. 435 436

Airborne Measurements
Two small unoccupied aircraft systems (sUAS) were flown to characterize surface and near-437 surface conditions (Fig. S1). During IOP1 (IOP2), a DJI S-1000 (e.g., Lee  In combination with an existing extensive database of foliar traits and image spectra (Wang et al. 471 in press), we will use the 122 foliar samples to develop and validate 1 m resolution maps for all 472 four dates of numerous foliar functional traits hypothesized to influence NPP, including LMA, 473 nitrogen concentration, chlorophyll and other pigments, phosphorus, non-structural 474 carbohydrates, fiber and lignin, and phenolics). From this, we will test the relationship between 475 functional traits and GPP (as derived from towers) and peak-season integrated NPP (early-July to 476 early-September, derived from the 41 plots). We will generate 1 m maps of NPP and GPP, and 477 identify the foliar factors that most influence each. 478 479 Additional plots were used to measure vegetation phenology as it changed through the season, 480 building upon several years of previous phenological observations collected in the domain. 481 Autumn tree leaf color and fall phenology levels were visually observed and recorded at least 482 twice weekly over six weeks during the senescence period (Sep 1 to Oct 25) for a group of 214 483 individual trees (at ten sites distributed over the 10 ×10 km area) that were representative of the 484 major species. 485 486 Forest canopy structure was characterized using an sUAS- can help attain the necessary information continuum from individual observation plots to model 558 grid scale. To achieve this, ERFs complement information across disciplines and observation 559 types by using a machine learning algorithm to find relationships between measured fluxes and 560 their meteorological and surface drivers within the flux footprints (Fig. 6A). This provides a 561 powerful approach not only for post-field data synthesis, but already in the experiment planning 562 stage e.g. in combination with Large Eddy Simulations (Fig. 6B). Over the course of the four-month study period the region exhibited light winds (diurnal means  597 from 1 -4 m s −1 ) from all directions, with the most prevalent direction being southwesterly (Fig.  598  7). Air and soil temperatures decreased over the period, while soil moisture increased (Fig. 7b,c). 599 Daily mean net radiation decreased over the course of the study, which showed a direct 600 relationship with ABL height (measured as the height of the inversion on the diurnal radiosonde 601 launches [ Fig. 7d]). One of the most relevant seasonal changes with respect to energy balance 602 was the change in the daytime Bowen Ratio (HS / HL) which averaged 0.5 in the summer and 1.0 603 in the fall, with the latter period having more variability than the former. Diurnal cycles of 604 sensible and latent heat flux show that latent heat flux is much larger in the summer when the 605 canopy is fully evapotranspiring compared to the fall, when senescence of broadleaf trees 606 reduces HL, allowing HS to comprise a larger share of the total heat flux over the region (Fig. 7f  607 -i). 608 609 610 Fig. 7. Daily mean (a) wind speed and direction, (b) temperature and relative humidity, (c) soil 611 moisture, (d) net radiation and ABL height, and (e) Bowen ratio averaged across all ISFS towers. 612 Aerial view of site NE2 on (f) July 12, 2020 and (g) October 9, 2020. Diurnal cycles of sensible 613 and latent heat averaged across all ISFS sites for the weeks of (h) Oct. 4 -11 and (i) July 7 -14. 614 Continuous data collection throughout the campaign linked the energy balance components to 616 the remotely sensed atmospheric environment (Fig. 8). As is typical for EC measurements, we 617 observed energy fluxes that were lower in magnitude than the net incoming energy (RN -G), 618 when averaged across all sites. The magnitudes of the energy balance residual (CEB) was largest 619 during the daytime, when incoming solar radiation was highest. The opposite sign of CEB from 620 day to night in part can be attributed to heat storage in the canopy. However, the magnitudes of 621 the daytime values are larger than the nighttime values, which results in a daily mean imbalance. The energy balance residual peaked under conditions of low turbulence (Fig. 9). It is during such 632 periods of calm wind and strongly unstable stratification in which thermally-induced mesoscale 633 eddies resulting from landscape-scale heterogeneity are expected (Steinfeld et al. 2007). This 634 lends support to the hypothesis that mesoscale eddies are responsible for the energy balance non-635 closure. Tower measurements, combined with in-situ measurements of air temperature and land surface 644 temperature from the DJI S-1000, were used to quantify variability in surface HS following Lee 645 et al. [2017], as shown in an example from 12 Jul 2019 (Fig. 10). On this day, as well as others, 646 there was significant temperature and HS variability; temperature (HS) differences were ~ 10°C 647 (100 W m −2 ) over the ~ 500 × 500 m area surrounding the SW2 tower. Fig. 10 illustrates even 648 finer scale resolution of surface temperature than the measures shown in Fig. 5. Such spatial 649 variation is directly related to underlying surface characteristics. Landscape heterogeneity was observed for a range of environmental variables, including 661 vegetation spectral characteristics and canopy height captured from downward-looking airborne 662 remote sensing instruments (Fig. 11). False color HySpex imagery is being used to differentiate 663 plant functional types at 1 m 2 resolution. Additional information on leaf-on canopy structure, 664 obtained from the Routescene LiDAR at 11 flux sites and across the entire domain from the State 665 of Wisconsin leaf-off LiDAR dataset, are being used to identify surface roughness in the flux 666 footprints of the EC towers. In addition, these spatial data are being used as input drivers within 667 the ERF-VCV machine learning approach. There is also spatial variation in the energy balance components across the domain on a typical 677 day (Fig. 12a). This variability includes the relative weighting of latent and sensible heat fluxes, 678 as well as the magnitude of the energy balance residual. Here we had the ability to calculate spatial fluxes from two different sources. First, the spatial 695 fluxes were calculated using a wavelet decomposition on the aircraft EC datasets. This dataset 696 has good spatial coverage but limited temporal resolution, though, with 72 flight hours spread 697 across 12 days, it is one of the largest airborne EC datasets ever collected. 698 699 The second data source for spatial EC was the set of 20 flux towers spread across the domain. 700 Calculations for flux footprints on Sep 26, 2019 (Fig. 13) show that spatial coverage of the 701 towers (including WLEF) covered roughly 8% of the domain (using Kljun et al. [2015]). This is 702 a significant increase compared to a single tower set up (typically <<1% of a 10 × 10 km area). 703 An additional benefit from the experiment design is that the towers cover a range of physical 704 environments. These data are being used to confirm the LES model results for improvements to 705 energy balance closure. By combining the tower and aircraft EC datasets we have excellent 706 coverage (~80%) of the study domain on flight days (Fig. 13). The characterization of the ABL and identification of mesoscale eddies will be performed using 717 lidar measurements of wind, water vapor, temperature, and backscatter. Figure 14 shows an 718 example of this on September 24. Increasing water vapor through the day is representative of a 719 large-scale warm, wet airmass entering the domain (Fig. 14c,d; Fig. S2a). This characterizes the 720 variation in water vapor throughout the collection of the morning UWKA CRL dataset (Fig.  721 14a). The afternoon CRL dataset (Fig. 14b) shows a more evenly mixed ABL, with variation in 722 water vapor due to local pockets of relatively moist and dry air. These two examples show the 723 varying applications of the CRL data depending on the atmospheric environment, with the 724 afternoon flight illustrating the potential of the dataset for determining the degree of ABL 725 heterogeneity arising from surface heterogeneity. Further analysis will investigate relationships 726 with underlying vegetation and LST. 727

728
Around 1200 UTC (7am local time) net radiation becomes positive (Fig. 8a) and soon after we 729 see the breakup of the surface inversion (Fig. 14d). Around 1400 -1500 UTC we see the ABL 730 grow (Fig. 8c) followed by development of large-scale structures revealed by strong oscillations 731 in vertical wind speed (±2 m s −1 ; Fig. 14e). During peak hours the angle of attack of the wind 732 vectors oscillate between roughly -30º to 50º degrees on time scales of 10 minutes to an hour. 733 These angles far exceed those of the underlying terrain, suggesting that these periodic updrafts 734 and downdrafts are the result of mesoscale eddies. 735 736 Around 1900 UTC the domain clouds over, seen in RN and backscatter ( Fig. 8a; Fig. 14f; Fig.  737 S2b). This causes the strength of the oscillation in vertical wind to decrease (Fig. 14e), which 738 coincides with a change in the relative weighting of the different energy balance components, 739 with both RN and HS decreasing strongly, while HL decreases only slightly (Fig. 8a). The data collected during the CHEESEHEAD19 field campaign show a distinct seasonal shift in 798 surface energy fluxes, as well as spatial patterning that appears to be directly related to the 799 characteristics of the underlying surface environment. Consequently, the imbalance in the energy 800 budget displays both temporal and spatial variability, with the imbalance becoming larger under 801 periods of low turbulence. The broad coverage of the measured fluxes using the 20-tower 802 network and airborne EC, combined with the collection of spatial data of surface characteristics 803 like LST, vegetation type, and canopy structure, will enable thorough investigation of the causes 804 of energy balance non-closure. Additionally, the suite of atmospheric profiling instrumentation 805 characterizes the mesoscale structure of atmospheric flows over the study domain to an 806 unprecedented degree, helping to determine how mesoscale eddies contribute to measured 807 imbalances. The observational dataset provided by CHEESEHEAD19 will also enable the use of 808 machine-learning approaches and LES for testing hypotheses on scaling and parameterization of 809 sub-grid processes in mesoscale meteorological models.