Spatial heterogeneity in CO₂, CH₄, and energy fluxes: insights from airborne eddy covariance measurements over the Mid-Atlantic region

Land cover information was taken from the National Land Cover Database (NLCD 2016), a high-resolution (30 m × 30 m) map based on Landsat imagery (Yang et al 2018). The CARAFE domain includes 14 of the 20 NLCD land classes. Dominant land classifications sampled during CARAFE include woody wetlands (45%), cultivated crops (22%), and dry forests (evergreen, deciduous, and mixed classes) (21%). The remaining types are developed land (open, low, medium, and high density), open water, emergent herbaceous wetlands (hereafter herbaceous wetlands), shrubs, pastures, and grasslands, which individually make up less than 5% of the cumulative footprint.

The flux footprint relates the spatial distribution of fluxes at the surface (${F}_{s}$) to the observed flux (${F}_{obs}$) measured at coordinates ${x}_{m},\,{y}_{m}$ and measurement height ${z}_{m}$ (Horst and Weil 1992, Schmid 1994):

$$\begin{eqnarray}&&{F}_{obs}\left({x}_{m},{y}_{m},{z}_{m}\right)=\displaystyle \int \displaystyle {\int }_{-\infty }^{\infty }f\left(x,y,{z}_{m}\right){F}_{s}\left(x,y,0\right)\,dx\,dy,\end{eqnarray} \tag{ 1 }$$

where, x and y are arbitrary horizontal coordinates and f is the flux footprint function, which expresses the contribution of each upwind unit surface element to ${F}_{obs}.$ We use the two-dimensional Flux Footprint Prediction (2D-FFP) developed by Kljun et al (2015), a parameterization based on a Lagrangian stochastic particle dispersion model (Kljun et al 2002) that is applicable to many turbulence regimes and measurement heights. The parameterization utilizes the following inputs: measurement height ${z}_{m},$ the mean horizontal wind speed U, the planetary boundary layer height ${z}_{bl},$ the Obukhov length ${L}_{OB},$ the standard deviation of the lateral (crosswind) velocity fluctuations ${\sigma }_{v},$ and the friction velocity u*. We derive ${\sigma }_{v}$ from the wavelet variances of the horizontal wind velocity vectors.

We calculate the 2D-FFP for all 1 Hz data points (∼75 m distance at typical flight speed) along all flux transects below 200 m. Each data point has an associated ${z}_{m},$ but leg-average values of the micro-meteorological variables (i.e. $U,\,u* ,\,{\sigma }_{v},\,{L}_{OB}$) are used as the FFP is based on a mean flow parameterization, and point-to-point momentum fluxes exhibit significant variability. We estimate the boundary layer height from vertical profiles before and after each set of flux transects as described in Wolfe et al (2018). Note that even a 20% error in ${z}_{{\rm{bl}}}$ has less than a 0.5% impact on the size and distribution of the footprint, except in highly stable conditions (Kljun et al 2015) atypical during the CARAFE flights. Once calculated, the 2D-FFP was rotated into the mean wind direction and transformed into the geographic coordinate space of the measurement, generating a gridded map of the footprint function associated with each flux observation. Figure 2 depicts an example of a single footprint for a flux measurement from the 18 May, 2017 flight to Choptank, MD superimposed on the NLCD 2016 land cover map. For all flux observations from the 2016 and 2017 campaigns below 200 m in altitude, the 90% upwind extent for calculated footprints ranged from 1.5 to 10 km.

Zoom In Zoom Out Reset image size

To derive fluxes representative of a single land class, we use the Disaggregation combining Footprint analysis and Multivariate Regression (DFMR) methodology described by Hutjes et al (2010). DFMR relies on the flux footprint and land cover to estimate a weighted contribution of each land class to the flux measurement. The observed flux, ${F}_{obs},$ can be written as a linear summation of component fluxes from each of n land classes within the footprint:

$$\begin{eqnarray}&&{F}_{obs}=\displaystyle \sum _{k=1}^{n}{C}_{k}{F}_{k},\end{eqnarray} \tag{ 2 }$$

where ${C}_{k}$ is the fractional area of the kth land class within the footprint and ${F}_{k}$ is the corresponding component flux, or the mean land-class flux for a given set of observations (i.e. a single flight). The values ${C}_{k}$ can be determined using the flux footprint function f to weight the relative contributions of land cover patches within the footprint, as patches closer to the sensor influence the measurement more heavily than patches farther away (equations (S1) and (S2)).

The multiple linear regression (equation (S3)) was performed on a flight-by-flight basis to derive land-class component fluxes for each flight. A grid of 2D-FFP values was superimposed onto the NLCD 2016 map to generate the weighted fractional area of each land class in every footprint (see figure 2). Although NLCD displayed 14 land classes in our sampling region, we down-selected for land classes that constituted more than 20% of the footprint-weighted area in at least 4 km of cumulative (but not necessarily consecutive) flux observations. This screening criterion, which was optimized via comparison with flux sub-samples from homogeneous footprints (figure S1), ensured that selected land types were sampled enough to provide a meaningful average.

In addition to residual error, random and systematic measurement errors (see Wolfe et al 2018) were propagated through the regression, and errors were summed in quadrature to yield the total uncertainty for each component flux. Uncertainties in a priori surface characterization and footprint extent are not included. While the footprint calculation should introduce minimal error (barring significant changes to the fractional areas), mischaracterization of the surface cover could introduce significant biases. Corrections for vertical flux divergence are likewise not included due to large uncertainties in the correction factors (see Wolfe et al 2018). The flights all took place near midday and targeted fair-weather conditions. However, the data are not screened for the presence of clouds, and such variations may contribute another source of flux variability in addition to those discussed below.

