Gas diffusivity-based characterization of aggregated soils linking to methane migration in shallow subsurface

Methane transport in soil is primarily affected by soil physical conditions such as soil texture and soil structure, soil moisture, soil-gas diffusivity, permeability, and soil temperature. Aggregated soils have distinct soil structure with two pore regions characteristics (i.e., interaggregate and intraaggregate regions) and therefore show bimodal behavior with respect to soil physical properties controlling gas migration. This study characterized an aggregated soil retrieved near a natural gas (NG) extraction site at Denver–Julesburg (D-J) basin in northeast Colorado (USA) with respect to soil-water characteristic (SWC), pore-size distribution, gas diffusivity and thermal conductivity. The investigated soil exhibited distinctive two-region characteristics, which were adequately parameterized with extended, existing, and newly developed bimodal functions. We carried out an analysis with integrated model parameters to obtain a graphical insight on the correlation of properties. In addition, CH 4 concentration profiles


INTRODUCTION
Natural gas leakage from buried transmission pipelines is the major source of greenhouse gases (GHGs) from the energy supply sector (EPA, United State). The energy supply sector (natural gas [NG] and oil extraction, conversion, storage, and transmission) constituted about 35% of the total anthropogenic emissions in 2010 (Edenhofer et al., 2014). Methane is the predominant component (>90%) of NG (U.S. EIA, 2019), and it is a high potent GHG with an average global warming potential 86 times that of CO 2 on a 20-yr basis and 25 times greater over a 100-yr time horizon (Jackson & Vengosh, 2013). In addition, CH 4 is a highly flammable (National Fire Protection Association hazard rating = 4; a severe flammability) and explosive gas, leading to NG pipeline incidents in recent years (PHMSA, 2020). Natural gas pipelines are generally buried within natural soil systems-for example, through agricultural fields where pipelines may potentially be subjected to damages due to aging infrastructure, excavation, and/or human error due to intensive agricultural operations. With increasing concerns over mitigating GHG emissions, understanding gas migration behavior within the soil environment in an effort to pinpoint leak locations and quantify emissions has become increasingly important.
Agricultural soils undergo frequent compaction (due to the use of agricultural implements) and tillage resulting in alterations to the soil structural arrangement, thus affecting the total porosity and pore structure. Well-structured agricultural or pasture soils are typically constitute aggregated pore regions with distinct interaggregate pores where pores associate between the aggregate and intraaggregate pores where pores associate inside the aggregate. This aggregated behavior in the agricultural and pasture soils results in a strong bimodal pore structure. Due to their unique pore structure in aggregated soils, the structure-dependent properties controlling gas migration (e.g., soil-water characteristics [SWCs], soil-gas diffusivity, and soil thermal conductivity) are contrastingly different from unimodal soils (Chamindu Deepagoda et al., 2019). Therefore, an adequate characterization of aggregated soils is an essential prerequisite to understand subsurface gas migration and thereby detect leaks early, as well as maintain safe underground infrastructure.
In the event of a failure in the pipeline carrying NG, predominantly composed of CH 4 , subsurface CH 4 flow is typically controlled by advective flow near the leak (due to the buoyancy effects of CH 4 , high pipeline pressure, and soil permeability) and becomes progressively diffusion controlled as the gas migrates far away from the gas leakage point. The transition from advective to diffusive flux domination is due to the available air-filled porosity and usually occurs within a few meters of the leak location (La Bolle et al., 1998). Soil air permeability (k a , μm 2 ) is the major parameter dominating the advective gas movement in soil and is markedly controlled by the soil physical characteristics, mainly soil texture and structure (Ball, 1981b;Unsal et al., 2005). Diffusion-controlled gas migration in soils is primarily described by the soil-gas diffusivity (D p /D o ; Buckingham, 1904), where the D p and D o are gas diffusing coefficients (m 2 s −1 ) in soil and free air, respectively. Notably, D p /D o is strongly linked to air-filled porosity (ε) and the tortuosity (τ) of the functional soil gaseous phase in the porous media. Since the volumetric water content and the air-filled porosity mutually share the soil total pore space, soil moisture causes a twofold effect on the diffusive flow of gases in soil (Moldrup et al., 2000). Moreover, soil moisture creates an additional tortuosity and discontinuity in the gaseous pore space in the soil (Thorbjørn et al., 2008), making D p /D o a strongly moisture-dependent parameter. The SWC can be expressed as a functional relationship between the soil matric suction (-ψ, kPa) and the corresponding volumetric moisture content (θ, cm 3 cm −3 ), which can provide a useful characterization of soil moisture status and, in turn, gas transport properties in soil.

