A multiphase flow and multispecies reactive transport model for DNAPL-involved Compound Specific Isotope Analysis

This study presents a multiphase flow and multispecies reactive transport model for the simultaneous simulation of NAPL and groundwater flow, dissolution, and reactive transport with isotope fractionation, which can be used for better interpretation of NAPL-involved Compound Specific Isotope Analysis in 3D heterogeneous hydrogeologic systems. The model was verified for NAPL-aqueous phase equilibrium partitioning, aqueous phase multi-chain and multi-component reactive transport, and aqueous phase multi-component transport with isotope fractionation. Several illustrative examples are presented to investigate the effect of DNAPL spill rates, degradation rate constants, and enrichment factors on the temporal and spatial distribution of the isotope signatures of chlorinated aliphatic hydrocarbon groundwater plumes. The results clearly indicate that isotope signatures can be significantly different when considering multiphase flow within the source zone. A series of simulations indicate that degradation and isotope enrichment compete with dissolution to determine the isotope signatures in the source zone: isotope ratios remain the same as those of the source if dissolution dominates the reaction, while heavy isotopes are enriched in reactants along groundwater plume flow paths when degradation becomes dominant. It is also shown that NAPL composition can change from that of the injected source due to the partitioning of components between the aqueous and NAPL phases even when degradation is not allowed in NAPL phase. The three-dimensional simulation is presented to mechanistically illustrate the complexities in determining and interpreting the isotopic signatures with evolving DNAPL source architecture.     

Impact of capillary hysteresis and trapping on vertically integrated models for CO2 storage

Vertically integrated models are frequently applied to study subsurface flow related to CO₂ storage scenarios in saline aquifers. In this paper, we study the impact of capillary-pressure hysteresis and CO₂ trapping on the integrated constitutive parameter functions. Our results show that for the initial drainage and a subsequent imbibition, trapping is the dominant contributor to hysteresis in integrated models. We also find that for advective processes like injection and plume migration in a sloped aquifer the correct treatment of the hysteretic nature of the capillary fringe is likely of secondary importance. However, for diffusive/dispersive processes such as a redistribution of the CO₂ plume due to buoyancy and capillary forces, the hysteretic nature of the capillary fringe may significantly impact the final distribution of the fluids and the timescale of the redistribution.     

A mixed-dimensional finite volume method for two-phase flow in fractured porous media

We present a vertex-centered finite volume method for the fully coupled, fully implicit discretization of two-phase flow in fractured porous media. Fractures are discretely modeled as lower dimensional elements. The method works on unstructured, locally refined grids and on parallel computers with distributed memory. An implicit time discretization is employed and the nonlinear systems of equations are solved with a parallel Newton-multigrid method. Results from two-dimensional and three-dimensional simulations are presented.     

油井多相流电磁成像测量敏感场仿真

为了考察油井多相流电磁成像测量的敏感场分布,本文根据阵列电极结构和测量电磁场特性,应用有限元求解流体截面的电磁场问题,对于不同介质分布模型,计算电势分布,进而模拟测量敏感场.仿真结果表明,测量敏感场呈马鞍面状分布,在靠近发射电极和测量电极的区域敏感性较强,在弧形区域的中间敏感性较小.     

A note on benchmarking of numerical models for density dependent flow in porous media

Verification of numerical models for density dependent flow in porous media (DDFPM) by the means of appropriate benchmark problems is a very important step in developing and using these models. Recently, Infinite Horizontal Box (IHB) problem was suggested as a possible benchmark problem for verification of DDFPM codes. IHB is based on Horton–Rogers–Lapwood (HRL) problem. Suitability of this problem for the benchmarking purpose has been investigated in this paper. It is shown that the wavelength of instabilities fails to be a proper criterion to be considered for this problem. However, the threshold of instability formation has been found to be appropriate for benchmarking purpose.     

A stochastic model for air injection into saturated porous media

