Detecting Hot Spots of Methane Flux Using Footprint-Weighted Flux Maps.

In this study, we propose a new technique for mapping the spatial heterogeneity in gas exchange around flux towers using flux footprint modeling and focusing on detecting hot spots of methane (CH4) flux. In the first part of the study, we used a CH4 release experiment to evaluate three common flux footprint models: the Hsieh model (Hsieh et al., 2000), the Kljun model (Kljun et al., 2015), and the K & M model (Kormann and Meixner, 2001), finding that the K & M model was the most accurate under these conditions. In the second part of the study, we introduce the Footprint-Weighted Flux Map, a new technique to map spatial heterogeneity in fluxes. Using artificial CH4 release experiments, natural tracer approaches and flux chambers we mapped the spatial flux heterogeneity, and detected and validated a hot spot of CH4 flux in a oligohaline restored marsh. Through chamber measurements during the months of April and May, we found that fluxes at the hot spot were on average as high as 6589 ± 7889 nmol m-2 s-1 whereas background flux from the open water were on average 15.2 ± 7.5 nmol m-2 s-1. This study provides a novel tool to evaluate the spatial heterogeneity of fluxes around eddy-covariance towers and creates important insights for the interpretation of hot spots of CH4 flux, paving the way for future studies aiming to understand subsurface biogeochemical processes and the microbiological conditions that lead to the occurrence of hot spots and hot moments of CH4 flux.

Introduction

Exploring the spatial heterogeneity in gas exchange around flux towers is a necessity for an increasing number of flux sites, particularly wetland sites, where the occurrence of “hot spots” of methane (CH4) flux adds a major uncertainty to the global CH4 budget (Morin, 2019). Wetlands are important carbon sinks in the biosphere (Mitsch et al., 2013), however, they also emit about a third of all the methane emissions to the atmosphere and the processes leading to these high methane emissions remains understudied (Kirschke et al., 2013). A large source of uncertainty in the methane budget is the existence of “hot spots”, defined as patches in the landscape with disproportionally higher rates of biogeochemical activity that can contribute disproportionally to the total fluxes in the ecosystem (McClain et al., 2003). Although the existence of hot spots can be associated with the underlying heterogeneity in land cover composition (Bansal et al., 2020; Fischer et al., 2010; Forbrich et al., 2011; Matthes et al., 2014; Morin et al., 2017; Sachs et al., 2010), there is still a gap in the understanding of when and where these hot spots occur and what are the subsurface processes that lead to these exceedingly high methane fluxes.

Traditionally, eddy‐covariance towers, which measure methane and other scalar fluxes, are expected to be located on sites with homogeneity in underlying sources/sinks, canopy height, and roughness (Lee, 2018). However, a great number of eddy‐covariance towers are being located on heterogeneous land and/or on sites where the underlying sources/sinks may not be completely homogeneous. In such cases, flux footprint analysis is instrumental to determine the representativeness of eddy‐covariance measurements (Chu et al., 2021), to quantify the sensor location bias over inhomogeneous surfaces (Schmid and Lloyd, 1999), and to upscale fluxes from the tower to the heterogeneous landscape (Chen et al., 2012; Kim et al., 2006).

The flux footprint can be defined as the source (or sink) area of a flux observation from eddy‐covariance measurements, and its location and extent depend primarily on tower height, wind direction, turbulence, atmospheric stability, and roughness length (Schmid, 2002). The footprint function or footprint model can also be defined mathematically as the transfer function between those sources/sinks and the flux measured by the tower (Kljun et al., 2015; Vesala et al., 2008). Footprint models can be derived through several approaches, including Lagrangian stochastic (LS) models, large eddy simulations (LES), or Eulerian analytical solutions. LS models have been applied in multiple studies (e.g., Baldocchi, 1997; Flesch et al., 1995; Kljun et al., 2002; Leclerc and Thurtell, 1990). They track particles back in time from the receptor location to the surface source or sink and are one of the most accurate approaches to represent the flow around eddy‐covariance towers, especially when there are obstacles or sharp changes in surface heterogeneity. . However, analytical solutions are often preferred for estimating average spatial fields because they simplify the physics of atmospheric transport, are less computational costly (Schmid, 2002) and because they have been largely incorporated into regular eddy‐covariance processing steps (Göckede et al., 2008; Leclerc and Foken, 2014). In recent decades, the evolution of analytical footprint models has been presented in multiple reviews (Leclerc and Foken, 2014; Schmid, 2002; Vesala et al., 2008). In short, it is worth noting that analytical approaches are based on the inverted plume assumption of Calder (1952) and the solution of the advection‐diffusion equation by van Ulden (1978), which have been progressively improved by the work of multiple authors (Gash, 1986; Haenel and Grünhage, 1999; Horst and Weil, 1992; Schmid and Oke, 1990; Schuepp et al., 1992). From this work stem some of the most commonly‐used analytical flux footprint models that are the focus of the present paper, namely: the Hsieh model (Hsieh et al., 2000), and the K&M model (Kormann and Meixner, 2001), and the Kljun model (Kljun et al., 2015), which is a dimensionalized parameterization of a backward Lagrangian stochastic dispersion model (Kljun et al., 2004).

With appropriate footprint models, data integration from multiple half‐hourly footprints can be used to resolve the source heterogeneity at a detailed level. Land cover heterogeneity can have a strong influence on the magnitude and interpretation of flux measurements (Chen et al., 2012; Griebel et al., 2016; Lloyd, 1995) and can systematically add variability to the upscaling of ecosystems fluxes measured by eddy covariance (Griebel et al., 2020). A way to estimate how strong is the spatial heterogeneity in the flux contributions of different areas around the tower is to aggregate footprints over a given period of time (i.e., create footprint climatologies (Amiro, 1998)) and overlap them against land cover maps (Amiro, 1998; Chen et al., 2009; Kim et al., 2006; Morin et al., 2017; Tuovinen et al., 2019). Footprint climatology studies can also be used to evaluate the spatial distribution of the quality of different fluxes and calculate spatially‐resolved values of roughness length (Göckede et al., 2006). When used at the half‐hourly level, footprints can also be used to create a time series of half‐hourly contributions from different land covers to the total flux (Chen et al., 2009; Giannico et al., 2018; Göckede et al., 2006; Morin et al., 2017; Rey‐Sanchez et al., 2018). However, most of these approaches derive the information on the spatial heterogeneity of fluxes from the features detected via remote sensors. Although this is a valid method, relying solely on this assumption ignores potential variability driven by subsurface processes. Amiro (1998) was one of the first researcher who first used a footprint aggregation approach to create a land‐cover‐independent surface flux map. Amiro (1998) achieved this by weighing the footprint function by the value of evapotranspiration (ET) and thus obtained a general representation of which regions around the tower had the highest and lowest ET. Similarly, Chen et al. (2009) obtained footprint‐weighted averages of carbon dioxide (CO2) flux by multiplying the flux at each half hour by the footprint and then normalizing it to a footprint climatology of interest. These authors have evaluated surface heterogeneity as a function of footprint representativeness but the potential variability in surface fluxes has not been captured fully. In principle, the mapping of the spatial variability of fluxes around eddy covariance towers can be done by using the superposition of flux footprints located on multiple grid elements on a surface grid, as applied for airborne measurements (Mauder et al., 2008; Schuepp et al., 1992). This approach calculates surface flux maps or “flux topographies” (Kohnert et al., 2017, 2018; Mauder et al., 2008) that can be applied to detect potentially high heterogeneity of surface fluxes around eddy covariance towers. In this study, we hypothesize that by slightly modifying the approach used in these airborne measurements, the spatial variability around stationary eddy‐covariance towers can be mapped, and we call this approach the Footprint‐weighted Flux Map.

