Ecosystem CO2 release driven by wind occurs in drylands at global scale

Abstract Subterranean ventilation is a non‐diffusive transport process that provokes the abrupt transfer of CO2‐rich air (previously stored) through water‐free soil pores and cracks from the vadose zone to the atmosphere, under high‐turbulence conditions. In dryland ecosystems, whose biological carbon exchanges are poorly characterized, it can strongly determine eddy‐covariance CO2 fluxes that are used to validate remote sensing products and constrain models of gross primary productivity. Although subterranean ventilation episodes (VE) may occur in arid and semi‐arid regions, which are unsung players in the global carbon cycle, little research has focused on the role of VE CO2 emissions in land–atmosphere CO2 exchange. This study shows clear empirical evidence of globally occurring VE. To identify VE, we used in situ quality‐controlled eddy‐covariance open data of carbon fluxes and ancillary variables from 145 sites in different open land covers (grassland, cropland, shrubland, savanna, and barren) across the globe. We selected the analyzed database from the FLUXNET2015, AmeriFlux, OzFlux, and AsiaFlux networks. To standardize the analysis, we designed an algorithm to detect CO2 emissions produced by VE at all sites considered in this study. Its main requirement is the presence of considerable and non‐spurious correlation between the friction velocity (i.e., turbulence) and CO2 emissions. Of the sites analyzed, 34% exhibited the occurrence of VE. This vented CO2 emerged mainly from arid ecosystems (84%) and sites with hot and dry periods. Despite some limitations in data availability, this research demonstrates that VE‐driven CO2 emissions occur globally. Future research should seek a better understanding of its drivers and the improvement of partitioning models, to reduce uncertainties in estimated biological CO2 exchanges and infer their contribution to the global net ecosystem carbon balance.


