Oxygen Propagation Fronts in Porous Media Under Evaporative Conditions at the Soil/Atmosphere Interface: Lab‐Scale Experiments and Model‐Based Interpretation

The interchange of gas components and volatile compounds between terrestrial and atmospheric compartments is critical for biogeochemical cycles and has important environmental and climate implications. In this study, we focus on oxygen and we explore the coupling between oxygen mass transfer and evaporation at the soil/atmosphere interface. We performed well‐controlled single‐phase and two‐phase laboratory experiments to determine the spatial and temporal evolution of oxygen fronts and to elucidate the coupling between mass and heat transfer in porous media with different grain sizes and under different evaporative conditions (i.e., no evaporation, natural, and enhanced evaporation). We also developed a non‐isothermal multiphase and multicomponent model to quantitatively interpret the experimental outcomes. The experiments and modeling allowed us to characterize the effects of external forcing (i.e., temperature gradients, humidity conditions) and internal factors (e.g., grain size) on the transport and distribution of oxygen in the different setups. Depth‐resolved spatial profiles and breakthrough curves of oxygen in the two‐phase experiments with evaporation are notably different from the single‐phase experiments due to the progressive gas invasion. The two‐phase experiments reveal a stepwise propagation pattern of oxygen that migrates considerably faster and penetrates deeper in the porous media in contrast to the relatively slow diffusion‐dominated transport regime in the absence of evaporation. The outcomes also show deeper and faster oxygen propagation in finer‐textured porous media under similar evaporative conditions, indicating the importance of internal factors for the distribution of the fluid phases and for the migration behavior of gas components in two‐phase systems.

The dynamics of evaporation in porous media and the associated multiphase processes are determined by external factors such as temperature/wind speed variations in the atmosphere and internal mechanisms like soil hydraulic properties Or et al., 2013;Shokri-Kuehni et al., 2020;Smits et al., 2011). There is a large body of literature characterizing the role of such mechanisms in different stages of evaporation from experimental and modeling perspectives at different scales (Mosthaf et al., 2014;Shahraeeni et al., 2012;Shokri et al., 2009;Trautz et al., 2014;Vanderborght et al., 2017;Weishaupt & Helmig, 2021). However, the investigation of component transport in porous media during evaporation has received less attention. The majority of studies in this direction focused primarily on saline water evaporation from porous media and described the salt precipitation patterns in drying soil (e.g., Shokri-Kuehni et al., 2018). Yet, a considerable gap remains in our understanding of gas component transport during the complex fluid displacement processes occurring due to evaporation. The detailed description of gas component transport in such systems requires the incorporation of the external and internal mechanisms affecting the distribution of liquid and gaseous phases because the transport time scale of gas components within these phases significantly varies. Despite the importance of such dynamics, a detailed analysis addressing the interplay between the external/internal factors and evolution of gas components during evaporation is still missing.
The objective of this study is to explore the impact of external atmospheric forcing (i.e., temperature gradients, humidity conditions) and soil intrinsic properties on the transport of gas components like oxygen in the shallow subsurface. The main hypothesis underlying this work is that the dynamic forcing in the atmospheric compartment influences the transport behavior of gas components in the subsurface by changing the spatiotemporal distribution of fluid phases in the pore space. We also hypothesize that internal factors such as the average grain size of the porous media impact the penetration depth of gas components in porous media under evaporative conditions. To this end, we performed a series of evaporation experiments in which soil columns with different grain sizes (i.e., coarse sand, medium glass beads, and fine sand) containing initially anoxic pore water were exposed to the atmosphere under two different evaporative conditions (i.e., natural and enhanced evaporation). In parallel, using the same porous media, we also performed a set of single-phase experiments in which evaporation was avoided and oxygen transport occurred in fully saturated porous media. High-resolution spatial and temporal profiles of oxygen and temperature were measured with non-invasive optode sensors, allowing us to capture the oxygen migration behavior with and without the influence of evaporation. We also developed a non-isothermal, multiphase, and multicomponent continuum model to quantitatively describe the experimental outcomes by considering the coupling and nonlinearities associated with the phase displacement process and component transport within and across phases.

Single-Phase and Two-Phase Experiments
We performed a series of single-phase and two-phase experiments to evaluate the oxygen migration in three porous media packed in glass columns. The schematic overview of the experimental systems is provided in Figure 1. In the single-phase experiments, evaporation in porous media was prevented and the diffusive transport of oxygen in a fully saturated porous medium was evaluated (Figure 1a), whereas in the two-phase experiments, we investigated the migration of oxygen under different evaporative conditions (i.e., natural and enhanced evaporation). The natural evaporation experiments were performed at room temperature (Figure 1b), whereas in the enhanced evaporation experiment, a heat source was applied on the top of the soil surface to increase the evaporation rate ( Figure 1c). This allowed us to investigate the influence of external conditions (i.e., temperature and humidity variations) on the transport of water and oxygen components. The columns used in all experiments were 16.85-cm long and had an inner diameter of 1.75 cm. The porous media used in the single-phase and two-phase experiments include fine sand (0.125-0.25 mm), medium glass beads (0.4-0.6 mm), and coarse sand (0.9-1.4 mm).