Air injection into porous media is investigated by laboratory experiments and numerical modelling. Typical applications of air injection into a granular bed are aerated bio-filters and air sparging of aquifers. The first stage of the dynamic process consists of air injection into a fixed or a quasi-fixed water-saturated granular bed. Later stages could include stages of movable beds as well, but are not further investigated here. A series of laboratory experiments were conducted in a two-dimensional box of the size 60 cm × 38 cm × 0.55 cm consisting of glass walls and using glass beads of diameter 0.4–0.6 mm as granular material. The development of the air flow pattern was optically observed and registered using a digital video camera. The resulting transient air flow pattern can be characterized as channelled flow in a fixed porous medium with dynamic tree-like evolution behaviour. Attempts are undertaken to model the air injection process. Multiphase pore-scale modelling is currently disregarded since it is restricted to very small scales. Invasion percolation models taking into account gravity effects are usually restricted to slow processes. On the other hand a continuum-type two-phase flow modelling approach is not able to simulate the observed air flow pattern. Instead a stochastic continuum-type approach is discussed here, which incorporates pore-scale features on a subscale, relevant for the immiscible processes involved. Consequently, the physical process can be modelled in a stochastic manner only, where the single experiment represents one of many possible realizations. However, the present procedure retains realistic water and air saturation patterns and therefore produces similar finger lengths and widths as observed in the experiments. Monte Carlo type modelling leads to ensemble mean water saturation and the related variance.     

Nonequilibrium effects in models of three-phase flow in porous media

In this paper we extend to three-phase flow the nonequilibrium formalism proposed by Barenblatt and co-workers for two-phase porous media flow. The underlying idea is to include nonequilibrium effects by introducing a pair of effective water and gas saturations, which are linked to the actual saturations by a local evolution equation. We illustrate and analyze how nonequilibrium effects lead to qualitative and quantitative differences in the solution of the three-phase flow equations.     

An accurate time integration method for simplified overland flow models

An accurate time integration method for the diffusion-wave and kinematic-wave approximated models for the overland flow is proposed. The discretization of the first- and second-order spatial derivatives in the basic equation is obtained by using the second-order Lax–Wendroff and the three-point centred finite difference schemes, respectively. For the solution in time, the system of ordinary differential equations, obtained by the linearization of the celerity and of the hydraulic diffusivity by Taylor series expansions, is integrated analytically. The stability and the numerical dissipation and dispersion are investigated by the Fourier analysis. A proper Courant number, and the corresponding time step for the numerical simulations can be established. In addition, the proposed diffusion-wave and kinematic-wave models are straightforwardly extended to the two-dimensional flow. Test cases for both one- and two-dimensional problems, compare the solutions of the diffusion-wave and kinematic-wave models with analytical solutions, with experimental results and with numerical solutions obtained by the Saint–Venant equations. These simulations show that the proposed numerical–analytical models accurately predict the overland flow for several situations, in particular for unsteady rainfall rate and for spatial variations of the surface roughness.     

Curvilinear immersed boundary method for simulating coupled flow and bed morphodynamic interactions due to sediment transport phenomena

The fluid-structure interaction curvilinear immersed boundary (FSI-CURVIB) numerical method of Borazjani et al. [3] is extended to simulate coupled flow and sediment transport phenomena in turbulent open-channel flows. The mobile channel bed is discretized with an unstructured triangular mesh and is treated as a sharp-interface immersed boundary embedded in a background curvilinear mesh used to discretize the general channel outline. The unsteady Reynolds-averaged Navier-Stokes (URANS) equations closed with the k − ω turbulence model are solved numerically on a hybrid staggered/non-staggered grid using a second-order accurate fractional step method. The bed deformation is calculated by solving the sediment continuity equation in the bed-load layer using an unstructured, finite-volume formulation that is consistent with the CURVIB framework. Both the first-order upwind and the higher-order hybrid GAMMA schemes [12] are implemented to discretize the bed-load flux gradients and their relative accuracy is evaluated through a systematic grid refinement study. The GAMMA scheme is employed in conjunction with a sand-slide algorithm for limiting the bed slope at locations where the material angle of repose condition is violated. The flow and bed deformation equations are coupled using the partitioned loose-coupling FSI-CURVIB approach [3]. The hydrodynamic module of the method is validated by applying it to simulate the flow in an 180° open-channel bend with fixed bed. To demonstrate the ability of the model to simulate bed morphodynamics and evaluate its accuracy, we apply it to calculate turbulent flow through two mobile-bed open channels, with 90° and 135° bends, respectively, for which experimental measurements are available.     

