**Research article**| 03 Jan 2023

# The paradox of assessing greenhouse gases from soils for nature-based solutions

Rodrigo Vargas and Van Huong Le

**Rodrigo Vargas and Van Huong Le**Rodrigo Vargas and Van Huong Le

- Department of Plant and Soil Sciences, University of Delaware, Newark, DE, USA

- Department of Plant and Soil Sciences, University of Delaware, Newark, DE, USA

**Correspondence**: Rodrigo Vargas (rvargas@udel.edu)

**Correspondence**: Rodrigo Vargas (rvargas@udel.edu)

Received: 12 Jul 2022 – Discussion started: 26 Jul 2022 – Revised: 02 Nov 2022 – Accepted: 14 Nov 2022 – Published: 03 Jan 2023

Quantifying the role of soils in nature-based solutions requires accurate estimates of soil greenhouse gas (GHG) fluxes. Technological advances
allow us to measure multiple GHGs simultaneously, and now it is possible to provide complete GHG budgets from soils (i.e., CO_{2}, CH_{4},
and N_{2}O fluxes). We propose that there is a conflict between the convenience of simultaneously measuring multiple soil GHG fluxes at fixed
time intervals (e.g., once or twice per month) and the intrinsic temporal variability in and patterns of different GHG fluxes. Information derived from
fixed time intervals – commonly done during manual field campaigns – had limitations to reproducing statistical properties, temporal dependence,
annual budgets, and associated uncertainty when compared with information derived from continuous measurements (i.e., automated hourly measurements)
for all soil GHG fluxes. We present a novel approach (i.e., temporal univariate Latin hypercube sampling) that can be applied to provide insights
and optimize monitoring efforts of GHG fluxes across time. We suggest that multiple GHG fluxes should not be simultaneously measured at a few fixed
time intervals (mainly when measurements are limited to once per month), but an optimized sampling approach can be used to reduce bias and
uncertainty. These results have implications for assessing GHG fluxes from soils and consequently reduce uncertainty in the role of soils in
nature-based solutions.

- Article
(1759 KB) -
Supplement
(1989 KB) - BibTeX
- EndNote

Soils are essential for nature-based solutions for their role in climate mitigation potential through implementing different natural pathways (Griscom et al., 2017; Bossio et al., 2020). The climate mitigation potential of soils is dependent on multiple factors such as weather variability (Kim et al., 2012), ecosystem type (Oertel et al., 2016), soil structure (Ball, 2013), management practices (Shakoor et al., 2021), or disturbances (Vargas, 2012), where soils can ultimately act as net sources or sinks of greenhouse gases (GHGs). Therefore, accurate quantification of the magnitudes and patterns of soil GHG fluxes is needed to understand the potential of soils to mitigate or contribute to global warming across ecosystems and different scenarios.

Most of our understanding of soil GHGs has come from manual measurements performed throughout labor-intensive field campaigns and experiments (Oertel
et al., 2016). While most studies around the world have focused on soil CO_{2} fluxes (Jian et al., 2020), early examples have reported coupled
measurements of soil CO_{2}, CH_{4}, and N_{2}O fluxes across tropical forests (Keller et al., 1986) and savannas (Hao et al., 1988),
temperate forests (Bowden et al., 1993), and peatlands (Freeman et al., 1993). These pioneer studies not only provided an early view of the importance of
integrated measurements of multiple soil GHG fluxes to understand the net global warming potential of soils but also demonstrated the technical
limitations and challenges associated with these efforts. For example, it is known that manual measurements have the strength of providing good
spatial coverage during field surveys but provide limited information about the temporal variability (Yao et al., 2009; Barba et al., 2021).

Technological advances have opened the opportunity to simultaneously measure multiple soil GHG fluxes (i.e., CO_{2}, CH_{4}, and
N_{2}O) at unprecedented temporal resolution (e.g., hourly). These efforts have demonstrated differences in diel patterns and pulse events (e.g.,
rewetting) due to wetting and drying cycles across tropical (Butterbach-Bahl et al., 2004; Werner et al., 2007), subtropical (Rowlings et al., 2012),
and temperate (Savage et al., 2014; Petrakis et al., 2017) ecosystems. These approaches provide more accurate information to calculate net GHG budgets
and the global warming potential of soils (Capooci et al., 2019). That said, performing automated measurements of multiple GHGs is expensive, and this
approach usually has a lower representation of the spatial heterogeneity within ecosystems (Yao et al., 2009; Barba et al., 2021).

Ideally, we would like to measure everything, everywhere, and all the time, but this is impossible due to logistical, technological, physical, and economic constraints. Lightweight and low-powered laser-based spectrometers have reduced technical barriers to simultaneously measuring multiple GHG fluxes from soils. It is now easier and faster to perform discrete manual surveys across time. This opportunity creates a paradox concerning when to measure different GHG fluxes from soils when performing manual measurements. Researchers generally tend to perform simultaneous measurements of multiple GHGs during manual surveys, but this convenience could result in biased information. We propose that there is a conflict between the convenience of measuring multiple GHGs at a few fixed time intervals and the intrinsic temporal variability in magnitudes and patterns of different GHG fluxes.

Here, we present a proof of concept and test how a subset of measurements derived from a fixed temporal stratification (FTS) for simultaneous
measurements (i.e., stratified sampling schedule) or using an optimized sampling (i.e., temporal univariate Latin hypercube
sampling, *tuLHs*) compared with automated measurements of soil CO_{2}
(*F*_{A}CO_{2}), CH_{4} (*F*_{A}CH_{4}), and N_{2}O (*F*_{A}N_{2}O) fluxes from a temperate forest
(Petrakis et al., 2018; Barba et al., 2021, 2019). The underlying assumption supporting any FTS approach is that a few measurements in time can
reproduce the statistical properties and temporal dependencies of soil CO_{2}, CH_{4}, and N_{2}O fluxes because these GHGs respond
similarly to biological and physical drivers. The tuLHs is a new statistical method for generating subsamples of parameter values (i.e., soil
GHG gas fluxes in this case study) to reproduce the probability distribution and the temporal dependence of each original time series of GHG
fluxes. We reveal that reporting GHG fluxes using an FTS for simultaneous measurements may result in biased information on temporal patterns and
magnitudes. This study shows how a biased sampling schedule could influence our understanding of GHG fluxes and, ultimately, the climate mitigation
potential of soils.

## 2.1 Study site

The experiment was performed in a temperate forest located at the St Jones Reserve (a component of the Delaware National Estuarine Research
Reserve (DNERR) in Delaware, USA). The site has a mean annual temperature of 13.3 ^{∘}C and mean annual precipitation of
1119 mm. Soils are classified as Othello silt loam with a texture of 40 % sand, 48 % silt, and 12 % clay within the first
10 cm (Petrakis et al., 2018). The dominant plant species are bitternut hickory (*Carya cordiformis*), eastern red cedar
(*Juniperus virginiana* L.), American holly (*Ilex opaca*), sweetgum (*Liquidambar styraciflua* L.), and black gum
(*Nyssa sylvatica* (Marshall)). The site has a mean tree density of 678 stems ha^{−1} and a diameter at breast height (DBH) of
25.7 ± 13.9 cm (mean ± SD) (Barba et al., 2021).

## 2.2 Automated measurements of soil GHG fluxes