Core Ideas
• Aggregated soils have bimodal properties that affect CH 4 migration in soil. • Gas diffusivity, water characteristic, and thermal conductivity are bimodal fingerprints. • Diffusive CH 4 migration is largely controlled by soil moisture than leakage rate.
Since CH 4 is essentially a temperature-sensitive gas, the soil thermal conductivity (λ,W m −1 Κ −1 ) becomes an important soil physical property affecting CH 4 migration in soil. Among the key temperature-dependent properties of CH 4 , density, viscosity, water solubility, and more importantly, the diffusion coefficient is of particular importance in relation to diffusive CH 4 migration in soils under different temperature conditions. As the heat transfer essentially takes place via grain-grain contacts in a dry porous system, λ is largely controlled by the solid content (or soil density) in dry soil. With increasing moisture, however, more bridges are progressively developed among otherwise unconnected soil grains, thereby facilitating heat transfer with increasing water content. Due to the presence of two distinct pore regions, aggregated media are presumed to exhibit bimodal behavior in the characteristic λ-θ relationship. Thus, understanding the behavior of soil thermal properties at varying moisture conditions is a prerequisite to describe its associated effects on gas transport properties. Only a handful of previous studies, however, have closely investigated the integrated link among soil thermal conductivity, soil moisture, and soil-gas diffusivity (Al-Nakshabandi & Kohnke, 1965;Bristow, 1998;Chamindu Deepagoda, Smits, Ramirez, & Moldrup, 2016).
Literature bears evidence for a wide array of experimental and numerical studies investigating the CH 4 leakage from the underground NG pipelines. Most of the studies involved one-dimensional column experiments (Hibi et al., 2009;Karl-Heinz et al., 2004) together with the numerical simulations based on either a simplified advection-dispersion (Fickian diffusion-based) model (ADM) or the dusty gas model (DGM). Okamoto and Gomi (2011) carried out a field-scale study to investigate CH 4 leakage from an underground pipeline placed in an engineered test bed resembling a typical road base. They found that the gas migration patterns are highly gas density dependent. In another field-scale study, Yan et al. (2015) suggested that the leak direction (i.e., up, down, left, or right) has a significant impact on the gas concentrations above the pipe, despite considerably little effect on the below pipe. Only a few studies specifically focused on the effect of subsurface soil conditions on CH 4 transport, most of which were performed using pretreated rather than natural soils. For example, Praagman and Rambags (2008) conducted an extensive study on the effects of soil conditions using a one-and two-dimensional benchtop experimental apparatus together with a numerical model, concluding that an increase in water saturation caused a decrease in spreading width and volumetric flux despite the increase in pipeline pressure, and adding soil layering increases the spreading width. They also found that increase in the leakage rate had no effect on the total spreading width. However, their numerical and experimental results did not have good agreement, since the actual soil parameters were not accounted for in their numerical study. Wakoh and Hirano (1991) used a simple analytical solution to the advection-diffusion equation, which predicted a spherical concentration distribution based on the leak rate, length of time leaking, porosity, diffusion coefficient, and distance from the leak. A limited number of laboratory-based experimental studies on heterogeneous porous media illustrated the differences in CH 4 transport behavior of different soil layers with varying saturation conditions (Chamindu Deepagoda, Smits, & Oldenburg, 2016;Okamoto et al., 2011;Mitton, 2018) in relation to advective and diffusive CH 4 migration. Despite many studies investigating the effect of soil conditions on subsurface gas transport, studies specifically focused on subsurface CH 4 transport characteristics in aggregated porous media are sparse in the literature.
In this study, we characterized aggregated soils sampled from the Denver-Julesburg basin (referred hereafter as D-J basin), in Colorado (USA) in relation to soil-water retention, soil-gas diffusivity, and thermal conductivity to simulate methane emissions from a leaky NG pipeline buried at 0.3m depth. The measured properties of SWC, D p /D o , and thermal conductivity were parameterized with existing and newly developed parametric functions, and thereby we investigated their integrated behavior linking to the subsurface methane migration originated from a leaky NG pipeline. Furthermore, we simulated CH 4 concentration profiles at dry (i.e., θ = θ r ) and different saturation (θ = 0.25 θ s, θ = 0.5 θ s , and θ = 0.75 θ s ) conditions of the soil at three different flow rates (6, 12, and 24 L min -1 ) using TOUGH2/EOS7CA numerical package for multiphase flow of multicomponent gas mixture (CH 4 , air and water vapor).

