From slots to tubes: The influence of dimensionality on fracture dissolution models

We briefly review the models of fracture dissolution process, discussing the experimental and numerical evidence showing that this phenomenon is inherently two-dimensional and hence cannot be accurately described by one-dimensional models. The physical reason for this incompatibility is that a dissolution front in a single rock fracture is potentially unstable to small variations in local permeability, leading to spontaneous formation of dissolution channels in the rock. This leads to a dramatic increase of fissure opening rates, which must be taken into account not only in the estimation of karstification times but also in the assessment of ground subsidence, dam collapse or toxic seepage risks.     

Effects of chemical reactions on density-dependent fluid flow: on the numerical formulation and the development of instabilities

A three-dimensional, reactive numerical flow model is developed that couples chemical reactions with density-dependent mass transport and fluid flow. The model includes equilibrium reactions for the aqueous species, kinetic reactions between the solid and aqueous phases, and full coupling of porosity and permeability changes that result from precipitation and dissolution reactions in porous media. A one-step, global implicit approach is used to solve the coupled flow, transport and reaction equations with a fully implicit upstream-weighted control volume discretization. The Newton–Raphson method is applied to the discretized non-linear equations and a block ILU-preconditioned CGSTAB method is used to solve the resulting Jacobian matrix equations. This approach permits the solution of the complete set of governing equations for both concentration and pressure simultaneously affected by chemical and physical processes. A series of chemical transport simulations are conducted to investigate coupled processes of reactive chemical transport and density-dependent flow and their subsequent impact on the development of preferential flow paths in porous media. The coupled effects of the processes driving flow and the chemical reactions occurring during solute transport is studied using a carbonate system in fully saturated porous media. Results demonstrate that instability development is sensitive to the initial perturbation caused by density differences between the solute plume and the ambient groundwater. If the initial perturbation is large, then it acts as a “trigger” in the flow system that causes instabilities to develop in a planar reaction front. When permeability changes occur due to dissolution reactions occurring in the porous media, a reactive feedback loop is created by calcite dissolution and the mixed convective transport of the system. Although the feedback loop does not have a significant impact on plume shape, complex concentration distributions develop as a result of the instabilities generated in the flow system.     

Simulation of steady-state flow in three-dimensional fracture networks using the boundary-element method

An efficient method for simulating steady-state flow in three-dimensional fracture networks is formulated with the use of the boundary-element method. The host rock is considered to be impervious, and the fractures can be of any orientation and areal extent. The fractures are treated as surfaces where fluid movement is essentially two-dimensional. Fracture intersections are regarded as one-dimensional fluid conduits. Hence, the three-dimensional geometric characteristics of the fracture geometry is retained in solutions of coupled sets of one- and two-dimentional equations. Use of the boundary-element method to evaluate the fluid responses in the fractures precludes the need to internally discretize the areal extent of the fractures.     

Matrix Diffusion Effects on Contaminant Migration from an Injection Well in Fractured Sandstone

Deep-well injection into fractured sandstone is an option for the disposal of contaminated mine dewatering discharge from an open pit uranium mine. As part of the assessment of potential contaminant migration from deep-well injection, the effect of matrix diffusion was evaluated. An analytical mathematical model was developed for the simulation of the radial movement of a contaminant front away from an injection point under steady flow conditions in a planar fracture with uniform properties. The model includes the effects of advection in the fracture, diffusion of contaminants from the fracture into the rock matrix, and equilibrium adsorption on the fracture surface as well as in the rock matrix. Effective diffusion coefficients obtained from laboratory experiments on 11 intact core samples varied from 3.4 × 10⁻⁸ to 3.2 × 10⁻⁷ cm²/s. Model simulations were made with diffusion coefficient values in this range and with single-fracture injection rates estimated from fracture frequencies in boreholes, and from bulk hydraulic conductivity values obtained from field tests. Because of matrix diffusion, the rate of outward movement of the front of the nonreactive contaminants from the injection well is much slower than the rate of water flow in the fractures. Simulations of the movement of contaminants that undergo adsorption indicate that even a small distribution coefficient for the rock matrix causes the contaminants to remain very close to the injection well during the one-year period. The results of the simplified model demonstrate that matrix diffusion is an important process that cannot be neglected in the assessment of a waste disposal scheme located in fractured porous rock. However, in order to make a definitive assessment of the capability of matrix diffusion and associated matrix adsorption to significantly limit the extent of contaminant migration around injection wells, it would be necessary to conduct field tests such as a preliminary or experimental injection.     