Single-Phase Diffusion Experiments
In the single-phase experiments, we studied the diffusive transport of O 2 from the atmosphere to anoxic pore water in three different porous media under fully saturated conditions. The glass columns were equipped with oxygen-sensitive sensor strips (10 × 0.5 cm, SP-PSt3-NAU, PreSens GmbH, Germany) glued onto the inner walls of the setups (Figure 1). This allowed us to non-invasively measure the spatial and temporal distribution of oxygen during the experiments using polymer optical fiber cables from the outside of the glass columns (Haberer et al., 2011).
The columns were wet packed and tapped several times while filling to minimize the air entrapment and layering during the packing procedure (Haberer et al., 2012). In order to prevent the ingress of oxygen, silicon septa together with a layer of Teflon were inserted at both ends of the columns, which were then tightly closed by plastic caps. A feed solution was prepared in a glass bottle filled with MilliQ water and sufficiently flushed with nitrogen to achieve oxygen-depleted water. The packed columns were flushed with the prepared anoxic solution from the bottom of the column with a constant flow rate using a high-precision multichannel peristaltic pump (IPC-N24, ColeParmer, United States). The flushing procedure continued until a low and uniform O 2 background concentration (<0.5 mg/L) was obtained at all depth-resolved measuring points. We then opened the caps at the top of the columns and placed the columns in a polystyrene box to maintain a high humidity and thus minimize evaporation by providing several water-filled beakers in the box. The oxygen diffusion from the atmosphere initiated after opening the caps and the oxygen concentration was continuously measured at different depths below the porous medium/atmosphere interface. In addition, O 2 spatial profiles were taken at different time intervals throughout the entire duration of the experiments (i.e., 24 hr). The temperature and relative humidity in the box were continuously monitored using a probe (Elitech GSP-6, USA) to ensure negligible temperature effects and high relative humidity (RH > 99%).

Two-Phase Diffusion Experiments
We performed two types of evaporation experiments, under natural and enhanced evaporative conditions, with the three different porous media to study the propagation of oxygen fronts in a two-phase system. In addition to the oxygen sensor foils, capable of measuring oxygen concentrations under both unsaturated and saturated conditions (Haberer et al., 2011;Huber & Krause, 2006), the columns were equipped with four temperature-sensitive spots (TPSP5, PyroScience, Aachen, Germany) that allowed us to non-invasively measure also the spatiotemporal evolution of temperature in the enhanced evaporation experiments. The measurement principle of such optical temperature sensors is based on thermal quenching of photoluminescence (Borisov et al., 2013). Note that for these sensors, as well as for the oxygen measurements, the assumption is that the non-invasive measurement at the wall is representative of the entire cross section of the thin columns used in the experiments.
In the natural evaporation experiments, after packing the setups and preparing the anoxic porous media (O 2 concentration <0.5 mg/L) as described above, the columns were opened at the top and placed under room temperature and humidity condition. This initiated the evaporation of pore water from the porous media to the atmosphere and the opposite migration of oxygen in the columns. The spatial and temporal dynamics of oxygen were monitored over the entire duration of the experiments (72 hr) to capture the oxygen migration in presence of evolving liquid/gas phase distributions. The temperature and relative humidity probe were placed at the top of the columns to continuously record the external conditions during the evaporation experiments. In parallel, we performed control experiments using the same columns and porous media to measure the water loss, determined by placing the entire column on a high-precision balance (KERN PCB, Germany) logging the weight at 10-min intervals.
In the enhanced evaporation experiments, we increased the evaporation demand by placing a heat source (i.e., 100-W lamp) at the soil surface ( Figure 1c). This led to heat propagation in the columns, which thus affected the associated multiphase and multicomponent transport processes. The packing and preparation of the columns were analogous to the natural evaporation experiments described above; however, in addition to the oxygen concentration, we measured the spatiotemporal variations of temperature using four sensor spots glued to the inner walls of the columns. The temperature rise and the relative humidity drop due to the presence of the heat source at the soil surface were continuously monitored during the enhanced evaporation experiments with probes at the top of the columns. Also, in this case, we conducted control experiments to measure the water loss during enhanced evaporation using the high-precision balance.