| INTRODUC TI ON
The eddy-covariance (EC) technique is a powerful tool (Aubinet et al., 2012;Baldocchi, 2014) that enables the continuous nondestructive monitoring of mass, energy and momentum exchange between the land and the atmosphere (Baldocchi et al., 2001). Its use helps the scientific community to reach a better understanding of how terrestrial ecosystems interact with the Earth's climate system, an outstanding scientific challenge in the global change context (Reich, 2010;Solomon et al., 2007). With a growing number of EC towers worldwide, a set of regional (e.g., Ameriflux, CarboEurope, AsiaFluz, OzFlux; Aubinet et al., 1999;Isaac et al., 2016;Novick et al., 2018;Reich, 2010;Yu et al., 2006) and global networks (FLUXNET; Pastorello et al., 2020) of EC flux measurement stations has been established. However, although arid and semiarid lands comprise approximately 40% of the global terrestrial surface (Reynolds et al., 2007) and play an important role in the global carbon (C) cycle (Ahlström et al., 2015;Lal, 2004), observational gaps (Baldocchi, 2008;Bell et al., 2012;Schimel et al., 2015;Xiao et al., 2008) have hampered a full understanding of their functional behavior.
As a result, the interpretation of net CO 2 fluxes derived from EC measurements has been constrained by available data and research, which has been extensively developed in temperate ecosystems where CO 2 fluxes have commonly been attributed purely biological processes (i.e., gross primary production, GPP, and ecosystem respiration, R eco ) that occur within the ecosystem boundaries (Chapin et al., 2006; i.e., from the upper soil layers to the close-to-surface atmosphere). However, in water-limited ecosystems, there is evidence of other processes, besides biological processes GPP and R eco , that can significantly contribute to net CO 2 flux observed by the EC technique, at least on short time scales (Emmerich, 2003;Mielnick et al., 2005;Serrano-Ortiz et al., 2010). These include photodegradation (Rutledge et al., 2010), geochemical weathering (Hamerlynck et al., 2013), and subterranean ventilation (Kowalski et al., 2008).
Neglecting these processes when partitioning net CO 2 fluxes can result in biased estimates of GPP and R eco fluxes (Ferlan et al., 2011;Inglima et al., 2009;Mielnick et al., 2005;Were et al., 2010;Xie et al., 2009), which are extensively utilized to calibrate models and satellite products (Jones & Cox, 2005;Ma et al., 2015;Verma et al., 2014) to assess current and future coupling between the terrestrial biosphere and the climate system (Naidu & Bagchi, 2021;Poulter et al., 2014). Here we pay special attention to subterranean ventilation, whose relevance is suggested to be notable compared to other contributing processes, particularly at those ecosystems where water scarcity over the drought season inhibits biological processes (López-Ballesteros et al., 2017;Pérez-Priego et al., 2013).
Subterranean ventilation is a non-diffusive transport process that provokes the abrupt transfer of CO 2 -rich air from the vadose zone to the atmosphere under drought and high-turbulence conditions (Kowalski et al., 2008). Recent studies developed in drylands have pointed out the importance of the vadose zone as a significant temporary CO 2 pool, given its capacity to store vast amounts of gaseous CO 2 in interconnected cracks, pores, and cavities belowground where volumetric CO 2 fractions can rise to two orders of magnitude greater than in the atmosphere (Benavente et al., 2010;Faimon et al., 2020;Fernandez-Cortes et al., 2015;Sánchez-Cañete et al., 2016). Regardless of its originating process, the enhanced density of CO 2 -rich air (Kowalski & Sánchez-Cañete, 2010)-due to its large molecular mass-may enable it to descend from shallow soil layers, where most root and microbial respiration takes place (Buchmann, 2000;Coleman & Crossley, 2004, down to the water table level, creating a concentration gradient that increases with depth. In order for subterranean ventilation episodes (VE) to happen, several conditions are required. First, low values of soil water content must be reached to allow gas flow across the porous medium between the vadose zone and the atmosphere (Cuezva et al., 2011;Risk et al., 2002;Roland et al., 2013;Serrano-Ortiz et al., 2010). This condition commonly occurs during daytime and especially during the dry season coinciding with the dormancy period, when water scarcity and high temperatures constrain biological activity (i.e., CO 2 emissions are likely inconsistent with simultaneous variations of organic carbon pools; Huxman et al., 2004). Conversely, during nighttime, atmospheric stability, water deposition, and vapor adsorption near the soil surface may limit this belowground and aboveground air interconnection. Second, high atmospheric turbulence conditions are indispensable to penetrate the vadose zone and produce the abrupt release of CO₂ from the vadose zone to the atmosphere via non-diffusive transport (Kowalski et al., 2008;Subke et al., 2003;Takle et al., 2004). For this reason, the wind velocity (or friction velocity) is considered the main trigger of this phenomena (Flechard et al., 2007;Maier et al., 2010;Redeker et al., 2015). Third, the air located in the vadose zone must be significantly CO 2 rich (a high storage term is needed) to make VE detectable for the EC system.
While the first and second conditions (i.e., low soil water content, and high turbulence) can be considered as predisposing factors enabling VE to occur, the third condition (i.e., vadose zone CO 2 concentration) directly affects the magnitude of VE-driven CO 2 emissions (Lang et al., 2017;Sanchez-Cañete et al., 2011).
These studies were mainly developed in Southwestern Europe and utilized specific experimental setups designed to carefully monitor the influencing factors mentioned above. However, the relevance of subterranean ventilation at global scale is still unknown. The general aim of this study is to test whether VE happen globally by analyzing flux data provided by FLUXNET and other related regional networks.
Our main hypothesis is that VE occur in water-limited ecosystems worldwide, and particularly during the dry season. To test our hypothesis, we designed an algorithm, based on previous scientific evidence, to systematically detect VE in the 145 open ecosystems selected for this study. Therefore, our specific objectives are as follows: (1) to detect daytime VE-driven CO 2 emissions in different ecosystems distributed around the world and (2) to investigate factors influencing VE.

| Database selection
To identify VE we used in situ, quality-controlled eddy-covariance observations of carbon dioxide fluxes and ancillary data from FLUXNET and main regional EC networks that are available as open data to the science community. The datasets used in this study include the FLUXNET2015 (https://fluxn et.fluxd ata.org), the AmeriFlux (http://ameri flux.lbl.gov), the OzFlux (http://www.ozflux. org.au), and the AsiaFlux (http://www.asiaf lux.net) databases. The basic instrumentation of EC technique consists of a sonic anemometer (to measure 3D wind velocity and virtual temperature at high frequency), an infrared gas analyzer (to measure CO 2 and water vapor densities at high frequency), and complementary sensors to measure other meteorological ancillary variables at the flux site (e.g., air temperature, relative humidity, net radiation, precipitation). The flux and meteorological half-hourly data used in this analysis corresponded to quality-controlled data directly measured (i.e., minimum processing level). Therefore, we excluded post-processed gap-filled data (Papale & Valentini, 2003;Reichstein et al., 2005) from our analysis.

| Data processing method
We designed a novel algorithm to automatically detect VE-driven CO 2 emissions among the analyzed flux data. Our algorithm aims at detecting VE over periods when subterranean ventilation may prevail over other processes. The algorithm was designed based on previous research reflected in the following assumptions: (1) atmospheric turbulence must be sufficient to pump CO 2 -rich air from the vadose zone to the atmosphere; (2) no nocturnal VE happens due to soil rehumidification; and (3) the required variability in air temperature and soil moisture must be insufficient to explain variations of net CO 2 emissions (i.e., incompatible with an ecophysiological interpretation of fluxes, such as the Birch effect at the end of the dry season).
A temporal window is needed to statistically discriminate subterranean CO 2 release due to changes in air temperature (e.g., with the passage of high-and low-pressure systems and fronts) or soil moisture. Although VE is an abrupt process (the smallest time window over which a VE could be identified is <1 s; Kimball & Lemon, 1970;Massman et al., 1997;Mohr et al., 2016;Takle et al., 2004), considering previous research (Sanchez-Cañete et al., 2011), a period of 5 days with sustained high turbulence, similar mean air temperature and dry conditions were considered as unequivocal to detect VE. Thus, we applied the algorithm to 5-day intervals over 1 year for each selected site.
Furthermore, to reduce the potential effect of data availability on the number of VE detected our algorithm, only 5-day intervals with a minimum number of data (N > 40) and CO 2 maximum quality (flag qc = 0) were used. Finally, we gathered an indicator of data availability for each site and year that were analyzed (see Table S1).
We computed Partial Spearman correlation coefficients over each 5-day period filtered to eliminate spurious correlation effects between F c (μmolCO 2 m −2 s −1 ) and ancillary data. Ancillary variables considered in partial correlation were as follows: air temperature, vapor pressure deficit, soil temperature, atmospheric pressure, soil water content, incoming photosynthetic photon flux density, incoming shortwave radiation and friction velocity. Only 5-day intervals with partial Spearman correlation coefficients between F c and u * above 0.2 and p < .05 (support the evidence of the alternative hypothesis, i.e., non-zero partial correlation, being significantly different from the null hypothesis, i.e., zero partial correlation) were considered as VE. The software Matlab was used for statistical analyses (Matlab R2017a).
To balance the different dataset duration among analyzed sites as well as to simplify our analysis, results shown in this study correspond to 1 year of data per site. We have considered the assumption that experimental sites where VE were not detected after analyzing 4 years of data are not VE predisposed sites (in accordance with our algorithm design). Thus, the algorithm was performed to the available datasets under the next conditions: (1) for site-specific databases lasting from 1 to 4 years, the algorithm was applied to the whole database; (2) for site-specific databases lasting from 5 to 6 years, the algorithm was applied to the first four consecutive years; and (3) for site-specific databases lasting more than 6 years, the algorithm was applied to the first four non-consecutive years.
In those sites where VE-driven CO 2 emissions were detected over several years, the year finally selected corresponded to the one with more VE detected. The list of FLUXNET and regional EC networks sites used in this study with respective basic information appears in the supporting information (Tables S1-S3).

| RE SULTS
Of the 145 sites analyzed in our study, we found that 34% (50 sites) The mean annual precipitation was below 1000 mm at all the experimental sites with VE detected ( Figure 1B). The mean annual temperature ranges between −1 and 26°C, although two ecosystems (Ru-Cok and US-A74 with −14.3°C and 33.9°C, respectively) far exceeded this range. The ecosystem types with more VE corresponded to grasslands, croplands, and open shrublands, in that order. Of the remaining ecosystems analyzed, only one (in savannas and barrens or sparsely vegetates sites) and two places (in woody savannas and closed shrublands sites) showed VE. The F I G U R E 1 (a) World map distribution of FLUXNET2015, AmeriFlux, OzFlux, and AsiaFlux sites analyzed in this study with distinction of experimental sites with (red marker) and without (white marker) VE-driven CO 2 emissions detected by our algorithm. Map colors indicate world dryland areas dataset classification (UNEP-WCMC; Sorensen, 2007) according to dryness index PET/P (potential evapotranspiration/ precipitation). Codes of sites with VE-driven CO 2 emissions have been included; (b) distribution of analyzed sites with (red) and without (white) VE emissions along mean annual temperature (MAT) and mean annual precipitation (MAP) gradients; (c) number of analyzed sites with (red) and without (white) VE emissions detected by IGBP land cover classification (grassland, GRA; savanna, SAV; woody savanna, WSA; open shrubland, OSH; closed shrubland, CSH; barren or sparsely vegetated, BSV; and cropland, CRO). CBD is the acronym for convention on biological diversity. supporting information (S1, S2, and S3) gives more details for each experimental site. percentiles. This SWC threshold was also evaluated in the other 46 experimental sites with VE and was fulfilled in sites classified as "dry subhumid," "semiarid," and "arid" in the dryland type classification (UNEP-WCMC). To point this assumption, the US-Ctn site (Figure 2a), which showed a progressive decrease in SWC values of one order of magnitude over a 15-day time window, has been graphed during the "pre-emission" and the beginning of the "emission" period. When the SWC reached its 30th percentile in US-Ctn, our filter detects the first 5-day periods of emissions . In contrast, ES-Amo, ZA-Kru and AU-DaP sites (Figure 2b-d), where SWC values were less variable, have been graphed during the "emission" period with more VE detected on previous days (not shown in Figure 2). VE detected in the four experimental sites selected were directly correlated with u * and the shortwave incoming radiation SW_IN ( Figure 3). The value of the F c -u * partial correlation or the F c -SW_IN partial correlation varies among the different VE detected; however, F I G U R E 3 Box-and-whisker plots of spearman coefficients obtained from partial correlation between CO 2 fluxes (F c ) and friction velocity (u * ), air temperature (TA), vapor pressure deficit (VPD), atmospheric pressure (PA), soil temperature (TS), soil water content (SWC) or shortwave radiation incoming (SW_IN) for the 5-day periods when VE-driven CO 2 emissions were detected at ZA-Kru, ES-Amo, US-Ctn, and AU-DaP sites. Each Spearman correlation coefficient obtained was not influenced by the indirect effect of the remaining variables considered in the analysis. Each data shown have an individual associated p-value <.05 (i.e., partial correlation coefficient different from zero). the median was similar among ecosystems (between 0.44 and 0.55 the F c -u * correlation and 0.45-0.67 the F c -SW_IN correlation).
During these periods, there were also positive and negative correlations with other variables measured in the ecosystems: air temperature (TA), vapor pressure deficit (VPD), atmospheric pressure (PA), soil temperature (TS), or/and soil water content (SWC). However, these correlations were more heterogeneous among ecosystems, weaker, and not always significant.
In the 50 sites with VE-driven CO 2 emissions, our algorithm detected 290 VE according to the established criteria. The number of VE detected were positively related to the aridity degree represented via the dryland type classification (UNEP-WCMC) and negatively related to the normalized SWC. The Spearman partial correlation coefficients obtained between F c and u * ranged from 0.2 to 0.8 with a median value around 0.5 (Figure 4) Table S1.

| DISCUSS ION
Our results suggest that CO 2 emissions driven by subterranean ventilation occur globally. After analyzing 145 open ecosystems (Bond, 2019) selected from main EC networks, our algorithm, designed upon previous scientific evidence, detected ventilation episodes (VE) in 34% of them ( Figure 1a, and Tables S1 and S2). Our main hypothesis was corroborated as 84% of the analyzed sites corresponded to open dryland ecosystems based on the IGBP classification and the UNCCD definition of Drylands (Sorensen, 2007). The detection of VE was quite balanced among the different kinds of drylands presented (50% in arid, 56% in semiarid, 57% in dry subhumid, and 64% in the additional dryland areas). During these VE, our results show that the relationship between net CO 2 fluxes and air and soil temperature, vapor pressure deficit, atmospheric pressure, and soil water content were more heterogeneous among ecosystems,

F I G U R E 4
Box-and-whisker plots of Spearman coefficients obtained from partial correlation between CO 2 fluxes (F c ) and friction velocity (u * ) for all the flux sites where 5-day VE were detected by our algorithm. Each Spearman correlation coefficient obtained was not influenced by the indirect effect of the rest of covariables considered: Air temperature, vapor pressure deficit, atmospheric pressure, soil temperature, soil water content, and shortwave radiation incoming; for the 5-day periods when VE-driven CO 2 emissions were detected at each experimental site. Sites cluster by ecosystem type [land cover (IGBP)] and Köppen climate classification (KGCC). Each data shown have an individual associated p-value <.05 (i.e., partial correlation coefficient different from zero). weaker, and not always significant, compared to incoming shortwave radiation and atmospheric turbulence (Figure 3). These results could be explained by the variability in the predisposing factors among VE occurred in a same site and suggest the predominance of subterranean ventilation over biological processes (i.e., GPP or R eco ) during drier periods.
The occurrence and magnitude of the VE-driven CO 2 release from the vadose zone to the atmosphere were directly related with the intensity of the atmospheric turbulence (u * value; Figures 2-4).
The Spearman coefficients obtained from partial correlation between F c and u * found among the climate-IGPB categories ranged from 0.2 and 0.8 with a median value around 0.5 (Figure 4). These results are in accordance with the assumption stating that eddies promoted by high atmospheric turbulence conditions can penetrate the soil matrix and produce an abrupt emission of CO 2 -rich air from the vadose zone to the atmosphere via non-diffusive transport (Kowalski et al., 2008;Subke et al., 2003;Takle et al., 2004).  Cuezva et al., 2011;Kosmas et al., 2001;Verhoef et al., 2006).
Most VE occurred during dry conditions, with SWC values below the 30th percentile of site-specific annual soil moisture in sites classified as "dry subhumid," "semiarid," and "arid" in the dryland type classification (UNEP-WCMC). Such behavior can be clearly observed in US-Ctn site (Figure 2a), where the first 5-day period of emissions  was detected when the SWC decreased down its 30th percentile during daytime and over the dry season. Due to the lack of water inputs, the interconnectivity between vadose zone and the atmosphere can allow for an effective gas transport across the land surface (Cuezva et al., 2011;Risk et al., 2002;Roland et al., 2013;Sanchez-Cañete et al., 2011). This dryness condition can correspond to an inherent feature of the ecosystem (i.e., US-Ctn, AU-DaP or ES-Amo) or to a climate extreme, such as a drought event (i.e., ZA-Kru).
It is also remarkable how the mean annual precipitation is a limiting factor, which was below 1000 mm in all the experimental sites with VE detected (Figure 1b).
The shortwave incoming radiation (SW_IN) also constitutes an important meteorological variable involved in the VE-driven CO 2 emissions detected (Figure 3) probably because solar heating is an important mechanism that triggers turbulence (López-Ballesteros et al., 2017;Stull, 1988). In relation with the VPD variable, López-Ballesteros et al. (2017) found that VE-driven CO 2 emissions coincided with high turbulence and high VPD conditions; however, our results do not clearly show such positive VPD-VE relationship.
In relation with the F c -pressure correlation, several processes would be involved. Pressure pumping (Massman et al., 1997;Mohr et al., 2016;Roland et al., 2015), pressure tides (Clements & Wilkening, 1974;Kimball & Lemon, 1970;le Blancq, 2011), and subterranean ventilation are non-diffusive transport processes affected by atmospheric pressure changes at different temporal scales. Enabling conditions for all of these processes are a high subsoil-atmosphere interconnectivity (López-Ballesteros, Oyonarte, et al., 2018;Pérez-Priego et al., 2013;, which, in turn, relies on low SWC values and high atmospheric turbulence, as these are indispensable to allow for the CO 2 -rich air transfer from the vadose zone to the atmosphere. For this reason, it is challenging to distinguish among pressure pumping, pressure tides, and subterranean ventilation. Some authors have shown a positive relationship between pressure pumping and ventilation (Mohr et al., 2016;Nachshon et al., 2012;Redeker et al., 2015), based on the correlation shown between pressure changes and wind perturbations at high frequencies (Maier et al., 2010;Redeker et al., 2015;Takle et al., 2004). However, other authors (Moya et al., 2019; found a temporal and spatial decoupling between VE and pressure tides. Previous research states that pressure tides were produced in the sub-surface at low frequencies by large-scale atmospheric dynamics (Elberling et al., 1998;Guoqing, 2005;Lindzen, 1979), while VE is produced at high frequencies and limited to the shallower horizon (Redeker et al., 2015;Takle et al., 2004). Regarding subsoil-air interconnectivity, we did not find any relation between site-specific soil type and occurrence of VE (data not shown); however, low SWC values are related to a minimum gas permeability of the soil matrix During growing periods, sufficient soil moisture enables autotrophic and heterotrophic biological activity and inhibits VE by lowering subsoil-atmosphere interconnectivity. However, during dry periods, the low soil moisture reduces biological processes involved in the absorption and production (and consequent emission) of CO 2 to a secondary role, allowing VE to become more relevant. Therefore, the high inter-annual variability of F c in ES-Amo and ZA-Kru could be a consequence of the frequency and magnitude of VE. ES-Amo acted as an annual net source of CO 2 ranging from +163 g C m −2 yr −1 to +324 g C m −2 yr −1 over (López-Ballesteros et al., 2017, whereas ZA-Kru varied from −138 to +155 g C m −2 yr −1 over 2000-2005(Archibald et al., 2009. Differences in the F c balance among years depended on the length and severity of the dry season, which was determined by the magnitude and timing of precipitation events (Gilmanov et al., 2006;Luo et al., 2007;S. Ma et al., 2007;Meyers, 2001;Pereira et al., 2007). Future climate projections suggest an increase in long-term droughts in the Mediterranean, west African, central Asian, and central American regions (Sheffield & Wood, 2008). Based on our results, this trend might be related to a higher occurrence of VE among these areas as minimum soil moisture is a prerequisite for VE to occur. However, it is not sure that drier drylands (in terms of drought duration and severity) will release more CO 2 to the atmosphere via subterranean ventilation, as one of the main factors that presumably determine the magnitude (and occurrence) of this ventilative CO 2 fluxes is the availability of CO 2 in the air stored in the vadose zone to be ventilated (Lang et al., 2017;Sanchez-Cañete et al., 2011).
Vadose zone can be considered as a spatial and temporal depot for CO 2 coming from different processes (Bourges et al., 2001;Chapin et al., 2006;Serrano-Ortiz et al., 2010). Depending on (1) the dimension and potential capacity of vadose zone subterranean CO 2 storage, (2) the occurrence of VE, (3) the simultaneous occurrence of other sink (e.g., dissolved CO 2 infiltration) and source mechanism (e.g., soil respiration), (4) changes in the drivers controlling the VE potential, and (5) changes in the degree of connection between the cavities and the aboveground system; the vadose zone net CO 2 balance could be negative and the subterranean systems could be degassed. In this situation, even if predisposing factors (low values of soil water content and high atmospheric turbulence conditions) are present, the lack of subterranean CO 2 available could disabled the occurrence of VE. Unfortunately, due to the absence of observations, we could not assess how vadose zone CO 2 concentration varies and affects VE-driven CO 2 emissions among the analyzed sites.
In general, vadose zone CO 2 may not necessarily be produced simultaneously and within the ecosystem where VE are detected.
Based on previous research, three hypotheses could explain potential origins of the ventilated CO 2 .
The first hypothesis (H1) states that CO 2 is produced in situ (i.e., within the ecosystem boundaries; Chapin et al., 2006) by microbial and root respiration, with delayed release to the atmosphere.
In this regard, apparent respiratory quotient (ARQ) measurements performed in two semiarid ecosystems revealed that actual soil CO 2 efflux only accounts for 1/3 and 1/4 of the total CO 2 emissions derived from soil respiration (Angert et al., 2015;Sánchez-Cañete et al., 2018). These studies demonstrated that not all CO 2 produced by microbial and root respiration is immediately released to the atmosphere. Potential processes involved in this discrepancy are suggested to be aqueous phase partitioning (Kern, 1960), calcite dissolution reactions , gravitational percolation due to a higher density , or CO 2 dissolution in xylem water (Aubrey & Teskey, 2009). However, in the ES-Amo site included in the present study, annual VE-driven CO 2 emissions measured systematically during six consecutive years were close to the estimated total soil organic carbon (SOC), suggesting that the origin of the ventilated CO 2 could not be explained by previous in situ biological activity (López-Ballesteros et al., 2017).
The second hypothesis (H2) states that CO 2 can be produced ex situ (i.e., out of the ecosystem boundaries) by respiration but transported laterally in aqueous and gaseous CO 2 forms. In this sense, Li et al. (2015) found in a closed arid basin that biological CO 2 produced at higher altitudes were leached as dissolved inorganic carbon (DIC) and accumulated in the groundwater beneath the basin center.
Finally, the third hypothesis (H3) relates the CO 2 production to non-biological processes. With regard, some estimations of CO 2 fluxes caused by photodegradation of senescent organic matter correspond to low CO 2 effluxes of <0.2 μmol m −2 s −1 (Austin & Ballaré, 2010;Austin & Vivanco, 2006;Brandt et al., 2009;King et al., 2012;Rutledge et al., 2010). Likewise, short-term estimates of geochemical weathering and calcite precipitation seemed to correspond to very low CO 2 effluxes of ca. 0.05 μmol m −2 s −1 (Goddéris et al., 2006;Hamerlynck et al., 2013;Roland et al., 2013;Serrano-Ortiz et al., 2010). Thus, these two processes are far from sufficient to explain the large magnitude of the daytime CO 2 emissions detected by our algorithm (Figure 2 and Table S1). However, karstic areas deserve special attention given their great capacity to temporally store large amounts of gaseous CO 2 belowground (Decarlo & Caylor, 2020;Emmerich, 2003;Serrano-Ortiz et al., 2010) within their carbonate vadose zone, which usually contains an interconnected system of caves, macropores, and fissures. According to Faimon et al. (2020), these places can be considered as "breathing spots" that may be relevant in the carbon cycle at basin scale. Finally, another non-biological process that might be involved in vadose zone CO 2 production is related to deep-seated CO 2 seepage from porous reservoirs usually related to volcanic or hydrothermal activity (Burton et al., 2013;Chiodini et al., 2008;Lewicki et al., 2003;Mörner & Etiope, 2002;Rey et al., 2012). However, to know how each mentioned process interact with VE, we need further investigation at site level.
The EC databases used in this study presented some other limitations to provide more in-depth information on the global VE relevance. These limitations include the following: (1) the presence of long-term gaps in many databases, which made it difficult to calculate the contribution of VE at the annual scale; (2) the limited length of the time series (the majority of EC flux sites analyzed are usually limited to only 3-5 years of measurements) could be an issue (Schimel et al., 2015;Sulkava et al., 2011;Sundareshwar et al., 2007) to determine the effect of VE in the interannual variation of the ecosystem C balance; (3) the heterogeneity among databases (absence of some variables collected and differences in data quality) and the lack of information and/ or parameter misspecification in the metadata made it difficult to establish relations between VE and other variables across the different types of ecosystems; and (4) the spatial representativeness in the EC databases (López-ballesteros, Beck, et al., 2018;Schimel et al., 2015;Sulkava et al., 2011;Sundareshwar et al., 2007) was not sufficient to capture the natural variability of climatological and biological conditions in some geographical regions. FLUXNET and regional network sites are mainly located in ecosystems with annual temperatures between 5 and 17°C, annual rainfall between 600 and 1250 mm, and latitudes above 30°N (Vargas et al., 2013).
Thus, regions more vulnerable to environmental or climate change are still underrepresented (Bell et al., 2012). Many of these regions are drylands or prone to desertification (Emmerich, 2003;Eswaran et al., 2000;Huenneke et al., 2002;Schlesinger, 1990), which makes VE potentially relevant now and even more so in the future (Ahlström et al., 2015). Furthermore, the lack of CO 2 concentration measurements within the vadose zone limited our analyses, as this is the main factor driving the magnitude of VEdriven CO 2 emissions (Sanchez-Cañete et al., 2011;Serrano-Ortiz et al., 2010).
These issues create uncertainties that severely limit our ability to better understand the occurrence and magnitude of VE and to make confident estimations of the contribution of VE in the global terrestrial carbon balance. Based on these uncertainties, the algorithm was designed with restrictive assumptions that may result in the underestimation of the total impact of CO 2 emissions released by VE, but that helped us test whether ventilation episodes occur globally (i.e., the main objective of our study). Furthermore, the incipient state of knowledge of this phenomenon hampers its detection when other processes simultaneously contribute to the net CO 2 flux observed by the EC system. Future research should involve a better characterization of belowground CO 2 production and transport processes by tracking the CO 2 concentration and isotopic signal within the vadose zone of ecosystems where VE are likely to happen. The improvement of partitioning models would also reduce uncertainties in estimated biological CO 2 exchanges.
Nevertheless, the results achieved in this work are very relevant because the presence of VE has been identified in many arid and semiarid lands, which constitute the largest biome in the world (Schimel, 2010).

| CON CLUS IONS
Due to the wide presence of arid and semiarid ecosystems, drylands' C balance strongly affects the inter-annual variability of C dynamics at a global scale. In this sense, a more detailed information of drylands behavior is necessary to predict how climate change will affect their ecological structure and functions. The present study is a step toward a better understanding of the processes involved in the global C cycle. After analyzing eddycovariance open data observations of C fluxes and ancillary data from 145 common semiarid ecosystem types from FLUXNET and its main regional networks, we can confirm that subterranean ventilation CO 2 emissions occur globally.
The occurrence and magnitude of CO 2 emissions from the vadose zone towards the atmosphere were directly related to changes in friction velocity (indicator of atmospheric turbulence), the triggering variable. The soil water content and annual rainfall were found to be limiting factors. Only when the soil water content declines to values below the 30 percentile of site-specific annual soil moisture, could the interconnectivity between vadose zone and the atmosphere allow for an effective gas transport across the land surface. These subterranean ventilated CO 2 emissions occurred in dry ecosystems, but also in ecosystems where dryness was not an inherent feature of the ecosystem. Subterranean ventilation episodes varied from days to weeks, and occurred only during daytime. The Spearman coefficients obtained from partial correlation between u * and CO 2 fluxes were found to be similar (with a median value around 0.5) across the different climatic zones analyzed. Finally, the CO 2 stored belowground available to be ventilated directly affects the occurrence and magnitude of CO 2 emissions induced by VE.
Our results suggest that our algorithm can efficiently detect CO 2 emissions produced by VE and constitutes a useful tool to its detection, for use in future research. Global warming and anthropogenic disturbance will increase the rate of drying in ecosystems leading to more frequent and intense drought events. This drying could provoke an ecosystem degradation favoring a greater exposure of subsoil CO 2 to ventilation and affecting negatively in the trend of global carbon sequestration. Improving our knowledge of the carbon exchanges between terrestrial ecosystems and the atmosphere is essential to better understand and model the Earth's climate system. Understanding the occurrence of processes, feedbacks, and driving factors that modulate the carbon (C) source capacity of natural ecosystems due to ventilation episodes is needed to advance towards more robust model projections for future climate in order to predict how climate change will affect biology structure and functions as well as more adequate design of mitigation policies.

ACK N OWLED G M ENTS
This work used eddy-covariance data acquired and shared by the EC global and regional networks FLUXNET, AsiaFlux, OzFlux, and AmeriFlux. The individual sites DOIs citation is available in Tables S1-S3. We acknowledge the support from Ministry of

CO N FLI C T O F I NTE R E S T
The authors declare no competing financial interests.

DATA AVA I L A B I L I T Y S TAT E M E N T S
The algorithm's Matlab code is available at https://figsh are.com/ s/1b670 62a4c c4b88 0b9a7. Data from experimental sites are freely available at FLUXNET2015 (https://fluxn et.fluxd ata.org), AmeriFlux