Analysis of a vertical dipole tracer test in highly fractured rock

The results of a vertical dipole tracer experiment performed in highly fractured rocks of the Clare Valley, South Australia, are presented. The injection and withdrawal piezometers were both screened over 3 m and were separated by 6 m (midpoint to midpoint). Due to the long screen length, several fracture sets were intersected, some of which do not connect the two piezometers. Dissolved helium and bromide were injected into the dipole flow field for 75 minutes, followed by an additional 510 minutes of flushing. The breakthrough of helium was retarded relative to bromide, as was expected due to the greater aqueous diffusion coefficient of helium. Also, only -25% of the total mass injected of both tracers was recovered. Modeling of the tracer transport was accomplished using an analytical one-dimensional flow and transport model for flow through a fracture with diffusion into the matrix. The assumptions made include: streamlines connecting the injection and withdrawal point can be modeled as a dipole of equal strength, flow along each streamline is one dimensional, and there is a constant Peclet number for each streamline. In contrast to many other field tracer studies performed in fractured rock, the actual travel length between piezometers was not known. Modeling was accomplished by fitting the characteristics of the tracer breakthrough curves (BTCs), such as arrival times of the peak concentration and the center of mass. The important steps were to determine the fracture aperture (240 microm) based on the parameters that influence the rate of matrix diffusion (this controls the arrival time of the peak concentration); estimating the travel distance (11 m) by fitting the time of arrival of the centers of mass of the tracers; and estimating fracture dispersivity (0.5 m) by fitting the times that the inflection points occurred on the front and back limbs of the BTCs. This method works even though there was dilution in the withdrawal well, the amount of which can be estimated by determining the value that the modeled concentrations need to be reduced to fit the data (approximately 50%). The use of two tracers with different diffusion coefficients was not necessary, but it provides important checks in the modeling process because the apparent retardation between the two tracers is evidence of matrix diffusion and the BTCs of both tracers need to be accurately modeled by the best fit parameters.     

Compaction of a Rock Fracture Moderated by Competing Roles of Stress Corrosion and Pressure Solution

Unusually rapid closure of stressed fractures, observed in the initial stages of loading and at low temperatures, is examined using models for subcritical crack growth and pressure solution. The model for stress corrosion examines tensile stress concentrations induced at the Hertzian contact of propping fracture asperities, and mediates fracture growth according to a kinetic rate law. Conversely, pressure solution is described by the rate-limiting process of dissolution, resulting from the elevated stresses realized at the propping asperity contact. Both models are capable of following the observed compaction of fractures in novaculite. However, closure rates predicted for stress corrosion cracking are orders of magnitudes faster than those predicted for pressure dissolution. For consistent kinetic parameters, predictions from stress corrosion better replicate experimental observations, especially in the short-term and at low temperature when mechanical effects are anticipated to dominate. Rates and magnitudes of both stress corrosion and pressure solution are dependent on stresses exerted over propping asperities. Rates of closure due to stress corrosion cracking are shown to be always higher than for pressure solution, except where stress corrosion ceases as contact areas grow, and local stresses drop below an activation threshold. A simple rate law is apparent for the progress of fracture closure, defined in terms of a constant and an exponent applied to the test duration. For current experimental observations, this rate law is shown to replicate early progress data, and shows promise to define the evolution of transport properties of fractures over extended durations.     

Catchment‐scale conceptual modelling of water and solute transport in the dual flow system of the karst critical zone