In this study, we test, validate and demonstrate the applicability of the footprint‐weighted flux maps using three methods. We start by validating three commonly used flux footprint models using a tracer release experiment. Although validation using trace gas release experiments like this exist, few have focused on daytime versus nighttime conditions and the detection of long‐range hot spots of methane flux. Using the model with the best performance, we then evaluate the use of footprint‐weighted flux maps to detect the spatial variability of fluxes in different sites, including (a) an alfalfa field with a methane gas release experiment, (b) a site with sharp transitions in land cover, and (c) wetlands with different mosaics of water and vegetation in the Sacramento‐San Joaquin River Delta, California.

Methods

Study Sites and Tracer Gas Release Experiments

Tracer Gas Release Experiment

To evaluate the accuracy of footprint models, we performed three tracer gas release campaigns in an alfalfa field located in the San‐Joaquin Sacramento River Delta in California (Ameriflux site: US‐Bi1). The site is ideal for validating footprint models because it is flat, homogeneous, with steady winds, long fetch, short canopy (Figure 1), and with a background CH4 flux near zero. The site presents prevailing westerly winds channeled into the Delta through the Carquinez Strait and driven by strong temperature gradients between the cooler San Francisco Bay and the warmer Delta. The tracer gas release experiments consisted of releasing 99.9% methane through a PVC pipe, 3.01 m long, 8.89 cm in diameter, and with 5 outlets regulated by a gas flow meter with apertures of 0.635 inches each. The pipe and outlets were positioned just above the canopy, as described below, and blew against the main wind to facilitate turbulent dispersion. The receptors for the released methane were two eddy‐covariance systems equipped with CH4 infrared gas analyzers (Li‐7700, Licor Biosciences, Lincoln, NE) located in the same tower but at two different heights, 3.9 and 2.1 m above ground (Figure 1c). More details on the data collection and processing from these towers are presented in the next section.

Figure 1

Setup of the tracer release experiment. (a) Location of the study site showing the location of the eddy covariance towers used in this study and their respective Ameriflux codes. (b) Aerial view of the daytime and nighttime locations, at 20 and 180 m, respectively, and the three exploratory locations for the initial daytime releases at 20, 40 and 60 m. (c) Photo of the eddy‐covariance setup with methane analyzers and sonic anemometers at 3.9 and 2.1 m. An additional sonic anemometer was located at 0.6 m.

We performed three CH4 release campaigns using different flow rates and distances between the source and the tower. In the first release, the exploratory campaign, the pipe was located at a 312° azimuth from the tower, which was the prevailing wind direction at the time, and was moved at distances of 20, 40, and 60 m and at variable flow rates (Figure 1b). This information allowed us to determine the use of a flow rate of 1.75 lpm, which was the optimal flow rate for the subsequent daytime experiment. For the second campaign, the release pipe was located at 20 m from the tower at an azimuth of 253° to try to match the forecast wind direction. This release point was activated during 3 days in roughly daytime hours (06:00–21:00 LT). The third release occurred three days after the end of the first release and was set at 180 m from the tower, at an azimuth of 253°, and activated only during nighttime hours (18:00–09:00 LT). For the second and third releases, a secondary tower was deployed at a lower height allowing us to assess flux divergence to bound the conservation of mass. We also obtained a second data set from this tower to validate the fluxes, doubling the number of points for the validation of footprints (Figure 1c). Additionally, we installed a sonic anemometer at a third level, 60 cm above the ground to obtain a three‐point wind profile and thus improve the current estimates of roughness length and displacement height necessary for these footprint models (Figure 1c).

Other Wetland and Crop Sites

We also analyzed data from three sites in the meso‐network of eddy‐covariance towers in the Sacramento San Joaquin River Delta, in California. The first site analyzed was the West Pond wetland (US‐Tw1), a restored wetland with a dense layer of cattails (Typha spp), and Tules (Schoenoplectus acutus), where micrometeorological measurements have been taken for almost 10 years. This site was selected because it presents a sharp land‐use transition between the nearby wetland in the predominant western wind direction and an alfalfa field toward the north. The main wetland site was the Sherman Wetland site (US‐Sne), a wetland restored in 2017 with large spatial variability in the proportion of vegetation patches to open‐water patches (Valach et al., 2021). The vegetation in this wetland is dominated by cattails and tules and common reeds to a lesser extent. An auxiliary site was the Mayberry wetland (US‐Myb), an 11‐year‐old wetland that sees a patch of cattails, tules, and common reeds (Phragmites spp) vegetation in the predominant westerly wind direction, but that in the north sees a channel of open water.

Data Handling and Processing

Release Experiment

Data for the release experiment was processed in 15 min intervals to maximize the number of points for validation. It was expected that low‐frequency losses from shorter averaging times of 15 or even 10 min would be negligible (Heidbach et al., 2017), and an analysis of the differences between 15 min and half‐hourly data demonstrated that there were no significant changes in the cospectra of methane and the vertical wind speed component between 15 and 30 min datasets (Figure S1 in Supporting Information S1). We also calculated the variance in methane concentration at each of the two‐time integrations. We did it as a way to evaluate the measurement frequency that would result in the least variability due to the meandering nature of the concentration from the trace release. Because the release was characterized by high fluctuations in the concentration of CH4, we increased the variance threshold allowed for CH4 concentrations (to 1.5 (μmol m−3)2), as well as the maximum CH4 flux threshold (from 800 to 4000 nmol m−2 s−1). Other authors have eliminated the despiking routines altogether (e.g., Heidbach et al., 2017). We applied a double rotation of the wind vectors to force the mean lateral and vertical wind components to zero. When the original vertical wind direction exceeded a threshold of 6°, the observation was discarded. We evaluated non‐stationarity in the time series using the approach of Foken and Wichura (1996). This stationarity test, commonly used for surface fluxes, needs to be revised for point‐source emissions (Dumortier et al., 2019), and since the fluxes were computed at 15‐min intervals, we eliminated the stationarity filter proposed by Foken & Wichura (1996) for CH4 concentration.

Processing of Other Flux Data

Data from the 3 sites described in Section 2.2.1 were analyzed in the standard 30‐min interval following common procedures described elsewhere (Detto et al., 2010; Hatala et al., 2012; Knox et al., 2015). In short, at each site, a 3‐D sonic anemometer (WindMaster Pro 1352 or 1590, Gill Instruments Ltd, Lymington, Hampshire, England) was used to measure high‐frequency fluctuations in the vertical, lateral, and along‐wind components. CO2 and H2O molar densities were measured with infrared gas analyzers (LI‐7500DS, LI‐COR, Lincoln, Nebraska, USA), while methane concentration was measured with a LI‐7700 (LI‐COR, Lincoln, Nebraska, USA). Each set of eddy covariance instrumentation collected data at 20 Hz. We employed a 3‐D coordinate rotation to force the main wind to the along‐wind coordinate (u), and force the mean lateral (v) and the mean vertical (w) wind to zero. Stationarity was evaluated using the approach of Foken and Wichura (1996), using a flag of 8 as the maximum threshold for data quality. Half‐hourly covariances were processed using the approach of Webb Pearman Leuning (Webb et al., 1980), to account for density fluctuations on the calculation of the true surface‐air flux. Data were filtered out based on a minimum friction velocity threshold of 0.2 m s−1 and on a minimum acceptable vertical rotation angle of −6 and 6 degrees.