Modeling Approach
Process-based models are useful tools to understand the interplay between the different nonlinear physical processes associated with the description of evaporation and component transport during soil drying and to interpret experimental observations. To this end, we developed a non-isothermal multiphase and multicomponent model to quantitatively assess the outcomes of the two-phase experiments. The model represents several mechanisms including heat transport, phase change, capillary flow, interphase mass transfer as well as component transport in two phases. In the conceptual model, we assumed two phases (i.e., liquid and gaseous phases) and two components (i.e., water and oxygen). The components can be present in both phases and transfer between them. As water is the dominant component in the liquid phase, the interphase mass transfer of this component (i.e., evaporation) leads to the change in the liquid-phase saturation, which is typically a non-isothermal process. The oxygen component also partitions between the gaseous and liquid phases in a dissolution process and is transported within these phases mainly via diffusion. The mechanisms incorporated into the conceptual model need to be mathematically described to simulate the evaporative drying in the porous media and to evaluate the propagation of water and oxygen components in the resulting two-phase system. A schematic representation of the conceptual model of the atmosphere/soil system, including the main governing processes, is provided in Figure 2. In saturated porous media, the model simplifies to a single-phase and single-component formulation that was used to simulate oxygen migration under fully saturated conditions. This model was used to evaluate O 2 transport in the single-phase experiments and facilitated the comparison with the complex oxygen propagation behavior in the two-phase systems under natural and forced evaporation.

Non-isothermal Compositional Two-Phase Flow
In order to quantitatively describe compositional two-phase flow in porous media, the mass balance equations for the liquid (wetting) and gaseous (nonwetting) phase need to be formulated. There are different approaches for the formulation of such mass balance equations depending on the choice of primary and secondary variables (Chen et al., 1999;Fetzer et al., 2017;Solovský et al., 2020). In this study, we used the saturation/pressure formulation, which considers the mass balance equation of one phase (either liquid or gas phase) and the total mass balance equation. The mass balance equation for the fluid phase α ∈ {l, g} in a porous medium can be written as: where ε is the porosity, s α is the saturation of the fluid phase α, ρ α is the density of the fluid phase α, u α is the Darcy flux of the fluid phase α, and R α is the source/sink term of the α phase. We can define u α using the multiphase formulation of Darcy's laws as follows: where K int is the intrinsic permeability of the porous medium, k r,α is the saturation-dependent relative permeability of phase α, g is the gravity vector, μ α is the α-phase dynamic viscosity, and p α is the pressure of phase α.
The total mass balance equation can be obtained by summing up the mass balance equations of the phases (Equation 1): In this formulation, the saturation of one phase (the gaseous phase in this study) and the pressure of the other phase (the liquid phase in this study) can be independently chosen as the primary variables in Equations 1 and 3, respectively. The gas phase pressure can be computed from the capillary pressure relation, which is defined as the difference between the pressure of the gaseous phase and liquid phase (p c = p g − p l ) and assumed to be saturation-dependent.
The system of equations is closed by a constraint denoting that the sum of the phase saturations is equal to one (∑ ∈{ } = 1 ) and thus the liquid-phase saturation can be determined as s l = 1 − s g . Additionally, the dependency of the capillary pressure and relative permeability on the phase saturations needs to be specified using the constitutive relations. In this study, we use the van Genuchten-Mualem model given by (Mualem, 1976;van Genuchten, 1980): where s e,l is the effective saturation of the liquid phase defined as s e,l = (s l − s r )/(1 − s r ) with s r being the residual saturation of the liquid phase, H c is the capillary pressure head calculated as H c = p c /ρ l g and α, n, l and m (=1 − 1/n) are empirical soil hydraulic properties.
In addition to the mass balance equations of the fluid phases, the mass balance equations for the components need to be considered. For the component k ∈ {w, o 2 } in phase α ∈ {l, g} we have: where , and are the mass fraction, diffusive flux and the sink/source term of the component k in the α phase, respectively. The diffusive flux of components can be described by Fick's law as: where is the effective diffusion coefficient of the component k in the α phase. The diffusion of components in each phase only occurs in the pore space containing the liquid and gaseous phases, and this is accounted for by an effective diffusion coefficient for each phase, parametrized as = with being the temperature-dependent free diffusion coefficient of the component k in the α phase given as (T/273.15) 2 (Campbell, 1985) and τ α being the tortuosity in the phase α given as τ α = ε −1/3 s α −7/3 (Millington, 1959).
To describe the evaporation experiments performed, we need to consider Equation 7: (a) for the water component in the gaseous phase (vapor) and (b) for the oxygen component in the liquid and gaseous phases. We additionally need to link such equations to account for the components interphase mass transfer. This can be incorporated in the system of equations by introducing a kinetic mass transfer expression in the sink/source terms as follows (e.g., Battistel et al., 2019;Solovský et al., 2020): where k eq,w and k eq,o2 are the mass transfer coefficients for the water and oxygen components, respectively, (= ) and 2 (= 2 ) are the density of water in the gaseous phase and density of oxygen in the liquid phase, respectively, ρ v,sat is the saturated vapor density in the gaseous phase and H is the dimensionless Henry's law coefficient. The saturated vapor density is a temperature-dependent function, which is given by Tetens' equation as (Monteith & Unsworth, 2014): where p 0 is 610.7 (Pa),  is the universal gas constant, T is the temperature in Kelvin, and M w is the molecular weight of water.
The Henry's law coefficient is also a function of temperature and specifies the solubility limit of oxygen at a given temperature. This can be expressed by Van't Hoff equation as: where  (Sander, 1999) as 1700 [K].
In addition to the mass balance equations for the phases and components, we need to incorporate the heat balance equation in the description of the non-isothermal multiphase and multicomponent model. If we assume a local thermal equilibrium between the porous medium, liquid, and gaseous phases, the heat balance equation can be formulated as (Whitaker, 1976): where λ eff is the effective thermal conductivity assuming local thermal equilibrium between the phases, ∆h vap is the enthalpy of vaporization, and (ρ b C p ) eff can be specified as: where C p,α is the heat capacity of the fluid phase α, C p,s is the heat capacity of the porous medium, and ρ s is the porous medium density. The effective thermal conductivity is a nonlinear function of the fluid phase saturations and can be computed according to the model proposed by (Johansen, 1977): where λ dry and λ sat are the effective thermal conductivities under dry and fully saturated conditions, respectively, and K e is a function that can be approximated as (Côté & Konrad, 2005): where κ is a dimensionless empirical fitting parameter and, for the simulation of the enhanced evaporation experiments (i.e., with heat source), we used the value of 25 and 15 for the fine and coarse sand, respectively, to reproduce the measured temperature data at four different locations along the porous media depth.
For the single-phase experiments, Equation 7 simplifies to a mass balance equation for a single component (i.e., O 2 ) in a single phase (i.e., liquid phase), which is the governing equation of oxygen diffusive transport in saturated porous media.