Hydrological and hydrochemical processes in the critical zone of karst environments are controlled by the fracture‐conduit network. Modelling hydrological and hydrochemical dynamics in such heterogeneous hydrogeological settings remains a research challenge. In this study, water and solute transport in the dual flow system of the karst critical zone were investigated in a 73.5‐km² catchment in southwest China. We developed a dual reservoir conceptual run‐off model combined with an autoregressive and moving average model with algorithms to assess dissolution rates in the “fast flow” and “slow flow” systems. This model was applied to 3 catchments with typical karst critical zone architectures, to show how flow exchange between fracture and conduit networks changes in relation to catchment storage dynamics. The flux of bidirectional water and solute exchange between the fissure and conduit system increases from the headwaters to the outfall due to the large area of the developed conduits and low hydraulic gradient in the lower catchment. Rainfall amounts have a significant influence on partitioning the relative proportions of flow and solutes derived from different sources reaching the underground outlet. The effect of rainfall on catchment function is modulated by the structure of the karst critical zone (e.g., epikarst and sinkholes). Thin epikarst and well‐developed sinkholes in the headwaters divert more surface water (younger water) into the underground channel network, leading to a higher fraction of rainfall recharge into the fast flow system and total outflow. Also, the contribution of carbonate weathering to mass export is also higher in the headwaters due to the infiltration of younger water with low solute concentrations through sinkholes.     

Dissolution in a variably confined carbonate platform: effects of allogenic runoff,hydraulic damming of groundwater inputs,and surface–groundwater exchange at the basin scale

In variably confined carbonate platforms, impermeable confining units collect rainfall over large areas and deliver runoff to rivers or conduits in unconfined portions of platforms. Runoff can increase river stage or conduit heads in unconfined portions of platforms faster than local infiltration of rainfall can increase groundwater heads, causing hydraulic gradients between rivers, conduits and the aquifer to reverse. Gradient reversals cause flood waters to flow from rivers and conduits into the aquifer where they can dissolve limestone. Previous work on impacts of gradient reversals on dissolution has primarily emphasized individual caves and little research has been conducted at basin scales. To address this gap in knowledge, we used legacy data to assess how a gradient of aquifer confinement across the Suwannee River Basin, north‐central Florida affected locations, magnitudes and processes of dissolution during 2005–2007, a period with extreme ranges of discharge. During intense rain events, runoff from the confining unit increased river stage above groundwater heads in unconfined portions of the platform, hydraulically damming inputs of groundwater along a 200 km reach of river. Hydraulic damming allowed allogenic runoff with SI_CAL < ?4 to fill the entire river channel and flow into the aquifer via reversing springs. Storage of runoff in the aquifer decreased peak river discharges downstream and contributed to dissolution within the aquifer. Temporary storage of allogenic runoff in karst aquifers represents hyporheic exchange at a scale that is larger than found in streams flowing over non‐karst aquifers because conduits in karst aquifers extend the area available for exchange beyond river beds deep into aquifers. Post‐depositional porosity in variably confined carbonate platforms should thus be enhanced along rivers that originate on confining units. This distribution should be considered in models of porosity distribution used to manage water and hydrocarbon resources in carbonate rocks. Copyright © 2013 John Wiley & Sons, Ltd.     

Challenges of using electrical resistivity method to locate karst conduits—A field case in the Inner Bluegrass Region, Kentucky

