Finite element modelling of recirculating density-driven saltwater intrusion processes in groundwater

This paper is concerned with an appropriate finite element model for solving transient and steady-state density-driven saltwater intrusion processes in groundwater with or without free surface. Alternative formulations and several numerical solution strategies are discussed. A comparative study is given for a cellular benchmark problem. An implicit predictor-corrector method with automatically adjusted time stepping involving a full Newton method is studied. Results on saltwater wedge displacement and gravity-driven sinking of saltwater from a deposit into the subsurface illustrate the efficacy of the model chosen. The results for the saltwater wedge problems are compared with analytical sharp-interface assessments.     

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.     

Hydraulic-Head Formulation for Density-Dependent Flow and Transport

Density-dependent flow and transport solutions for coastal saltwater intrusion investigations, analyses of fluid injection into deep brines, and studies of convective fingering and instabilities of denser fluids moving through less dense fluids typically formulate the groundwater flow equation in terms of pressure or equivalent freshwater head. A formulation of the flow equation in terms of hydraulic head is presented here as an alternative. The hydraulic-head formulation can facilitate adaptation of existing constant-density groundwater flow codes to include density-driven flow by avoiding the need to convert between freshwater head and hydraulic head within the code and by incorporating density-dependent terms as a compartmentalized “correction” to constant-density calculations already performed by the code. The hydraulic-head formulation also accommodates complexities such as unconfined groundwater flow and Newton-Raphson solution schemes more readily than the freshwater-head formulation. Simulation results are presented for four example problems solved using an implementation of the hydraulic-head formulation in MODFLOW.     

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 general approach to advective–dispersive transport with multirate mass transfer

A numerical solution that is significantly more general than other semi-analytical solutions is presented for governing equations describing advective–dispersive transport with multirate mass transfer between mobile and immobile domains. The new solution approach is general in the sense that it does not impose any restrictive assumption on the spatial or temporal variability of advective and dispersive processes in the mobile domain. A single integro-differential equation (IDE) is developed for the concentration in the mobile domain by separating the concentration in the immobile domain from the set of two partial differential equations. The solution to the IDE requires the evaluation of a temporal integral of the concentration in the mobile domain, which is a function of the Laplace transform of the distribution of the mass transfer rate coefficient. The Laplace transform is not limited to flow fields with known constant velocities. The solutions for one- and two-dimensional examples obtained using the new approach agree with those obtained by existing semi-analytical and numerical approaches.     

A new solution and data analysis for gas flow to a barometric pumping well

When an open well is installed in an unsaturated zone, gas will flow to/from the well depending on the pressure difference between the well and the surrounding media. This process is called barometric pumping and the well is called a barometric pumping well (BPW). Understanding subsurface gas pressure distribution and gas flow rate to/from a BPW is indispensable to optimize passive soil vapor extraction. This study presents a 2-D semi-analytical solution to calculate the subsurface gas pressure and gas flow rate to/from a BPW with and without a check valve. The problem is conceptualized as a mixed-type boundary value problem. The solution for pumping without a check valve is used to analyze the behavior of the radius of influence (ROI). Results show that ROI is time-dependent. It increases with radial gas permeability and decreases with vertical gas permeability. Field application of the solution without a check valve demonstrates the high accuracy of the developed solution.     

Semi-analytical Solution for Viscoelastic Relaxation in Eccentrically-nested Spheres Induced by Surface Toroidal Traction

A semi-analytical solution to the 2-D forward modelling of viscoelastic relaxation in a heterogeneous sphere induced by a surface toroidal force is derived. The model consists of a concentrically-nested elastic lithosphere, a viscoelastic mantle, and an eccentrically-nested viscoelastic core. Since numerical codes based on finite-element or spectral-finite-difference techniques for modelling viscoelastic relaxation in a spherical geometry in the presence of lateral viscosity variations are becoming more popular, reliable examples for testing and validating such codes are essential. The eccentrically-nested sphere solution has been tested by comparing it with two distinct results: The analytical solution for viscoelastic relaxation in concentrically-nested spheres and the time domain, spectral finite-element numerical solution for viscoelastic relaxation in eccentrically-nested spheres, with excellent agreement being obtained.     

Variable-density groundwater flow and solute transport in irregular 2D fracture networks