We analyzed data from automated measurements (1 h time intervals) of soil emissions of three GHGs (i.e., CO_{2}, CH_{4}, and
N_{2}O) between January and December 2015. This was a typical year with a mean annual temperature of 13.4 ^{∘}C and annual
precipitation of 1232 mm. Continuous measurements of soil GHGs were taken by coupling a closed-path infrared gas analyzer (LI-COR LI-8100A,
Lincoln, Nebraska) and nine dynamic soil chambers (LI-COR 8100-104) controlled by a multiplexer (LI-COR 8150) with a cavity ring-down spectrometer (Picarro G2508, Santa Clara, California). A detailed description of the
experimental design and measurement protocol has been described in previous studies (Petrakis et al., 2018; Barba et al., 2021, 2019). Briefly, for each
flux observation, we measured CO_{2}, CH_{4}, and N_{2}O concentrations every second with the Picarro G2508 for 300 s and
calculated fluxes (at 1 h time intervals) from the mole dry fraction of each gas (i.e., corrected for water vapor dilution) using the
SoilFluxPro software (v4.0; LI-COR, Lincoln, Nebraska, USA). Fluxes were estimated using linear and exponential fits, and we kept the flux calculation
with the highest *R*^{2}. We applied quality assurance and quality control protocols using information from all three GHGs as established in previous
studies (Petrakis et al., 2018, 2017; Barba et al., 2021, 2019; Capooci et al., 2019). Using these time series, we extracted values to represent
discrete temporal measurements based on FTS and used the optimization approach described below.

## 2.3 Temporal subsampling of time series

Subsampling of time series was performed using FTS and a temporal optimization following a univariate Latin hypercube (tuLHs) approach. The difference between FTS and temporal optimization is that the first approach is focused on a fixed schedule (e.g., sampling once per month) and the second is focused on reproducing the statistical properties and temporal dependence relationship of the original GHG time series with a subset of measurements. This means optimized subsamples may not be spaced systematically (e.g., every 15 d), and selected dates may vary for each GHG flux due to their specific statistical properties and temporal dependence.

FTS represents a traditional schedule for performing manual measurements of GHG fluxes from soils. The FTS is usually performed with manual
measurements because they require extensive logistical coordination due to travel time and costs, availability of instrumentation (e.g., gas
analyzers), personnel to perform the measurements, and weather conditions. During these scheduled visits, researchers usually collect fluxes from all
three GHGs and analyze them systematically to calculate magnitudes and patterns throughout the experiment. Usually, researchers perform manual samples
during the early hours of the day (between 09:00 and 12:00 LT) to avoid confounding effects due to large changes in temperature and moisture, as
demonstrated by information summarized by the soil respiration global database (Cueva et al., 2017; Jian et al., 2020). Consequently, we selected
subsamples from each original GHG time series (derived from automated measurements) using flux measurements from 10:00 LT at fixed intervals of once per
month (*n* = 12), twice per month (*n* = 24), or four times per month (*n* = 48) starting from the first week of available data from
automated measurements.

We applied tuLHs as an alternative subsampling approach to obtain an optimized subsample with the same univariate statistical properties and
temporal dependence relationship as the original GHG time series. Optimization was performed to select subsamples for each GHG flux using the same
number of samples as for FTS: 12 (*k* = 12), 24 (*k* = 24), or 48 (*k* = 48) measurements throughout the year of
available data from automated measurements.

## 2.4 Temporal univariate Latin hypercube sampling (tuLHs)

Let $S=\mathit{\left\{}\right({x}_{\mathrm{1}},{y}_{\mathrm{1}},{z}_{\mathrm{1}}),({x}_{\mathrm{2}},{y}_{\mathrm{2}},{z}_{\mathrm{2}}),\mathrm{\dots},({x}_{n},{y}_{n},{z}_{n}\left)\mathit{\right\}}$ be observations of the variables *X*, *Y*, and *Z* in a time
series, where *X*, *Y*, and *Z* are soil GHGs (i.e., CO_{2}, CH_{4}, and N_{2}O). Each measured variable is characterized by the
univariate probability distribution function and the temporal dependency function. Once these two functions are known, then the behaviors of the
variable can be reproduced (Le et al., 2020; Chilès and Delfiner, 2009; Trangmar et al., 1986; Pyrcz and Deutsch, 2014). The tuLHs
consists of three steps: (1) modeling the univariate behavior of the variable using the empirical cumulative univariate probability distribution function;
(2) modeling the temporal dependence using the empirical variogram function; and (3) optimizing a subsample applying a global optimization method,
differential evolution, using the previously obtained variogram function as an objective function.

First, to model the univariate behavior of the variables from the observations of *S*, the empirical univariate cumulative distribution function
${F}_{n}^{\ast}\left(x\right)$ of *X* is estimated by

where *I* represents an indicator function equal to 1 when its argument is true and is 0 otherwise. Similarly, the empirical univariate distribution
function of the variables *Y* and *Z* can be derived.

Second, to model the temporal dependence of the variables from the observations of *S*, the empirical temporal correlation function (i.e., temporal
variogram function) *γ*^{∗}(*t*) of *X* is estimated by

where *N*(*t*) is the number of pairs and *X*(*t*_{i}+*t*) and *X*(*t*_{i}) are separated by a time *t*. The variogram functions of the variables *Y* and *Z* are
analogous. Third, to optimize the subsample, it is required to choose the “optimal” data points with the selected sample size (i.e., *k* = 12,
24, or 48, where *k*≪*n*) that will have the same behavior as the original observations of *S* (i.e., GHG fluxes derived from automated
measurements). To achieve this objective, we use differential evolution, a global optimization method (Storn and Price, 1997), using the variogram
function as an objective function. The procedure consists of dividing the univariate empirical probability distribution in Eq. (1) into
*k* equiprobable strata, which is equivalent to *k*-ordered data subsets. From each subset, only one value must be chosen to satisfy the condition of a
univariate Latin hypercube. The differential evolution method is applied to find the optimal points that minimize the difference between the subsample
variogram *γ*(*t*) and the data variogram *γ*^{∗}(*t*) in Eq. (3).

where OF is the objective function and the variograms *γ*(*t*) and *γ*^{∗}(*t*) are calculated using Eq. (2).

## 2.5 Statistical analyses

The *t* test was used to compare the means, and the Kolmogorov–Smirnov test was used to compare the probability distribution of measurements derived
from each sampling protocol. All tests were performed with a 95 % confidence level. In addition, their statistical properties, such as the mean, median,
standard deviation, and first and third quartiles, were compared. The differences in the experimental semivariograms were calculated as a comparison
measure for the temporal dependence of the samples and the original time series of GHG fluxes. For cumulative sums of GHG flux, their mean is
calculated as the most likely value, and their quantile difference between 97.5 % and 2.5 % is used to quantify the range of uncertainty. All analyzes
were performed using the R program (R Code Team, 2013).

## 3.1 Relationships among GHG fluxes from soils

Justification in support of FTS for simultaneous measurements of GHG fluxes would require evidence of strong linear correlations between magnitudes
and temporal dependence among soil GHG fluxes. First, we did not find strong linear relationships between any combination of GHG fluxes from soils
derived from automated measurements (Fig. S1 in the Supplement). Therefore, our data did not support
the assumption that the magnitude of one GHG flux was associated with a linear increase or decrease in another GHG flux. Second, semivariogram models
demonstrated differences in the temporal dependence for each GHG flux. Automated measurements of soil CO_{2} fluxes (*F*_{A}CO_{2})
showed a temporal dependence following a Gaussian variogram model, with a nugget of 4, a sill plus nugget of 28, and a correlation range
of 80 d (Fig. S2a). Automated measurements of soil CH_{4} fluxes (*F*_{A}CH_{4}) also showed a temporal
dependence but followed a spherical variogram model, with a nugget of 7 × 10^{−8}, a sill plus nugget of 1.5 × 10^{−7}, and a
correlation range of 110 d (Fig. S2b). In contrast, automated measurements of soil N_{2}O fluxes (*F*_{A}N_{2}O) did not show
a temporal dependence, where a pure nugget effect was present and with a correlation range of 0 d (Fig. S2c). Consequently, these GHG fluxes'
magnitudes and temporal patterns were different and did not support FTS for simultaneous measurements of GHG from soils.