Conduits serve as major pathways for groundwater flow in karst aquifers. Locating them from the surface, however, is one of the most challenging tasks in karst research. Geophysical methods are often deployed to help locate voids by mapping variations of physical properties of the subsurface. Conduits can cause significant contrasts of some physical properties that can be detected; other subsurface features such as water-bearing fractures often yield similar contrasts, which are difficult to distinguish from the effects of the conduits. This study used electrical resistivity method to search for an unmapped karst conduit that recharges Royal Spring in the Inner Bluegrass karst region, Kentucky, USA. Three types of resistivity techniques (surface 2D survey, quasi-3D survey, and time-lapse survey) were used to map and characterize resistivity anomalies. Some of the major anomalies were selected as drilling targets to verify the existence of the conduits. Drilling near an anomaly identified by an electrical resistivity profile resulted in successful penetration of a major water-filled conduit. The drilling results also suggest that, in this study area, low resistivity anomalies in general are associated with water-bearing features. However, differences in the anomaly signals between the water-filled conduit and other water-bearing features such as water-filled fracture zones were undistinguishable. The electrical resistivity method is useful in conduit detection by providing potential drilling targets. Knowledge of geology and hydrogeology about the site and professional judgment also played important roles in locating the major conduit.     

Rates of dissolution of aluminosilicates in seawater

Dissolution of eight clay minerals, four zeolites, and quartz in seawater has been monitored for81/2 years. For most of the minerals, dissolution can be described as a first-order reaction in which dissolved silica approaches from undersaturation steady concentration values with time. Characteristic reaction rate constants (k₁) are of the order of 10^?7 sec^?1. One of the zeolites, clinoptilolite, shows a different dissolution behavior: SiO₂ concentration in solution reaches a high value within one year, followed by a decline to a lower value, suggestive of precipitation of another silicate phase (possibly sepiolite).A mathematical solution is given for a kinetic equation combining the parabolic-rate and first-order rate processes. It is shown that in a wide range of silicate dissolution reactions taking place over long periods of time, the presence of the parabolic-rate dissolution processes cannot be detected, thereby making its inclusion in the kinetic equations unnecessary. The experimental rates of dissolution are comparable to the SiO₂^? dissolution rates in oceanic sediments near the sediment/water interface. But deeper in the sediment, the calculated dissolution rates are significantly lower than the near-interface and experimental values.     

An analysis of crystal dissolution fronts in flows through porous media part 2: incompatible boundary conditions

A model for transport of solutes in a porous medium participating in a dissolution–precipitation reaction, in general not in equilibrium, is studied. Ignoring diffusion–dispersion the initial value problem for piecewise constant initial states is studied, which e.g. for ionic species include a change of the ionic composition of the solution. The mathematical solution, nearly explicitly found by the method of characteristics up to the (numerical) solution of an integral equation for the position of the dissolution front, exhibits a generalized expanding plateau-structure determined by the dissolution front and the water flow (or salinity) front.     

CO2-induced dissolution of low permeability carbonates. Part II: Numerical modeling of experiments

We used the 3D continuum-scale reactive transport models to simulate eight core flood experiments for two different carbonate rocks. In these experiments the core samples were reacted with brines equilibrated with pCO₂ = 3, 2, 1, 0.5 MPa (Smith et al., 2013 [27]). The carbonate rocks were from specific Marly dolostone and Vuggy limestone flow units at the IEAGHG Weyburn-Midale CO₂ Monitoring and Storage Project in south-eastern Saskatchewan, Canada. Initial model porosity, permeability, mineral, and surface area distributions were constructed from micro tomography and microscopy characterization data. We constrained model reaction kinetics and porosity–permeability equations with the experimental data. The experimental data included time-dependent solution chemistry and differential pressure measured across the core, and the initial and final pore space and mineral distribution. Calibration of the model with the experimental data allowed investigation of effects of carbonate reactivity, flow velocity, effective permeability, and time on the development and consequences of stable and unstable dissolution fronts.The continuum scale model captured the evolution of distinct dissolution fronts that developed as a consequence of carbonate mineral dissolution and pore scale transport properties. The results show that initial heterogeneity and porosity contrast control the development of the dissolution fronts in these highly reactive systems. This finding is consistent with linear stability analysis and the known positive feedback between mineral dissolution and fluid flow in carbonate formations. Differences in the carbonate kinetic drivers resulting from the range of pCO₂ used in the experiments and the different proportions of more reactive calcite and less reactive dolomite contributed to the development of new pore space, but not to the type of dissolution fronts observed for the two different rock types. The development of the dissolution front was much more dependent on the physical heterogeneity of the carbonate rock. The observed stable dissolution fronts with small but visible dissolution fingers were a consequence of the clustering of a small percentage of larger pores in an otherwise homogeneous Marly dolostone. The observed wormholes in the heterogeneous Vuggy limestone initiated and developed in areas of greater porosity and permeability contrast, following pre-existing preferential flow paths.Model calibration of core flood experiments is one way to specifically constrain parameter input used for specific sites for larger scale simulations. Calibration of the governing rate equations and constants for Vuggy limestones showed that dissolution rate constants reasonably agree with published values. However the calcite dissolution rate constants fitted to the Marly dolostone experiments are much lower than those suggested by literature. The differences in fitted calcite rate constants between the two rock types reflect uncertainty associated with measured reactive surface area and appropriately scaling heterogeneous distribution of less abundant reactive minerals. Calibration of the power-law based porosity–permeability equations was sensitive to the overall heterogeneity of the cores. Stable dissolution fronts of the more homogeneous Marly dolostone could be fit with the exponent n = 3 consistent with the traditional Kozeny–Carman equation developed for porous sandstones. More impermeable and heterogeneous cores required larger n values (n = 6–8).     