Numerical simulations of variable-density flow and solute transport have been conducted to investigate dense plume migration for various configurations of 2D fracture networks. For orthogonal fractures, simulations demonstrate that dispersive mixing in fractures with small aperture does not stabilize vertical plume migration in fractures with large aperture. Simulations in non-orthogonal 2D fracture networks indicate that convection cells form and that they overlap both the porous matrix and fractures. Thus, transport rates in convection cells depend on matrix and fracture flow properties. A series of simulations in statistically equivalent networks of fractures with irregular orientation show that the migration of a dense plume is highly sensitive to the geometry of the network. If fractures in a random network are connected equidistantly to the solute source, few equidistantly distributed fractures favor density-driven transport. On the other hand, numerous fractures have a stabilizing effect, especially if diffusive transport rates are high. A sensitivity analysis for a network with few equidistantly distributed fractures shows that low fracture aperture, low matrix permeability and high matrix porosity impede density-driven transport because these parameters reduce groundwater flow velocities in both the matrix and the fractures. Enhanced molecular diffusion slows down density-driven transport because it favors solute diffusion from the fractures into the low-permeability porous matrix where groundwater velocities are smaller. For the configurations tested, variable-density flow and solute transport are most sensitive to the permeability and porosity of the matrix, which are properties that can be determined more accurately than the geometry and hydraulic properties of the fracture network, which have a smaller impact on density-driven transport.     

正交加肋地基板在谐荷载作用下的动力响应

将加肋地基板作为板-梁耦合系统，考虑弹性地基模型，建立土—基础耦合系统的运动微分方程，通过富氏积分变换及其逆变换，演绎出双向正交加肋钢筋混凝土地基板在设备谐激力作用下的振动响应的半解析解法。结果说明，在土建设计中，为了减少振动对环境的有害影响，应考虑动力设备基础的隔振和减振问题。数值算例也表明，本方法有较好的计算精度。     

Balancing the source terms in a SPH model for solving the shallow water equations

A shallow flow generally features complex hydrodynamics induced by complicated domain topography and geometry. A numerical scheme with well-balanced flux and source term gradients is therefore essential before a shallow flow model can be applied to simulate real-world problems. The issue of source term balancing has been exhaustively investigated in grid-based numerical approaches, e.g. discontinuous Galerkin finite element methods and finite volume Godunov-type methods. In recent years, a relatively new computational method, smooth particle hydrodynamics (SPH), has started to gain popularity in solving the shallow water equations (SWEs). However, the well-balanced problem has not been fully investigated and resolved in the context of SPH. This work aims to discuss the well-balanced problem caused by a standard SPH discretization to the SWEs with slope source terms and derive a corrected SPH algorithm that is able to preserve the solution of lake at rest. In order to enhance the shock capturing capability of the resulting SPH model, the Monotone Upwind-centered Scheme for Conservation Laws (MUSCL) is also explored and applied to enable Riemann solver based artificial viscosity. The new SPH model is validated against several idealized benchmark tests and a real-world dam-break case and promising results are obtained.     

Vertical structure of residual slope circulation driven by JEBAR and tides: an idealised model

The residual circulation over the continental slope, and in particular, its vertical structure, is analysed by means of an idealised hydrodynamic model. The model is based on the depth-dependent shallow-water equations under uniform along-isobath conditions and is forced by a prescribed meridional density gradient and tidal velocities. By means of expansion in the small Rossby number solutions are analysed for conditions representative for the continental slopes off the Hebrides and in the Bay of Biscay. The steady solution at zeroth order consists of a linear density-driven flow. At order a tidally rectified flow is found and a stationary flow due to self-interaction of the zeroth-order density-driven flow. At order ² the leading-order effect of the interaction between the zeroth-order density-driven flow and the tides is found: the ‘interaction current’. The solutions up to and including order ² constitute an along-isobath steady slope current which is comparable to field data. The slope current and the accompanying cross-shelf circulation depend strongly on the shelf and flow characteristics. For the Hebridean case the density forcing predominates, but for the Biscay case the tidal effects are of the same order of magnitude as the density effects. Under those conditions the interaction current is significant which implies that linear superposition of density and tidal effects differs from the non-linear combination of both. It is also shown that the depth-average of the interaction current differs essentially from the solution obtained from a depth-averaged model.     

Theoretical analysis of the worthiness of Henry and Elder problems as benchmarks of density-dependent groundwater flow models

Computer models must be tested to ensure that the mathematical statements and solution schemes accurately represent the physical processes of interest. Because the availability of benchmark problems for testing density-dependent groundwater models is limited, one should be careful in using these problems appropriately. Details of a Galerkin finite-element model for the simulation of density-dependent, variably saturated flow processes are presented here. The model is tested using the Henry salt-water intrusion problem and Elder salt convection problem. The quality of these benchmark problems is then evaluated by solving the problems in the standard density-coupled mode and in a new density-uncoupled mode. The differences between the solutions indicate that the Henry salt-water intrusion problem has limited usefulness in benchmarking density-dependent flow models because the internal flow dynamics are largely determined by the boundary forcing. Alternatively, the Elder salt-convection problem is more suited to the model testing process because the flow patterns are completely determined by the internal balance of pressure and gravity forces.     

