Fluid dispersion effects on density-driven thermohaline flow and transport in porous media

This study introduces the dispersive fluid flux of total fluid mass to the density-driven flow equation to improve thermohaline modeling of salt and heat transports in porous media. The dispersive fluid flux in the flow equation is derived to account for an additional fluid flux driven by the density gradient and mechanical dispersion. The coupled flow, salt transport and heat transport governing equations are numerically solved by a fully implicit finite difference method to investigate solution changes due to the dispersive fluid flux. The numerical solutions are verified by the Henry problem and the thermal Elder problem under a moderate density effect and by the brine Elder problem under a strong density effect. It is found that increment of the maximum ratio of the dispersive fluid flux to the advective fluid flux results in increasing dispersivity for the Henry problem and the brine Elder problem. The effects of the dispersive fluid flux on salt and heat transports under high density differences and high dispersivities are more noticeable than under low density differences and low dispersivities. Values of quantitative indicators such as the Nusselt number, mass flux, salt mass stored and maximum penetration depth in the brine Elder problem show noticeable changes by the dispersive fluid flux. In the thermohaline Elder problem, the dispersive fluid flux shows a considerable effect on the shape and the number of developed fingers and makes either an upwelling or a downwelling flow in the center of the domain. In conclusion, for the general case that involves strong density-driven flow and transport modeling in porous media, the dispersive fluid flux should be considered in the flow equation.     

Coupled groundwater flow and transport: 1. Verification of variable density flow and transport models

This work examines variable density flow and corresponding solute transport in groundwater systems. Fluid dynamics of salty solutions with significant density variations are of increasing interest in many problems of subsurface hydrology. The mathematical model comprises a set of non-linear, coupled, partial differential equations to be solved for pressure/hydraulic head and mass fraction/concentration of the solute component. The governing equations and underlying assumptions are developed and discussed. The equation of solute mass conservation is formulated in terms of mass fraction and mass concentration. Different levels of the approximation of density variations in the mass balance equations are used for convection problems (e.g. the Boussinesq approximation and its extension, fully density approximation). The impact of these simplifications is studied by use of numerical modelling.Numerical models for nonlinear problems, such as density-driven convection, must be carefully verified in a particular series of tests. Standard benchmarks for proving variable density flow models are the Henry, Elder, and salt dome (HYDROCOIN level 1 case 5) problems. We studied these benchmarks using two finite element simulators - ROCKFLOW, which was developed at the Institute of Fluid Mechanics and Computer Applications in Civil Engineering and FEFLOW, which was developed at the Institute for Water Resources Planning and Systems Research Ltd. Although both simulators are based on the Galerkin finite element method, they differ in many approximation details such as temporal discretization (Crank-Nicolson vs predictor-corrector schemes), spatial discretization (triangular and quadrilateral elements), finite element basis functions (linear, bilinear, biquadratic), iteration schemes (Newton, Picard) and solvers (direct, iterative). The numerical analysis illustrates discretization effects and defects arising from the different levels of the density of approximation. We contribute new results for the salt dome problem, for which inconsistent findings exist in literature. Applications of the verified numerical models to more complex problems, such as thermohaline and three-dimensional convection systems, will be presented in the second part of this paper.     

Modeling Transport in Transient Ground-Water Flow: An Unacknowledged Approximation

Abstract. During unsteady or transient ground-water flow, the fluid mass per unit volume of aquifer changes as the potentiometric head changes, and solute transport is affected by this change in fluid storage. Three widely applied numerical models of two-dimensional transport partially account for the effects of transient flow by removing terms corresponding to the fluid continuity equation from the transport equation, resulting in a simpler governing equation. However, fluid-storage terms remaining in the transport equation that change during transient flow are, in certain cases, held constant in time in these models. For the case of increasing heads, this approximation, which is unacknowledged in these models'documentation, leads to transport velocities that are too high, and increased concentration at fluid and solute sources. If heads are dropping in time, computed transport velocities are too low. Using parameters that somewhat exaggerate the effects of this approximation, an example numerical simulation indicates solute travel time error of about 14 percent but only minor errors due to incorrect dilution volume. For horizontal flow and transport models that assume fluid density is constant, the product of porosity and aquifer thickness changes in time: initial porosity times initial thickness plus the change in head times the storage coefficient. This formula reduces to the saturated thickness in unconfined aquifers if porosity is assumed to be constant and equal to specific yield. The computational cost of this more accurate representation is insignificant and is easily incorporated in numerical models of solute transport.     

