Disentangling the drivers of soil CO₂ ventilation in a Mediterranean dryland using in situ and remote sensing techniques

Subterranean CO₂ concentrations are driven by complex interactions between biological and physical processes. In semiarid ecosystems, atmospheric processes can play a relevant role in modulating soil CO₂ storage and release. In the current study, a multi-instrumental dataset, collected in a Mediterranean shrubland in southern Spain, was analyzed, and the main atmospheric drivers controlling soil CO₂ and radon (Rn) dynamics were investigated. Based on a precise methodology, 10 significant ventilation events were detected, and the Spearman correlation coefficients between the soil CO₂ and Rn concentrations and the different atmospheric variables were calculated.

The results identified surface atmospheric pressure as the most consistent and independent driver across the events, exhibiting strong negative correlations with the subterranean CO₂ and Rn concentrations. Surface-level friction velocity (u_∗), boundary-layer turbulent kinetic energy dissipation rate (ϵ) and wind shear (sh) showed significant positive correlations. However, their independence was not consistent when compared with diluting ventilation events, when u_∗ was more relevant, with enriching ventilation periods, that were more influenced by boundary-layer ϵ and sh. In contrast, at lower altitudes ϵ, sh, atmospheric boundary layer height and mixing layer height were less strongly correlated with soil CO₂ and Rn concentration changes.

These findings provide new insights into the mechanisms that promote soil-atmosphere transport in drylands, especially those regarding the carbon cycle, and highlight the need to incorporate such mechanisms into Earth system models to improve carbon cycle predictions under future climate scenarios.

1 Introduction

The current context of climate change and increasing CO₂ concentrations in the atmosphere (IPCC, 2021; WMO, 2024) highlights the importance of characterizing CO₂ emission and sequestration mechanisms between ecosystems and the atmosphere, since CO₂ is the second most important greenhouse gas after water vapor (GHG; Schimel, 1995). Terrestrial ecosystems were found to be significant CO₂ net sinks, with increasing CO₂ uptake in recent decades (e.g., Piao et al., 2018). Additionally, terrestrial uptake is the main driver of interannual variability of atmospheric CO₂ (IPCC, 2021), especially in semiarid ecosystems (e.g., Ahlström et al., 2015) where the contribution to interannual variability is larger. The eddy covariance (EC) technique has emerged as a powerful tool to track CO₂ exchanges, identifying and quantifying the mechanisms involved (e.g., Baldocchi, 2014). CO₂ fluxes obtained with the EC technique have been traditionally associated with biological processes, such as respiration and photosynthesis (e.g., Pastorello et al., 2020), while abiotic geochemical and mechanical processes (e.g., Brantley et al., 2008; Serrano-Ortiz et al., 2010; Moya et al., 2022) were traditionally neglected. However, the importance of these mechanisms in the carbon cycle has been proven over carbonate terrains (Sánchez-Cañete et al., 2011; López-Ballesteros et al., 2017).

Carbonate-rich substrates are commonly found in dryland regions and represent the world's largest carbon reservoir (Liu and Zhao, 2000), while they are filled with pores and even cavities across the subterrain. These cavities are formed by the infiltration of water from the surface through weathering processes (Brantley et al., 2008), particularly after rain events, which dissolves inorganic carbon present in the surface, transporting it deeper. At depth, the dissolved inorganic carbon may precipitate, releasing CO₂ into the soil atmosphere (Spötl et al., 2005), and thereby increasing CO₂ levels and creating significant underground soil carbon pools in the vadose zone. In addition, biotic processes such as root respiration and microbial activity, along with abiotic influences such as geological CO₂ sources and the downward movement of dense, CO₂-rich air masses (Sánchez-Cañete et al., 2013), contribute to elevated subsurface CO₂ levels. The combined effect of these processes produces increases of underground CO₂ storage at concentrations 10 to 100 times greater than in the atmosphere (e.g., Faimon et al., 2020).

Once stored, CO₂ is transported within the soil or emitted to the atmosphere through two main processes: diffusion and ventilation. Mixing by random motions causes CO₂ transport down the gradient from enriched soil to the more diluted atmosphere through molecular and turbulent diffusion. However, CO₂ can be aerodynamically controlled and vented through mechanical airflow, normally during dry periods (Kowalski et al., 2008, Sánchez-Cañete et al., 2011), when the soil pores are not saturated by water. Subterranean ventilation is known to occur mainly over drylands (Moya et al., 2022), impacting CO₂ levels in the atmosphere and modifying the net ecosystem carbon balance (NECB) at significant levels both in the short term (e.g., Serrano-Ortiz et al., 2010; Sánchez-Cañete et al., 2011, 2013) and in the long term (e.g., López-Ballesteros et al., 2017). These processes can pump fresh atmospheric air into the pores and cavities, displacing the subterrain CO₂-rich air masses and producing subsequent release of CO₂ to the atmosphere, and consequently an increase in atmospheric CO₂ concentrations (Sánchez-Cañete et al., 2011). The emission is measurable with the EC technique and has been related at different sites to high friction velocities (Sánchez-Cañete et al., 2011) and atmospheric pressure differences (Sánchez-Cañete et al., 2013; Moya et al., 2019). Soil and cavities structures, and consequently soil ventilation, are highly heterogeneous (Sánchez-Cañete et al., 2016) making the interaction between the soil and the atmosphere very complex to characterize. Additionally, diurnal fluctuations in subsurface CO₂ concentrations have been identified (e.g., Moya et al., 2019). The amplitude of these fluctuations decreases with depth, indicating the diminishing influence of atmospheric processes (Sánchez-Cañete et al., 2013).

Ventilation events affect concentrations of numerous species. While molecules like CO₂ or CH₄ exhibit high interactivity with biological processes, others like ²²²Rn (hereinafter Rn) are mostly inert. Thus, CO₂ levels are highly affected by biological processes, while Rn natural production process is exclusively radioactive decay from ²²⁶Ra present in the Earth's crust (Moed et al., 1988). However, the amount of Rn released into the pore air (emanation) can be influenced by soil moisture content (Zhuo et al. 2006). Thus, in poorly vented pores and cavities, Rn levels are expected to accumulate at a constant rate, while a forced change in its concentration will be caused by ventilation of the pores and cavities (Kowalski et al., 2008; Pla et al., 2023). Consequently, Rn is an adequate proxy to track natural CO₂ ventilation events (e.g.,Sánchez-Cañete et al., 2011). In atmospheric studies, the Rn concentration is typically tracked by its radioactivity level (usually in Bq m⁻³) since they are directly proportional.