One-dimensional uniform and time varying solute dispersion along transient groundwater flow in a semi-infinite aquifer

An analytical solution for the space-time variation of contaminant concentration in one-dimensional transient groundwater flow in a homogenous semi-infinite aquifer, subjected to time-dependent source contamination, is derived. The uniform and time varying dispersion along transient groundwater flow is investigated under two conditions. First, the flow velocity distribution in the aquifer is considered as a sinusoidally varying function, and second, the flow velocity distribution is treated as an exponentially increasing function of time. It is assumed that initially the aquifer is not solute free, so the initial background concentration is considered as an exponentially decreasing function of the space variable which is tending to zero at infinity. It is assumed that dispersion is directly proportional to the square of the velocity, noting that experimental observations indicate that dispersion is directly proportional to the velocity with a power ranging from 1 to 2. The analytical solution is illustrated using an example and may help benchmark numerical codes and solutions.     

Eulerian spatial moments for solute transport in three-dimensional heterogeneous,dual-permeability media

A Eulerian analytical method is developed for nonreactive solute transport in heterogeneous, dual-permeability media where the hydraulic conductivities in fracture and matrix domains are both assumed to be stochastic processes. The analytical solution for the mean concentration is given explicitly in Fourier and Laplace transforms. Instead of using the fast fourier transform method to numerically invert the solution to real space (Hu et al., 2002), we apply the general relationship between spatial moments and concentration (Naff, 1990; Hu et al., 1997) to obtain the analytical solutions for the spatial moments up to the second for a pulse input of the solute. Owing to its accuracy and efficiency, the analytical method can be used to check the semi-analytical and Monte Carlo numerical methods before they are applied to more complicated studies. The analytical method can be also used during screening studies to identify the most significant transport parameters for further analysis. In this study, the analytical results have been compared with those obtained from the semi-analytical method (Hu et al., 2002) and the comparison shows that the semi-analytical method is robust. It is clearly shown from the analytical solution that the three factors, local dispersion, conductivity variation in each domain and velocity convection flow difference in the two domains, play different roles on the solute plume spreading in longitudinal and transverse directions. The calculation results also indicate that when the log-conductivity variance in matrix is 10 times less than its counterpart in fractures, it will hardly influence the solute transport, whether the conductivity field is matrix is treated as a homogeneous or random field.     

Variable-density groundwater flow and solute transport in porous media containing nonuniform discrete fractures

Variations in fluid density can greatly affect fluid flow and solute transport in the subsurface. Heterogeneities such as fractures play a major role for the migration of variable-density fluids. Earlier modeling studies of density effects in fractured media were restricted to orthogonal fracture networks, consisting of only vertical and horizontal fractures. The present study addresses the phenomenon of 3D variable-density flow and transport in fractured porous media, where fractures of an arbitrary incline can occur. A general formulation of the body force vector is derived, which accounts for variable-density flow and transport in fractures of any orientation. Simulation results are presented that show the verification of the new model formulation, for the porous matrix and for inclined fractures. Simulations of variable-density flow and solute transport are then conducted for a single fracture, embedded in a porous matrix. The simulations show that density-driven flow in the fracture causes convective flow within the porous matrix and that the high-permeability fracture acts as a barrier for convection. Other simulations were run to investigate the influence of fracture incline on plume migration. Finally, tabular data of the tracer breakthrough curve in the inclined fracture is given to facilitate the verification of other codes.     

A semi-analytical solution based on a numerical solution of the solute transport of a conservative and nonreactive tracer

In the evaluation of potential risk from ingestion of groundwater near an impacted site, numerical simulation of fate and transport processes of chemicals of concern is often required. If there is potential concern about multiple chemicals, numerical simulation of each chemical separately is often needed. In this paper, a semi-analytical solution is presented based on a numerical solution of the transport of a conservative and nonreactive tracer. When multiple chemicals undergoing sorption and first-order degradation need to be modeled, we can avoid performing individual numerical simulations for each chemical by applying the semi-analytical solution. Numerical test runs were conducted to verify the semi-analytical solution; simulation results reveal that the concentrations derived from the semi-analytical solution are identical to those derived from the individual numerical fate and transport model simulations. The semi-analytical solution requires steady-state flow conditions, no continuing contaminant source, and similar initial source concentration distributions.     