## 3.2 Optimization of GHG sampling protocols

We applied a tuLHs approach to identify subsamples with the same statistical properties and temporal dependence for each of the original GHG
time series from automated measurements. Subsamples were identified for 12 (*k* = 12), 24 (*k* = 24), or 48
(*k* = 48) measurements throughout the year for each GHG time series. Our results show that the optimized measurement dates were different for
each GHG flux (Fig. 1), and we provide explicit examples for *k* = 24 (Fig. 1) and *k* = 12 and 48 (Figs. S3 and S4).

The optimized CO_{2} subsamples were well distributed throughout the year for all sampling scenarios (i.e., *k* from 12 to 48) because
*F*_{A}CO_{2} had a strong temporal dependence and a small nugget effect with respect to the sill (Fig. S2a). The optimized CH_{4}
subsamples were also relatively well distributed throughout the year, especially for scenarios of *k* = 24 and *k* = 48, as
*F*_{A}CH_{4} also had a temporal dependence but with a higher nugget effect with respect to the sill (Fig. S2b). Finally, the optimized
N_{2}O subsamples were more challenging to define, especially with a small sample size (i.e., *k* = 12; Fig. S3c) because
*F*_{A}N_{2}O did not have a temporal dependence (Fig. S2c).

## 3.3 Differences in statistical properties and temporal dependency of subsamples

Overall, there were no statistically significant differences between the mean values derived from automated measurements and those from FTS or the
tuLHs approach (Fig. 2 for *k* = 24; Fig. S5 for *k* = 12; Fig. S6 for *k* = 48; Tables S1
and S2 in the Supplement). Although this appears promising, more than a simple comparison of the
means is needed to evaluate the information derived from different sampling approaches. In other words, it is possible to have a similar mean value
without reproducing the probability distribution or the temporal dependence of the original time series (i.e., correct answer but for the wrong
reasons). Here, we present results based on comparing the means, standard deviation, probability distributions, and semivariograms derived from
automated measurements and the different sampling scenarios for all GHG fluxes.

The mean of *F*_{A}CO_{2} was 5.9 $\mathrm{\mu}\mathrm{mol}\phantom{\rule{0.125em}{0ex}}{\mathrm{CO}}_{\mathrm{2}}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{-\mathrm{2}}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{-\mathrm{1}}$, while the mean for FTS was 5.5 and
5.9 $\mathrm{\mu}\mathrm{mol}\phantom{\rule{0.125em}{0ex}}{\mathrm{CO}}_{\mathrm{2}}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{-\mathrm{2}}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{-\mathrm{1}}$ for the tuLHs approach with *k* = 24 (Fig. 3a–c). These results were comparable with the
means derived from FTS (5.4 and 5.4 $\mathrm{\mu}\mathrm{mol}\phantom{\rule{0.125em}{0ex}}{\mathrm{CO}}_{\mathrm{2}}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{-\mathrm{2}}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{-\mathrm{1}}$) and the tuLHs approach (6.2 and
5.9 $\mathrm{\mu}\mathrm{mol}\phantom{\rule{0.125em}{0ex}}{\mathrm{CO}}_{\mathrm{2}}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{-\mathrm{2}}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{-\mathrm{1}}$) using *k* = 12 and *k* = 48, respectively (Figs. S5 and S6; Table S1). The standard deviation of
*F*_{A}CO_{2} was 3.9 and 3.2 $\mathrm{\mu}\mathrm{mol}\phantom{\rule{0.125em}{0ex}}{\mathrm{CO}}_{\mathrm{2}}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{-\mathrm{2}}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{-\mathrm{1}}$ for FTS and 3.9 $\mathrm{\mu}\mathrm{mol}\phantom{\rule{0.125em}{0ex}}{\mathrm{CO}}_{\mathrm{2}}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{-\mathrm{2}}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{-\mathrm{1}}$ for the
tuLHs approach with *k* = 24 (Fig. 3a–c). These results were comparable with the standard deviations derived from FTS (3.1 and
3.3 $\mathrm{\mu}\mathrm{mol}\phantom{\rule{0.125em}{0ex}}{\mathrm{CO}}_{\mathrm{2}}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{-\mathrm{2}}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{-\mathrm{1}}$) and the tuLHs approach (4.1 and 3.9 $\mathrm{\mu}\mathrm{mol}\phantom{\rule{0.125em}{0ex}}{\mathrm{CO}}_{\mathrm{2}}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{-\mathrm{2}}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{-\mathrm{1}}$) using
*k* = 12 and *k* = 48, respectively (Figs. S5 and S6; Table S1). Our results show that the semivariograms of optimized samples using the
tuLHs approach closely approximate the semivariograms of automated measurements for *k* = 24 (Fig. 4a) and *k* = 12 and 48 (Figs. S7a
and S8a). These results are consistent with the sums of absolute differences between the semivariograms of the samples and the
semivariogram of *F*_{A}CO_{2} with differences of 69.31, 54.39, and 49.42 for FTS and of 5.69, 1.99, and 1.39 for the tuLHs approach for
*k* = 12, 24, and 48, respectively (Table S2).

The mean of *F*_{A}CH_{4} was −0.93 $\mathrm{nmol}\phantom{\rule{0.125em}{0ex}}{\mathrm{CH}}_{\mathrm{4}}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{-\mathrm{2}}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{-\mathrm{1}}$, while it was −0.86 $\mathrm{nmol}\phantom{\rule{0.125em}{0ex}}{\mathrm{CH}}_{\mathrm{4}}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{-\mathrm{2}}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{-\mathrm{1}}$ for FTS and
−0.94 $\mathrm{nmol}\phantom{\rule{0.125em}{0ex}}{\mathrm{CH}}_{\mathrm{4}}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{-\mathrm{2}}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{-\mathrm{1}}$ for the tuLHs approach with *k* = 24 (Fig. 3d–f). These results were also comparable with
the means derived from FTS (−0.83 and −0.88 $\mathrm{nmol}\phantom{\rule{0.125em}{0ex}}{\mathrm{CH}}_{\mathrm{4}}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{-\mathrm{2}}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{-\mathrm{1}}$) and the tuLHs approach (−0.87 and
−0.92 $\mathrm{nmol}\phantom{\rule{0.125em}{0ex}}{\mathrm{CH}}_{\mathrm{4}}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{-\mathrm{2}}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{-\mathrm{1}}$) using *k* = 12 and 48, respectively (Figs. S5 and S6; Table S1). The standard deviation of
*F*_{A}CH_{4} was 0.36 and 0.26 $\mathrm{nmol}\phantom{\rule{0.125em}{0ex}}{\mathrm{CH}}_{\mathrm{4}}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{-\mathrm{2}}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{-\mathrm{1}}$ for FTS and 0.34 $\mathrm{nmol}\phantom{\rule{0.125em}{0ex}}{\mathrm{CH}}_{\mathrm{4}}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{-\mathrm{2}}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{-\mathrm{1}}$ for the
tuLHs approach with *k* = 24. These results were comparable with the standard deviations derived from FTS (0.27 and
0.29 $\mathrm{nmol}\phantom{\rule{0.125em}{0ex}}{\mathrm{CH}}_{\mathrm{4}}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{-\mathrm{2}}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{-\mathrm{1}}$) and the tuLHs approach (0.33 and 0.35 $\mathrm{nmol}\phantom{\rule{0.125em}{0ex}}{\mathrm{CH}}_{\mathrm{4}}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{-\mathrm{2}}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{-\mathrm{1}}$) using *k* = 12 and
*k* = 48, respectively (Figs. S5 and S6; Table S1). The semivariograms of optimized samples using the tuLHs approach closely approximate
the semivariogram of automated measurements for *k* = 24 (Fig. 4b) and *k* = 12 and 48 (Figs. S7b and S8b). Consequently, the sums of absolute
differences between the semivariograms of the samples and the semivariogram of *F*_{A}CH_{4} were 0.63, 0.48, and 0.49 for FTS and 0.06, 0.04, and
0.02 for the tuLHs approach with *k* = 12, 24, and 48, respectively (Table S2).