Study area
We selected a location (104.8772˚N, 40.2951˚E) near a NG extraction well at D-J basin in northeast Colorado, USA. The D-J basin subsurface is characterized as Poleogene sandstone and conglomerate, a layer now known as Denver formation. The Denver basin is bounded with an area of about 180,000 km 2 in eastern Colorado, southeastern Wyoming, and western Nebraska. More than ∼10 11 m 3 of NG have been produced from more than 52,000 wells across the basin. Depths of production across the basin ranges from less than 270 m (900 ft) to about 2,700 m (9,000 ft) (Higley & Cox, 2007).
There are several studies focusing on the atmospheric surface emission of CH 4 , propane, and other NGs from the Denver basin (Levi, 2012;Peischl et al., 2018;Pétron et al., 2014), but with limited attention on accurate soil characterization.
In the neighborhood of the sampling location, a near-surface CH 4 concentration of 2.578 mg L −1 was measured.
The soil samples were retrieved in close proximity to an extraction well and near the main distribution pipeline laid at ∼0.3-m depth (pipe diameter = 15-20 cm). The soils were transferred to the laboratory and air dried prior to measurements. The samples were prepared in triplicates for the determination of all the required properties. Table 1 shows the measured soil-physical properties of the investigated soil. Note the relatively high saturated hydraulic conductivity due to the low bulk density (1.01 g cm −3 ) and aggregated nature of the soil.

SWC and thermal conductivity
We used Tempe cell apparatus (Smits et al., 2013) installed with a porous cup (7.2-cm length, 1.0-cm diameter, 51-kPa air-entry value) to measure the SWC and thermal conductivity for the matric suction above −1,000 kPa (pF < 3). The cell was hydraulically connected to a 200-cm high hanging water reservoir for capillary pressure adjustments. Two moisture sensors were installed into the cell for the continuous measurement of soil moisture (ECH2O EC-5, Decagon Devices) and thermal conductivity (SH-1, Decagon Devices) at different drainage levels. A detailed explanation of the sensor calibration, experimental procedure, and schematic of the experimental setup are presented in Smits et al. (2013). The soil was first wet-packed into the Tempe cell to the intended bulk density. Next, the water level of the reservoir was set to the top surface of the soil packed and the ADC (dimensionless analog-to-digital converter) count at fully saturated condition (ADC sat ) was recorded. The water level of the reservoir was then systematically lowered up to ψ = −10 kPa and the ADC count was recorded at each corresponding suction level. The volumetric water content belonging to each matric potential was calculated using the two-point α mixing model developed by Sakaki et al. (2008) (Equation 3). The WP4-T Dewpoint Potentiometer (Decagon Devices) was used to measure the matric potential-moisture status at low matric potentials (ψ < -1,000 kPa or pF > 3). The WP4-T performs based on the chilled-mirror dew point technique and provides a rapid determination of substrate water potential in equilibrium with water vapor (Gee et al., 1992;Resurreccion et al., 2011). Moisture content was adjusted by adding water to an ovendried soil sample (∼10 g); the sample was then sealed for 3 d until equilibrium was reached. At equilibrium, the sample was packed in the sample container and placed in the cell for measurements (in triplicates). Relative humidity in the chamber was used to determine the sample water potential at each moisture stage.