Footprint Models and Their Parameters

Under the assumption of turbulence being horizontally homogeneous, the flux density measured via eddy‐covariance above the origin of the coordinate system, and at a height z = z , can be defined as (Lee, 2018; Schmid, 2002): where is the vertical flux of the scalar of interest in kg m−2 s−1, Q is the surface flux in kg m−2 s−1, and f is the footprint function in m−2. We tested three commonly used footprint functions (models), the Hsieh model (Hsieh et al., 2000), the Kljun model (Kljun et al., 2015), and the K&M model (Kormann and Meixner, 2001). In addition to their widespread use, these models were chosen because they can produce 2‐dimensional footprint fields and consider the stability effects on the plume width. The total flux footprint in 2‐dimensions was defined as the product of the cross‐integrated footprint function () and a crosswind dispersion function (Amiro, 1998; Detto et al., 2006): where the footprint function, , is a function of the effective height (, tower height minus displacement height), the crosswind distance from the tower (y), and the streamwise distance (x). is the cross‐integrated footprint, and is the crosswind dispersion, which is usually calculated as a Gaussian function: where σ is the standard deviation of the plume in the y dimension in response to distance in the streamwise direction, x, and can be calculated from the standard deviation in lateral wind fluctuations as in (Detto et al., 2006): where and are similarity parameters equal to 0.3 and 0.86, respectively, is the friction velocity and is the roughness length for momentum. In the case of the KL model, which is based on a scaling approach, crosswind dispersion is also Gaussian, and is equal to: where is a non‐dimensional standard deviation of the crosswind distance, and is a proportionality factor (Kljun et al., 2015).

The cross‐integrated footprint function for the Hsieh model is equal to: where, k is the von Karman constant (0.4), x is the streamwise distance from the tower, z is a length scale (Hsieh et al., 2000) and D and P are similarity constants whose values change for different regimes of atmospheric stability. For the K&M model the cross‐integrated footprint is equal to: where is a flux length scale, and is equal to (1 + m)/r. The variable m is the exponent of the wind velocity power law, and n is the exponent of the eddy diffusivity power law. The Kljun model stems from a parameterization of a backward Lagrangian particle model tested for specific stability conditions and is applicable for different regimes of atmospheric stability and includes multiple fitting functions described in (Kljun et al., 2015). The Hsieh and K&M models were coded in Matlab, while the Kljun code in Matlab was downloaded directly from the author's website (https://footprint.kljun.net/).

After applying the 2‐D dispersion model to the cross‐integrated footprint model, the footprint models are run on grids of 3 × 3m (5 × 5 m for higher z e.g., during nighttime conditions) with the tower in the center of a grid, and the weight of the footprint distributed for each of the cells in the grid. The sum of the footprint weights in all the cells of the grid approximates to unity, more commonly to 0.9 given the incapability of the model to capture 100% of the sources. Due to this impossibility to calculate 100% of the sources, the footprint weights were normalized by diving them by the sum of the original footprint weights (i.e., 0.9 or 0.8). Once the gridded footprint weights are created, footprint contours are delimited using the function “contour” in Matlab. The maximum footprint contour was equal to 80% given the last 20% of the cumulative probability distribution has an exceeding long tail, but other contours, including 50%, 60%, and 70%, were also calculated.

Calculation of Roughness Length, Canopy Height, and Displacement Height in the Release Experiment

In the alfalfa field, canopy height was measured directly and estimated as the aerodynamic canopy height using the algorithm of Pennypacker and Baldocchi (2016). The difference between aerodynamic canopy height and measured height was minimal, so for the final footprint calculations, we used the average aerodynamic canopy height, expecting it to be more representative of the area that the tower observes and captures small variability in the height of alfalfa (Chu et al., 2018).

We used a 3‐level wind profile to calculate roughness length (z ). The wind profiler, which was available for the second and third campaigns, consisted of wind velocity measurements at three heights, Z = 3.9, 2.1, and 0.6 m. To calculate roughness length, we fitted the logarithmic wind profile on multiple observations under neutral conditions and selected the median of all the z calculated. For the first release, where the wind profiler was not available, we used the approach of (Maurer et al., 2013) to calculate roughness length. The K&M and Kljun calculate a z internally, while z needs to be prescribed for the Hsieh model. However, if z is known, the Kljun model would preferentially use this value in replacement of wind speed for computations. Therefore, we decided to use our internal calculation of roughness length for both Kljun and Hsieh model.

Calculation of Roughness Length, Canopy Height, and Displacement Height in the Other Wetland Sites

For calculating the canopy height in other sites, including wetland sites where direct measurements of canopy height are not technically feasible, we exclusively used the algorithm of Pennypacker and Baldocchi (2016). In non‐inundated fields such as the alfalfa crops, the tower height was constant, while in wetland sites, the tower height can change in response to increases in water level above the surface. However, it is important to note, that since the plants in our wetlands are rooted in the sediments, their height does not change with water level, so the distance between the instruments and the top of the vegetation is independent of water height. One could then propose eliminating the need to look at the water level to evaluate canopy height. However, because the Pennypacker and Baldocchi (2016) algorithm evaluates the effect of both the canopy and the water level on the shape of the wind profile, it is important to consider both effects. Given this situation, it is necessary to have a dynamic tower height when using the algorithm to better account for the dynamic canopy changes (Chu et al., 2018). Roughness length for the wetland sites was estimated as 0.1 of canopy height while displacement height was equal to 0.6 of canopy height.

Footprint Validation

Considering a grid cell array of m by n cells, a numerical discretization of Equation 1 can be given as where F is the flux calculated based on footprint modeling, which should be close to the flux measured by the eddy‐covariance, F () if the footprint model is accurate and if all the sources are estimated precisely. f are the values of the footprint function for cell i,j, and F is the source/sink flux emitted by the tile i,j. Finally, is the size of the grid cell used in the analysis.

The approach to validate the three footprint models consisted of overlaying the values of the footprint function on grids of m by n cells with a size m2 (where in our case and varied from 3 to 5 m). Such grid would include the 2‐D dispersion model such that the sum of all the grids within the 80% footprint would be equal to 0.8 (Morin et al., 2014).

If we imagine a field of homogeneous sources in methane flux, such as is the case for the alfalfa field, all the tiles within the grid emit a near‐zero flux resulting in an effective zero‐flux measured at the tower. Now, if we add a known source of gas in one of these tiles, i = r, j = t, the flux measured by the tower will be modified and will be equal to:

Given that the background flux is zero, we can assume that the flux measured by the EC system would be equal to: where is the value of the footprint function at cell i = r, j = t, and F is the flux from the release at cell i = r, j = t. The flow of CH4 in L s−1 calculated based on the flow rate of the release can be converted into a flux using ideal gas law and normalizing the flow to a source area of m2: where is the flow rate of pure CH4 in L s−1, R is the gas constant (8.314462 Pa L K−1 mol−1), T is the temperature in K, P is the atmospheric pressure in Pa, and A is the area () in m2, to which the flux is being normalized, and F is the resulting flux in mol m−2 s−1. In our case, the line source of the release had a length of 3 m, which we then normalized to an area of 3 m by 3 m or 5 m by 5 m, depending on the size of the grid. Based on this principle, a footprint validation model can be created, where different source functions will create different values of F that can be compared against the flux measured by the eddy‐covariance system (F ) to find the best model.

Footprint‐Weighted Flux Maps