Finally, the mean of *F*_{A}N_{2}O was 0.45 and 0.61 $\mathrm{nmol}\phantom{\rule{0.125em}{0ex}}{\mathrm{N}}_{\mathrm{2}}\mathrm{O}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{-\mathrm{2}}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{-\mathrm{1}}$ for FTS and
0.51 $\mathrm{nmol}\phantom{\rule{0.125em}{0ex}}{\mathrm{N}}_{\mathrm{2}}\mathrm{O}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{-\mathrm{2}}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{-\mathrm{1}}$ for the tuLHs approach with *k* = 24 (Fig. 3g–i). These results were also comparable with the
means derived from FTS (0.59 and 0.25 $\mathrm{nmol}\phantom{\rule{0.125em}{0ex}}{\mathrm{N}}_{\mathrm{2}}\mathrm{O}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{-\mathrm{2}}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{-\mathrm{1}}$) and the tuLHs approach (0.58 and
0.49 $\mathrm{nmol}\phantom{\rule{0.125em}{0ex}}{\mathrm{N}}_{\mathrm{2}}\mathrm{O}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{-\mathrm{2}}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{-\mathrm{1}}$) using *k* = 12 and 48, respectively (Figs. S5 and S6; Table S1). The standard deviation of
*F*_{A}N_{2}O was 1.62 and 1.97 $\mathrm{nmol}\phantom{\rule{0.125em}{0ex}}{\mathrm{N}}_{\mathrm{2}}\mathrm{O}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{-\mathrm{2}}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{-\mathrm{1}}$ for FTS and 1.54 $\mathrm{nmol}\phantom{\rule{0.125em}{0ex}}{\mathrm{N}}_{\mathrm{2}}\mathrm{O}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{-\mathrm{2}}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{-\mathrm{1}}$ for the
tuLHs approach with *k* = 24. These results were comparable with the standard deviations derived from FTS (1.38 and
0.91 $\mathrm{nmol}\phantom{\rule{0.125em}{0ex}}{\mathrm{N}}_{\mathrm{2}}\mathrm{O}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{-\mathrm{2}}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{-\mathrm{1}}$) and the tuLHs approach (1.58 and 1.54 $\mathrm{nmol}\phantom{\rule{0.125em}{0ex}}{\mathrm{N}}_{\mathrm{2}}\mathrm{O}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{-\mathrm{2}}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{-\mathrm{1}}$) using *k* = 12 and
*k* = 48, respectively (Figs. S5 and S6; Table S1). Our results show no temporal dependence for N_{2}O fluxes, but the semivariograms of
optimized samples using the tuLHs approach closely approximate the semivariogram of automated measurements for *k* = 24 (Fig. 4c) and
*k* = 12 and 48 (Figs. S7c and S8c). Consistently, the sum of absolute differences between the semivariograms of the samples and the semivariogram
of *F*_{A}N_{2}O was 10.01, 12.25, and 16.75 for FTS and 0.82, 1.13, and 3.57 for the tuLHs approach with *k* = 12, 24, and 48,
respectively (Table S2).

These results show that the tuLHs approach reproduced the probability distribution and the temporal dependence of the time series derived from automated measurements with more precision than FTS for all GHGs. In the next section, we explore the implications of these differences for calculating cumulative GHG fluxes.

## 3.4 Calculation of cumulative GHG fluxes

We calculated the cumulative flux for all GHGs using available information from automated measurements (Fig. 2; Table S3). The
cumulative sum for available measurements for *F*_{A}CO_{2} was 5758.5 g CO_{2} m^{−2} ([893.9, 13860.8], 95 % CI), for
*F*_{A}CH_{4} was −0.47 g CH_{4} m^{−2} ([−0.81, −0.19], 95 % CI), and for *F*_{A}N_{2}O was 0.63 g N_{2}O m^{−2} ([−0.75, 5.19],
95 % CI).

We used the mean for each GHG flux derived from the tuLHs approach or the FTS to calculate the cumulative sum (Table S3). We found that the
FTS underestimated the cumulative flux (−8.4 %, −6.2 %, −7.1 %) and the uncertainty (−32.6 %, −21.6 %, −19.3 %) of
*F*_{A}CO_{2} for *k* = 12, 24, and 48, respectively (Fig. 5a). In contrast, the tuLHs approach slightly overestimated the
cumulative flux (6.5 %, 1.1 %, 0.1 %) and slightly underestimated the uncertainty (−9.1 %, −4.4 %, −3.7 %) for
*k* = 12, 24, and 48, respectively (Fig. 5a).

The FTS underestimated the cumulative flux (−9.1 %, −6.1 %, −3.1 %) and the uncertainty (−31.8 %, −27.3 %,
−15.9 %) of *F*_{A}CH_{4} for *k* = 12, 24, and 48, respectively (Fig. 5b). In contrast, the tuLHs approach underestimated
the cumulative flux (−6.1 %) only for *k* = 12 but slightly underestimated the uncertainty (−15.9 %, −6.8 %, −4.5 %) for
*k* = 12, 24, and 48, respectively (Fig. 5b).

The FTS substantially underestimated the cumulative flux (−168 %, −170 %, −173 %) of *F*_{A}N_{2}O for *k* = 12, 24, and
48, respectively. Uncertainty was overestimated for *k* = 12 and 24 (3.6 % and 26 %) and underestimated for *k* = 48 (−31 %;
Fig. 5c). In contrast, the tuLHs approach overestimated the cumulative flux less (29.5 %, 13.4 %, 9.1 %) for *k* = 12, 24, and
48, respectively (Fig. 5c). This approach underestimated the uncertainty for *k* = 12 (−11.2 %) and *k* = 24 (−13.8 %) but
overestimated the uncertainty by 2.9 % for *k* = 48 (Fig. 5c). These results show that the tuLHs approach consistently provided
closer estimates for cumulative sums and uncertainty ranges than an FTS for all GHG fluxes.

Applied challenges, such as quantifying the role of soils in nature-based solutions, require accurate estimates of GHG fluxes. To do this, two fundamental problems exist for designing environmental monitoring protocols: where and when to measure. Ultimately a monitoring protocol aims to quantify the attributes of an ecosystem so that it can be compared in time within that ecosystem or with other ecosystems. Because we cannot measure everything, everywhere, and all the time, we can argue that any monitoring protocol has assumptions based on physical, economic, social, and practical reasons to address a specific scientific question. These assumptions for designing monitoring protocols could result in misleading, biased, or wrong conclusions, and therefore it is critical to assess the consequences of different monitoring efforts. As Hutchinson (1953) described in “The Concept of Pattern in Ecology”, we do not always know if a given pattern is extraordinary or a simple expression of something which we may learn to expect all the time.

Automated measurements have revolutionized our understanding of the temporal patterns and magnitudes of soil GHG fluxes (Savage et al., 2014;
Bond-Lamberty et al., 2020; Tang et al., 2006; Capooci and Vargas, 2022b). These measurements have limitations in representing spatial variability and
have higher equipment costs that limit their broad applicability across study sites (Vargas et al., 2011). Consequently, discrete manual measurements
are a common approach to simultaneously measure multiple GHG fluxes and report patterns, budgets, and information to parameterize empirical and
process-based models (Phillips et al., 2017; Wang and Chen, 2012). In this study, we argue that the convenience of simultaneously measuring multiple
GHGs using FTS may result in biased estimates. Therefore, optimization of sampling protocols is needed to provide insights to improve measurement
protocols when there is a limited number of measurements in time (i.e., *k* = 12, 24, 48).