Worldwide networks have been established to track CO₂ exchanges, such as FLUXNET (Pastorello et al., 2020), which also tracks water vapor and energy exchanges. In the case of CO₂, the technique and instrumentation have been improved with time, while CO₂ exchanges have been linked exclusively to biological processes. However, the coverage of this network is not homogeneous across the world, and some regions lack representativity (Villarreal et al., 2018). This might have influenced the fact that communities like FLUXNET have underestimated abiotic processes specific to carbonate terrains, such as carbonaceous rock dissolution and soil ventilation, as these processes might not be significant in the majority of the stations but should be considered for others. A proper characterization of the processes involved in soil CO₂ ventilation should be conducted, in order to improve CO₂ exchange models and produce reliable products.

The surface and the atmosphere interact constantly with each other, from daily boundary layer development due to surface irradiation to the exchanges of energy, molecules, or particle fluxes (Stull, 1988). Atmospheric dynamics play a crucial role in the dispersion of atmospheric constituents and the importance of its dynamics has been proven to be a driver of soil ventilation. The characterization of the atmospheric dynamics for turbulent exchange investigation has been traditionally measured with in-situ instrumentation, with sonic anemometers installed atop towers (e.g., Jones and Smith et al., 1977). Recently, new instrumentation based on the lidar technique has opened opportunities to characterize the upper atmosphere and the dynamics of its constituents (e.g., Ortiz-Amezcua et al., 2022a; Abril-Gago et al., 2023, Andújar-Maqueda et al., 2025). However, their implementation to the study of turbulent exchanges between the surface and the atmosphere is yet to be developed and standardized (e.g., Gibert et al., 2011).

Ventilation of underground CO₂ is a well-known phenomenon occurring in Mediterranean drylands, specifically in Cabo de Gata-Níjar Natural Park, in Southeastern Spain, where atmospheric pressure and turbulence have been found to be its main drivers (e.g., Sánchez-Cañete et al., 2013; López-Ballesteros et al., 2017, 2018; Moya et al., 2019), typically during the long dry summer period or sufficiently after precipitation pulses. In the current study, carried out in Cabo de Gata-Níjar Natural Park, turbulence within the boundary layer is investigated for the first time using a range-resolved Doppler lidar as a potential driver of soil ventilation processes, together with the assessment of the importance of variables such as atmospheric pressure and friction velocity, whose impact is already known.

2 Experimental setup

2.1 Station

This study focuses on the field campaign SCARCE (Synchronized Characterization of Aerosol, Radon and Carbon dioxide Exchanges in drylands), carried out between 18 July 2023 and 15 January 2024 at Balsa Blanca (36.94° N, 2.03° W, 205 m a.s.l., GHG-Europe site code: ES-Agu, referred to as Aguamarga in some studies) in Southeastern Spain. The station is located within the Cabo de Gata-Níjar Natural Park, a Mediterranean semiarid grassland of significant ecological importance at 6 km from the western coast of the Mediterranean Sea. The climate of the region is classified as dry subtropical semiarid, or BSh Köppen classification (El-Kenawy et al., 2022), coincident with the most arid region of Spain with a mean annual precipitation of around 220 mm and temperature of around 18 °C. Additionally, precipitation occurs in the form of rain pulses followed by long dry periods, with the most significant being the dry season from June to September/October (López-Ballesteros et al., 2016). These precipitation pulses trigger significant CO₂ emissions mainly due to soil respiration caused by rapid reactivation of microbial and root activity in the ecosystem (López-Ballesteros et al., 2016). Winds exhibit two dominant patterns, northeast and southwest, aligned with the coastline and the local orography, corresponding to the easterly and westerly wind regimes of the Alborán Sea (Abril-Gago et al., 2025a).

The landscape is characterized by patches of vegetation and bare soil. The dominant vegetation species is Macrochloa tenacissima, strongly adapted to severe water stress, with other less common species including Chamaerops humilis, Rhamnus lycioides, Asparagus horridus, Olea europaea var. sylvestris, Pistacia lentiscus, and Rubia peregrina (Rey et al., 2011). The growing season occurs from late autumn to early spring, when vegetation is more active (Serrano-Ortiz et al., 2014) due to the availability of water resources (López-Ballesteros et al., 2016). This period is followed by a non-growing period, characterized by strong water stress, higher temperatures, and a lack of precipitation.

The station is located over quaternary conglomerates and Neogene-Quaternary sediments cemented by lime (caliche), in a flat alluvial fan of gentle slopes (2 % to 6 %; Rey et al., 2017). Typical soils correspond to Calcaric Lithic Leptosol soils, shallow (10 cm maximum) and alkaline (pH above 8), with petrocalcic horizons and saturated in carbonates. The texture of the soil is sandy loam, with sand (61 %), silt (23 %), and clay (16 %), with a bulk density of 1.25 g cm⁻³. Groundcover consists of gravel and rock bare soil (49 %) and vegetation (51 %; Rey et al., 2011).

2.2 Atmospheric instrumentation

The station centers on a flux tower equipped for eddy covariance measurements. The 3.5 m tower was equipped with an enclosed-path LI-7200RS (LI-COR, Lincoln, NE, USA) gas analyzer that measures CO₂ and H₂O densities and atmospheric pressure, synchronized with a Gill HS-100-3 (Gill Instruments Ltd, Lymington, UK) three-axis sonic anemometer, together with a HMP45C thermohygrometer (Campbell scientific, Logan, UT, USA; hereafter CSI). These data were recorded by a CR3000 (CSI) data logger at 10 Hz. Micrometeorological data were later processed with EddyPro (LI-COR) version 7.0.8 over 30 min intervals. Finally, quality flags provided by EddyPro (Mauder and Foken, 2006) were considered and only data with 0 or 1 flag values were labeled as valid and consequently included in the analysis. Additionally, an ARG100 (CSI) pluviometer provided 30 min precipitation data registered on a CR1000X (CSI) with a resolution of 0.2 mm per tip.