A test case for the simulation of three-dimensional variable-density flow and solute transport in discretely-fractured porous media

A test case has been developed for three-dimensional simulations of variable-density flow and solute transport in discretely-fractured porous media. The simulation domain is a low-permeability porous matrix cube containing a single non-planar fracture. The initial solute concentration is zero everywhere. A constant solute concentration is assigned to the top of the domain, which increases near-top fluid density and induces downward density-driven flow. The test case is therefore comparable to downwelling of a dense brine below a saline disposal basin or a waste repository. Numerous fingers and distinct convection cells develop early in the fracture but the fingers later coalesce and convection becomes less apparent. To help test other variable-density flow and transport models, results of the test case are presented both qualitatively (concentration contours and velocity fields) and quantitatively (penetration depth, mass flux, total mass stored, maximum fracture and matrix velocity).     

Immiscible two-phase fluid flows in deformable porous media

Macroscopic differential equations of mass and momentum balance for two immiscible fluids in a deformable porous medium are derived in an Eulerian framework using the continuum theory of mixtures. After inclusion of constitutive relationships, the resulting momentum balance equations feature terms characterizing the coupling among the fluid phases and the solid matrix caused by their relative accelerations. These terms, which imply a number of interesting phenomena, do not appear in current hydrologic models of subsurface multiphase flow. Our equations of momentum balance are shown to reduce to the Berryman–Thigpen–Chen model of bulk elastic wave propagation through unsaturated porous media after simplification (e.g., isothermal conditions, neglect of gravity, etc.) and under the assumption of constant volume fractions and material densities. When specialized to the case of a porous medium containing a single fluid and an elastic solid, our momentum balance equations reduce to the well-known Biot model of poroelasticity. We also show that mass balance alone is sufficient to derive the Biot model stress–strain relations, provided that a closure condition for porosity change suggested by de la Cruz and Spanos is invoked. Finally, a relation between elastic parameters and inertial coupling coefficients is derived that permits the partial differential equations of the Biot model to be decoupled into a telegraph equation and a wave equation whose respective dependent variables are two different linear combinations of the dilatations of the solid and the fluid.     

Modeling changes in permeability characteristics of gypsified rocks accompanying brine flow

The study is aimed at changes in the permeability of rocks, partially sealed by mineral filling, under the effect of brine, which infiltrates from a salt dumb and is diluted by groundwater. Changes in the ionic strength of solution determine saturation deficiency with respect to minerals, resulting in the dissolution of the filling and changes in the permeability of fractures. Flow in porous medium, mass transport, and thermodynamic equilibria where jointly simulated using PHT3D program with additional revaluation of hydraulic characteristics. A narrow domain with higher permeability, which was found to form at the upper boundary of the gypsified zone, has almost no effect on the waterproof properties of the above-salt stratum in a potassium salt deposit during the operation of a salt dumb.     

Hindered settling with an apparent particle diameter concept

In this work a new expression has been developed to predict the settling velocity of a sediment particle which is dispersed in a sediment-fluid mixture during a turbulent flow. A concept of apparent particle diameter has been introduced and is defined by the diameter of the spherical volume in which the particle can move randomly after collision with other particles in suspension. The effect of suspension concentration is studied on the mass density of the sediment-fluid mixture. It has been shown that the settling velocity of sediment particle in a sediment-fluid mixture is a function of different characteristics of the sediment particle such as settling velocity in clear fluid, suspension concentration, relative mass density and Reynolds number. The model has shown good agreement when compared with previously published experimental data and it’s prediction accuracy is superior than the other existing models.     