CARAFE deployments included three flights to the Choptank agricultural area (12 September 2016, 4 May 2017, and 18 May 2017). Flux transects spanned the Delmarva peninsula from the Chesapeake Bay to the Atlantic Ocean (typical length 60 km) and included overflights of the USDA-Chop and US-StJ towers (figure 1). This region had mixed terrain, with six land classes meeting the down-selection criteria described in section 2.3.3. Table 2 summarizes footprint-weighted contributions of each land class.

Disaggregated fluxes highlight the variability in carbon dynamics between land classes and over time (figure 3). Of particular interest, cultivated crops (e.g. annual crops such as soybean or corn) and forested lands (e.g. woody wetlands and deciduous forest) display substantial differences in F_CO2 for the sampling periods. The CO₂ uptake from cultivated crops ranged from −3.4 ± 0.7 to −11.5 ± 1.6 μmol m⁻² s⁻¹, whereas deciduous and wetland forests display a much larger uptake, ranging from −12.1 ± 3.9 to −22.7 ± 3.9 μmol m⁻² s⁻¹ and −9.1 ± 1.5 to −22.7 ± 3.2 μmol m⁻² s⁻¹, respectively. Forest uptake of CO₂ also dominates that by croplands in other regions with substantial cropland fraction (figures S2, S3). While the difference in CO₂ uptake between croplands and forest will be strongly dependent on crop type and phenology (Lokupitiya et al 2009), crops in the CARAFE region are typically in their early growth stages in May and undergoing senescence in September. Developed open lands, which comprise mostly lawn grasses and vegetation with <20% impervious surface area, also draw down substantial CO₂ (∼−13 to −30 μmol m⁻² s⁻¹) during the May sampling period.

Zoom In Zoom Out Reset image size

Choptank data also exhibit a general anticorrelation between F_CO2 and LE, expected for vegetated land surfaces where transpiration and stomatal control is a major contributor to evapotranspiration. The sampling is too limited to infer much about seasonal flux response for the various land types, but forested lands are comparably photosynthetically active during the growing season between May and September.

The disaggregation methodology also illustrates F_CH4 variability with land type. F_CH4 observations were at or below the detection limit for most CARAFE flights, and uncertainties are large due to poorly constrained regression results. However, it is known that soils from forested ecosystems represent a weak CH₄ sink (Subke et al 2018), whereas tree stems represent a weak CH₄ source (Vargas and Barba 2019) that may counterbalance ecosystem scale CH₄ fluxes in upland forested ecosystems. In contrast, herbaceous wetlands, located primarily on the Eastern end of the track near the St. Jones tower, exhibit relatively strong CH₄ emissions of 76.1 ± 29.4 nmol m⁻² s⁻¹ on Sep-12 and 32.3 ± 17.0 nmol m⁻² s⁻¹ on May-04. This region has a mix of herbaceous wetlands that extend across a salinity gradient, where lower CH₄ emissions may be associated with wetlands in brackish waters and larger CH₄ emissions with freshwater wetlands (Poffenbarger et al 2011, Capooci et al 2019).

Three flights over the Alligator River region took place during the CARAFE deployments on 24 September 2016, 15 May 2017, and 26 May 2017. Flux transects spanned the Alligator River National Wildlife Refuge in the N–S direction, with the US-NC4 tower located near the middle of the flight transects (see figure 1). The dominant land cover contributions included woody wetlands, open water, and some minor areas of cultivated crops and herbaceous wetlands. Table 2 contains a summary of the land-cover contributions to the footprint for Alligator River region.

The component fluxes from Alligator River display significant variability with land type (figure 4). For example, the open water component of H is at or near zero for all flights, and evaporation dominates the surface energy fluxes for this class, as displayed by LE values generally greater than 200 W m⁻². Note that although classified as open water in the NLCD, the coastal waters sampled near Choptank and Alligator River are actually comprised of estuarine waters and tidal mudflats. We observed occasional CO₂ emissions from these waters of 6.2 ± 2.2 μmol m⁻² s⁻¹ over the Alligator River (Sep-24, figure 4) and 6.4 ± 3.7 μmol m⁻² s⁻¹ in the Choptank region (May-18, figure 3). Both regions also display positive fluxes of CH₄, with means of 40.5 ± 12.2 nmol m⁻² s⁻¹ over the Alligator River and 20.5 ± 10.1 nmol m⁻² s⁻¹ in Choptank. These values are within the range of prior flux estimates, which can be up to ∼10 μmol CO₂ m⁻² s⁻¹ in low salinity estuarine waters, and ∼30–35 nmol CH₄ m⁻² s⁻¹ in tidal mudflats (Abril and Borges 2005).

Zoom In Zoom Out Reset image size

The woody wetlands land class, a freshwater forested swamp in the Alligator River region, also displays persistent, large CH₄ emissions ranging from 58.4 ± 12.0 to 181.2 ± 36.8 nmol m⁻² s⁻¹. These remarkably high values are recurrent and consistent over long periods at this site (Mitra et al 2019a), which has recorded among the highest CH₄ emissions from wetlands globally, including tropical wetlands. Other temperate swamplands exhibit mean CH₄ emissions of ∼35 nmol m⁻² s⁻¹ (Turetsky et al 2014), but few global observational records exist for wetland ecosystems.