Dynamics of lava flow fronts, Arenal Volcano, Costa Rica

Arenal Volcano has effused basaltic andesite lava flows nearly continuously since September, 1968. The two different kinds of material in flows, lava and lava debris, have different rheologic properties and dynamic behavior. Flow morphology depends on the relationship between the amount and distribution of the lava and the debris, and to a lesser extent the ground morphology.Two main units characterize the flows: the channel zone and the frontal zone. The channel zone consists of two different units, the levées and the channel proper. A velocity profile in the channel shows a maximum value at the plug where the rate of shear is zero, and a velocity gradient increasing outward until, at the levées, the velocity becomes zero. Cooling produces a marked temperature gradient in the flow, leading to the formation of debris by brittle fracture when a critical value of shear rate to viscosity is reached. When the lava supply ceases, much of this debris and part of the lava is left behind after the flow nucleus drains out, forming a collapsed channel.Processes at the frontal zone include levée formation, debris formation, the change in shape of the front, and the choice of the flow path. These processes are controlled primarily by the rheological properties of the lava.Frontal zone dynamics can be understood by fixing the flow front as the point of reference. The lava flows through the channel into the front where it flows out into the levées, thereby increasing the length of the channel and permitting the front to advance. The front shows a relationship of critical height to the yield strength (τ₀) surface tension, and slope; its continued movement is activated by the pressure of the advancing lava in the channel behind. For an ideal flow (isothermal, homogeneous, and isotropic) the ratio of the section of channel proper to the section of levées is calculated and the distance the front will have moved at any time t_x can be determined once the amount of lava available to the front is known. Assuming that the velocity function of the front {G(t)} during the collapsing stage is proportional to the entrance pressure of the lava at the channel-front boundary, an exponential decrease of velocity through time is predicted, which shows good agreement with actual frontal velocity measurements taken on two flows. Local variations in slope have a secondary effect on frontal velocities.Under conditions of constant volume the frontal zone can be considered as a machine that consumes energy brought in by the lava to perform work (front advancement). While the front will use its potential energy to run the process, the velocity at which it occurs is controlled by the activation energy that enters the system as the kinetic energy of the lava flowing into the front. A relation for the energy contribution due to frontal acceleration is also derived. Finally the entrance pressure, that permits the front to deform, is calculated. Its small value confirms that the lava behaves very much like a Bingham plastic.     

Exceedance probability distributions of steady conduit flow in karst aquifers