We show that the magnitude of one GHG flux is not associated with a linear increase or decrease in another GHG flux, and the temporal dependencies of each GHG flux are different (Fig. S1). Therefore, it is not possible to infer the dynamics of one GHG flux based solely on information from another under the assumption that they share similar (or autocorrelated) biophysical drivers. These results imply that the magnitudes and temporal patterns of GHGs are different and therefore do not support an FTS approach for simultaneous measurements of GHG fluxes from soils.

Multiple studies have shown that the relevance of different biophysical drivers (e.g., temperature, moisture, light) is different for soil
CO_{2}, CH_{4}, or N_{2}O fluxes (Luo et al., 2013; Tang et al., 2006; Ojanen et al., 2010). Our results show that soil
CO_{2} fluxes have a strong temporal dependence (Fig. S2a), likely due to the strong relationship between these fluxes and soil temperature in
this and other temperate mesic ecosystems (Hill et al., 2021; Bahn et al., 2010; Barba et al., 2019). The temporal dependence decreased for soil
CH_{4} fluxes (Fig. S2b), where there is less evidence for such a strong correlation with soil temperature in this and other temperate mesic
ecosystems (Bowden et al., 1998; Castro et al., 1995; Warner et al., 2019; Barba et al., 2019). It has been reported that multiple variables and
complex relationships are usually needed to explain the variability in soil CH_{4} fluxes in forest soils, as there is a delicate balance between
methanogenesis and methanotrophy (Luo et al., 2013; Castro et al., 1994; Murguia-Flores et al., 2018). In contrast, soil N_{2}O fluxes had no
temporal dependence (Fig. S2c), showing decoupling from the observed patterns of soil CO_{2} and CH_{4} fluxes (Wu et al., 2010), likely as
a result of independent biophysical drivers regulating soil N_{2}O fluxes (Luo et al., 2013; Bowden et al., 1993; Ullah and Moore, 2011).

To address the limitations of an FTS protocol, we propose a novel optimization approach (i.e., tuLHs) to reproduce the probability distribution and the temporal dependence of each original time series of GHG fluxes. Traditional methods usually optimize subsamples by either individually focusing on reproducing the probability distribution of the original information (Huntington and Lyrintzis, 1998) or reproducing the temporal dependence of the original information (Gunawardana et al., 2011). The tuLHs is a simple approach that uses the univariate probability distribution function and the temporal correlation function (i.e., variogram) as objective functions for each GHG flux. Our results show that optimized subsamples do not coincide in time for the three GHGs, suggesting that information should be collected based on each GHG flux's specific statistical and temporal characteristics (Fig. 1). This study provides a proof of concept for the application of the tuLHs. It demonstrates how an optimization approach provides insights to design monitoring protocols and improve soil GHG flux estimates.

The more temporal data we can collect, the better, but in many cases, measurement protocols are limited to a few measurements per year (i.e.,
*k* = 12 to 48). Our results demonstrate that for a small sample size (i.e., *k* = 12), the optimized measurements for soil CO_{2} fluxes
are consistently spread across the year, and for soil CH_{4} fluxes are centered within the growing season because of their strong temporal
dependence. For the case of soil N_{2}O fluxes, the variogram shows a constant temporal variability, meaning there is no temporal
dependence. Therefore, the optimized measurements are concentrated within the fall season due to their distribution probability (Fig. 1a). Our
optimization approach shows how measurements can be distributed across time as more samples are available (i.e., *k* = 24 to 48; Fig. 1b and c)
and demonstrates that optimization is critical when a limited number of measurements are available. In other words, a few measurements properly
distributed across time provide better agreement with information derived from automated measurements. A similar conclusion was proposed for the
spatial distribution of environmental observatory networks, where a network of a few sites properly distributed (e.g., across a country) improves our
understanding of the target variable more than a spatially biased network (Villarreal et al., 2019). Thus, a representativeness assessment of
information collected across time and space is needed to evaluate environmental measurements and quantify nature-based solutions accurately.

We highlight that this optimization approach should be implemented across different ecosystems as it will result in site-specific recommendations. The
tuLHs can be applied to any time series length and with any time step (e.g., hours, days), but specific results will be representative of the
probability distribution and the temporal dependence of the selected time series. What is essential is to question if a few measurements from an
experiment represent the *reality* of the physical world because if limited information is available, then the actual probability distribution
and temporal dependence of the phenomena could be an unknown unknown. In other words, with few measurements, we may not be aware of and we will not be
able to know which is the actual probability distribution and temporal dependence of the studied phenomena. To address this challenge, we tested the
tuLHs approach with high-temporal-frequency information representing the probability distribution of multiple soil GHG fluxes at the daily time step
across a calendar year.

In this case study, the year chosen had typical climatological conditions, and we demonstrated that the statistical properties of the different GHG fluxes differ. Consequently, this study questions the application of the FTS approach to measuring multiple GHGs simultaneously with a limited number of sampling dates (mainly once a month). We recognize that longer time series (e.g., multi-year) could provide more robust optimizations that can be applied to monitoring efforts. We recommend co-locating automated measurements with manual survey efforts to adequately capture the temporal and spatial variability in soil GHG fluxes at study sites.

There are several implications of biased monitoring protocols for understanding soil GHG fluxes and nature-based solutions. First, temporal patterns
and temporal dependency may need to be revisited for studies using an FTS approach. Soil GHG fluxes have complex temporal dynamics that vary from
diurnal to seasonal and annual scales (Vargas et al., 2010) that a few measurements following an FTS approach cannot reproduce (Barba et al., 2019;
Bréchet et al., 2021). Second, soil GHG fluxes could present hot moments, which are transient events with disproportionately high values that are
often missed with an FTS approach (Vargas et al., 2018; Butterbach-Bahl et al., 2004). Third, cumulative sums and uncertainty ranges are biased or
misleading when derived using an FTS approach (Tallec et al., 2019; Lucas-Moffat et al., 2018; Capooci and Vargas, 2022b). Our study demonstrates that
an optimized approach consistently provided closer estimates for cumulative sums and uncertainty ranges when compared with automated measurements
(Fig. 5). We postulate that representing the variability in soil N_{2}O fluxes is more sensitive to the FTS approach (> 170 %
and > 30 % for cumulative sums and uncertainty ranges, respectively) than for soil CH_{4} and CO_{2} fluxes. Fourth, it is possible
that if the information derived from an FTS approach is biased, then functional relationships could also be different from those derived from
automated measurements (Capooci and Vargas, 2022a). It has been argued that hypothesis testing and our capability of forecasting responses of soil GHG
fluxes to changing climate conditions is also biased with information from the FTS approach (Vicca et al., 2014). Finally, because soils have a
central role in nature-based solutions within countries and across the world (Griscom et al., 2017; Bossio et al., 2020), accurate measurements are
required to assess management practices, environmental variability, and the contribution of GHGs from soils.

We highlight that we only sometimes know if a given pattern is extraordinary or a simple expression of something which we may learn to expect all the time. Arguably, there is bias in our understanding of the probability distribution and temporal dependency of soil GHG fluxes across the world because most results are based on a few manual measurements (e.g., once a month) following an FTS approach. Currently, it is unknown how large such bias could be across studies and ecosystems, but because most studies lack high-temporal-frequency information, the real probability distribution and temporal dependency of soil GHG fluxes may remain unknown in most study sites. What is essential is to question if the observed patterns, derived from an FTS approach, are enough for improving our understanding of soil processes or are results that we have learned to expect.