By observing multiple footprints in a cell array with the flux tower in its center, we can calculate a footprint‐weighted flux following a similar approach to that of Schuepp et al. (1992) and Mauder et al. (2008) as applied to airborne measurements. This footprint‐weighted flux for cell i,j (Fw ) can be calculated as: where F is the half‐hourly eddy‐covariance flux for a given scalar (e.g., CH4, CO2, H2O, sensible heat flux) at time m. M is the number of half hours for the period of interest for which data are available, and f is the value of the footprint function at cell i,j, for the m run. With this approach, a footprint‐weighted flux map can be created that indicates the average flux that the flux tower measures when the footprint intersects a given grid cell. This value can give us an indication of the spatial heterogeneity of fluxes around the EC tower. We call this approach the footprint‐weighted flux map.

Use of Footprint‐Weighted Flux Maps

To test the ability of footprint‐weighted flux maps to capture the spatial heterogeneity in fluxes around flux towers, we tested the technique on several sites differing in landcover and landscape configuration. We first used data from the alfalfa field, where CH4 fluxes were near zero and where the tracer gas release of CH4 was deployed. Here, we evaluated the capacity of the technique to detect a short‐range (20 m from the tower) and a long‐range (170 m from the tower) hot spot of CH4 flux. We then tested the technique on West Pond, a small 3 ha wetland located between an alfalfa and a corn field in the NW and SW, respectively, thus showing land‐use transitions across the main wind direction. Finally, footprint‐weighted relative flux maps were applied to Mayberry and Sherman wetland, large wetland sites of unknown heterogeneity in methane flux.

To create a footprint‐weighted flux map it is important to minimize the temporal variability in fluxes within the analysis timeframe. In line with this,we ran footprints analysis for day and night times separately to minimize diurnal variations associated with temperature responses and over short periods of time (30 days approximately), to minimize temporal variations associated with seasonal changes. For the detection of hot spots of methane flux, we also evaluated the time series of CH4 fluxes and identified those hot moments or clusters of high CH4 flux that deviated from the normal diurnal cycle in CH4 emissions and calculated a footprint‐weighted flux map for a period of 3–4 weeks centered around that event. Another consideration of the application of this tool is the use of 80% footprints in the creation of the surface flux maps. Although 80% could neglect the contributions from distant sources, it accomplishes a good compromise in the signal‐noise ratio at a spatial level. Those distant sources have a larger uncertainty associated with a much larger area responsible for a minimum amount of the signal. Footprint‐weighted flux maps were created using heat maps in Matlab, and color maps were obtained from https://github.com/bids/colormap.

Estimation of the Magnitude of the Hot Spot Flux

A footprint‐weighted flux map only indicates where the fluxes are higher or lower relative to other positions in a map. To estimate the actual flux from a potential hot spot, it is necessary to go a step further and create an additional map, the footprint‐derived hot spot map. Considering a hypothetical scenario where a hot spot of methane flux is identified, the area of this hot spot can be estimated. Using the estimate of the background methane flux, we can then provide a first‐order estimate of the flux from the hot spot. The flux from the hot spot grid cell or cells () can be calculated based on the cumulative footprint values as: where is the footprint‐weighted methane flux in the area of interest (i.e., the average flux measured by the flux tower when the footprint intersects the area of the hotspot), is the sum of the background footprint contributions, which should be closer to 1 for small hotspot areas (such as the release experiment in this study, where F only covers 1 cell), is the background methane flux, which in the case of the alfalfa field can be approximated to zero, but in other cases can be taken as the average of the footprint‐weighted flux for the cells outside of the hotspot. Finally, is the footprint contribution from the cell or cells where the hot spot is located, which in the case of the release experiment should be much smaller than 1 given that these are contributions only from one grid cell.

Chamber Fluxes

Non‐steady‐state static chambers (Subke et al., 2021) were used to measure fluxes of CH4 and CO2 at several locations in the Sherman wetland site to corroborate the existence of a CH4 hotspot detected using footprint‐weighted flux maps. We performed chamber measurements at the hot spot location and contrasted them against measurements at the point of the highest footprint weight (i.e., the spot that has the highest contribution to the total flux), which was an open water area. We performed two chamber campaigns, one in April and one in May 2021. In April, 13 individual observations were taken in the open water, and 22 in the hot spot location. In May, 9 and 21 individual observations were taken, respectively. To prevent investigator disturbance, chamber collars and an elevated boardwalk were installed 1 week before the first sampling event and left in place at the hot spot location, and a small boat was used to measure fluxes at open water areas. Measurements were taken between 10 a.m. and 4 p.m.

In the open water, we used closed floating chambers made of acrylic with dimensions of 28 cm width by 28 cm length and a height of 11.4 cm above the water level. At the hot spot location, floating chambers were used where standing water was higher than ∼8 cm. At sites where standing water was at or close to the surface, collars (28 × 28 × 15 cm depth) were inserted ∼5 cm into the soil and left to settle for 1 week before the first set of measurements were made and left in place for subsequent measurements. The upper end of the chamber collars was outfitted with a groove that held water during sampling and thus created a seal from the outside air when chamber tops were inserted. All chambers were equipped with a mixing fan to aid mixing, and a thermocouple for air temperature measurements inside the chamber. The chambers were connected to an infrared gas analyzer (LI‐7810, LI‐COR Biosciences, Lincoln, NE) in a closed‐path system which sampled the chamber concentration of CO2 and CH4 at a rate of 1 Hz. Once the chambers were settled on the surface, we waited for about 1 min for the entire system to equilibrate. After recording the start time, gas flux measurements were made for approximately 2 min, then the end time was recorded. Diffusive fluxes of CH4 and CO2 were calculated based on the slope of the rate of change in gas concentrations with time. Gas molar mixing ratios were converted to molar densities using the ideal gas law. Where present, we accounted for ebullitive fluxes of CH4 using the approach of Villa et al. (2021). In short, an ebullition event is identified based on an excessively high increase in the concentration of the gas. The beginning and end times for the ebullition event were identified, and a new flux was calculated based on the slope of gas concentrations with time during that period. This new flux is normalized based on the proportion of the time of the ebullition event to the period of time of the whole chamber measurement. Total CH4 flux was calculated as the sum of the ebullitive and diffusive fluxes.

Results

Foundational Data Set in the Methane Release Experiment

A total of 266 events were observed during the daytime CH4 releases, 221 of which were finally selected after filtering out for low turbulence and/or non‐stationarity. For the nighttime release, the number of 15‐min periods was equal to 410, and as expected, there was a higher number of filtered events due to lower turbulence, which resulted in 173 final observations (Table 1). From the nighttime events, only 20 observations were categorized as unstable, and only a very small percentage of stable events were included as most of them were filtered out based on low turbulence (u* < 0.2 m s−1). The number of highly daytime convective events (zL −15.5) was equal to zero. This was important as the Kljun model is restricted to operate only above that stability threshold.

Table 1

Number of 15‐Min Events and Quality Controls for the CH4 Releases During the Daytime and the Nighttime

	Daytime release	Nighttime release
Number of 15‐min flux periods evaluated	266	410
Filtered by wind direction	0	6
Filtered by low turbulence and stationarity	45	234
Total number of observations	221	173
Unstable	154	20
Neutral‐stable	67	150