Discharge exceedance probabilities are calculated for a simple model karst aquifer composed of a few multilevel conduits with recharge from a single sinking stream with an exponential flow exceedance distribution. It is assumed that outflow instantaneously matches inflow, so that the conduit volume is constant but discharge is governed by the head in a volumeless shaft at the top end of the system. It is shown that small single conduit aquifers will frequently overflow at the surface during floods and the exceedance probability of flow through the aquifer and over the surface can be defined as a function of the inflow distribution and the form of the aquifer. Systems with multiple conduits will overflow less frequently, but each conduit will exhibit a flow distribution characteristic of its form and position in the vertical hierarchy. Comparison of these findings with actual flow data from a conduit aquifer shows that the approach is valid, although imprecise. The model is unlikely to be applied directly, as it requires unusually detailed data. However, it provides for the first time an indication of typical flow durations for surface overflows and individual conduits in a karst aquifer. Contrasts in flow duration will have a profound influence on solutional and sedimentary processes in the karst system.     

Simulation of dam-break flow with grid adaptation

The unsteady free surface flow caused by sudden collapse of a dam produces discontinuities in the flow variables. As the flow surges downstream, it forms a moving bore front with steep gradients of water height and velocity. In the numerical simulation of this flow, proper grid distribution can play a crucial part in the prediction and resolution of the solutions. The use of presently available numerical schemes to solve this problem on a uniform course grid system fails to resolve the characteristic flow features and hence do a poor job in simulating this flow. In this paper, an adaptive grid which adjusts itself as the solution evolves is used for a better resolution of the flow properties. Rai and Anderson's¹² method is used to determine the grid speed; however, a different partial differential equation based on the conservative principle of grid arc lengths for clustering grids in one-dimensional flow is used along with the St. Venant equations to numerically simulate the flow. Both the subcritical and the supercritical flows under extreme boundary conditions are solved using this technique. With a specified number of grid points, this provides better quality solutions as compared to those obtained with uniformly distributed grids.     

Nonequilibrium flow in volcanic conduits and application to the eruptions of Mt. St. Helens on May 18, 1980, and Vesuvius in AD 79

A steady-state, one-dimensional, and nonhomogeneous two-phase flow model was developed for the prediction of local flow properties in volcanic conduits. The model incorporates the effects of relative velocity between the phases and for the variable magma viscosity. The resulting set of nonlinear differential equations was solved by a stiff numerical solver and the results were verified with the results of basaltic fissure eruptions obtained by a homogeneous two-phase flow model, before applying the model to the eruptions of Mt. St. Helens and Vesuvius volcanoes. This verification, and a study of the sensitivity of several modeling parameters, proved effective in establishing the confidence in the predicted nonequilibrium results of flow distribution in the conduits when the mass flow rate is critical or maximum. The application of the model to the plinian eruptions of Mt. St. Helens on May 18, 1980, and Vesuvius in AD 79, demonstrates the sensitivity of the magma discharge rate and distributions of pressure, volumetric fraction, and velocities of phases, on the hydrous magma viscosity feeding the volcanic conduits. Larger magma viscosities produce smaller mass discharge rates (or greater conduit diameters), smaller exit pressures, larger disequilibrium between the phases, and larger difference between the local lithostatic and fluid pressures in the conduit. This large pressure difference occurs when magma fragments and may cause a rupture of the conduit wall rocks, producing a closure of the conduit and cessation of the volcanic eruption, or water pouring into the conduit from underground aquifers leading to phreatomagmatic explosions. The motion of the magma fragmentation zone along a conduit during an eruption can be caused by the varying viscosity of magma feeding the volcanic conduit and may cause intermittent phreatomagmatic explosions during the plinian phases as different underground aquifers are activated at different depths. The variation of magma viscosity during the eruptions of Mt. St. Helens in 1980 and Vesuvius in AD 79 is normally associated with the tapping of magmas from different depths of the magma chambers. This variation of viscosity, which can include different crystal and dissolved water contents, can also produce conduit wall erosion, the onset and collapse of volcanic columns above the vent, and the onset and cessation of pyroclastic flows and surges.     

Analytical Solution of Advective Mixing in a Conduit