Evaporites: Relative humidity control of primary mineral facies revisited

It is widely accepted that the deposition of mineral facies of evaporite basins is controlled by the average annual relative humidity of the contiguous atmosphere, which dictates the equilibrium activity of the evaporating brine. This concept has far reaching implications in salt works and for the investigation of paleoenvironmental settings affecting depositional sequences within evaporite basins. The above concept, which dominated the scientific thought of evaporite basin investigations, suffers from two serious flaws: (a) the assumption of a static decoupled atmosphere and (b) the total neglect of energy input and thermodynamic feedbacks resulting from evaporation suppression. The present investigation will resolve the underlying mechanisms controlling the equilibrium activity of hypersaline solutions using a theoretical framework that combines energy and mass transport across the surface–atmosphere boundary. Calculations of the equilibrium activity of hypersaline solutions under isothermal conditions, as implied in the original concept, are not in line with the basic physical principles defining heat and mass exchange across the brine–atmosphere boundary and lead to substantial overestimation of actual evaporation and the activity itself. It is demonstrated that in addition to atmospheric relative humidity, the activity of hypersaline solutions is determined by numerous meteorological forcings along with hydrological, geochemical, and thermodynamic feedback mechanisms. Evaporation suppression resulting from a drop in brine activity causes substantial increase in brine temperature, which enhances vapour pressure differential across the interface, leading to more evaporation. This negative feedback shifts the brine activity downward for equilibrium to be attained. It is also demonstrated that evaporation from a brine surface usually proceeds when the relative humidity of the contiguous atmosphere is similar or even higher than that of the brine due to energy input and the strong negative feedback caused by evaporation suppression. The present investigation re‐establishes a new paradigm concerning the processes controlling evaporite basin sedimentation and palaeoclimate reconstruction as deduced from evaporite/hypersaline basin deposits. Findings have operational ramifications in the industrial applications of dissolved salt mineral extraction.     

Modeling interaction of fluid and salt in an aquifer/lagoon system

To simulate the dynamic interaction between a saline lagoon and a ground water system, a numerical model for two-dimensional, variable-density, saturated-unsaturated, and coupled flow and solute transport (saltwater intrusion by finite elements and characteristics [SIFEC]) was modified to allow the volume of water and mass of salt in the lagoon to vary with each time step. The modified SIFEC allows the stage of a lagoon to vary in accordance with a functional relation between the stage and water volume of the lagoon, and also allows the salt concentration of the lagoon to vary in accordance with the salt budget of the lagoon including chemical precipitation and dissolution of salt. The updated stage and salt concentration of the lagoon are in turn used as transient boundary conditions for the coupled flow and solute transport model. The utility of the modified model was demonstrated by applying it to the eastern Mediterranean coastal region of Turkey for assessing impacts of climate change on the subsurface environment under scenarios of sea level rise, increased evaporation, and decreased precipitation.     

Effects of fluid displacement pattern on complex electrical impedance in Berea sandstone over frequency range 104–106 Hz

To better understand the effect of fluid distribution on the electric response of rocks saturated with oil and brine, we conducted experimental studies on the complex electrical impedance in a Berea sandstone, together with in situ acquisitions of oil distribution images employing a high‐resolution medical X‐ray computed tomography. We performed two tests of brine displacement by oil under high (10 MPa) and low (5 MPa) pressures, which were accompanied by fingering and stable displacement patterns, respectively. The measured complex impedance data were fitted to the Cole model to obtain the resistance, capacitance, peak frequency of the imaginary impedance, and the exponent α of the rock–fluid system. With increasing oil saturation, the resistance showed an increasing trend, whereas the other three parameters decreased. The fingering displacement exhibited lower resistance and capacitance than the stable displacement. The analysis of the resistance changes using a simple parallel connection model indicates that there are more components of residual brine in parallel connections in the fingering pattern than in the stable displacement pattern at the same saturation. We also interpreted the normalised changes in the capacitance (or apparent dielectric constant) with respect to the oil saturation via an analysis of the shape factor of fluid distribution based on the Maxwell–Wagner–Brugermann–Hanai model. The changes in the shape factor suggest that the pinch‐off of the brine in parallel connection by the oil is a dominant mechanism reducing the capacitance. In the stable displacement, most of the connections in the brine phase are immediately pinched off by oil displacement front at a local oil saturation of 65%. Conversely, in the fingering displacement, there is a transition from the bulk or layered brine to the pinched‐off at a local oil saturation below 60%. The analyses indicate that the difference in the fluid distribution under different fluid conditions is responsible for the non‐Archie behaviour.     