We postulate that with emergent technologies, there is convenience in measuring multiple GHGs from soils; however, few measurements collected at fixed time intervals result in biased estimates. We recognize that potential measurement bias depends on each GHG flux's magnitudes and temporal patterns and could be site-specific. Nevertheless, evaluations are needed to quantify potential bias in estimates of GHG budgets and information used for model parameterization and environmental assessments. Furthermore, the underlying assumption that each GHG flux responds similarly to biophysical drivers may need to be tested across multiple ecosystems to quantify how few measurements influence our understanding of magnitudes and temporal patterns of soil GHG fluxes.

In this study, we present a proof of concept and propose a novel approach (i.e., temporal univariate Latin hypercube sampling) that can be applied with site-specific information on different ecosystems to improve monitoring efforts and reduce the bias in GHG flux measurements across time. We highlight that constant biased environmental monitoring may provide confirmatory information, which we have learned to expect, but modifications of monitoring protocols could shed light on new or unexpected patterns. These new patterns are the ones that will test paradigms and push science frontiers.

All data used for this analysis are available at https://doi.org/10.6084/m9.figshare.19536004.v1 (Petrakis et al., 2022). The R code used in this study is available at https://doi.org/10.5281/zenodo.7323913 (Le and Vargas, 2022).

The supplement related to this article is available online at: https://doi.org/10.5194/bg-20-15-2023-supplement.

RV conceived this study, and VHL designed and performed the primary analysis with input from RV in all phases. RV wrote the manuscript with contributions from VHL.

The contact author has declared that neither of the authors has any competing interests.

Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

The authors thank the Delaware National Estuarine Research Reserve (DNERR) and the St Jones Reserve personnel for their support throughout this study. The authors acknowledge the land on which they realized this study as the traditional home of the Lenni-Lenape tribal nation (Delaware Nation).

This study was funded by a grant from the National Science Foundation (grant no. 1652594) and NASA Carbon Monitoring System (grant no. 80NSSC21K0964).

This paper was edited by Sara Vicca and reviewed by two anonymous referees.

Bahn, M., Reichstein, M., Davidson, E. A., Grünzweig, J., Jung, M., Carbone, M. S., Epron, D., Misson, L., Nouvellon, Y., Roupsard, O., Savage, K., Trumbore, S. E., Gimeno, C., Curiel Yuste, J., Tang, J., Vargas, R., and Janssens, I. A.: Soil respiration at mean annual temperature predicts annual total across vegetation types and biomes, Biogeosciences, 7, 2147–2157, https://doi.org/10.5194/bg-7-2147-2010, 2010.

Ball, B. C.: Soil structure and greenhouse gas emissions: a synthesis of 20 years of experimentation, Eur. J. Soil Sci., 64, 357–373, https://doi.org/10.1111/ejss.12013, 2013.

Barba, J., Poyatos, R., and Vargas, R.: Automated measurements of greenhouse gases fluxes from tree stems and soils: magnitudes, patterns and drivers, Sci. Rep.-UK, 9, 4005, https://doi.org/10.1038/s41598-019-39663-8, 2019.

Barba, J., Poyatos, R., Capooci, M., and Vargas, R.:
Spatiotemporal variability and origin of CO_{2} and CH_{4} tree stem fluxes in an upland forest, Glob. Change Biol., 27, 4879–4893, https://doi.org/10.1111/gcb.15783, 2021.

Bond-Lamberty, B., Christianson, D. S., Malhotra, A., Pennington, S. C., Sihi, D., AghaKouchak, A., Anjileli, H., Altaf Arain, M., Armesto, J. J., Ashraf, S., Ataka, M., Baldocchi, D., Andrew Black, T., Buchmann, N., Carbone, M. S., Chang, S., Crill, P., Curtis, P. S., Davidson, E. A., Desai, A. R., Drake, J. E., El-Madany, T. S., Gavazzi, M., Görres, C., Gough, C. M., Goulden, M., Gregg, J., Gutiérrez del Arroyo, O., He, J., Hirano, T., Hopple, A., Hughes, H., Järveoja, J., Jassal, R., Jian, J., Kan, H., Kaye, J., Kominami, Y., Liang, N., Lipson, D., Macdonald, C. A., Maseyk, K., Mathes, K., Mauritz, M., Mayes, M. A., McNulty, S., Miao, G., Migliavacca, M., Miller, S., Miniat, C. F., Nietz, J. G., Nilsson, M. B., Noormets, A., Norouzi, H., O'Connell, C. S., Osborne, B., Oyonarte, C., Pang, Z., Peichl, M., Pendall, E., Perez-Quezada, J. F., Phillips, C. L., Phillips, R. P., Raich, J. W., Renchon, A. A., Ruehr, N. K., Sánchez-Cañete, E. P., Saunders, M., Savage, K. E., Schrumpf, M., Scott, R. L., Seibt, U., Silver, W. L., Sun, W., Szutu, D., Takagi, K., Takagi, M., Teramoto, M., Tjoelker, M. G., Trumbore, S., Ueyama, M., Vargas, R., Varner, R. K., Verfaillie, J., Vogel, C., Wang, J., Winston, G., Wood, T. E., Wu, J., Wutzler, T., Zeng, J., Zha, T., Zhang, Q., and Zou, J.: COSORE: A community database for continuous soil respiration and other soil-atmosphere greenhouse gas flux data, Glob. Change Biol., 249, 434, https://doi.org/10.1111/gcb.15353, 2020.

Bossio, D. A., Cook-Patton, S. C., Ellis, P. W., Fargione, J., Sanderman, J., Smith, P., Wood, S., Zomer, R. J., von Unger, M., Emmer, I. M., and Griscom, B. W.: The role of soil carbon in natural climate solutions, Nature Sustainability, 3, 391–398, https://doi.org/10.1038/s41893-020-0491-z, 2020.

Bowden, R. D., Castro, M. S., Melillo, J. M., Steudler, P. A., and Aber, J. D.: Fluxes of greenhouse gases between soils and the atmosphere in a temperate forest following a simulated hurricane blowdown, Biogeochemistry, 21, 61–71, 1993.

Bowden, R. D., Newkirk, K. M., and Rullo, G. M.: Carbon dioxide and methane fluxes by a forest soil under laboratory-controlled moisture and temperature conditions, Soil Biol. Biochem., 30, 1591–1597, https://doi.org/10.1016/S0038-0717(97)00228-9, 1998.

Bréchet, L. M., Daniel, W., Stahl, C., Burban, B., Goret, J.-Y., Salomn, R. L., and Janssens, I. A.:
Simultaneous tree stem and soil greenhouse gas (CO_{2}, CH_{4}, N_{2}O) flux measurements: a novel design for continuous monitoring towards improving flux estimates and temporal resolution, New Phytol., 230, 2487–2500, https://doi.org/10.1111/nph.17352, 2021.

Butterbach-Bahl, K., Kock, M., Willibald, G., Hewett, B., Buhagiar, S., Papen, H., and Kiese, R.:
Temporal variations of fluxes of NO, NO_{2}, N_{2}O, CO_{2}, and CH_{4} in a tropical rain forest ecosystem, Global Biogeochem. Cy., 18, GB3012, https://doi.org/10.1029/2004gb002243, 2004.

Capooci, M. and Vargas, R.:
Diel and seasonal patterns of soil CO_{2} efflux in a temperate tidal marsh, Sci. Total Environ., 802, 149715, https://doi.org/10.1016/j.scitotenv.2021.149715, 2022a.

Capooci, M. and Vargas, R.: Trace gas fluxes from tidal salt marsh soils: implications for carbon–sulfur biogeochemistry, Biogeosciences, 19, 4655–4670, https://doi.org/10.5194/bg-19-4655-2022, 2022b.