For the present analysis, atmospheric CO₂ concentration (CO${}_{\mathrm{2}}^{\mathrm{air}}$ in ppm), CO₂ flux ($F^{{\mathrm{CO}}_{\mathrm{2}}}$, in µmol m⁻² s⁻¹), atmospheric pressure (P^air, in kPa), relative humidity (RH^air, in %) and temperature (T^air, in °C), horizontal wind speed (V^h, in m s⁻¹) and direction (in degrees), and friction velocity (u_∗, in m s⁻¹) information was taken from the flux tower upon processing and filtering whereas precipitation data (in mm) were directly obtained from the pluviometer.

A HALO Photonics StreamLine XR Doppler lidar was deployed at Balsa Blanca to characterize boundary-layer winds and turbulence during the SCARCE campaign, from 18 July 2023 to 15 January 2024. The system emitted laser radiation at 1565 nm with a pulse repetition rate of 10 kHz and detected the backscattered signal using a heterodyne detector (Pentikäinen et al., 2020). The Doppler shift introduced by atmospheric constituents was then retrieved from the signal. For this purpose, the standardized software package “HALO lidar toolbox” (Manninen, 2019) was used, producing multiple products including 3D wind fields and turbulence parameters. The system operated at a 1 s pulse rate, although processed products were generated with various temporal resolutions, with a default vertical resolution of 30 m. The system's effective full-overlap height was estimated at 90 m above ground level, defining the lowest reliable altitude (the altitude of 105 m was used as a reference, corresponding to the center of the 30 m interval starting at 90 m). The maximum detectable altitude varied depending on atmospheric turbidity, found within the atmospheric boundary layer (Ortiz-Amezcua et al., 2022b). Details of the deployment during SCARCE were provided by Abril-Gago et al. (2025a).

In the current study, two different products from the Doppler lidar were used: the wind shear (sh), representing the vertical gradient of horizontal wind velocity and direction, used as an indicator of mechanical atmospheric mixing (Ortiz-Amezcua et al., 2022a; Andújar-Maqueda et al., 2025); and the turbulent kinetic energy (TKE) dissipation rate (ϵ), quantifying the rate at which the TKE was dissipated into internal thermal energy (O'Connor et al., 2010). Atmospheric turbulence was assumed to be wind-shear-driven when sh exceeded 0.03 s⁻¹ (Manninen et al., 2018). Values of ϵ greater than 10⁻⁴ m² s⁻³ were indicative of intense turbulent activity (e.g., O'Connor et al., 2010), often associated with convection (Manninen et al., 2018), while lower values corresponded to atmospheric stability and stratification. Values of ϵ were filtered based on the signal-to-noise ratio (SNR). Vertical intervals with SNR lower than 0.006 (or equivalently −22.2 dB) were excluded from the analysis (Manninen et al., 2018). Finally, the original 3 min resolution of the sh and ϵ profiles was averaged into 30 min intervals, while the 30 m vertical resolution was aggregated into three broader layers, encompassing the ranges up to 195 m (sh^195 m, ϵ^195 m), up to 495 m (sh^495 m, ϵ^495 m) and up to 1005 m (sh^1005 m, ϵ^1005 m). Additionally, the particle transport velocities (v_t), directly related to upward and downward particle transport, calculated by Abril-Gago et al. (2025a) for the same location and period were incorporated into the statistical analysis.

A CHM15k-Nimbus ceilometer (Lufft, Fellbach, Germany), part of the ICENET network (Cazorla et al., 2017), was deployed alongside the Doppler lidar. This instrument operates by emitting laser pulses at a wavelength of 1064 nm with repetition rates between 5 and 7 kHz, and detecting the atmospheric backscatter using an avalanche photodiode. The device records the range-corrected signal (RCS) at a 1 min interval and 15 m vertical resolution. While full overlap of the laser beam with the receiver field of view occurs around 1500 m above ground level, a significant portion of the signal (approximately 90 %) is captured within the 555–885 m range. The RCS data were processed using the STRATfinder algorithm (Kotthaus et al., 2020) to retrieve the Mixing Layer Height (MLH), the convective mixing layer that develops during daytime due to surface heating, and the Atmospheric Boundary Layer Height (ABLH), which represents the top of the atmospheric boundary layer including the mixing layer and residual layers of previous days, marking the transition to the free troposphere (Kotthaus et al., 2023). To mitigate the impact of incomplete overlap at low altitudes, the RCS data were corrected using the overlap function provided by the manufacturer, providing reliable ABLH and MLH retrievals down to approximately 210–255 m (Kotthaus et al., 2020). Further methodological details related to the derivation of these layers during the present campaign are available in Abril-Gago et al. (2025a).

2.3 Vadose-zone instrumentation

The vadose zone of the station was monitored using a set of CO₂ sensors, Rn radioactivity detectors, and volumetric water content (θ), temperature, and pressure sensors installed in a similar configuration at three depths: 15, 50 and 150 cm. A similar setup was installed in a borehole approximately 18 m deep, which did not intersect the groundwater table. The borehole was sealed at the surface with a cap and equipped with a casing 15 cm in diameter. The first meter of the borehole casing closer to the surface was ungrooved, followed by a slotted section extending to the bottom, allowing interaction between the surrounding soil and the borehole interior. The additional borehole setup provides an alternative perspective on subsurface gas dynamics, as the instruments sample air circulating within the borehole, which is less affected by soil porosity than sensors installed directly within the soil profile. Three GMP252 (Vaisala Inc., Vantaa, Finland) CO₂ mole fraction probes were installed at the depths of 15, 50 and 150 cm (CO${}_{\mathrm{2}}^{\mathrm{15}\mathrm{cm}}$, CO${}_{\mathrm{2}}^{\mathrm{50}\mathrm{cm}}$ and CO${}_{\mathrm{2}}^{\mathrm{150}\mathrm{cm}}$, respectively, in ppm), each paired with a CS655 (Campbell Scientific) water content reflectometer, providing temperature (T^15 cm, T^50 cm and T^150 cm, in °C) and θ (θ^15 cm, θ^50 cm and θ^150 cm, in m³ m⁻³) measurements. Additionally, the 50 cm depth included an AlphaE (Bertin Technologies, Montigny-le-Bretonneux, France) Rn radioactivity detector (Rn^50 cm, in Bq m⁻³) which also provided pressure (P^50 cm, in kPa) and humidity (RH^50 cm, in %) measurements. Finally, the same configuration was applied in the borehole, with a GMP252, a HygroVUE5 (Campbell Scientific) temperature and relative humidity sensor and an AlphaE sensor were located within the first meter, providing CO₂ concentration (CO${}_{\mathrm{2}}^{\mathrm{bh}}$), temperature (T^bh), relative humidity (RH^bh), Rn activity (Rn^bh) and pressure (P^bh). Every sensor was connected to a CR1000X (CSI) data logger and data were averaged and recorded with a 30 min temporal resolution.