Numerical modeling of lock-exchange gravity/turbidity currents by a high-order upwinding combined compact difference scheme

This study presents two-dimensional direct numerical simulations for sediment-laden current with higher density propagating forward through a lighter ambient water.The incompressible NavierStokes equations including the buoyancy force for the density difference between the light and heavy fluids are solved by a finite difference scheme based on a structured mesh.The concentration transport equations are used to explore such rich transport phenomena as gravity and turbidity currents.Within the framework of an Upwinding Combined Compact finite Difference(UCCD)scheme,rigorous determination of weighting coefficients underlies the modified equation analysis and the minimization of the numerical modified wavenumber.This sixth-order UCCD scheme is implemented in a four-point grid stencil to approximate advection and diffusion terms in the concentration transport equations and the first-order derivative terms in the Navier-Stokes equations,which can greatly enhance convective stability and increase dispersive accuracy at the same time.The initial discontinuous concentration field is smoothed by solving a newly proposed Heaviside function to prevent numerical instabilities and unreasonable concentration values.A two-step projection method is then applied to obtain the velocity field.The numerical algorithm shows a satisfying ability to capture the generation,development,and dissipation of the Kelvin-Helmholz instabilities and turbulent billows at the interface between the current and the ambient fluid.The simulation results also are compared with the data in published literatures and good agreements are found to prove that the present numerical model can well reproduce the propagation,particle deposition,and mixing processes of lock-exchange gravity and turbidity currents.     

Discharge estimation of submarine springs in the Dead Sea based on velocity or density measurements in proximity to the water surface

The Dead Sea is a closed lake, the water level of which is lowering at an alarming rate of about 1 m/year. Factors difficult to determine in its water balance are evaporation and groundwater inflow, some of which emanate as submarine groundwater discharge. A vertical buoyant jet generated by the difference in densities between the groundwater and the Dead Sea brine forms at submarine spring outlets. To characterize this flow field and to determine its volumetric discharge, a system was developed to measure the velocity and density of the ascending submarine groundwater across the center of the stream along several horizontal sections and equidistant depths while divers sampled the spring. This was also undertaken on an artificial submarine spring with a known discharge to determine the quality of the measurements and the accuracy of the method. The underwater widening of the flow is linear and independent of the volumetric spring discharge. The temperature of the Dead Sea brine at lower layers primarily determines the temperature of the surface of the upwelling, produced above the jet flow, as the origin of the main mass of water in the submarine jet flow is Dead Sea brine. Based on the measurements, a model is presented to evaluate the distribution of velocity and solute density in the flow field of an emanating buoyant jet. This model allows the calculation of the volumetric submarine discharge, merely requiring either the maximum flow velocity or the minimal density at a given depth.     

A numerical model for simulation of fluid mud with different rheological behaviors