Capooci, M., Barba, J., Seyfferth, A. L., and Vargas, R.: Experimental influence of storm-surge salinity on soil greenhouse gas emissions from a tidal salt marsh, Sci. Total Environ., 686, 1164–1172, https://doi.org/10.1016/j.scitotenv.2019.06.032, 2019.

Castro, M. S., Melillo, J. M., Steudler, P. A., and Chapman, J. W.: Soil moisture as a predictor of methane uptake by temperate forest soils, Can. J. Forest Res., 24, 1805–1810, https://doi.org/10.1139/x94-233, 1994.

Castro, M. S., Steudler, P. A., Melillo, J. M., Aber, J. D., and Bowden, R. D.: Factors controlling atmospheric methane consumption by temperate forest soils, Global Biogeochem. Cy., 9, 1–10, https://doi.org/10.1029/94gb02651, 1995.

Chilès, J.-P. and Delfiner, P.: Geostatistics: Modeling Spatial Uncertainty, John Wiley & Sons, 720 pp., 2009.

Cueva, A., Bullock, S. H., López-Reyes, E., and Vargas, R.:
Potential bias of daily soil CO_{2} efflux estimates due to sampling time, Sci. Rep.-UK, 7, 11925, https://doi.org/10.1038/s41598-017-11849-y, 2017.

Freeman, C., Lock, M. A., and Reynolds, B.:
Fluxes of CO_{2}, CH_{4} and N_{2}O from a Welsh peatland following simulation of water table draw-down: Potential feedback to climatic change, Biogeochemistry, 19, 51–60, https://doi.org/10.1007/bf00000574, 1993.

Griscom, B. W., Adams, J., Ellis, P. W., Houghton, R. A., Lomax, G., Miteva, D. A., Schlesinger, W. H., Shoch, D., Siikamäki, J. V., Smith, P., Woodbury, P., Zganjar, C., Blackman, A., Campari, J., Conant, R. T., Delgado, C., Elias, P., Gopalakrishna, T., Hamsik, M. R., Herrero, M., Kiesecker, J., Landis, E., Laestadius, L., Leavitt, S. M., Minnemeyer, S., Polasky, S., Potapov, P., Putz, F. E., Sanderman, J., Silvius, M., Wollenberg, E., and Fargione, J.: Natural climate solutions, P. Natl. Acad. Sci. USA, 114, 11645–11650, https://doi.org/10.1073/pnas.1710465114, 2017.

Gunawardana, A., Meek, C., and Xu, P.: A model for temporal dependencies in event streams, Advances in Neural Information Processing Systems 24 (NIPS 2011), edited by: Shawe-Taylor, J., Zemel, R., Bartlett, P., Pereira, F., and Weinberger, K. Q., Purchase Printed Proceeding, 1962–1970, ISBN: 9781618395993, 2011.

Hao, W. M., Scharffe, D., Crutzen, P. J., and Sanhueza, E.:
Production of N_{2}O, CH_{4}, and CO_{2} from soils in the tropical savanna during the dry season, J. Atmos. Chem., 7, 93–105, https://doi.org/10.1007/bf00048256, 1988.

Hill, A. C., Barba, J., Hom, J., and Vargas, R.:
Patterns and drivers of multi-annual CO_{2} emissions within a temperate suburban neighborhood, Biogeochemistry, 152, 35–50, https://doi.org/10.1007/s10533-020-00731-1, 2021.

Huntington, D. E. and Lyrintzis, C. S.: Improvements to and limitations of Latin hypercube sampling, Probab. Eng. Mech., 13, 245–253, https://doi.org/10.1016/S0266-8920(97)00013-1, 1998.

Hutchinson, G. E.: The Concept of Pattern in Ecology, P. Acad. Nat. Sci. Phila., 105, 1–12, 1953.

Jian, J., Vargas, R., Anderson-Teixeira, K., Stell, E., Herrmann, V., Horn, M., Kholod, N., Manzon, J., Marchesi, R., Paredes, D., and Bond-Lamberty, B.: A restructured and updated global soil respiration database (SRDB-V5), Earth Syst. Sci. Data, 13, 255–267, https://doi.org/10.5194/essd-13-255-2021, 2021.

Keller, M., Kaplan, W. A., and Wofsy, S. C.:
Emissions of N_{2}O, CH_{4} and CO_{2} from tropical forest soils, J. Geophys. Res., 91, 11791, https://doi.org/10.1029/jd091id11p11791, 1986.

Kim, D.-G., Vargas, R., Bond-Lamberty, B., and Turetsky, M. R.: Effects of soil rewetting and thawing on soil gas fluxes: a review of current literature and suggestions for future research, Biogeosciences, 9, 2459–2483, https://doi.org/10.5194/bg-9-2459-2012, 2012.

Le, V. H. and Vargas, R.: A workflow to apply temporal univariate Latin Hypercube (tuLH), Zenodo [code], https://doi.org/10.5281/zenodo.7323913, 2022.

Le, V. H., Díaz-Viera, M. A., Vázquez-Ramírez, D., del Valle-García, R., Erdely, A., and Grana, D.: Bernstein copula-based spatial cosimulation for petrophysical property prediction conditioned to elastic attributes, J. Petol. Sci. Eng., 193, 107382, https://doi.org/10.1016/j.petrol.2020.107382, 2020.

Lucas-Moffat, A. M., Huth, V., Augustin, J., Brümmer, C., Herbst, M., and Kutsch, W. L.:
Towards pairing plot and field scale measurements in managed ecosystems: Using eddy covariance to cross-validate CO_{2} fluxes modeled from manual chamber campaigns, Agr. Forest Meteorol., 256–257, 362–378, https://doi.org/10.1016/j.agrformet.2018.01.023, 2018.

Luo, G. J., Kiese, R., Wolf, B., and Butterbach-Bahl, K.: Effects of soil temperature and moisture on methane uptake and nitrous oxide emissions across three different ecosystem types, Biogeosciences, 10, 3205–3219, https://doi.org/10.5194/bg-10-3205-2013, 2013.

Murguia-Flores, F., Arndt, S., Ganesan, A. L., Murray-Tortarolo, G., and Hornibrook, E. R. C.: Soil Methanotrophy Model (MeMo v1.0): a process-based model to quantify global uptake of atmospheric methane by soil, Geosci. Model Dev., 11, 2009–2032, https://doi.org/10.5194/gmd-11-2009-2018, 2018.

Oertel, C., Matschullat, J., Zurba, K., Zimmermann, F., and Erasmi, S.: Greenhouse gas emissions from soils—A review, Geochem.-Explor. Env. A., 76, 327–352, https://doi.org/10.1016/j.chemer.2016.04.002, 2016.

Ojanen, P., Minkkinen, K., Alm, J., and Penttilä, T.:
Soil–atmosphere CO_{2}, CH_{4} and N_{2}O fluxes in boreal forestry-drained peatlands, Forest Ecol. Manag., 260, 411–421, https://doi.org/10.1016/j.foreco.2010.04.036, 2010.

Petrakis, S., Seyfferth, A., Kan, J., Inamdar, S., and Vargas, R.: Influence of experimental extreme water pulses on greenhouse gas emissions from soils, Biogeochemistry, 133, 147–164, https://doi.org/10.1007/s10533-017-0320-2, 2017.

Petrakis, S., Barba, J., Bond-Lamberty, B., and Vargas, R.: Using greenhouse gas fluxes to define soil functional types, Plant Soil, 423, 285–294, https://doi.org/10.1007/s11104-017-3506-4, 2018.

Petrakis, S., Barba, J., and Vargas, R.: Soil respiration dataset from a temperate forest (St Jones Reserve), figshare [data set], https://doi.org/10.6084/m9.figshare.19536004.v1, 2022.