Model Implementation
The model description for the non-isothermal compositional two-phase flow in porous media includes mass conservation equations for the phases and components, constitutive relations, as well as heat balance equations. The expressions representing the interphase mass transfer processes (Equations 9 and 10) link these equations using the thermodynamic relations (Equations 11 and 12) and induce a strong interplay between the temporal and spatial distribution of temperature, component concentrations, and fluid phase saturations and pressures.
The incorporation of such feedback effects in the model description requires solving the entire set of equations governing the non-isothermal compositional two-phase flow (Equations 1-16) in a fully coupled manner. The model formulation was completed by prescribing relevant initial and boundary conditions and then solved numerically in COMSOL Multiphysics (version 5.6). The finite element scheme implemented in this simulator allows the spatial discretization of the 2-D computational domain considered in this work. The domain includes a rectangular geometry representing a cross section of the column setups used in the experiments. The 2-D representation of the setups in the simulations allowed us to account for the heat loss from the side walls and at the bottom of the columns, but the properties across the cross section were assumed to be uniform. The details about the geometry of the setup as well as the initial and boundary conditions are provided in Figure S1 and Table S1 in Supporting Information S1. The time discretization was performed by an implicit time-marching scheme (backward differentiation formula) implemented in the simulator. A direct solver (PARDISO) combined with a nonlinear Newton solver was used to solve the coupled set of equations in the considered domain.

Single-Phase Diffusion Experiments
In the single-phase experiments, the diffusion of oxygen from the atmosphere into different porous media (fine sand, medium glass beads, and coarse sand) was investigated under fully saturated conditions. The entire period of the experiments was 24 hr, during which the spatial and temporal profiles of oxygen were regularly measured. The spatial distribution of normalized oxygen measured in the columns at different time intervals (i.e., in 1-12-hr time intervals) is presented in Figure 3. The normalized concentration is defined as (C − C bg )/(C 0 − C bg ) with C bg being the background oxygen concentration and C 0 the oxygen concentration at the top boundary.
The spatial profiles show that the oxygen penetrates from the top boundary in the initially anoxic porous media. The oxygen front progressively advances in the media and reaches a depth of ∼3.5 cm after 24 hr. The observed spatial profiles show smooth diffusive patterns and a similar behavior in the different columns. This indicates that the propagation of the oxygen fronts was not significantly affected by the grain size in the single-phase experiments. This can be explained by the relatively similar porosity (0.40, 0.42, and 0.41 for the fine sand, medium glass beads, and coarse sand, respectively), which results in similar tortuosity values. The tortuosity controls the magnitude of the effective diffusion coefficient and, fitting the single-phase diffusive transport model to the experimental observations, resulted in similar values for the considered porous media (i.e., 1.54, 1.43, and 1.72 for the fine sand, medium glass beads, and coarse sand, respectively).