3 Methodology

3.1 Preprocessing

The raw 30 min series of every variable were detrended in order to remove the dominant diurnal cycle (24 h variations; Sánchez-Cañete et al., 2013). The resulting series reflected the underlying trends of the changes, which would highlight ventilation events not caused by daily patterns. A simple 24 h rolling average was applied to the series. Afterwards, the normality of the detrended series was refuted by Shapiro-Wilk and Anderson-Darling normality tests. Thus, non-parametric statistics had to be performed on the series. Additionally, an augmented Dickey-Fuller test confirmed the stationarity of the 30-min series, so no abrupt changes were present.

3.2 Event detection

The study's statistical analysis focused on specific selected periods (ventilation events) and not on the whole available data series. Therefore, a precise identification of the events was first performed. Each event would exhibit a sudden increase in Rn activity (from now on understood as Rn concentration) and CO₂ concentration (enriching ventilation), followed by an equally fast decrease in the concentrations (diluting ventilation). However, some of these phenomena might be caused by precipitation, which triggers significant CO₂ emissions from the soil (Vargas et al., 2018) to the ecosystem (López-Ballesteros et al., 2016) in the experimental site. Thus, periods during and up to one day after precipitation were not considered in this study, thereby reducing potential biases associated with soil moisture variability and its effect on Rn emanation (Zhuo et al., 2006). Additionally, based on previous studies (i.e., Sánchez-Cañete et al., 2013; Moya et al., 2019) the concentration increase had to coincide with an atmospheric pressure decrease, and equivalently a concentration decrease coincided with a pressure increase. Atmospheric pressure was used here exclusively as a physical criterion for event detection, while the subsequent analysis focused on the role of boundary layer dynamics in modulating the events through different meteorological variables. Furthermore, the soil humidity has also been identified as an important factor for CO₂ ventilation, so a θ^15 cm below 0.07 m³ m⁻³ was required, ensuring that the soil humidity was sufficiently low (Moya et al., 2022), allowing airflow through soil pores. Additionally, high-turbulence conditions were imposed, with a daily average u_∗ threshold of 0.2 m s⁻¹ at surface level, as established for different sites across the world by Moya et al. (2022). For this study, the high-turbulence condition was applied also to the Doppler lidar turbulence indicator, and the maximum ϵ value during the event should exceed 10⁻⁴ m² s⁻³ (O'Connor et al., 2010). The high-turbulence conditions were met easily, with daily average u_∗ exceeding 0.2 m s⁻¹ for every potential event, and maximum ϵ values exceeding 10⁻² m² s⁻³ for the selected events. Furthermore, the whole event had to last more than one day.

The beginning of each event was defined as the onset of a sustained (minimum 8 h) and rapid increase in tracer concentrations, exceeding 90 Bq m⁻³ h⁻¹ for Rn (Rn^50 cm or Rn^bh) or 80 ppm h⁻¹ for CO₂ (CO${}_{\mathrm{2}}^{\mathrm{150}\mathrm{cm}}$), the latter being the most variable CO₂ series and most sensitive to the concentration changes. The end of the event corresponded to the completion of a similarly sustained and rapid decrease in concentrations during the diluting phase. Sensitivity tests using minimum event durations between 6 and 12 h showed no significant impact on the statistical results. Therefore, given the use of 30 min averaged data and events lasting longer than one day, small uncertainties in the exact start and end times are not expected to significantly affect the correlation analyses. The inter-event periods were not included in the analysis, since they are dominated by slow, steady CO₂ production driven by biogeochemical processes (recharge, Sánchez-Cañete et al., 2013) rather than by mechanical processes linked to atmospheric phenomena. In the end, a total of 10 events were identified, encompassing 661 intervals of 30 min resolution, entailing 7.6 % of the available database.

3.3 Statistical analysis

Following standardization of the series and identification of enriching and diluting ventilation periods, bivariate Spearman correlation coefficients, ρ_S, were computed to assess the strength and direction of monotonic relationships between the study variables (Y variables) and the potential drivers (X variables). The study variables included CO₂ and Rn concentrations (CO${}_{\mathrm{2}}^{\mathrm{air}}$, CO${}_{\mathrm{2}}^{\mathrm{15}\mathrm{cm}}$, CO${}_{\mathrm{2}}^{\mathrm{50}\mathrm{cm}}$, CO${}_{\mathrm{2}}^{\mathrm{150}\mathrm{cm}}$, CO${}_{\mathrm{2}}^{\mathrm{bh}}$, Rn^50 cm and Rn^bh), as well as CO₂ flux ($F^{{\mathrm{CO}}_{\mathrm{2}}}$). CO${}_{\mathrm{2}}^{\mathrm{air}}$ and $F^{{\mathrm{CO}}_{\mathrm{2}}}$ were included among the study variables to initially assess their relationship with the potential drivers, but later the analysis focused on the soil CO₂ and Rn concentrations. The potential drivers were surface atmospheric pressure (P^air) and the surface friction velocity (u_∗) from the flux tower, the Doppler lidar's wind shear and turbulent kinetic energy dissipation rate at different height ranges (sh^195 m, sh^495 m, sh^1005 m, ϵ^195 m, ϵ^495 m and ϵ^1005 m) and the ceilometer's ABLH and MLH estimations. Rn concentration and CO₂ fluxes were also included among the X variables to confirm their known relationship with CO₂ concentrations, but they were not considered as potential drivers.