A simple overland flow calculation method for distributed groundwater recharge models

A method to improve the calculation of overland flow in distributed groundwater recharge models is presented and applied to two sub‐catchments in the Thames Basin, UK. Recharge calculation studies tend to simulate the runoff flow component of river flow in a simplistic way, often as a fraction of rainfall over a particular period. The method outlined in this study intends to improve the calculation of groundwater recharge estimates in distributed recharge models but does not present an alternative to complex overland flow simulators. This method uses seasonally varying coefficients to calculate runoff for specified hydrogeological classes or runoff zones, which are used to model baseflow index variations across the basin. It employs a transfer function model to represent catchment storage. Monte Carlo simulation was applied to refine the runoff values. Decoupling the runoff zones between the two sub‐catchments produces a better match between the simulated and observed values; however, the difference between observed runoff and the simulated output indicates other factors, such as landuse and topographical characteristics that affect the generation of runoff flow, need to be taken into account when classifying runoff zones. British Geological Survey © NERC 2011. Hydrological Processes © 2011 John Wiley & Sons, Ltd.     

A numerical method for two-phase flow in fractured porous media with non-matching grids

We propose a novel computational method for the efficient simulation of two-phase flow in fractured porous media. Instead of refining the grid to capture the flow along the faults or fractures, we represent the latter as immersed interfaces, using a reduced model for the flow and suitable coupling conditions. We allow for non matching grids between the porous matrix and the fractures to increase the flexibility of the method in realistic cases. We employ the extended finite element method for the Darcy problem and a finite volume method that is able to handle cut cells and matrix-fracture interactions for the saturation equation. Moreover, we address through numerical experiments the problem of the choice of a suitable numerical flux in the case of a discontinuous flux function at the interface between the fracture and the porous matrix. A wrong approximate solution of the Riemann problem can yield unphysical solutions even in simple cases.     

A contribution to risk analysis for leakage through abandoned wells in geological CO2 storage

The selection and the subsequent design of a subsurface CO₂ storage system are subject to considerable uncertainty. It is therefore important to assess the potential risks for health, safety and environment. This study contributes to the development of methods for quantitative risk assessment of CO₂ leakage from subsurface reservoirs. The amounts of leaking CO₂ are estimated by evaluating the extent of CO₂ plumes after numerically simulating a large number of reservoir realizations with a radially symmetric, homogeneous model. To conduct the computationally very expensive simulations, the ‘CO₂ Community Grid’ was used, which allows the execution of many parallel simulations simultaneously. The individual realizations are set up by randomly choosing reservoir properties from statistical distributions. The statistical characteristics of these distributions have been calculated from a large reservoir database, holding data from over 1200 reservoirs. An analytical risk equation is given, allowing the calculation of average risk due to multiple leaky wells with varying distance in the surrounding of the injection well. The reservoir parameters most affecting risk are identified. Using these results, the placement of an injection well can be optimized with respect to risk and uncertainty of leakage. The risk and uncertainty assessment can be used to determine whether a site, compared to others, should be considered for further investigations or rejected for CO₂ storage.     

Approaches to modelling coupled flow and reaction in a 2D cementation experiment