Two-Phase Diffusion Experiments
In the two-phase experiments, we performed two sets of evaporation experiments with the three porous media under different evaporative conditions. Figure 4 shows the measured external conditions and evaporation rates for the considered porous media during the natural and enhanced evaporation experiments, the latter induced by a heat source above the columns.
The average relative humidity and temperature in the natural evaporation experiments are around 45% and 21°C, respectively, which ensures a consistent and stable external condition during the experiments performed with the different grain sizes. The enhanced evaporation experiments were also performed under stable conditions, but due to the heat source at the top of the setups the temperature reached almost 30°C and the relative humidity dropped to ∼35%. In the two sets of evaporation experiments, the cumulative water loss shows a linear trend within the time frame of the experiments, indicating the first stage of evaporation in which the capillary-driven liquid flow continuously provides water to the porous medium surface and supplies the evaporation demand of the atmosphere (Mosthaf et al., 2014;Or et al., 2013). The evaporation rate in this stage is typically independent of the grain size and is controlled by the external conditions (e.g., ambient temperature and relative humidity) as well as the evaporating fluid properties (Shokri-Kuehni et al., 2020). The porous media are exposed to similar external conditions in each group of evaporation experiments, which explains the reasonably similar slope of the cumulative water loss profiles and therefore the similar average evaporation rate of 3.36 and 5.05 mm/day determined for the three porous media in the natural and enhanced evaporation experiments, respectively. The higher evaporation rates (i.e., steeper slope of cumulative water loss) in the enhanced evaporation experiments is induced by the different external forcing (i.e., lower relative humidity and higher temperature), which increases the evaporation demand.
Despite the similar cumulative water loss in the evaporation experiments, the experimental outcomes show different patterns of oxygen spatial and temporal distribution in the porous media. Figure 5 shows the results of the natural evaporation experiments.
As it can be observed in the spatial profiles, at the beginning of the experiments (t = 0 hr), the pore water was anoxic along the depth (C norm = 0), and the oxygen diffusion from the top boundary initiated as the setup is exposed to the atmosphere. The O 2 diffusive front propagated in the porous media for the first few hours of the experimental period (e.g., until t = 3 hr for the fine sand), but the transport regime notably changed at the next time step in which an abrupt increase in the oxygen concentration occurred (e.g., at t = 6 hr for the fine sand), leading to a fast penetration of high oxygen concentrations in the porous media. Such behavior was followed by a diffusive trend (at t = 9 hr for the fine sand), and afterward by another fast propagation step (at t = 12 hr for the fine sand). Similar patterns were also observed for the fine sand at later times and deeper locations in the column (e.g., from 24 to 54 hr), and were consistently observed in the other porous media. In other words, the O 2 transport regime alternated between the slow diffusion-dominated stage and the fast migration step at different periods of the two-phase experiments. This represents a markedly different O 2 transport behavior in comparison to the diffusive transport observed in the single-phase experiments, suggesting the contribution of different mechanisms. Such mechanisms are mainly associated with the dynamic phase distribution and the progressive drying front in the evaporation experiments. In fact, in these experiments, the evaporation of pore water led to the formation of a drying front and created a partially saturated zone that gradually advanced along the column depth. The appearance of the gaseous phase in this zone significantly increased oxygen diffusion (the diffusivity of oxygen is 4 orders of magnitude larger in the gas phase), enhanced the oxygen mass transfer between the gaseous and liquid phases and thus allowed an increased oxygen penetration along the depth. Below the drying front (i.e., in the single-phase region), however, diffusion in the liquid phase was the dominant O 2 transport mechanism, thus resulting in the observed diffusive profiles until the drying front migrated further downwards and facilitated the interphase oxygen mass transfer (i.e., abrupt increase in the oxygen concentration). This can explain the alternating behavior between the slow and fast oxygen transport regimes observed in the experiments.
These key processes in the two-phase experiments thus collectively induced the oxygen front propagation that was considerably faster and penetrated deeper in the different porous media compared to the relatively slow diffusion-dominated transport regime in the single-phase experiments.