Once the importance of potential drivers was confirmed using ρ_S, different statistical procedures were conducted to investigate their role as real drivers of the soil CO₂ and Rn dynamics.

A partial correlation analysis was performed. The Spearman correlation coefficients between the potential drivers and the soil CO₂ and Rn concentrations were calculated, but separately controlling for each potential driver to be constant for the calculation. The way the ρ_S varies highlights the importance of each potential driver over soil CO₂ and Rn levels and among them. For example, let us take A as a potential driver which shows a strong Spearman correlation with CO₂ and Rn levels. If the ρ_S for A are invariant when the rest of drivers are forced to be constant, this would indicate that A is a real driver whose effect over the concentrations is independent of the other drivers. Similarly, if the ρ_S of other drivers were to reduce to zero when A is held constant, this would suggest that A could rule the other drivers. Equivalently, if the ρ_S of the other drivers were not to change when A is forced to be constant, this would indicate that A is not a real driver, but it is rather ruled by other real drivers. This partial correlation analysis was conducted taking as potential drivers: P^air, u_∗, ϵ^195 m, ϵ^495 m, ϵ^1005 m, sh^195 m, sh^495 m, sh^1005 m, ABLH and MLH.

To justify the use of bivariate Spearman correlations, a multicollinearity assessment was performed by calculating the Variance Inflation Factors (VIFs) of the potential drivers, yielding values below the common threshold of 5. This result suggests that multicollinearity was not significant and that the partial correlation results reliably reflect the independent effects of the drivers. Notably, significant collinearity was observed only between ϵ and sh at different altitudes, as expected, whereas their collinearity with the remaining variables remained below the threshold.

Finally, temporal shifts of each detrended potential driver series were performed and the ρ_S were recalculated. Since drivers experience changes before the variables they govern, if correlations with the concentration improve by shifting the driver series a given interval, this suggests which drivers actuate the underlying mechanisms. Thus, time series of a subset of potential drivers were shifted by a time lag between 180 min before and 180 min after the origin, in 30 min intervals. For these 12 shifted series the ρ_S with the soil CO₂ and Rn concentrations were recalculated. If the ρ_S increases for a series shifted by any interval after the origin (delayed in time), this suggests that the driver changes precede the concentration changes, and possibly its changes were affecting the concentrations. On the other hand, if the ρ_S increases for a series shifted by any interval before the origin (earlier in time), then the potential driver is experiencing changes later than the concentrations, suggesting that this is not a primary driver of the concentration's changes.

4 Results and discussion

4.1 Case study

Figure 1 displays an example of two consecutive ventilation events that occurred between 1 and 5 November 2023. The diluting ventilation period of a previous event is also visible on 30 October 2023.

Figure 1a shows the atmospheric and vadose-zone CO₂ concentration series at different depths. Shaded red areas indicate enriching ventilation periods, while blue areas indicate diluting ventilation periods. Note that the transition periods (no shading in Fig. 1) between enriching (red shaded) and diluting (blue shaded) ventilation periods, typically lasting around 12 h, were excluded from the statistical analysis, to ensure consistent tendencies among variables, since each concentration may respond to atmospheric changes at slightly different times. It can be observed that the deeper the sensor, the higher the measured concentration, especially during the maximum values recorded, with averages values ± standard deviation of 430 ± 4 ppm, 950 ± 180 ppm, 1800 ± 450 ppm, 1840 ± 760 ppm and 1100 ± 550 ppm for CO${}_{\mathrm{2}}^{\mathrm{air}}$, CO${}_{\mathrm{2}}^{\mathrm{15}\mathrm{cm}}$, CO${}_{\mathrm{2}}^{\mathrm{50}\mathrm{cm}}$, CO${}_{\mathrm{2}}^{\mathrm{150}\mathrm{cm}}$, CO${}_{\mathrm{2}}^{\mathrm{bh}}$, respectively.

Figure 1b exhibits how CO₂ and Rn concentrations behaved similarly throughout the process, confirming the usefulness of soil Rn as a proxy for soil CO₂ concentration under mechanically driven, non-biological influences. Unfortunately, the Rn sensor at 50 cm depth failed during this interval. During enriching ventilation periods, soil concentrations increased significantly, indicating that CO₂ and Rn from deeper layers were being transported upward, coinciding with a marked drop in P^air. In contrast, during diluting ventilation periods, concentrations decreased significantly, coinciding with a significant P^air increase, suggesting that atmospheric air (low in CO₂ and Rn) was infiltrating soil pores. On 5 November 2023, P^air began a steady increase, remaining high until 8 November. During this period, soil concentrations remained low and nearly constant, indicating that negligible recharge of CO₂ or Rn was possible during this period.

The borehole behaved differently, facilitating mixing from 0 to 18 m depth, although the CO${}_{\mathrm{2}}^{\mathrm{bh}}$ instruments were located closer to the surface. After the enriching phase, CO${}_{\mathrm{2}}^{\mathrm{bh}}$ approximated CO${}_{\mathrm{2}}^{\mathrm{50}\mathrm{cm}}$, while after the dilution it approached CO${}_{\mathrm{2}}^{\mathrm{15}\mathrm{cm}}$ or even CO${}_{\mathrm{2}}^{\mathrm{air}}$, indicating high permeability in the borehole that allowed both deep CO₂-rich air and atmospheric air to move freely. Rn^bh dropped to near-zero (atmospheric level) during diluting ventilation events, confirming that atmospheric air with low Rn concentrations was easily entering the borehole. In fact, the average Rn over this period was 750 ± 750 Bq m⁻³, with a variance coefficient (ratio between the standard deviation and the mean value) of 100 %, indicating that the Rn level varies between maxima values after enrichment and near zero concentrations after dilution.