Phillips, C. L., Bond-Lamberty, B., Desai, A. R., Lavoie, M., Risk, D., Tang, J. W., Todd-Brown, K., and Vargas, R.: The value of soil respiration measurements for interpreting and modeling terrestrial carbon cycling, Plant Soil, 413, 1–25, https://doi.org/10.1007/s11104-016-3084-x, 2017.

Pyrcz, M. J. and Deutsch, C. V.: Geostatistical Reservoir Modeling, Oxford University Press, 433 pp., ISBN: 978-0-19-973144, 2014.

R Code Team: R: A language and environment for statistical computing, R Foundation for Statistical Computing, Vienna, Austria, https://www.R-project.org/ (last access: 1 September 2022), 2021.

Rowlings, D. W., Grace, P. R., Kiese, R., and Weier, K. L.:
Environmental factors controlling temporal and spatial variability in the soil-atmosphere exchange of CO_{2}, CH_{4} and N_{2}O from an Australian subtropical rainforest, Glob. Change Biol., 18, 726–738, https://doi.org/10.1111/j.1365-2486.2011.02563.x, 2012.

Savage, K., Phillips, R., and Davidson, E.: High temporal frequency measurements of greenhouse gas emissions from soils, Biogeosciences, 11, 2709–2720, https://doi.org/10.5194/bg-11-2709-2014, 2014.

Shakoor, A., Shahbaz, M., Farooq, T. H., Sahar, N. E., Shahzad, S. M., Altaf, M. M., and Ashraf, M.: A global meta-analysis of greenhouse gases emission and crop yield under no-tillage as compared to conventional tillage, Sci. Total Environ., 750, 142299, https://doi.org/10.1016/j.scitotenv.2020.142299, 2021.

Storn, R. and Price, K.: Differential Evolution – A Simple and Efficient Heuristic for global Optimization over Continuous Spaces, J. Global Optim., 11, 341–359, https://doi.org/10.1023/A:1008202821328, 1997.

Tallec, T., Brut, A., Joly, L., Dumelié, N., Serça, D., Mordelet, P., Claverie, N., Legain, D., Barrié, J., Decarpenterie, T., Cousin, J., Zawilski, B., Ceschia, E., Guérin, F., and Le Dantec, V.:
N_{2}O flux measurements over an irrigated maize crop: A comparison of three methods, Agr. Forest Meteorol., 264, 56–72, https://doi.org/10.1016/j.agrformet.2018.09.017, 2019.

Tang, X., Liu, S., Zhou, G., Zhang, D., and Zhou, C.:
Soil-atmospheric exchange of CO_{2}, CH_{4}, and N_{2}O in three subtropical forest ecosystems in southern China, Glob. Change Biol., 12, 546–560, https://doi.org/10.1111/j.1365-2486.2006.01109.x, 2006.

Trangmar, B. B., Yost, R. S., and Uehara, G.: Application of Geostatistics to Spatial Studies of Soil Properties, in: Advances in Agronomy, Vol. 38, edited by: Brady, N. C., Academic Press, https://doi.org/10.1016/S0065-2113(08)60673-2, 45–94, 1986.

Ullah, S. and Moore, T. R.: Biogeochemical controls on methane, nitrous oxide, and carbon dioxide fluxes from deciduous forest soils in eastern Canada, J. Geophys. Res., 116, G03010, https://doi.org/10.1029/2010jg001525, 2011.

Vargas, R.:
How a hurricane disturbance influences extreme CO_{2} fluxes and variance in a tropical forest, Environ. Res. Lett., https://doi.org/10.1088/1748-9326/7/3/035704, 2012.

Vargas, R., Detto, M., Baldocchi, D. D., and Allen, M. F.:
Multiscale analysis of temporal variability of soil CO_{2} production as influenced by weather and vegetation, Glob. Change Biol., 16, 1589–1605, https://doi.org/10.1111/j.1365-2486.2009.02111.x, 2010.

Vargas, R., Carbone, M. S., Reichstein, M., and Baldocchi, D. D.: Frontiers and challenges in soil respiration research: from measurements to model-data integration, Biogeochemistry, 102, 1–13, https://doi.org/10.1007/s10533-010-9462-1, 2011.

Vargas, R., Sánchez-Cañete P., E., Serrano-Ortiz, P., Curiel Yuste, J., Domingo, F., López-Ballesteros, A., and Oyonarte, C.:
Hot-Moments of Soil CO_{2} Efflux in a Water-Limited Grassland, Soil Systems, 2, 47, https://doi.org/10.3390/soilsystems2030047, 2018.

Vicca, S., Bahn, M., Estiarte, M., van Loon, E. E., Vargas, R., Alberti, G., Ambus, P., Arain, M. A., Beier, C., Bentley, L. P., Borken, W., Buchmann, N., Collins, S. L., de Dato, G., Dukes, J. S., Escolar, C., Fay, P., Guidolotti, G., Hanson, P. J., Kahmen, A., Kröel-Dulay, G., Ladreiter-Knauss, T., Larsen, K. S., Lellei-Kovacs, E., Lebrija-Trejos, E., Maestre, F. T., Marhan, S., Marshall, M., Meir, P., Miao, Y., Muhr, J., Niklaus, P. A., Ogaya, R., Peñuelas, J., Poll, C., Rustad, L. E., Savage, K., Schindlbacher, A., Schmidt, I. K., Smith, A. R., Sotta, E. D., Suseela, V., Tietema, A., van Gestel, N., van Straaten, O., Wan, S., Weber, U., and Janssens, I. A.:
Can current moisture responses predict soil CO_{2} efflux under altered precipitation regimes? A synthesis of manipulation experiments, Biogeosciences, 11, 2991–3013, https://doi.org/10.5194/bg-11-2991-2014, 2014.

Villarreal, S., Guevara, M., Alcaraz-Segura, D., and Vargas, R.: Optimizing an Environmental Observatory Network Design Using Publicly Available Data, J. Geophys. Res.-Biogeo., 124, 1812–1826, https://doi.org/10.1029/2018JG004714, 2019.

Wang, G. and Chen, S.: A review on parameterization and uncertainty in modeling greenhouse gas emissions from soil, Geoderma, 170, 206–216, https://doi.org/10.1016/j.geoderma.2011.11.009, 2012.

Warner, D. L., Guevara, M., Inamdar, S., and Vargas, R.:
Upscaling soil-atmosphere CO_{2} and CH_{4} fluxes across a topographically complex forested landscape, Agr. Forest Meteorol., 264, 80–91, https://doi.org/10.1016/j.agrformet.2018.09.020, 2019.

Werner, C., Kiese, R., and Butterbach-Bahl, K.:
Soil-atmosphere exchange of N_{2}O, CH_{4}, and CO_{2}and controlling environmental factors for tropical rain forest sites in western Kenya, J. Geophys. Res., 112, D03308, https://doi.org/10.1029/2006jd007388, 2007.

Wu, X., Brüggemann, N., Gasche, R., Shen, Z., Wolf, B., and Butterbach-Bahl, K.:
Environmental controls over soil-atmosphere exchange of N_{2}O, NO, and CO_{2} in a temperate Norway spruce forest, Global Biogeochem. Cy., 24, https://doi.org/10.1029/2009gb003616, 2010.

Yao, Z., Zheng, X., Xie, B., Liu, C., Mei, B., Dong, H., Butterbach-Bahl, K., and Zhu, J.:
Comparison of manual and automated chambers for field measurements of N_{2}O, CH_{4}, CO_{2} fluxes from cultivated land, Atmos. Environ., 43, 1888–1896, https://doi.org/10.1016/j.atmosenv.2008.12.031, 2009.