Figure 5.
Measured spatial profiles (a-c) and temporal evolution (d-f) of oxygen in the three porous media during the natural evaporation experiments in which the columns were exposed to ambient conditions. Another interesting observation in the O 2 spatial profiles is related to the impact of the grain size on the transport behavior. Oxygen transport in the considered porous media shows a similar trend in the two-phase experiments, but the oxygen penetration depth is larger for the fine sand and the medium glass beads compared to the coarse sand. This is attributed to the dynamics of the drying front depth (i.e., length of the partially saturated zone) in the fine and coarse sand, which influences the distribution of the gaseous/liquid phases above the partially saturated zone and thus changes the oxygen transport. During soil water evaporation, the drying front depth and the fluid phase distribution in the pores are controlled by a balance between capillary forces, viscous losses and gravity, as well as liquid phase properties (Lehmann et al., 2008;Or et al., 2013). Therefore, the differences in the drying front depth and the oxygen transport behavior observed for the fine and coarse sand can be explained by the different capillary forces in the fine-textured and coarse-textured porous media as reflected in their corresponding van Genuchten parameters. This observation shows that unlike the single-phase experiments, the spatial distribution of oxygen in the two-phase system is affected by the porous media grain size.
In addition to the spatial profiles, the temporal evolution of oxygen was monitored at three different locations (1, 3, and 5 cm below the porous media/atmosphere interface). Figures 5d-5f show the oxygen breakthrough curves at these locations for the natural evaporation experiments. The outcomes demonstrate three different stages in the evolution of the oxygen concentration in the porous media. The first stage includes a gradual increase in the oxygen concentration (i.e., slow diffusive regime) under saturated conditions, but in the second stage, this behavior shifts to an abrupt increase (i.e., fast migration regime) under partially saturated conditions, and eventually the oxygen concentration reaches the maximum value (C norm = 1) and remains constant in the third stage. Such pattern in the oxygen temporal profiles can be observed at the three locations in which the oxygen breakthrough curves were measured. However, with increasing depth, the duration of the first stage (i.e., slow diffusive regime) is longer and therefore the abrupt increase in the oxygen concentration occurs at later times. This is due to the drying front that first invades the shallower depth of the soil (i.e., 1 cm) and thus initiates the oxygen transport regime under partially saturated condition at earlier times in comparison to the deeper locations (i.e., 3 and 5 cm). Due to such dynamics of the drying front, the duration of the slow diffusive transport regime in the temporal evolution of oxygen in the porous media (i.e., first stage) is prolonged with increasing depth. Figures 5d-5f also confirm the impact of the grain size on the O 2 temporal evolution during the natural evaporation experiments. The measured oxygen concentrations at 1, 3, and 5 cm below the atmosphere/soil interface show that with increasing grain sizes, the slow O 2 diffusive transport (i.e., first stage) occurs over a longer period and therefore the abrupt increase in the oxygen concentration (i.e., the second stage) initiates at later times. This indicates that the drying front in the coarse sand arrives at later times at the respective locations (i.e., 1, 3, and 5 cm) in comparison to the fine sand. Figure 6 shows the observed spatiotemporal distribution of oxygen in the three porous media in the experiments performed under enhanced evaporation conditions. Overall, the spatial profiles show a similar trend in the oxygen transport behavior observed in the natural evaporation experiments with diffusive fronts that alternate with abrupt changes over the course of the experiments. However, in these experiments, the oxygen fronts advanced faster and deeper in the porous media in comparison to the natural evaporation experiments. In fact, the high demand for evaporation imposed by the heat source at the top of the porous media increased the evaporation rate, which accelerated the advancement of the drying front and, thus, enhanced the invasion of the gaseous phase and the oxygen migration. Deeper and faster oxygen penetration is also reflected in the oxygen breakthrough curves measured at different locations (Figures 6d-6f). It can be observed that the abrupt increase in the oxygen concentration (i.e., second stage in the temporal evolution of oxygen), associated with the arrival of the drying front at the respective locations, takes place at earlier times as evaporation is enhanced by the heat source.
In these experiments, we also measured the temporal evolution of temperature at 1, 3, 5, and 7 cm below the porous media surface (Figures 6g-6i). Overall, such profiles show a similar trend in the temporal evolution of temperature in the three porous media. The trends observed at all locations show an increase of temperature due to heat conduction in the first 5 hr of the experiments after which the temperature stabilizes and remains almost constant, approaching a pseudo steady-state condition. The temperature closer to the porous media surface (i.e., sensor measurements at 1-cm depth) is higher than the values recorded by the sensors at deeper locations, after reaching the steady-state condition. This indicates that, despite the insulation used in the experiments, some heat loss occurred in the columns. Such distribution of temperature affected the propagation of the oxygen fronts in the porous media. This influence is reflected in the O 2 spatial profiles particularly in the unsaturated zone where oxygen dissolution occurs. It is interesting to note that the measured oxygen concentrations in this zone are higher at the deeper locations (i.e., cooler bottom of the columns) reflecting the fact that the vertical thermal gradients imply spatially variable solubility (i.e., the solubility of oxygen increases at lower temperature).

Figure 6.
Measured spatial profiles (a-c), temporal evolution of oxygen (d-f), and temporal evolution of temperature (g-i) in the three porous media during the enhanced evaporation experiments in which the columns are exposed to a heat source.