Figure 1c shows u_∗ from the sonic anemometer and the Doppler-lidar-derived wind shear at different atmospheric height ranges. Their behavior resembled the concentration series, with both turbulence and wind shear increasing during periods of decreasing P^air, and decreasing when P^air rose. Both u_∗ and sh in the whole column remained low throughout the high-pressure period at the end of the series. Finally, Fig. 1d presents the Doppler lidar's ϵ. Two distinct behaviors were observed. ϵ^195 m appeared unrelated to the other variables, while ϵ^495 m and ϵ^1005 m were more clearly connected. Specially, ϵ^495 m and ϵ^1005 m exhibited significant increases (note the logarithmic scale on the vertical axis) during enriching ventilation periods, and significant decreases during diluting ventilation, similar to the trends observed in Fig. 1a, b and c. It is noteworthy that the timing of u_∗, sh and ϵ maxima and minima generally coincided with those of the CO₂ and Rn concentrations, reinforcing the fact that these variables were synchronized by specific atmospheric drivers.

Figure 1 Rolling-averaged temporal series of (a) CO₂ concentration in the atmosphere (air), at different vadose-zone depths (15, 50 and 150 cm) and in the borehole (bh); (b) Rn concentration in the borehole (left axis) and atmospheric pressure (right axis); (c) surface friction velocity (left axis) and wind shear (right axis) at the selected atmospheric intervals; and (d) ϵ at the selected atmospheric height ranges. Blue and red shading indicate diluting and enriching ventilation periods, respectively.

4.2 Statistical analysis

This section gathers the calculation of the ρ_S for the different studied variables. A total of 10 events, encompassing 661 intervals of 30 min resolution, were detected and analyzed. Initially, all the variables mentioned in Sect. 2 were considered. However, some were omitted from the final analysis shown here to avoid repetition for the sake of readability. Temperature (T^air, T^15 cm, T^50 cm, T^150 cm and T^bh), relative humidity (RH^air, RH^50 cm, and RH^bh) and soil water content (θ^15 cm, θ^50 cm and θ^150 cm) were correlated with the CO₂ and Rn concentrations. While some of these correlations were statistically significant, they are not presented here, as none of these state variables could be considered dynamic drivers of the ventilation events; rather, they seem to vary as a consequence of ventilation or because of common external forcing. Additionally, P^50 cm, and P^bh exhibited behavior similar to P^air and were also excluded from the correlation matrix figures.

Figure 2 displays the ρ_S matrix for the study variables (CO₂ and Rn concentrations and CO₂ flux) and the selected atmospheric variables, calculated with the selected 661 intervals of 30 min resolution. All ρ_S greater than 0.4 or less than −0.4 were statistically significant, with p-values below 0.05. Most ρ_S outside this range were also statistically significant. Figure 2a shows the ρ_S for the enriching ventilation periods, while Fig. 2b exhibits those for the diluting ventilation intervals. In both cases, soil CO₂ and Rn concentrations showed a strong positive correlation, with an average ρ_S of 0.92 ± 0.04, confirming that these variables were closely linked. Since Rn is an inert species, this strong correlation confirmed that Rn is an ideal proxy for identifying soil CO₂ ventilation events driven by atmospheric mechanisms (e.g., Kowalski et al., 2008; Sánchez-Cañete et al., 2011; Pla et al., 2023). In contrast, CO${}_{\mathrm{2}}^{\mathrm{air}}$ showed a lower but still significant correlation with soil CO₂ and Rn concentrations (0.67 ± 0.20), suggesting that, although related, their drivers did not actuate equally. On the other hand, $F^{{\mathrm{CO}}_{\mathrm{2}}}$ showed no correlation during enriching ventilation periods (−0.23 ± 0.11), but a significant negative correlation during diluting ventilation (−0.72 ± 0.04). This suggests that $F^{{\mathrm{CO}}_{\mathrm{2}}}$ was not influenced by enriching mechanisms, while dilution of subterranean CO₂ happens to the atmosphere, producing a significant emission of CO₂ to the air (e.g., Sánchez-Cañete et al., 2011, 2013, 2016). As preconditioned, P^air showed a strong negative correlation in both periods (−0.94 ± 0.04); however, the importance of P^air as a potential driver is already known (e.g., Sánchez-Cañete et al., 2013; Moya et al., 2019). Additionally, u_∗ correlated significantly and positively, especially during diluting ventilation (0.92 ± 0.03), supporting previous reports of its role in controlling CO₂ and Rn dynamics in soils at nearby sites (e.g., Sánchez-Cañete et al., 2016). As observed in Fig. 1, P^air increases correspond to CO₂ and Rn concentration and u_∗ decreases, and vice versa.

Figure 2 also shows the ρ_S between sh and ϵ, and the soil CO₂ and Rn concentrations. A consistent and significant positive correlation was observed for the wind shear and turbulence in the higher height ranges (0.85 ± 0.03, 0.87 ± 0.03, 0.84 ± 0.04, 0.83 ± 0.05, for sh^495 m, sh^1005 m, ϵ^495 m and ϵ^1005 m, respectively). However, low-altitude wind shear and turbulence (sh^195 m, ϵ^195 m) were only weakly correlated during enriching ventilation periods (0.64 ± 0.07 and 0.74 ± 0.06, respectively). sh^195 m showed a significant positive correlation during diluting ventilation (0.83 ± 0.04), while ϵ^195 m showed no correlation in that period (0.17 ± 0.15). Several considerations should be considered before interpreting these correlations. First, sh^195 m and ϵ^195 m were more affected by noise, due to the smaller vertical integration of the signal. Additionally, turbulence at these ranges was more sensitive to local heterogeneities and short-lived eddies, which may have different effects on the soil and atmosphere. In contrast, sh^495 m, sh^1005 m, ϵ^495 m and ϵ^1005 m represented averages over much deeper height ranges, typically covering the whole ABL in this location (Abril-Gago et al., 2025a), reducing noise and emphasizing large-scale eddies and synoptic-scale mechanisms. Therefore, boundary layer averaged wind shear and turbulence appeared to be better predictors of soil CO₂ and Rn dynamics, versus the Doppler-derived turbulence at lower altitudes. As observed in Fig. 1, increases in CO₂ and Rn concentrations were accompanied by increases in sh^195 m, sh^495 m, sh^1005 m, ϵ^495 m and ϵ^1005 m. Although not shown here, a strong correlation was observed between u_∗, sh and ϵ, only for the higher-altitude Doppler lidar ranges.