Porosity evolution at reactive interfaces is a key process that governs the evolution and performances of many engineered systems that have important applications in earth and environmental sciences. This is the case, for example, at the interface between cement structures and clays in deep geological nuclear waste disposals. Although in a different transport regime, similar questions arise for permeable reactive barriers used for biogeochemical remediation in surface environments.The COMEDIE project aims at investigating the coupling between transport, hydrodynamics and chemistry when significant variations of porosity occur. The present work focuses on a numerical benchmark used as a design exercise for the future COMEDIE-2D experiment. The use of reactive transport simulation tools like Hytec and Crunch provides predictions of the physico-chemical evolutions that are expected during the future experiments in laboratory. Focus is given in this paper on the evolution during the simulated experiment of precipitate, permeability and porosity fields.A first case is considered in which the porosity is constant. Results obtained with Crunch and Hytec are in relatively good agreement. Differences are attributable to the models of reactive surface area taken into account for dissolution/precipitation processes. Crunch and Hytec simulations taking into account porosity variations are then presented and compared. Results given by the two codes are in qualitative agreement, with differences attributable in part to the models of reactive surface area for dissolution/precipitation processes. As a consequence, the localization of secondary precipitates predicted by Crunch leads to lower local porosities than for predictions obtained by Hytec and thus to a stronger coupling between flow and chemistry. This benchmark highlights the importance of the surface area model employed to describe systems in which strong porosity variations occur as a result of dissolution/precipitation. The simulation of highly non-linear reactive transport systems is also shown to be partly dependent on specific numerical approaches.     

2DV modelling of sediment transport processes over full-scale ripples in regular asymmetric oscillatory flow

Wave-induced, steep vortex ripples are ubiquitous features in shallow coastal seas and it is therefore important to fully understand and model the sediment transport processes that occur over them. To this end, two two-dimensional vertical (2DV) models have been critically tested against detailed velocity and sediment concentration measurements above mobile ripples in regular asymmetric oscillatory flow. The two models are a k–ω turbulence-closure model and a discrete-vortex, particle-tracking (DVPT) model, while the data are obtained in the Aberdeen oscillatory flow tunnel (AOFT). The models and the data demonstrate that the time-dependent velocity and suspended sediment concentration above the ripple are dominated by the generation of lee-side vortices and their subsequent ejection at flow reversal. The DVPT model predicts the positions and strengths of the vortices reasonably well, but tends to overpredict the velocity close to the ripple surface. The k–ω model, on the other hand, underpredicts the height to which the vortices are lifted, but is better able to predict the velocity close to the bed. In terms of the cycle- and ripple-averaged horizontal velocity, both models are able to reproduce the observed offshore flow close to and below the ripple crest and the DVPT model is able to produce the onshore flow higher up. In the vicinity of the vortices, the DVPT model better represents the concentration (because of its better prediction of vorticity). The k–ω model, on the other hand, better represents the concentration close to the ripple surface and higher up in the flow (because of the better representation of the near-bed flow and background turbulence). The measured and predicted cycle- and ripple-averaged suspended sediment concentrations are in reasonable agreement and demonstrate the expected region of exponential decay. The models are able to reproduce the observed offshore cycle- and ripple-averaged suspended sediment flux from the ripple troughs upwards, and as a result, produce net offshore suspended sediment transport rates that are in reasonable agreement. The net measured offshore suspended transport rate, based on the integration of fluxes, was found to be consistent with the total net offshore transport measured in the tunnel as a whole once the onshore transport resulting from ripple migration was taken into account, as would be expected. This demonstrates the importance of models being able to predict ripple-migration rates. However, at present neither of the models is able to do so.     

A new lattice Boltzmann model for interface reactions between immiscible fluids