Model-Based Interpretation
We used the non-isothermal multiphase and multicomponent model to interpret the outcomes of the two-phase experiments performed with the fine and coarse sand, to characterize the external/internal factors controlling the oxygen transport, and to understand the differences observed in the oxygen transport behavior during the single-phase and two-phase experiments. The model inputs including the initial and boundary conditions used for the simulation of the natural and enhanced evaporation experiments are presented in Table 1 and Table S1 in Supporting Information S1.
The simulation outcomes for the natural evaporation experiments are presented in Figure 7 and show that the model is capable of reproducing the spatial distribution of oxygen and the dynamics of the water component in the fine and coarse sand by considering the multiple feedback processes involved in the migration of the components. The interpretation of the oxygen transport behavior in the two-phase experiments requires information about the dynamics of the water component, the distribution of the liquid and gaseous phases, and the thickness of the unsaturated zone. The simulated water loss agrees with the experimental observations and shows a linear trend with a similar slope for the fine and coarse sand, indicating that the water loss and thus the evaporation rates within the period of the experiments are independent of the soil hydraulic properties (i.e., internal factors). The   reorganization of the remaining liquid phase in the porous media, however, is influenced by the hydraulic properties of the porous media as shown by the simulated liquid-phase saturation (Figures 7c and 7d). Such results indicate differences in the spatial distribution of the liquid/gaseous phase saturation for the fine and coarse sand at similar time steps. The major distinctions are reflected in the presence of the gaseous phase at deeper locations in the fine sand (i.e., larger unsaturated zone) and larger proportion of the gaseous phase in the upper part of the coarse sand column. Such differences in the distribution of the fluid phases change the dynamics of the key processes associated with the oxygen transport in the two-phase system and can explain the influence of mean grain size on the oxygen transport observed in the evaporation experiments (Figures 7e and 7f). In fact, the deeper penetration of the gaseous phase in the fine sand allows faster oxygen diffusion and facilitates the oxygen mass transfer between the gaseous and liquid phases in a larger portion of the fine porous medium, leading to deeper penetration of oxygen in the fine sand. Figure 8 summarizes the model results for the enhanced evaporation experiments. Due to the presence of the heat source above the porous media, heat transport, and its coupling with the multiphase and multicomponent transport processes was included in these simulations. The good agreement between the simulated and measured temperature values indicates that the model is able to capture the temporal evolution of temperature recorded by the non-invasive sensors located at different depths (Figures 8a and 8b). The advancement of the heat front in the porous media induces complex non-isothermal processes, which affect the dynamics of phase change, the displacement of fluids, and the component transport during the experiments. As shown in Figure 8, the model can reproduce the cumulative water loss in both porous media (panels c and d) and, also for these enhanced evaporation experiments, the simulated spatial distribution of the liquid phase shows a deeper invasion of the gaseous phase in the fine sand (panels e and f). Finally, the spatial distribution of oxygen and its evolution at different times is also correctly captured by the simulations (panels g and h). These results indicate the capability of the model to reproduce the dynamics of the water and oxygen components by incorporating the feedback effects induced by the heat transport in the porous media. In particular, the model is able to capture the impact of the temperature gradients on the oxygen concentrations in the unsaturated zone (i.e., higher concentration at deeper locations). The results indicate a further invasion of the drying fronts and thus of the gaseous phase into the porous media relative to the natural evaporation experiments discussed above. This quantitatively explains the faster advancement of the O 2 front in the enhanced evaporation experiments and highlights the impact of the external forcing (i.e., temperature gradients) on oxygen transport in the porous media.
The model also helps explaining the differences observed in the oxygen transport behavior during the single-phase and two-phase experiments. Figure 9 provides a comparison between the measured and simulated oxygen breakthrough curves in the single-phase and two-phase experiments under natural evaporative conditions for the case of the fine sand.
The duration of the single-phase experiments was 24 hr, but the simulations were performed for 60 hr to allow a comparison of the simulation outcomes with the evaporation experiments. Overall, the model is capable of capturing the trend and concentration values measured at the three locations both in the single-phase and two-phase experiments. The simulation outcomes highlight a significant distinction between the single-phase and two-phase experiments with respect to the propagation of the oxygen front. In the single-phase experiment, the temporal evolution of oxygen in the fine sand shows a slow diffusive front with decreasing velocity under fully saturated conditions, whereas in the two-phase experiment, the O 2 front initiates with a slow diffusive front followed by a significantly fast front as evaporation occurs and unsaturated conditions prevail in the fine sand. Such distinction is particularly reflected in the simulated oxygen concentrations at 1 and 3 cm in the single-phase experiment, showing the penetration of almost 65% and 20% of the maximum oxygen concentration (C norm = 1) after 60 hr, respectively. However, in the two-phase experiments, the simulation results at the same locations demonstrate the arrival of the maximum oxygen concentration (C norm = 1) at much earlier times (i.e., 10 hr at 1 cm and 35 hr at 3 cm) when the drying front invades these locations. These outcomes indicate the significant influence of evaporation as well as the interphase mass transfer processes on gas component transport in porous media, and highlight the importance of incorporating soil/atmosphere interactions in the description of gas component migration in the subsurface. Such interactions can be particularly important for the dynamics of biogeochemical processes that rely on the delivery of oxygen from the atmospheric compartment to the subsurface. For instance, oxygen-dependent transformation kinetics of organic compounds and the spatial location of biodegradation hotspots could be impacted by the fast supply of high oxygen concentrations in porous media during evaporation.