Finally, the relevance of ceilometer-derived ABLH and MLH was also explored. Negative correlations were found with CO₂ and Rn concentrations, being the strongest for ABLH during diluting ventilation events (−0.77 ± 0.08). This indicates that a significant increase in ABLH coincided with periods of reduced CO₂ and Rn levels.

Figure 2 ρ_S between CO₂ and Rn concentrations and CO₂ flux, and the selected atmospheric variables during (a) enriching and (b) diluting ventilation periods. Matching variable pairs have been indicated as blank cells.

Wind direction was not included in the correlation due to its circular nature. Nevertheless, it was consistently observed in every event, both enriching and diluting, that the wind blew from E or NE, indicating a preferred direction for their occurrence. This direction aligned with one of the two prevailing wind directions observed at the station (Abril-Gago et al., 2025a). Enriching ventilation periods began in the evening, between 18 and 00 UTC, and lasted until midday of the following day. Diluting ventilation periods started in the evening and continued until the afternoon of the next day. This pattern contrasts with observations at Llano de los Juanes (36.93° N, 2.75° W, 1600 m), located about 65 km west of Balsa Blanca, where nighttime ventilation is absent (Sánchez-Cañete et al., 2011), possibly due to site-specific differences in the dominant ventilation driver: P^air at Balsa Blanca versus u_∗ at Llano de los Juanes. Ventilation periods also showed marked meteorological differences. Enriching periods were characterized by lower surface turbulence and ABLH, with an average ± standard deviation values of $u_{\ast }=\mathrm{0.31}$ ± 0.13 m s⁻¹ and ABLH = 490 ± 200 m. In contrast, diluting periods exhibited higher surface turbulence and ABLH, with $u_{\ast }=\mathrm{0.46}$ ± 0.17 m s⁻¹ and ABLH = 700 ± 470 m. Other meteorological variables did not show significant differences between these periods.

Enriching ventilation periods were also associated with reduced vertical particle transport in the atmosphere, with average $v_{\mathrm{t}}^{\mathrm{495}\mathrm{m}}=\mathrm{0.99}$ ± 0.87 cm s⁻¹, similar to days with low convective mixing (Abril-Gago et al., 2025a). Diluting ventilation periods corresponded to enhanced vertical particle transport, with $v_{\mathrm{t}}^{\mathrm{495}\mathrm{m}}=\mathrm{2.0}$ ± 1.9 cm s⁻¹, comparable to days characterized by strong turbulence, whether from convective or mechanical mixing. This value was higher than those reported by Abril-Gago et al. (2025a) in Tables 1 and 2 for this location (0.90 cm s⁻¹ for the whole day and 1.56 cm s⁻¹ for the daytime period), suggesting that vertical transport, and consequently convective or mechanical mixing, during the identified diluting ventilation periods was exceptionally intense, even during nighttime. These results further highlight the connection between soil CO₂ and Rn dynamics and atmospheric particle exchanges, indicating that both were likely driven by similar atmospheric processes.

4.3 Driver analysis

Several atmospheric mechanisms have been proposed as potential drivers of soil CO₂ and Rn concentrations. In this section, different statistical analyses were conducted to determine which of these variables function as real drivers, and which were merely correlated due to being influenced by common underlying processes. To that end, partial correlation coefficients were recomputed to evaluate how the ρ_S changed when a specific variable (i.e., potential driver) is held constant.

Figure 3 gathers the average ρ_S of the potential drivers with the soil CO₂ and Rn concentrations, when each potential driver was independently held constant. A group of potential drivers was easily discarded, namely ϵ^195 m, sh^195 m, ABLH and MLH, as holding them constant did not meaningfully affect the ρ_S of the remaining variables. This suggested that these were not real drivers of the soil CO₂ and Rn dynamics and they were consequently omitted from Fig. 3. The same applies to u_∗, but only during enriching ventilation periods. During dilution ventilation events, average ρ_S were slightly reduced (in absolute terms) when ϵ^495 m, ϵ^1005 m, sh^495 m and sh^1005 m were individually held constant, implying that while these variables may influence the dynamics, they were unlikely to be the main drivers.

Among all potential drivers, P^air emerged as the most significant driver for both enriching and diluting ventilation periods. When held constant, it caused a notable decrease in the ρ_S of the other variables (0.20 ± 0.14), indicating its fundamental role as a determinant of ventilation. Moreover, its ρ_S changed the least when other variables were controlled (−0.75 ± 0.08), suggesting a strong independence. u_∗ also played a significant role during diluting ventilation periods, while ϵ^495 m, ϵ^1005 m, sh^495 m and sh^1005 m appeared to be more relevant as an enriching ventilation mechanism, with ϵ^495 m and ϵ^1005 m showing the most consistent influence. The results were almost identical when calculating the average ρ_S separately for the correlations between the potential drivers and the soil CO₂ at different depths, and for those between the potential drivers and soil Rn concentrations at different depths.

Figure 3 Average ρ_S between the potential drivers (X axis) and the soil CO₂ and Rn concentrations at different depths when a specific potential driver is held constant (Y axis), during (a) enriching and (b) diluting ventilation periods. Matching variable pairs have been indicated as blank cells.

The assessment of multicollinearity among the potential drivers was conducted, and the VIFs were all below the threshold of 5 (ranging from 1.9 to 4.2), indicating that each variable exhibited an acceptable level of independence. This suggests that the relationships identified between the drivers and the soil concentrations were not artificially influenced by redundant collinearity, and that the previous bivariate partial correlation results were reliable. As expected, ϵ^495 m, ϵ^1005 m, sh^495 m and sh^1005 m displayed a high degree of collinearity among themselves. This behavior is anticipated, as both ϵ and sh are derived from the same Doppler lidar signal, and measurements extending to higher altitudes also encapsulate information from lower height ranges.

A temporal shift of the series P^air, u_∗, ϵ^495 m, ϵ^1005 m, sh^495 m and sh^1005 m was performed, and the ρ_S were recalculated to identify the primary driver of the soil CO₂ and Rn concentrations.