Canopy height for the first release was equal to 0.75 m, as measured directly. For the second and third releases, when the alfalfa had been recently cut, the canopy heights were equal to 0.18 and 0.19 m, respectively, based on calculations of aerodynamic canopy height using the algorithm of Pennypacker & Baldocchi (2016). Roughness length was calculated using the algorithm of (Maurer et al., 2013) and was equal to 0.052 m for the first campaign, 0.017 m for the second, and 0.019 m for the third. For the second and third campaigns, which were done within a couple of days from each other, we took information from the sonic‐anemometer wind profile during neutral conditions to graphically calculate the value of roughness length (Figure S2 in Supporting Information S1). Based on the sonic anemometer wind profile under neutral conditions, the median value of roughness length from the 22 profiles under neutral conditions was equal to 0.022 m, which was very close to the value calculated using the approach of Maurer et al. (2013). This provided us with more confidence in the roughness length values used in the study. From the sonic anemometer profile, we observed that the wind shear between lower and top sonic anemometers was higher during stable conditions and under low wind speeds and low u*. In fact, no significant wind shear was encountered when u* was higher than 0.2 m s−1 (Figure S3 in Supporting Information S1), which gave us more confidence on the selection of 0.2 m s−1 as a threshold for turbulence.

Footprint Model Validation

The average CH4 flux of all the events considered for the validation of the models was equal to 329.4 ± 427.0 nmol m−2 s−1 (mean ± standard deviation) for the daytime release and 4.12 ± 11.28 nmol m−2 s−1 for the nighttime release. A histogram of the CH4 fluxes measured during the release and used for the validation is shown in Figure S4 in Supporting Information S1. Overall, the footprints were shorter using the Kljun model and longer using the K&M model, particularly during the nighttime (Figure 2). There were no strong differences in daytime footprint between the Hsieh and the K&M model, while during the night, the footprints from the Hsieh model tended to be longer. As expected, footprints under unstable conditions were much shorter than those under neutral or stable conditions (Figure S5 in Supporting Information S1). Footprints from the Kljun model tended to be shorter, especially during neutral conditions while the Hsieh model produced much longer footprints.

Figure 2

Cumulative footprints during (a) daytime and (b) nighttime releases using three footprint models: Hsieh: (Hsieh et al., 2000), Klujn: (Kljun et al., 2015), and K&M: (Kormann and Meixner, 2001). The outer contours represent the 80% and 70% footprint contours for daytime and nighttime conditions, respectively, whereas the inner contours represent 50% contours for both daytime and nighttime conditions.

The comparison of observed methane fluxes (F ) against calculated methane fluxes with the footprint model (F ) (Figure 3) suggested that the three models performed satisfactorily. However, the Kljun and Hsieh models tended to heavily overestimate the predictions of CH4 flux at the lower range, given they tend to predict peak weight closer to the tower than the K&M model. The overestimation was particularly high in the Kljun model, which produced much shorter footprints, and thus the daytime fluxes tended to be overestimated (i.e., higher expected contribution from the 20 m hot spot). During the daytime, the normalized mean bias estimate was equal to 45% for the Hsieh model, 113% for the Kljun model, and 8.9% for the K&M model. For the nighttime, the normalized mean bias estimate was 267% for the Hsieh model, −0.08% for the Kljun model, and 162% for the K&M model. It was the K&M model that showed the best performance for the daytime release, as it showed the lower RMSE and the lowest MBE (Figure 3).

Figure 3