基于时空守恒元和解元(CE/SE)方法的孔隙介质多相流动计算

将时空守恒元/解元(CE/SE)方法推广到二维孔隙介质多相流问题的数值计算中,采用人工压缩法耦合速度和压力,同时结合杂交粒子水平集方法捕捉物质界面.提出一套完整的二维欧拉型孔隙介质非稳态多相不可压缩黏性流动计算方案.通过对溃坝和液滴在重力作用下的运动和变形问题的数值模拟,验证了方法的精度和有效性.在此基础上,提出了一个新的孔隙介质两相流物理模型——双层流体顶盖驱动方腔流.     

Wave-induced dynamic response of saturated multi-layer porous media: Analytical solutions and validity regions of various formulations in non-dimensional parametric space

In this paper, dynamic response of saturated-layered porous media under harmonic waves is evaluated through a semi-analytical solution. The coupled differential equations governing the dynamics of saturated or nearly saturated porous media such as soils containing all the inertial terms of solid and fluid phases are presented for a multi-layer system. Possible simplifications of the equations which are called formulations are introduced based upon the presence of inertial terms associated with the phases. The semi-analytical solutions to the response of multiple layers for all the formulations are presented in terms of pore water pressure and stress variations considering a set of non-dimensional parameters and their respective ratios. Validity of the formulations is presented in a non-dimensional parametric space. The maximum discrepancies in the pore pressure response of the formulations leading to validity regions are illustrated for typical dynamic problems. Subsequently, the effects of layering and drainage conditions on these regions are also presented. The proposed semi-analytical solution may be served as a benchmark one for validating the coupled numerical solutions, which can be used to deal with real scientific and geo-engineering problems in the emerging field of computational geomechanics.     

3D finite volume groundwater and heat transport modeling with non-orthogonal grids,using a coordinate transformation method

Many popular groundwater modeling codes are based on the finite differences or finite volume method for orthogonal grids. In cases of complex subsurface geometries this type of grid either leads to coarse geometric representations or to extremely fine meshes. We use a coordinate transformation method (CTM) to circumvent this shortcoming. In computational fluid dynamics (CFD), this method has been applied successfully to the general Navier–Stokes equation. The method is based on tensor analysis and performs a transformation of a curvilinear into a rectangular unit grid, on which a modified formulation of the differential equations is applied. Therefore, it is not necessary to reformulate the code in total. We applied the CTM to an existing three-dimensional code (SHEMAT), a simulator for heat conduction and advection in porous media. The finite volume discretization scheme for the non-orthogonal, structured, hexahedral grid leads to a 19-point stencil and a correspondingly banded system matrix. The implementation is straightforward and it is possible to use some existing routines without modification. The accuracy of the modified code is demonstrated for single phase flow on a two-dimensional analytical solution for flow and heat transport. Additionally, a simple case of potential flow is shown for a two-dimensional grid which is increasingly deformed. The result reveals that the corresponding error increases only slightly. Finally, a thermal free-convection benchmark is discussed. The result shows, that the solution obtained with the new code is in good agreement with the ones obtained by other codes.     

Wind effects on the lateral structure of density-driven circulation in Chesapeake Bay

The response of the density-driven circulation in the Chesapeake Bay to wind forcing was studied with numerical experiments. A model of the bay with realistic bathymetry was first applied to produce the density-driven flow under average river discharge and tidal forcing. Subsequently, four spatially uniform wind fields (northeasterly, northwesterly, southwesterly, and southeasterly) were imposed to examine the resulting cross-estuary structure of salinity and flow fields. In general, northeasterly and northwesterly winds intensified the density-driven circulation in the upper and middle reaches of the bay, whereas southeasterly and southwesterly winds weakened it. The response was different in the lower bay, where downwind flow from the upper and middle reaches of the bay competed with onshore/offshore coastal flows. Wind remote effects were dominant, over local effects, on volume transports through the bay entrance. However, local effects were more influential in establishing the sea-level slopes that drove subtidal flows and salinity fields in most of the bay. The effect of vertical stratification on wind-induced flows was also investigated by switching it off. The absence of stratification allowed development of Ekman layers that reached depths of the same order as the water depth. Consequently, bathymetric effects became influential on the homogeneous flow structure causing the wind-induced flow inside the bay to show a marked transverse structure: downwind over the shallow areas and upwind in the channels. In the presence of stratification, Ekman layers became shallower and the wind-induced currents showed weaker transverse structure than those that developed in the absence of stratification. In essence, the wind-driven flows were horizontally sheared under weak stratification and vertically sheared under stratified conditions.