For the periods and potential drivers considered in this study, the shifted correlations revealed that P^air maximized its ρ_S (increasing slightly by 2.2 %) during enriching ventilation intervals when shifted 30 min earlier, indicating that changes in P^air occurred after changes in CO₂ and Rn concentrations. Conversely, during diluting ventilation events, the maximum ρ_S were observed with no temporal shift, suggesting that P_air changes occurred simultaneously with soil concentration changes. The opposite pattern was observed for u_∗, with ρ_S maximizing (increasing slightly by 2.1 %) when shifted 30 min earlier for diluting ventilation events and maximizing during enriching ventilation periods with the unshifted series.

ϵ^1005 m, sh^495 m and sh^1005 m reached their maximum ρ_S with the original series during diluting ventilation, while ϵ^495 m maximized ρ_S (increasing significantly by 3.9 %) when the series was shifted 30 min later, indicating that changes in ϵ^495 m occurred before dilution and suggesting its importance as primary driver. During enriching ventilation periods, ϵ^495 m, ϵ^1005 m, sh^495 m and sh^1005 m showed maximum ρ_S (increasing by 1.2 %, 0.6 %, 1.9 % and 1.6 %, respectively) when shifted 30, 30, 60 and 30 min later, respectively, suggesting their roles as primary drivers.

4.4 Phenomenological description

To summarize the previous findings, the following conceptual model of the ventilation phenomenon is proposed.

At the start of an enriching ventilation period (Fig. 4a), boundary layer ϵ and sh start to increase anomalously (in the afternoon, even after daytime convection has ceased). Shortly thereafter, soil CO₂ and Rn concentrations begin to increase, accompanied by a drop in surface pressure and a rise in surface friction velocity. This state persists until noon of the following day, coinciding with a shallow ABLH. The combination of a poorly developed boundary layer, and elevated turbulence and wind shear suggests a situation of limited convection but enhanced mechanical mixing within the boundary layer. Consistently, CO₂ fluxes remain relatively low during this phase (0.05 ± 0.11 µmol m⁻² s⁻¹), approximately half of those observed under non-ventilation conditions, supporting the interpretation that CO₂ is being temporarily stored within the soil.

Then, in the afternoon, boundary-layer ϵ and sh begin to decrease notably (Fig. 4b), followed by a similar decrease in surface friction velocity and a rise in surface pressure, coinciding with the dilution of soil CO₂ and Rn. This period extends into the next afternoon and coincides with a higher ABLH and enhanced vertical particle transport within the boundary layer. The increase in surface pressure suggests the pumping of atmospheric air poor in CO₂ and Rn into the soil pores, compressing air within the soil and displacing CO₂-rich soil air to the atmosphere producing a significant CO₂ flux (0.30 ± 0.13 µmol m⁻² s⁻¹) up to three times that observed under non-ventilation conditions.

Figure 4 Graphical depiction of the ventilation model proposed, considering the (a) enriching and (b) diluting ventilation phases. The color gradient from black to gold indicates the transition from high to low CO₂ and Rn concentrations. Blue represents atmospheric concentration.

5 Conclusions

The present study advances our understanding of the CO₂ cycle in semiarid ecosystems, improving the characterization of CO₂ and Rn dynamics within the upper soil layers at a Mediterranean shrubland with carbonate soils, clarifying how their subsurface storage and release are governed by atmospheric processes. This work presents the first assessment of the role of above-surface atmospheric turbulence, across multiple atmospheric ranges, via remote sensing techniques, in driving subsurface CO₂ and Rn ventilation. Over the SCARCE campaign, 10 ventilation events of subsurface CO₂ and Rn, driven exclusively by abiotic mechanisms, were identified.

A set of meteorological variables were identified as potential drivers of subsurface ventilation, namely atmospheric pressure (P^air) and friction velocity (u_∗) at surface level, along with turbulent kinetic energy dissipation rate (ϵ) and wind shear (sh) within the boundary layer (ϵ^495 m, ϵ^1005 m, sh^495 m and sh^1005 m). The Spearman correlation coefficients (ρ_S) revealed significantly strong correlations between these variables and soil CO₂ and Rn concentrations at different depths, being negative for P^air and positive for the others, with slight differences between the absolute values of their ρ_S. Lower-height ϵ and sh (ϵ^195 m, sh^195 m), and also the atmospheric boundary layer height (ABLH) and mixing layer height (MLH) were discarded due to their weaker ρ_S.

Partial Spearman correlations confirmed that P^air was the most influential and independent factor during both enriching and diluting ventilation periods at the location. u_∗ was also significantly relevant during diluting ventilation, and ϵ^495 m, ϵ^1005 m, sh^495 m and sh^1005 m were notably influential during enriching ventilation periods, particularly ϵ^1005 m. This confirmed that these atmospheric parameters are actual drivers of air movement within the soil, thereby altering soil gas concentrations.

The temporal shifting analysis further supported these relationships. During diluting ventilation, most drivers maximized their ρ_S without temporal offset, except for ϵ^495 m, whose variations preceded those of the soil concentrations, identifying it as a primary driver. In contrast, during enriching ventilation, ϵ^495 m, ϵ^1005 m, sh^495 m and sh^1005 m reached their maximum ρ_S when shifted forwards in time, suggesting that ϵ and sh function as primary drivers during these intervals, whereas changes in P^air were observed after those of concentration.

Turbulent fluxes show that the ecosystem is a net source of CO₂. During enriching ventilation, CO₂ emissions drop to about half the background non-ventilation fluxes, suggesting temporary underground storage of CO₂. In contrast, diluting ventilation emissions increase by three times the background non-ventilation fluxes, suggesting the release to the atmosphere of stored CO₂. These results emphasize the role of the soil ventilation as a transient mechanism, modulating CO₂ emissions in semiarid ecosystems.

Overall, these findings enhance our understanding of the physical processes controlling soil-atmosphere gas exchanges in semiarid environments. In addition, this study reveals the potential of Doppler lidar observations, typically used for atmospheric characterization, for investigating complex soil–atmosphere interactions. Accounting for the impact of soil ventilation processes can therefore improve the representation of these exchanges in carbon cycle models, particularly in ecosystems where underground storage and release of CO₂ occur under specific atmospheric conditions, reducing uncertainties in climate change projections.