Soil-gas diffusivity (D p /D o )
Samples were repacked to stainless steel sample cores (100 cm 3 ) at three layers while maintaining the packing energy constant as five tamps per each layer (where "tamp" refers to the number of taps with equal energy around the core during packing). Soil cores were saturated for 72 h prior to draining to achieve different moisture conditions. Saturation was achieved by systematically increasing the depth of soaking water by 1 cm in 24-h intervals. Different moisture conditions were achieved by stepwise evaporation of saturated samples (moisture reduction in each step is 5 g). At each evaporation step, samples were kept in the air sealed bags for uniform distribution of the moisture across the sample. Diffusivity was measured following the method outlined by Moldrup et al. (2000), using a one-chamber diffusion apparatus (Currie, 1960

Numerical simulations
Experimental determination of CH 4 migration in soil under a wide range of potential variables (e.g., soil physical conditions, leakage rages, moisture status) is challenging; it is a common practice to conduct numerical simulations for a range of possible scenarios while validating a selected cases with experimental results. Subsurface migration and fate of CH 4 originated from the leaky NG pipe was simulated using a numerical simulator, TOUGH2, together with the equationof-state module (EOS7CA) (Oldenburg, 2015;Pruess et al., 1999). The modeling program uses a cubic equation-of-state with a multiphase version of Darcy's law, solved by an integral finite difference method, to model subsurface flow and transport of gas phases containing five components (i.e., H 2 O, brine, noncondensable gas, gas tracer, and air) under isothermal or nonisothermal conditions. TOUGH2/EOS7CA uses thermodynamic properties in real gas mixtures calculated by following the Peng and Robinson (1976) equation-of-state model. Soil-gas diffusion for the simulation is modeled with the temperature-dependent Fickian molecular diffusion coefficient. The solubility of methane in the aqueous phase is modeled by Henry's law, with Henry's coefficients calculated by Cramer (1982)'s method. Since the Henry's law approach formulated in the model is applicable for low pressure (less than ∼1 MPa) scenarios, TOUGH2/EOS7CA is essentially valid for the shallow subsurface (Oldenburg, 2015) as considered in this study. Note here that homogeneous and isotropic conditions across the entire domain with respect to transport parameters and no (or small) uptake of methane due to oxidation were assumed for the conservative simulation. A two-dimensional Cartesian numerical domain (width = 5 m, height = 0.3 m) was discretized into 336 soil elements for the subsurface CH 4 simulation originated from leaky NG pipe (Figure 1). Typical burial depth (0.3 m) of the NG distribution pipes near the location was set as the height of the domain. The top boundary was set as an open boundary for both gas (Dirichlet-type) and heat transport whereas the left, right, and the bottom boundaries of the domain were set as no-gas flow condition (Neumann-type) and adiabatic for heat flow (Figure 1a). The bottom mid-element was set as a CH 4 point source (width = 1 cm, length = 1 cm) to imitate a 1-cm × 1-cm hole in the top surface of the NG pipe. Simulations were performed at near-dry (i.e., θ = θ r ) and different saturated (θ = 0.25θ s , θ = 0.5θ s , and θ = 0.75θ s ) conditions for three leakages (6, 12, and 24 L min -1 ); these leakage rates were selected based on the historical fingerprints for the NG leakage incidents as reported in previous studies (Lamb et al., 2015;von Fischer et al., 2017;Weller et al., 2018). The Vadose Zone Journal F I G U R E 1 Discretized two-dimensional computational domains for CH 4 migration simulations in (a) completely dry soil where w is the water phase saturation, P c is the capillary pressure, and Γ is the relative permeability, and Panels b, c, and d are partially saturated (25, 50, and 75%) soil conditions with a known soil moisture gradient. Location of the leakage point as a diffusive CH 4 source (width [1 cm] × length [1 cm] × depth [0.05 cm]) also shown gas permeability is determined using intrinsic permeability (see Table 1) and gas properties. To simulate gas flow in the partially saturated condition, a pre-simulation was performed before injecting the CH 4 flow to achieve the gravity-capillary equilibrium throughout the domain, which has resulted in a narrow moisture gradient across the domain as illustrated in Figure 1b.

SWC
For a unimodal medium, van Genuchten (1980) proposed a multiparametric function for parameterizing water characteristic as follows: where θ (cm 3 cm −3 ) is the volumetric water content, θ s (cm 3 cm −3 ) and θ r (cm 3 cm −3 ) are the saturated and residual water contents, respectively, ψ (cm) is the matric potential, α (cm −1 ) is the model scaling factors associated to the airentry potential, and n and m (m = 1 -1/n) are the model shape factors. Durner (1994) modified the van Genuchten (1980) function by algebraically superimposing two pore regions to parameterize continuous relationship between capillary suction (ψ) and the volumetric water content (θ) in two-pore structured porous medias (interaggregate region and intraaggregate region), which can be written as follows: where w is the weighing factor (model parameter to numerically distinguish the two-pore regions), and the subscripts 1 and 2 indicate the two regions of the pore size distribution. The measured ADC (dimensionless analog-to-digital converter) counts using EC H 2 O EC-5 soil moisture sensors were converted into volumetric water content (θ, cm 3 cm −3 ) by an empirical two-point α mixing model (Sakaki et al., 2008). The two-point α mixing model can the take the form of where α is the model coefficient, ADC dry and ADC sat are the ADC counts corresponding to the soil at dry and saturation conditions, respectively, and Φ (cm 3 cm −3 ) is the total porosity. The first derivative of the van Genuchten (1980) model provides the equivalent pore size distribution of the soil medium (Durner, 1994), which can be written as where r (cm) is the pore radius and pF = log|ψ, 10 −1 kPa| (Schofield, 1935). Here dθ/dψ can be calculated as

Soil-gas diffusivity and tortuosity
We modified a pervious bimodal function (Chamindu Deepagoda et al., 2012) to parameterize the measured bimodal D p /D o data. The modified diffusivity function can be written as: Region 2 and S is the water saturation (= θ/Φ, where Φ [cm 3 cm −3 ] is the soil total porosity), w is the weighing factor (fitted model parameter) to numerically distinguish the two regions, and n is the shape factor. The subscripts 1 and 2 represent Region 1 and Region 2, respectively. Pore network tortuosity (τ) can be derived from the measured gas diffusivity following the method developed by Ball (1981a). Derivation of the tortuosity model is as follows:

Thermal conductivity
Predictive models expressing the soil thermal conductivity (λ, W m −1 K −1 ) as a function of water saturation (S) are abundant in the literature (Campbell et al., 1994;Côté & Konrad, 2005;de Vries, 1963;Haigh, 2012). Among all, the Côté and Konrad (2005) model (named hereafter as C-K model) is simple and adaptable to distinguish different thermal regimes in the characteristic λ-S relationship. The original C-K model is presented below: where λ sat and λ dry are thermal conductivities at completely saturated and dry conditions of the soil, respectively, K e (referred as the Kersten number) is a function of saturation (S), k is a dimensionless fitting parameter. Chamindu Deepagoda, Smits, Ramirez, and Moldrup (2016) modified the Kersten number in the original C-K model to a multiparameter multiregion function, with each parameter characterizing the distinct regions observed in the characteristic λ-S curve. In this study, since the measured λ-S curve exhibited the tworegion behavior, we used modified C-K model developed by Chamindu Deepagoda, Smits, Ramirez, and Moldrup (2016), and the measured λ dry and λ sat (see Table 1) values were used as the model-constraining parameters in λ-S relations. The model for two-region pore distribution can be written as for Region 1 and λ ( for Region 2, where the w is the weighing factor. Note that Φ inter−agg = and Φ intra−agg = (1 − ) .

F I G U R E 2
Volumetric soil-water content (θ, cm 3 cm −3 ) as a function of soil metric potential (pF) (denoted by solid circles) together with parameterized Durner (1994) bimodal function, Equation 4 (denoted by solid lines). (b) Pore density as a function of pore radius (cm) derived from the parameterized Durner (1994) model (Equation 4). Note that the pore radius is plotted on a log scale The parameter λ I denotes the thermal conductivity corresponding to Φ inter-agg.

Simulation of methane transport in subsurface
We used measured soil physical properties in TOUGH2/EOS7CA, a multiphase transport simulator that can simulate density-dependent flow of multicomponent gas mixture (CH 4 , water vapor, and air) in partially saturated porous media. The model has been successfully applied in many gas transport studies (Chamindu Deepagoda, Smits, & Oldenburg, 2016;Cho et al., 2020;Oldenburg et al., 2004;Mitton, 2018). The equations involved in this simulation of CH 4 migration can be found in Shanujah et al. (2020).

SWC and pore-size distribution
The SWC for the investigated soil is shown in Figure 2a together with the parameterized Durner (1994) model (Equation 2). Notably, the investigated soil exhibited a bimodal behavior with distinct interaggregate and intraaggregate pore regions, which is common for aggregated media. This can be more evidently seen in the pore-size distribution (Figure 2b), which shows two distinct pore regions characterized by two peaks. It can also be seen that the soil has a small air-entry pressure (P b = −1.2 kPa) to initiate desaturation, indicating a rapid gravity drainage of the interaggregate pores within a narrow matric potential range. At -1.2 kPa, the moisture content drops from saturated moisture content (θ = 0.62 cm 3 cm −3 ) to a moisture content of 0.40 cm 3 cm −3 , which apparently represents the interaggregate porosity of the soil. There seems to be a secondary, albeit small, air-entry pressure at θ = 0. 39 cm 3 cm −3 , as the intraaggregate pores start to drain. The subsequent draining, however, seems to be steady as the intraaggregate pore space consists of a wide pore size distribution (see Figure 2b). Towards the dry end (pF > 3.5), moisture is predominantly held in microcapillary pores by adsorptive forces, which require higher suction for draining. Clearly, the presence of two-pore spaces as well as the distinct SWCs may have notable implications on gas and CH 4 transport in the soil. Until the air-entry pressure is exceeded, there is virtually no distinct gas phase in soil; therefore, the CH 4 migration will occur as dissolved CH 4 transport through the aqueous phase and/or as CH 4 bubbles. Upon exceeding the air-entry pressure, nearly 20% of the pore space drains immediately, allowing CH 4 to rapidly access the pore domain. At the neighborhood of air-entry pressure, however, random connection or disconnection of variably sized and shaped pores in the heterogeneous pore space may potentially create additional tortuosity to the gaseous phase, yielding a discrepancy in expected emissions. Further draining above the air-entry exposes the intraaggregate pore space for gas transport, resulting in a systemic increase in gas migration.

Soil-gas diffusivity and gas-phase tortuosity
Soil-gas diffusivity with air-filled porosity is illustrated in Figure 3, together with the descriptive newly developed gas diffusivity function (Equations 6 and 7) performance. Data show mean diffusivity values (in solid symbols) together with error bars (denoted by solid lines) showing one standard deviation around the mean (n = 3). Notably, D p /D o also showed strong bimodal characteristics with distinct interaggregate and intraaggregate pore regions. It is worth noting that D p /D o showed nonlinear behavior in both the regions in contrast with some literature that showed linear variation of diffusivity with air-filled porosity within the intraaggregate region (Resurreccion et al., 2007). The similar variation of diffusivity within the two regions suggests a strong analogy of the two-pore spaces with respect to gas diffusion, although their size distributions were observed to be markedly different. The parametric diffusivity model gave a very good fit (R 2 = .98) to the measured data. The variation of gaseous phase tortuosity with air-filled porosity is also illustrated in Figure 3b. The simulated result from Equation 8 is also shown (solid line). Noticeably, high tortuosity near the water saturation is due to the pronounced water-induced tortuosity. The measured data, however, showed a markedly lower tortuosity than the modelpredicted values. In fact, below and around the air-entry pressure, the tortuosity values are generally found to be largely scattered due to two types of uncertainties: the intrinsic uncertainty of soil, as a result of randomly connected and disconnected water-filled pore space yielding considerably different values even among replicate samples, and measurement uncertainty which is due to the practical difficulty in handling the samples in an identical manner at high moisture conditions. With further draining, however, the tortuosity reduced sharply towards a minimum of 1.7 when the interaggregate pores were completely drained. Interestingly, a slight increase in tortuosity was observed at the onset drainage of intraaggregate pores as the gas is attempting to overcome the secondary air entry pressure, followed by a decrease as the draining of intraaggregate pores continued.

4.3
Thermal conductivity Figure 4 shows the measured thermal conductivity (λ) against volumetric water content for the investigated soil together with the improved C-K model (Equations 10 and 11). The original C-K (2005) model is also shown for reference. Importantly, the effect of water on λ is not equally weighed across the entire saturation range; two distinct regions in λ were particularly evident with increasing degree of saturation, which also implies a strong bimodal characteristic. As can be seen in Figure 4, the improved C-K (2005) model accurately represented the two different regions, whereas the original C-K model mischaracterized the two region behavior in the λ−θ relationship. The randomly connected or disconnected water bridges in between soil particles indicate the large scatter in λ at the low moisture content region. Evidently, the thermal conductivity in a porous system increases with increasing moisture content as the soil particles are bridged by more thermally conductive water (0.580 W m −1 K −1 ) than less thermally conductive air (λ = 0.024 W m −1 K −1 ). The scatter reduces and λ increases as more and more water is present in the pore space until all interaggregate pores are filled with water. The demarcation between the two regions occurred at θ = 0.33 cm 3 cm −3 , a slightly lower value than we observed for D p /D o and SWC. It is also worth noting that in both regions, the rate of increase of λ at low moisture contents is higher as the water is present as bound water to the grain surfaces, which contributes preferentially to the increase in λ. Further increasing Vadose Zone Journal F I G U R E 4 Thermal conductivity (λ, W m −1 K −1 ) as a function of volumetric soil-water content (θ, cm 3 cm -3 ) (denoted by circles) together with modified and parameterized Côté and Konrad (2005) function (Equations 10 and 11) (denoted by solid line). The original Côté and Konrad (2005) model (Equation 9) is also shown (in dotted line) moisture will merely contribute to a steady increase in λ as air is replaced with water.

Inter-parameter relations (pF, D p /D o , S, λ)
By combining the above observations on gas and thermal properties with varying moisture, Figure 5 shows a twodimensional representation of four parameters: soil-gas diffusivity (y axis), thermal conductivity (in color contours) and matric suction (pF, in solid line) across saturation S (= θ/θ s , x axis). The figure clearly demonstrates the two-region behavior in the observed data, suggesting SWC, D p /D o , and λ are strong fingerprints of a bimodal soil. In particular, D p /D o , and λ are correlative transport parameters in porous media and are highly moisture-dependent, though they respond oppositely to varying moisture status of the porous system. Note also that the pF contours and λ color contours are not congruent, as these parameters are nonlinearly related. The figure thus gives a useful graphical insight to the observed parametric relations for the direct measurement of D p /D o and λ measurements.

4.5
Methane concentration profiles Figure 6a presents the steady-state CH 4 concentration profiles for three NG leakage rates at the dry condition. Note here that the TOUGH2/ESO7CA accounts for a temperature correction on density, viscosity, solubility, and gas diffusivity of CH 4 within an advection-diffusion modeling framework. Though TOUGH2/ESO7CA modeling framework is not provisioned to incorporate soil bimodality, the aggregated soil is adequately represented by measured total porosity, gas diffusivity and tortuosity, intrinsic permeability, and thermal conductivity within the simulated moisture regimes.
Notably, due to the low specific gravity of CH 4 gas, the concentration contours show the tendency to move upward by buoyancy. At close proximity to the point of leakage, both pressure-driven advective and concentration-driven diffusive flow dominates the CH 4 transport and the flow becomes more and more diffusion-controlled as the CH 4 moves away from the source (Okamoto & Gomi, 2011). Diffusivity-based lateral and downward movements are virtually the same in all directions for a given leakage rate as can be seen in Figure 6a. Notably, the increase in leakage rate generally increased the surface CH 4 concentration due to the pressure-induced advective flow, while diffusion-controlled lateral movement is not particularly affected by leakage rate as also reported by Praagman and Rambags (2008). Simulated steady-state CH 4 concentration profiles for the several saturated (25, 50, and 75%) condition are shown in Figure 6b, 6c, and 6d. The saturation contours, after simulation, are also shown in white color. Methane concentration contours are distinctly different from those in dry soil systems due to the marked moisture-induced tortuosity for CH 4 diffusion in all directions. This is also partly due to the changes in the thermal conductivity of the soil at dry and partially saturated conditions (see Figures 2 and 5), causing increased density and decreased diffusivity of methane, resulting in F I G U R E 6 Simulated steady-state CH 4 concentration profiles for three leakage rates (6, 12, and 24 L min -1 ) in (a) dry condition and (b) 25% (c) 50%, and (c) 75% of saturated conditions of the soil. Colors denote the mass fraction of gaseous CH 4 . White colored horizontal lines indicate the saturation contours constrained movement of CH 4 through the gaseous phase with the increase of saturation level. In fact, the diffusioncontrolled lateral movement of CH 4 is largely different across all leakage rates under different statured condition (Figure 6b, 6c, and 6d) as compared with buoyancy (advection)dependant vertical movement, suggesting a clear effect of soil moisture on the diffusive movement of the CH 4 in the subsurface, as also observed in many past studies (Chamindu Deepagoda, Smits, & Oldenburg, 2016;Moldrup et al., 2000). Noticeably, a slight anomaly in moisture contours in partially saturated conditions near the source due to the displacement of moisture by leaking methane. Evidently, subsurface methane concentration profiles depicted pronounced effects of moisture compared with leakage rate, suggesting that moisture has made noticeable control on the subsurface migration of CH 4 . This is primality due to the diffusion-dominated movement of methane in subsurface which is predominantly controlled by soil moisture rather than the leakage rate. We further note that in all cases, the surface soils constitute <5% of CH 4 (by volume), which generally sharply decrease within the atmospheric boundary layer (Chamindu Deepagoda et al., 2017).
It also worthy to be noted herein that underground methane sources are broadly of two types: biogenic (i.e., CH 4 emitted from CH 4 -producing microorganisms, e.g., in landfills, and agricultural systems; Nielson et al., 2017) and thermogenic (e.g., fossil fuel production; and leakages from NG). A leaky pipeline spanning across an agricultural site, as the case in the present study, could have CH 4 emitted from both thermogenic and biogenic sources, yielding an overprint of CH 4 . To demonstrate that the measured CH 4 is predominantly thermogenically mediated, it is a common practice to trace another gaseous component: for example, ethane (C 2 H 6 ), which has no biogenic origin regardless of the methenogenic pathway. Literature reported that the ethane-tomethane (C 2 H 6 /CH 4 ) ratio for thermogenic origin methane is >20% Rowe & Muehlenbachs, 1999;Schoell, 1980). The measured average concentrations of nearsurface CH 4 (2.578 mg L −1 ) and C 2 H 6 (99.77 mg L −1 ) concentrations in this study clearly demonstrated that the measured CH 4 has a predominant thermogenic origin.
It is also worth mentioning here that the effects of near-surface atmospheric dynamics (e.g., wind-induced nearsurface pressure fluctuations) were not accounted for in the above simulations, which should also be considered for more realistic simulations. Subsurface movement of CH 4 originated from the leaky pipe can be further affected by leakage direction, and depth of pipe laying. Furthermore, diffusivity calculations ignored the effect of microbial O 2 consumption, which should also be needed for more realistic analysis. Nevertheless, Schjønning et al. (1999) reported that the O 2 depletion due to microbial consumption can be negligible within the experimental period. Therefore, care needs to be taken when comparing the results of this study with field-based CH 4 concentration measurements.

CONCLUSIONS
The SWC, pore-size distribution, gas-diffusivity (D p /D o ), and thermal conductivity (λ) of an aggregated soil retrieved at D-J basin in northeast Colorado (USA) were examined to investigate subsurface CH 4 transport originated from a buried leaky NG pipeline. Measured D p /D o , SWC, and λ showed strong bimodal behavior due to the presence of inter-and intraaggregated pore regions and thus demonstrated their applicability to obtain soil structural fingerprint in aggregated porosity media. The measured D p /D o and λ were found to be strongly moisture dependent, which markedly affected the subsurface gas transport behavior. Measured SWC, D p /D o , pore-size distribution, tortuosity , and λ were adequately described with extended, existing and newly developed two-region parametric functions. An integrated analysis on SWC, D p /D o , and thermal conductivity was also conducted to obtain a graphical insight on correlation of properties. We further simulated subsurface concentration profiles of CH 4 emitted from a point source representing a buried pipeline leaking at three different flow rates (6, 12, and 24 L min -1 ) using the multiphase transport simulator TOUGH2-EOS7CA under dry and different saturation conditions (25, 50, and 75% of saturation levels). Results highlighted pronounced effects of soil moisture and, to a lesser degree, of gas leakage rate on subsurface CH 4 concentrations.
Although the results provided a valuable insight to understand the soil properties and leakage rate on subsurface gas transport characteristic, further experimental and numerical studies that account additional soil and atmospheric complexities (e.g., wind-induced atmospheric dynamics, CH 4 oxidation, atmospheric temperature, and porous media anisotropy), are needed to conduct more realistic simulations.

A C K N O W L E D G M E N T S
Financial support from Indian-Sri Lankan Inter Governmental Science and Technology Co-operation Programme (MTR/TRD/AGR /3/2/17) is gratefully acknowledged. Any opinion, findings, and conclusions or recommendations expressed herein are those of the authors and do not necessarily reflect the views of those providing technical input or financial support. Trade names are mentioned only for identification purpose and do not constitute endorsement from any party involved in this study.

C O N F L I C T O F I N T E R E S T
The authors declare no conflict of interest.