Conclusions
The coupled mass transfer of oxygen and water in the shallow subsurface is essential for a wide range of biotic and abiotic processes. The mechanisms controlling the migration of these important components in the subsurface are influenced by interchange processes and physical and chemical gradients at the soil/atmosphere interface. We carried out a detailed experimental and model-based investigation of oxygen transport in porous media under fully saturated conditions (i.e., without evaporation) and under natural and enhanced evaporation (i.e., without and with heat source). The outcomes of our investigation show a markedly different oxygen transport behavior in single-phase (i.e., under water-saturated conditions) and two-phase systems (i.e., under partially saturated condition) and corroborate our initial hypothesis that the atmospheric forcing greatly impacts the transport of gas components by controlling the fluid phase distribution in porous media. In particular, the results quantitively show that the water evaporation and the invasion of the gaseous phase in the two-phase experiments lead to larger O 2 penetration and increase the velocity of the oxygen propagation front in the initially anoxic porous media in comparison to the single-phase experiments. The observed propagation patterns were found to occur as sequential contributions of multiple transport mechanisms, including the O 2 diffusive migration in the gaseous and liquid phases as well as the interphase mass transfer processes. The dynamics of such transport mechanisms are strongly linked to the evolution of phases, and this leads to a highly complex and nonlinear behavior in the system that is controlled by external mechanisms (i.e., temperature and humidity conditions) and internal factors (i.e., grain size). The experiments performed using porous media with different mean grain sizes demonstrate that this internal factor does not greatly affect oxygen transport under fully saturated conditions, but becomes important during evaporation. We developed and applied a non-isothermal multiphase and multicomponent model to quantitatively understand the coupling and feedback between the governing nonlinear processes and to characterize the contribution of individual mechanisms. The proposed model was tested with the series of column experiments performed and was capable of reproducing the experimental observations, including evaporative water losses and oxygen spatiotemporal profiles, as well as the distinct behavior observed in the fine and coarse sand, leading to different propagation of the oxygen fronts in the two porous media. Finally, the model also allowed capturing oxygen transport in presence of forced evaporation induced by a heat source, thus highlighting the coupling of heat propagation with the distribution of phases and the transport of components in the porous media.
In this study, we focused on oxygen transport, but the approach can be used to explore the displacement of other components, characterized by different physico-chemical properties, at the soil/atmosphere interface. In particular, investigating the exchange and fluxes of greenhouse gases (e.g., CH 4 , CO 2 , and N 2 O) and volatile organic contaminants in multiphase systems and under evaporative atmospheric forcing would be of timely scientific and practical relevance. Furthermore, the investigation of oxygen transport could be extended to chemical and/ or biological reactive systems to explore the interactions between the multiphysics dynamics highlighted in this study and the kinetics of reactive processes. Future investigation should also address the processes investigated in Figure 9. Comparison between measured (symbols) and simulated (lines) oxygen breakthrough curves in the fine sand during the single-phase experiment (a) and two-phase experiments under natural evaporation (b). The inset shows the progressive invasion of the gaseous phase during the two-phase experiment, which strongly affects the oxygen transport regime.
this work, considering both chemical heterogeneity, such as the presence of reactive minerals in the solid matrix (Battistel et al., 2021;Erfani et al., 2021;Salehikhoo & Li, 2015;Tartakovsky et al., 2008) and physical heterogeneity, with spatially variable properties such as porosity and permeability (Heidari & Li, 2014;Ye et al., 2015). Finally, we also envision future research focusing on the exchange of gas components across the soil/atmosphere interface in the presence of variable dynamic forcing in the atmosphere by applying different wind speed in the atmospheric compartment (Bahlmann et al., 2020) and considering pressure fluctuations.

Data Availability Statement
The experimental data sets for the single-phase and two-phase (i.e., natural and enhanced evaporation) experiments are available in the data repository: https://doi.org/10.11583/DTU.16988899.