In this paper, a three-dimensional isopycnal approach is presented to simulate the dynamics of fluid mud covering the formation, development, transport, and disappearance of fluid mud. The basic assumption is the assignment of the fluid’s density as the indicating parameter for the rheological behavior. Considering stable stratification, as is usually the case for fluid mud, layers of constant density discretize the vertical domain. The non-Newtonian dynamics of fluid mud is simulated by solving the Cauchy equations for general continuum dynamics. Instead of using a turbulent viscosity approach, the viscosity is allowed to vary according to the rheological behavior of mud suspensions. This apparent viscosity can be determined for different rheological formulations in dependence of the volume solid fraction and the shear rate. An existing three-dimensional isopycnal hydrodynamic model was extended for vertical mass transport processes and was applied on a schematic system with hindered settling. For including the rheological behavior of fluid mud, the Worrall–Tuliani approach was parameterized and implemented. The resulting flow behavior is shown on a model application of fluid mud layers moving down an inclined plane. With these changes, it is demonstrated that the isopycnal model is capable of simulating fluid mud dynamics.     

Seismic time‐lapse effects of solution salt mining – a feasibility study

This article addresses the question whether time‐lapse seismic reflection techniques can be used to follow and quantify the effects of solution salt mining. Specifically, the production of magnesium salts as mined in the north of the Netherlands is considered. The use of seismic time‐lapse techniques to follow such a production has not previously been investigated. For hydrocarbon production and CO₂ storage, time‐lapse seismics are used to look at reservoir changes mainly caused by pressure and saturation changes in large reservoirs, while for solution mining salt is produced from caverns with a limited lateral extent, with much smaller production volumes and a fluid (brine) replacing a solid (magnesium salt). In our approach we start from the present situation of the mine and then study three different production scenarios, representing salt production both in vertical and lateral directions of the mine. The present situation and future scenarios have been transformed into subsurface models that were input to an elastic finite‐difference scheme to create synthetic seismic data. These data have been analysed and processed up to migrated seismic images, such that time‐lapse analyses of intermediate and final results could be done. From the analyses, it is found that both vertical and lateral production is visible well above the detection threshold in difference data, both at pre‐imaging and post‐imaging stages. In quantitative terms, an additional production of the mine of 6 m causes time‐shifts in the order of 2 ms (pre‐imaging) and 4 ms (post‐imaging) and amplitude changes of above 20% in the imaged sections. A laterally oriented production causes even larger amplitude changes at the edge of the cavern due to replacement of solid magnesium salt with brine introducing a large seismic contrast. Overall, our pre‐imaging and post‐imaging time‐lapse analysis indicates that the effects of solution salt mining can be observed and quantified on seismic data. The effects seem large enough to be observable in real seismic data containing noise.     

Modelling evaporation processes in soils from the Huasco salt flat basin,Chile

The need to understand and simulate hydrological phenomena and their interactions, and the impact of anthropogenic and climate changes on natural environments have promoted the study of evaporation from bare soils in arid climates. In closed Altiplano basins, such as those encountered in arid and hyper arid basins in northern Chile, evaporation from shallow groundwater is the main source of aquifer depletion, and thus, its study is crucial for water resources management. The objective of this work is to understand the mechanisms of evaporation in saline soils with shallow water tables, in order to better quantify evaporation fluxes and improve our understanding of the water balance in these regions. To achieve this objective, a model that couples fluid flow with heat transfer was developed and calibrated using column experiments with saline soils from the Huasco salt flat basin, Chile. The model enables determination of both liquid and water vapour fluxes, as well as the location of the evaporation front. Experimental results showed that salt transport inside the soil profile modified the water retention curve, highlighting the importance of including salt transport when modelling the evaporation processes in these soils. Indeed, model simulations only agreed with the experimental data when the effect of salt transport on water retention curves was taken into account. Model results also showed that the evaporation front is closer to the soil surface as the water table depth reduces. Therefore, the model allows determining the groundwater level depth that results in disconnection of liquid fluxes in the vadose zone. A sensitivity analysis allowed understanding the effect of water‐flux enhancements mechanisms on soil evaporation. The results presented in this study are important as they allow quantifying the evaporation that occurs in bare soils from Altiplano basins, which is typically the main water discharge in these closed basins. Copyright © 2016 John Wiley & Sons, Ltd.