The permeable conduit wall in a karst aquifer allows for water and solute to be exchanged between conduits and the limestone matrix. Contaminant sequestered in the limestone matrix is flushed into conduits following flood events. The contaminant released from the permeable wall will then mix with conduit water and will be transported downgradient in the conduit. A one-dimensional advection-dispersion equation is presented to describe this mixing-transport incorporating water flow and solute flux through the conduit wall. An analytical solution ignoring conduit dispersion is derived using the method of characteristics. Scale analysis is performed to provide a general guideline to estimate when conduit dispersion can be neglected. The solution also can be used to compute the distribution of solute in the matrix before flushing.     

A low‐dimensional model of bedrock weathering and lateral flow coevolution in hillslopes: 2. Controls on weathering and permeability profiles,drainage hydraulics,and solute export pathways

The advance of a chemical weathering front into the bedrock of a hillslope is often limited by the rate weathering products that can be carried away, maintaining chemical disequilibrium. If the weathering front is within the saturated zone, groundwater flow downslope may affect the rate of transport and weathering—however, weathering also modifies the rock permeability and the subsurface potential gradient that drives lateral groundwater flow. This feedback may help explain why there tends to be neither “runaway weathering” to great depth nor exposed bedrock covering much of the earth and may provide a mechanism for weathering front advance to keep pace with incision of adjacent streams into bedrock. This is the second of a two‐part paper exploring the coevolution of bedrock weathering and lateral flow in hillslopes using a simple low‐dimensional model based on hydraulic groundwater theory. Here, we show how a simplified kinetic model of 1‐D rock weathering can be extended to consider lateral flow in a 2‐D hillslope. Exact and approximate analytical solutions for the location and thickness of weathering within the hillslope are obtained for a number of cases. A location for the weathering front can be found such that lateral flow is able to export weathering products at the rate required to keep pace with stream incision at steady state. Three pathways of solute export are identified: “diffusing up,” where solutes diffuse up and away from the weathering front into the laterally flowing aquifer; “draining down,” where solutes are advected primarily downward into the unweathered bedrock; and “draining along,” where solutes travel laterally within the weathering zone. For each pathway, a different subsurface topography and overall relief of unweathered bedrock within the hillslope is needed to remove solutes at steady state. The relief each pathway requires depends on the rate of stream incision raised to a different power, such that at a given incision rate, one pathway requires minimal relief and, therefore, likely determines the steady‐state hillslope profile.     

Patterns of Entrapped Air Dissolution in a Two‐Dimensional Pilot‐Scale Synthetic Aquifer

Past studies of entrapped air dissolution have focused on one‐dimensional laboratory columns. Here the multidimensional nature of entrapped air dissolution was investigated using an indoor tank (180 × 240 × 600 cm³) simulating an unconfined sand aquifer with horizontal flow. Time domain reflectometry (TDR) probes directly measured entrapped air contents, while dissolved gas conditions were monitored with total dissolved gas pressure (P_TDG) probes. Dissolution occurred as a diffuse wedge‐shaped front from the inlet downgradient, with preferential dissolution at depth. This pattern was mainly attributed to increased gas solubility, as shown by P_TDG measurements. However, compression of entrapped air at greater depths, captured by TDR and leading to lower quasi‐saturated hydraulic conductivities and thus greater velocities, also played a small role. Linear propagation of the dissolution front downgradient was observed at each depth, with both TDR and P_TDG, with increasing rates with depth (e.g, 4.1 to 5.7× slower at 15 cm vs. 165 cm depth). P_TDG values revealed equilibrium with the entrapped gas initially, being higher at greater depth and fluctuating with the barometric pressure, before declining concurrently with entrapped air contents to the lower P_TDG of the source water. The observed dissolution pattern has long‐term implications for a wide variety of groundwater management issues, from recharge to contaminant transport and remediation strategies, due to the persistence of entrapped air near the water table (potential timescale of years). This study also demonstrated the utility of P_TDG probes for simple in situ measurements to detect entrapped air and monitor its dissolution.