In this paper, we describe a lattice Boltzmann model to simulate chemical reactions taking place at the interface between two immiscible fluids. The phase-field approach is used to identify the interface and its orientation, the concentration of reactant at the interface is then calculated iteratively to impose the correct reactive flux condition. The main advantages of the model is that interfaces are considered part of the bulk dynamics with the corrective reactive flux introduced as a source/sink term in the collision step, and, as a consequence, the model’s implementation and performance is independent of the interface geometry and orientation. Results obtained with the proposed model are compared to analytical solution for three different benchmark tests (stationary flat boundary, moving flat boundary and dissolving droplet). We find an excellent agreement between analytical and numerical solutions in all cases. Finally, we present a simulation coupling the Shan Chen multiphase model and the interface reactive model to simulate the dissolution of a collection of immiscible droplets with different sizes rising by buoyancy in a stagnant fluid.     

Soil pipe flow tracer experiments: 2. Application of a streamflow transient storage zone model

Soil pipes are important subsurface flow pathways in many soil erosion phenomena. However, limited research has been performed on quantifying and characterizing their flow and transport characteristics. The objectives of this research were to determine the applicability of a streamflow model with transient storage in deriving flow and transport characteristics of soil pipes. Tracer data from pulse inputs were collected in four different soil pipes after a fluorescein dye was injected in the upstream end of each soil pipe network in three branches (west, middle, and east) of a main catchment and a back catchment in Goodwin Creek Experimental Watershed in Mississippi. Multiple sampling stations were positioned along each soil pipe network. The transient storage zone model OTIS‐P was executed inversely to estimate transport parameters by soil pipe reach such as the soil pipe cross‐sectional area (A), soil storage zone cross‐sectional area (A_s), and exchange rate between the soil pipe and the soil storage zone (α_s). Model convergence was achieved, and simulated breakthrough curves of the reaches were in good agreement with actual tracer data for eight of the nine reaches of the three branches of the Main Catchment and five of the seven reaches of the Back Catchment soil pipe. Simulation parameters for the soil pipe networks were similar to the range of values reported for flow and transport characteristics commonly observed in streams. Inversely, estimated soil pipe flow velocities were higher with increased tortuosity, which led to a smaller cross‐sectional areas predicted for the soil pipe flowpaths, while other parameters were not sensitive to tortuosity. In general, application of One‐Dimensional Transport with Inflow and Storage‐P to this unique soil pipe condition suggested larger transient storage (A_s and α_s) compared with most stream systems. This was hypothesized to be because of relatively higher ratio of the wetted perimeter to flow area in the soil pipe, the hydraulic roughness of the soil pipe, potential retention in collapsed portions of the pipe, and interaction with smaller preferential flow systems. Copyright © 2015 John Wiley & Sons, Ltd.     

A generalized transformation approach for simulating steady-state variably-saturated subsurface flow

A numerical method is proposed to accurately and efficiently compute a direct steady-state solution of the nonlinear Richards equation. In the proposed method, the Kirchhoff integral transformation and a complementary transformation are applied to the governing equation in order to separate the nonlinear hyperbolic characteristic from the linear parabolic part. The separation allows the transformed governing equation to be applied to partially- to fully-saturated systems with arbitrary constitutive relations between primary (pressure head) and secondary variables (relative permeability). The transformed governing equation is then discretized with control volume finite difference/finite element approximations, followed by inverse transformation. The approach is compared to analytical and other numerical approaches for variably-saturated flow in 1-D and 3-D domains. The results clearly demonstrate that the approach is not only more computationally efficient but also more accurate than traditional numerical solutions. The approach is also applied to an example flow problem involving a regional-scale variably-saturated heterogeneous system, where the vadose zone is up to 1 km thick. The performance, stability, and effectiveness of the transform approach is exemplified for this complex heterogeneous example, which is typical of many problems encountered in the field. It is shown that computational performance can be enhanced by several orders of magnitude with the described integral transformation approach.     

Efficient approximations for the simulation of density driven flow in porous media

Due to the non-linear coupling between flow and transport equations the simulation of real density driven flow problems requires a lot of computational time and/or heavy equipments. We suggest some approximations and numerical recipes to reduce the CPU costs for these strongly non-linear coupled equations without loss in accuracy.