Comparison of methane flux measured through eddy covariance against the methane flux calculated with the release flow and the footprint model. Three different models are evaluated (Hsieh, Kljun, and Korman & Meixner (K&M) for daytime (upper panel) and nighttime (lower panel).

For the nighttime release, the performance of the models was more favorable for the K&M model and the Kljun model. The Hsieh model heavily overestimated the fluxes, and this was due to the long footprints produced by this model. The K&M model had a better performance with a slope of 1.1 and a lower RMSE, but its MBE was high. The Kljun model, on the other hand, had a lower MBE and a lower RMSE, but it tended to underestimate fluxes at the higher end. This was due to the short footprint extension of the Kljun model in comparison to the other two models. The evaluation with the 30‐min data set offered very similar results (Figure S6 in Supporting Information S1). Given these evaluations, we decided to use the K&M model for daytime and nighttime conditions for the subsequent analysis in this manuscript.

Different variables driving the error in the models were evaluated for daytime (Figure S9 in Supporting Information S1) and nighttime (Figure S10 in Supporting Information S1) conditions. We found that the standard deviation of the vertical velocity component () was correlated with the error in the Hsieh model during daytime conditions but not during nighttime, stable conditions (Figures S9 and S10 in Supporting Information S1). We found that wind direction tended to drive some of the model error, particularly during the nighttime, but we did not find particularly strong relationship between model error and variables such as wind shear, friction velocity, and sensible heat flux (data not shown).

Testing Footprint‐Weighted Flux Maps With Artificial Hot Spots of Methane

The analysis of the methane release experiment using footprint‐weighted flux maps is shown in Figure 4. Figures 4a and 4d show the cumulative footprint for the ∼3 days of the continuous release, while Figures 4b and 4e, show the relative flux maps for daytime and nighttime conditions, respectively. These footprint‐weighted maps are used to locate areas with the highest flux rather than for flux quantification. The fluxes from this footprint‐weighted flux map can be interpreted as the average magnitude of the flux measured by the eddy covariance tower when the footprint intersects that cell. In the footprint‐weighted flux maps, it can be appreciated that the artificial hot spot created with the release experiment was located accurately in terms of its angle of orientation toward the tower, and to a minor extent, to its proximity, with the daytime release being closer than the nighttime point. Both footprint‐weighted flux maps, however, contain undesired gradients in flux that do not represent the true fluxes of the site (Figures 4b and 4e), what we call the “blurring effect”, which is an artifact of the technique given that (a) the footprint is a probability function and (b) no prior information about the location of the hot spot is given. However, because we know the location of the hot spot we can now try to back‐calculate the flux using the approach from Equation 13, which results in the plots shown in Figures 4c and 4f. For this calculation, we used the assumption of a near‐zero flux in the alfalfa field and the known size of the hot spots (equal to one grid cell of 3 × 3 or 5 × 5 m during daytime and nighttime, respectively). Finally, because we know the true flux based on the flow rate and the dimensions of the pipe used for the release, we can compare estimated fluxes against real fluxes as shown in Figures 4c and 4f. These results show that the magnitude of the hot spots calculated by implementing Equation 13 matches closely the hot spot's prescribed magnitude from the release after normalizing it by the area of one grid cell (Figures 4c and 4f).

Figure 4

Footprint‐weighted flux maps to detect the artificial hot spots of methane from the methane release experiment during the daytime (top row) and nighttime (bottom row) conditions. A cumulative footprint is obtained for the times during which the tracer release was active (a and d), then the footprint‐weighted flux map is calculated (b and e), and finally, an absolute flux from the hot spot can be calculated if the locations of the hotspot is known and the background methane flux can be estimated (c and f). Both daytime and nighttime footprints are calculated using the K&M model.

Evaluating Sharp Transitions in Land Cover With Footprint‐Weighted Flux Maps

Another way to test the applicability of footprint‐weighted flux maps is to test them on sites of known heterogeneity or flux discontinuity (e.g., Hsieh et al., 2000). The maps in Figures 5a and 5b show an eddy covariance tower located at the border between two contrasting land covers: a wetland in the center, a highly productive alfalfa crop toward the northwest and a fallow corn field toward the south‐southwest. In this scenario, which occurred during March 2013, one could expect much higher CH4 emissions from the wetlands than from the alfalfa or corn fields. Indeed, this flux pattern was nicely captured by the footprint‐weighted flux map both during the daytime (Figure 5b) and during the nighttime (Figure 5g). As a way to corroborate the validity of this method, CO2 and H2O flux measurements were also evaluated using the footprint‐weighted flux maps (Figures 5c–5g and  5h). Since the alfalfa is photosynthetically active during spring daytime conditions while the wetland plants are still brown and inactive during these months, we can observe in Figure 5c much higher photosynthesis from the northern alfalfa than from the western wetland. Similarly, the flux maps also show higher evapotranspiration from the alfalfa than from the inactive blanket of tules and cattails from the wetlands (Figure 5d). This pattern also occurs during the nighttime (Figure 5i), which is surprising because transpiration from alfalfa is restricted during the nighttime. This behavior could be a consequence of the wetlands receiving colder waters that arrive from the Sacramento River during the period of April, where water is still very cold from the snowmelt in the mountains. It is worth remembering that no assumption about the underlying surface variability is taken in the application of footprint‐weighted relative flux maps, and thus the technique alone can pick up these landscape features. Figures 5e and 5j also show that the mean wind came consistently from the west and that there is not a strong pattern of wind direction with time of day, which guarantees that the observed spatial patterns are not being biased by temporal patterns.

Figure 5

Footprint‐weighted flux maps from an eddy covariance tower located at the border between three contrasting land covers: a wetland in the center, a highly productive alfalfa crop toward the northwest and a fallow corn field toward the south. Plots (a and f) show the 50%–80% footprint climatology for DOY 090–120, 2013, for daytime and nighttime values, respectively. Plots (b and g); (c and h); and (d and i) show the corresponding footprint‐weighted flux map for methane flux, water vapor flux, and CO2 flux, respectively. The dashed red line indicates the section where the wetland ends, and there is a sharp transition to a fallow corn field. The colorbar units are footprint‐weighted flux. Plots (e and j) show the average diel pattern in methane emissions and the wind direction during the aggregation time. Error bars are the standard deviation. The spatial resolution for the footprint‐weighted flux maps is 4 and 6 m for daytime and nighttime conditions, respectively.

During the nighttime, the expected behavior in surface fluxes was also captured (Figures 5f–5h, and  5i) and in addition, the nighttime data provided an opportunity to evaluate the ability to detect transitions along the main axis of the footprint. In this case, it can be noted in Figure 5f that about 200 m upwind from the tower, the wetland ends, and there is a sharp transition to a fallow corn field. One could expect much lower methane emissions, lower transpiration, and higher soil respiration in this corn field in line with what the footprint‐weighted flux maps were able to detect (dashed line in Figures 5f–5h and  5i). This transition is more evident in the case of CH4 emissions, given that CH4 fluxes from the fallow corn field are close to zero, thus the difference in surface fluxes between the wetland and the fallow corn is contrasting. This finding is relevant because it is evidence that the proposed technique not only maps areas with different fluxes at different wind directions around the tower, but i.e., also capable of detecting differences determined by distance from the sensor.

Hot Spots of CH4 Flux in a Recently Restored Wetland

The next wetland we evaluated was Sherman wetland, where irregular patterns in CH4 flux were identified. As expected, CH4 fluxes were low during the winter and high during the summer, which is a usual response to seasonal temperature variation. In addition, the fluxes in the second part of the summer were characterized by regular patterns that closely track the air temperature. However, during the first part of the summer, the patterns were more irregular, with frequent spikes. To better understand these irregular patterns, we evaluated two periods from DOY 155 to DOY 189 in 2019, and in DOY 213–243, 2020, where important hot moments of methane flux occurred. Figure 6 shows one hot spot of CH4 flux in the western side of the wetland, both in the selected period in 2019 (Figure 6b) and in the selected period of 2020 (Figure 6e). Figure 6 only shows nighttime footprints to reduce the variability associated with the diurnal patterns of CH4 flux driven by plant photosynthesis. However, this hotspot location was also detected during the daytime for DOY 155–185 2019 (Figure S7 in Supporting Information S1). Analysis of hot spots in other wetlands in the Delta (US‐Myb) also showed the existence of hot spots in the same location for daytime and nigh time locations (Figure S8 in Supporting Information S1).

Figure 6

Footprint‐weighted flux maps in Sherman Wetland from DOY 155–185 2019 (top) and DOY 213–243 (bottom) 2020 for nighttime conditions. Plots (a and c) show the cumulative footprints for the respective times. Plots (b and e) show the corresponding footprint‐weighted flux map for methane flux, and plots (c and f) show the estimation of the size of the hot spot and their magnitude. Plots (d and h) show the average diel pattern in methane emissions and the wind direction during the aggregation time. Error bars are the standard deviation. The spatial resolution of the footprint‐weighted flux maps is 4 m.

The hot spot of methane flux was found consistently in the far west side of the Sherman wetland (red area in Figures 6a and 6e). The footprint‐weighted map indicated that the general location of this hotspot was about 300 m away from the tower, by the road berm that separates two wetland parcels. Through visual inspection at the site and analysis of satellite imagery, we were able to attribute this hot spot location to an area of lower water tables running along the levee for about 50 m as shown in Figure 6. With this knowledge, the hotspot area was estimated, its relative contribution to the total area of the footprint‐weighted map was calculated, and using Equation 13 the flux from the hot spot was estimated. For the period of DOY 155–185 2019, the magnitude of the hot spot was estimated as 15717 nmol m−2 s−1, whereas for the period of DOY 213–243 2020 the magnitude of the hot spot was estimated to be 7166 nmol m−2 s−1.

Comparison Against Chamber Measurements

Chamber fluxes from the hot spot were significantly higher than the fluxes from the open water (Figure 7). The mean CH4 flux from the open water was equal to 9.35 ± 0.99 nmol m−2 s−1 (mean +S.D.) in April (n = 13) and 15.19 ± 7.5 nmol m−2 s−1 in May (n = 9). In contrast, the mean flux from the hot spot was equal to 679 ± 688 nmol m−2 s−1 in April (n = 22) and 6589 ± 7889 nmol m−2 s−1 in May (n = 21). These mean values included both ebullitive and diffusive fluxes. The diffusive fluxes were on average equal to 509 ± 507 nmol m−2 s−1 in April and 5879 ± 7208 nmol m−2 s−1 in May.

Figure 7

Methane fluxes as measured with chambers for the area of highest contribution to the footprint (open water) and the hot spot. The chamber measurements were performed on 06 April 2021 (left), and 16 May 2021 (right). Horizontal line is median, and whiskers extend from first and third quartiles to minimum and maximum chamber methane (CH4) flux rates.

The mean chamber flux results from May (6589 ± 7889 nmol m−2 s−1) were of the same order of magnitude to the flux estimated from the footprint‐weighted flux maps of Figure 7 (15717 nmol m−2 s−1 for a time around June 2019 and 7166 nmol m−2 s−1 for a time around August 2020).

Discussion

Validation of Footprint Models

As the availability of remote sensing products at high spatial resolution increases, and the number of eddy covariance sites with heterogeneous land covers also increases, we need to rely more on flux footprint models for data analysis. In this study, we validated three of the most commonly‐used footprint models using an artificial tracer in contrast to the more‐commonly used natural tracer approach. Natural tracer approaches take advantage of well‐defined heterogeneity (e.g., contrasting land uses) that result in sharp changes on the surface fluxes with different flux footprints (Arriga et al., 2017; Cooper et al., 2003; Göckede et al., 2005). On the other hand, artificial tracer release experiments like the one used in this study, have the advantage of a predetermined flow rate at a specific distance from the tower, thus offering more precise calibrations. Traditionally, sulfur hexafluoride (SF6) has been used to validate both analytical and Lagrangian models, using cross‐wind line sources on homogeneous areas (Finn et al., 1996), over tall canopies under different regimes of stability (Leclerc, Karipot, et al., 2003), and for canopies of intermediate roughness (Leclerc, Meskhidze, et al., 2003). Other tracer gases used for validation include CO2 (Arriga et al., 2017; Kumari et al., 2020) and CH4 (Dumortier et al., 2019; Heidbach et al., 2017). The uniqueness of this study is the use of the artificial tracer release over a wide range of stable and unstable conditions, and its subsequent applications for testing landscape discontinuities and detecting hot spots of methane flux.

We found that the K&M model had the best performance for daytime and nighttime conditions using methane as a tracer. Footprint validation work from other authors corroborates the validation results from this study. For example, Kumari et al. (2020) found that the K&M performed better than the Hsieh and Kljun models under different source‐receptor configurations. Dumortier et al. (2019), who also used CH4 for their release experiment, also found that the K&M model had better performance than the Kljun model, and as in this study, the Kljun model tended to produce shorter footprints, with most of the sources being projected at a distance very close to the tower. This feature in the Kljun model resulted in a large over‐estimation of CH4 fluxes when the release was at 20 m, and in an underestimation when the release was located at 180 m. However, not all footprint validation studies agree on the K&M model having the best performance, and it is worth noting that footprint model performance may also be affected by specific local conditions and likely on the tracer being avaluated. Using a similar approach to the one in this study, Heidbach et al. (2017) evaluated the K&M, Hsieh and Kljun models focusing on evaluating upwind flux contributions from a sharp change in roughness length nearby the eddy covariance tower. Unlike the results from this study, they found that the Hsieh model underestimated the flux. This discrepancy could have arisen partly because Heidbach et al. (2017) focused on daytime and unstable conditions. They also found that the Kljun model tended to shorten peak distances, as we found in our study. However, in their case, this resulted in the Kljun model having better performance, as opposed to our study where the Kljun model heavily overestimated fluxes under unstable conditions.

Footprint model validations tend to perform better under near‐neutral conditions (e.g., Finn et al., 1996). Validations in stable conditions are complicated by the difficulty of obtaining accurate small fluxes in only slightly stationary conditions (Finn et al., 1996). Our study is one of the few that includes validations across a range of thermal stratification conditions, which is especially important when there is interest in evaluating surface fluxes at longer distances from the tower. Different variables driving the error in the models were evaluated for daytime (Figure S9 in Supporting Information S1) and nighttime (Figure S10 in Supporting Information S1) conditions but none of these variables was a key driver of the models error. Other variables might be driving this variability, such as advection. Indeed, Leclerc, Karipot, et al. (2003) point out the importance of non‐local advection outside of the footprint in driving some of this error. Thus, wind direction‐specific advection may explain the relationship between model error and wind direction, but further research in local advection patterns is needed to make a stronger conclusion about this. These footprints models are built under the assumption of horizontal homogeneity (Lee, 2018). However, it is has been found that in the atmospheric surface layer Taylor's frozen hypothesis fails (Cheng et al., 2017), thus the size and properties of eddies may change horizontally and a new assumption for the creation of footprint models may need to be implemented in the future.

On the Use of Footprint‐Weighted Flux Maps

Most flux aggregation methods have been developed for airborne flux measurements (Kohnert et al., 2017; Mauder et al., 2008; Metzger et al., 2013), and fewer for stationary towers (Xu et al., 2017, 2018). One of the first approaches to mapping the spatial variability of fluxes around the eddy covariance towers by using the superposition of flux footprints from different wind directions was done by Schuepp et al. (1992), who developed the concept of a two‐dimensional surface flux map (or later called a flux density map) constructed from flux footprints. Airborne measurements include a larger representation of a surface grid, and thus provide a good representation of the spatial heterogeneity of fluxes at a given site (Desjardins et al., 1989). In addition, the development of environmental response functions (ERF), when applied to aircraft measurements, can provide a surface map of sensible and latent heat flux and their spatially‐resolved drivers (Metzger et al., 2013). However, airborne measurements do not have the continuous sampling at the smaller scale that might be needed to detect hotspots of CH4 flux. It was Xu et al. (2017) who adated the ERF to stationary towers, thus allowing to create one of the first flux maps of CO2, latent heat flux and sensible flux at a regional scale (20 by 20 km), by relating turbulent flux observations to meteorological forcings and surface properties across the area encompassed by the footprint. This concept was further extended (Xu et al., 2018) to evaluate eddy covariance fluxes in a box, and considering the effects of the underlying spatial heterogeneity on storage and vertical advection fluxes to thus create more robust estimates of net vertical flux. Other tower surface flux maps to detect leakages of CO2 using footprint modeling are reported in the literature (Lewicki et al., 2009), they are not widely used, and the application in this study is novel. Tower surface flux maps have the disadvantage of their representation being biased toward the point of the maximum footprint contribution (i.e., most of the signal comes from this area). However, this situation can be overcome if multiple footprints encompass different wind directions and extensions. For example, if a hot spot is located relatively far from the tower, a proper scenario to detect that hot spot would be to have a combination of shorter footprints that do not encompass the hot spot, or have a very minor contribution, and several other footprints that extend longer and do encompass the hot spot and thus its contribution to the footprint is higher.

One concern with the flux footprint aggregation is that aggregating footprints at different times of day may introduce bias due to diurnal and/or seasonal patterns of fluxes. In the case of CH4 flux, there is higher confidence in using the footprint‐weighted flux map because the spatial variability in methane emissions is potentially much higher than the temporal variability (Peltola et al., 2015). This is the opposite of other trace gases such as CO2 or H2O, where the temporal variability is usually higher than the spatial variability. In addition, we have selected aggregations for daytime and nighttime conditions separately and minimized seasonal variation by maintaining a limit of aggregation of 2 months. In addition, the lack of a pattern of wind direction changes with time of day in our study sites (Figures 5e, 5j and 6d, 6h) also aided in the reduction of potential biases in the spatial patterns detected.

We have seen how the footprint‐weighted flux maps are extremely useful delineating sharp transitions along a radial gradient from the tower (Figure 5) as well as detecting gradients in surface fluxes along the mean wind direction axis, as is the case for changes in land‐uses (e.g., from wetland to agricultural land) (Figure 5f). The footprint maps also perform well at locating hot spots of CH4, particularly, its direction relative to the tower, but also approximating its location downstream from the tower. For example, we observed from the release experiment that a hot spot near the tower is correctly mapped close to it rather than away from it (Figure 4b). The same was the case for the hot spot far away from the tower (Figure 4e), although in this case, the region of location of the hot spot was much wider.

Recommendations for the Use of Footprint‐Weighted Flux Maps

The proposed footprint‐weighted flux maps can aid in understanding the spatial variability of fluxes of different scalars, but it is important to understand the limitations of this tool. We have seen how to detect hotspots of CH4 flux or other trace gases using footprint‐weighted flux maps. These maps work well for methane because methane emissions can be highly heterogeneous and the difference between hot spots and median fluxes can be of several orders of magnitude. Identifying hot spots of other trace gases that are more spatially homogeneous, such as water vapor or carbon dioxide may be more challenging. Nevertheless, the tool can still be applied to detect broad patterns in water vapor or carbon dioxide flux in naturally heterogeneous ecosystems such as savannas or mixed forests (Griebel et al., 2016, 2020). In addition, this tool can be highly relevant for the new developments in mapping greenhouse gases and atmospheric pollutants around urban environments, which are highly heterogeneous and require the use of tall towers to measure them (e.g., Davis et al., 2017).

It is important to remember that the footprint‐weighted flux maps only identify areas with relatively higher fluxes, and the magnitude indicates the average flux that the eddy‐covariance tower detects when the footprint intersects a given location. In addition, the location of a potential hot spot is only an approximation, therefore to determine its exact location, it is necessary to combine these flux maps with a land cover map or a satellite image of the site. An example of this is shown in Figures 6c and 6f, where it was clear that the hot spot of methane flux coincided with the a strip of land close to the levee with different properties than the vegetation around it. This is important because without this information we could not have constrained the size of the hot spot. It can be noted that there is a difference between the initially suggested and the final hot spot location. This indicates that the hot spot location could have been located even outside of the footprint map. However, because we were able to inspect the conditions of the field and measure the hot spots we were confident in assigning this size and location.

Footprint‐weighted flux maps can also be used in conjunction with high‐resolution remote sensing products to better interpret differences in spectral indices around eddy covariance towers. We have seen in Figure 5 how strong transitions in land cover, which can be obtained at high spatial and temporal resolution with satellites such as Planet Labs, can be matched to footprint‐weighted flux maps. Moreover, we have seen how patterns in evaporation from images at 60 m resolution from the ECOSTRESS mission (Fisher et al., 2020; Hook and Fisher, 2019) match the patterns in footprint‐weighted flux maps of evaporation (data not shown). Matching these remote sensing products with these tested models has the potential to facilitate the understanding of the spatial heterogeneity of water and carbon fluxes as driven by species diversity or vegetation composition.

The use of the 80% footprints has been preferred by different authors (e.g., Chu et al., 2021) mainly for two reasons. First, it is often difficult to reproduce 90% of the sources especially during neutral and stable conditions that generate long footprints. Second, the area between the 80% and the 90% contours tends to be excessively large, often larger than the area within the 80% footprint. At the same time, the flux contributions from this area are minimal. In Fact, Arriga et al. (2017) suggested that 75% of the sources are confined to an area 10–20 times the effective measurement height. Thus, one could argue that the noise from this 80%–90% area might be too large and thus could be ignored for most practical cases.

Understanding Hot Spots of CH4 in Wetlands

With a larger representation of methane fluxes in flux networks (Delwiche et al., 2021; Irvin et al., 2021; Knox et al., 2019), there is an increasing need to study the spatial heterogeneity of methane emissions from wetlands. In this study, we focus on evaluating hot spots at scales from tens to hundreds of square meters around the footprint of the eddy‐covariance tower. However, there are different scales of heterogeneity that are important for evaluating hotspots. For example, in peatlands, the surface area of biogeochemical hotspots could be accurately characterized with a pixel size of 25 cm to encompass all the hummocks and hollows (Becker et al., 2008). Such high resolution is likely not required in freshwater wetlands, where changes in microscale topography are not as significant. In ecosystem such as the Tundra, relevant scales occur at the scale of the footprint, where variations in land cover and soil conditions can drive some of the methane flux variability (Rößger et al., 2019). Methane from wetlands can be transported to the atmosphere via three mechanisms: diffusion from saturated surface soils and water, ebullition or through plant‐mediated transport via the aerenchyma tissue. Plant‐mediated transport has been found to be one of the main drivers of CH4 fluxes in freshwater wetlands, both via aerenchyma transport (Villa et al., 2020) and through the enhancement of acetoclastic methanogenesis through root exudates (Angle et al., 2017; Matthes et al., 2014; Knox et al., 2021; Mitra et al., 2020), lagging ecosystem photosynthesis by 1–4h (Knox et al., 2021). Thus, the spatial heterogeneity in methane emissions can be linked to the spatial heterogeneity in emergent vegetation cover. In open water areas and mudflats, the high spatial variability and the occurrence of ebullition events (Villa et al., 2021) can result in large intermittency and non‐stationarity that may bias the budgets of methane flux (Göckede et al., 2019). Ebullitive effects may be associated with hot spots of methane flux as they drive high variability in the measured methane fluxes (Zorzetto et al., 2021). Thus, there is an increased interest in partitioning methane fluxes between background and ebullitive fluxes using wavelet transforms (Iwata et al., 2018; Schaller et al., 2019), and surface‐renewal techniques (Zorzetto et al., 2021). We believe there is an opportunity to further understand hot spots and hot moments of methane flux by combining these novel techniques with the proposed footprint‐weighted flux maps. For example, footprint‐weighted flux maps can be created for only those half‐hours with a high proportion of ebullition flux and thus the location of a hot spot that is driven by ebullition events can be determined.

In this study, we found that the location of the hot spots was associated with the existence of mudflats or shallow pans where the water level is close to the surface and conditions are less favorable for emergent vegetation to survive, resulting in areas that remain at early stages of ecological succession after wetland restoration. It has been found that mudflats can act as hot spots of methane, contributing disproportionally to ecosystem scale CH4 flux despite occupying a small fraction of the ecosystem area. For example, Rey‐Sanchez et al. (2018) found that mudflats occupying only 1.5% of a wetland can contribute about 6.8% of the total annual methane flux. Similar situations are found in pastures. Moreover, Baldocchi et al. (2012), and Teh et al. (2011) estimated that background emissions from the drained areas in the pasture were relatively low (0.96 ± 0.19 nmol m−2 s−1), while ditches and shallow water pans created CH4 flux hotspots on the order of 449 ± 75 nmol m−2 s−1. Through chamber measurements, Rey‐Sanchez et al. (2019) showed variations of up 3 orders of magnitude in methane fluxes from different transects within the same vegetation cover type. However, this type of spatial heterogeneity in chamber measurements might not be related to land cover changes, but to subsurface processes and microbial composition, particularly the ratio of methanogens to methanotrophs. Further understanding of why the hot spot of methane flux detected at Sherman wetland occurs requires microbiological and biogeochemical analysis that are outside the scope of this study but will be addressed in a future study. We hope that this study can serve as a useful tool to help researchers understand the spatial heterogeneity in their fluxes and thus tackle specific question about microbial and biogeochemical processes around flux towers.

Conclusions

In the first part of this study, we have tested three of the most common flux footprint models using artificial tracer gas experiments. We have used the validation results to find the best model to create footprint‐weighted flux maps, a new method to map the spatial variability of scalar fluxes around eddy covariance towers. These maps, which are developed based on cumulative footprints for periods of several weeks to months, indicate where any scalar fluxes may be higher, although they are particularly suited for identifying hot spots of methane flux. We validated footprint‐weighted flux maps using artificial methane release experiments, a natural tracer approach taking advantage of well‐defined source/sink changes associated with different land‐uses, and gas flux chambers. Using footprint‐weighted flux maps we were able to detect spatial flux heterogeneity accurately across ecosystems, and most importantly, hotspots of methane flux in a freshwater wetland in California. The hot spot consisted of a shallow pan or mudflat area with ∼103 times larger methane fluxes than mean fluxes observed elsewhere within the wetland. We were able to corroborate the existence of the hot spot using chamber measurements and found that by knowing the area of the hot spot and a background flux, the flux of the hot spot could be approximated from the footprint‐weighted flux maps. This new method has the potential to be used with remote sensing techniques to constrain the spatial heterogeneity of fluxes around eddy‐covariance towers and thus create better models of carbon, water, and methane flux across large areas to upscale ecosystem processes. Additionally, this method can inform wetland management for climate change mitigation by helping identify areas that require contrasting management actions.

Supporting Information S1

Click here for additional data file.