Development and evaluation of a local grid refinement method for block-centered finite-difference groundwater models using shared nodes

A new method of local grid refinement for two-dimensional block-centered finite-difference meshes is presented in the context of steady-state groundwater-flow modeling. The method uses an iteration-based feedback with shared nodes to couple two separate grids. The new method is evaluated by comparison with results using a uniform fine mesh, a variably spaced mesh, and a traditional method of local grid refinement without a feedback.Results indicate: (1) The new method exhibits quadratic convergence for homogenous systems and convergence equivalent to uniform-grid refinement for heterogeneous systems. (2) Coupling the coarse grid with the refined grid in a numerically rigorous way allowed for improvement in the coarse-grid results. (3) For heterogeneous systems, commonly used linear interpolation of heads from the large model onto the boundary of the refined model produced heads that are inconsistent with the physics of the flow field. (4) The traditional method works well in situations where the better resolution of the locally refined grid has little influence on the overall flow-system dynamics, but if this is not true, lack of a feedback mechanism produced errors in head up to 3.6% and errors in cell-to-cell flows up to 25%.     

Evaluation of MODFLOW-LGR in connection with a synthetic regional-scale model

This work studies costs and benefits of utilizing local-grid refinement (LGR) as implemented in MODFLOW-LGR to simulate groundwater flow in a buried tunnel valley interacting with a regional aquifer. Two alternative LGR methods were used: the shared-node (SN) method and the ghost-node (GN) method. To conserve flows the SN method requires correction of sources and sinks in cells at the refined/coarse-grid interface. We found that the optimal correction method is case dependent and difficult to identify in practice. However, the results showed little difference and suggest that identifying the optimal method was of minor importance in our case. The GN method does not require corrections at the models' interface, and it uses a simpler head interpolation scheme than the SN method. The simpler scheme is faster but less accurate so that more iterations may be necessary. However, the GN method solved our flow problem more efficiently than the SN method. The MODFLOW-LGR results were compared with the results obtained using a globally coarse (GC) grid. The LGR simulations required one to two orders of magnitude longer run times than the GC model. However, the improvements of the numerical resolution around the buried valley substantially increased the accuracy of simulated heads and flows compared with the GC simulation. Accuracy further increased locally around the valley flanks when improving the geological resolution using the refined grid. Finally, comparing MODFLOW-LGR simulation with a globally refined (GR) grid showed that the refinement proportion of the model should not exceed 10% to 15% in order to secure method efficiency.     

Three-dimensional local grid refinement for block-centered finite-difference groundwater models using iteratively coupled shared nodes: a new method of interpolation and analysis of errors

This paper describes work that extends to three dimensions the two-dimensional local-grid refinement method for block-centered finite-difference groundwater models of Mehl and Hill [Development and evaluation of a local grid refinement method for block-centered finite-difference groundwater models using shared nodes. Adv Water Resour 2002;25(5):497–511]. In this approach, the (parent) finite-difference grid is discretized more finely within a (child) sub-region. The grid refinement method sequentially solves each grid and uses specified flux (parent) and specified head (child) boundary conditions to couple the grids. Iteration achieves convergence between heads and fluxes of both grids. Of most concern is how to interpolate heads onto the boundary of the child grid such that the physics of the parent-grid flow is retained in three dimensions. We develop a new two-step, “cage-shell” interpolation method based on the solution of the flow equation on the boundary of the child between nodes shared with the parent grid. Error analysis using a test case indicates that the shared-node local grid refinement method with cage-shell boundary head interpolation is accurate and robust, and the resulting code is used to investigate three-dimensional local grid refinement of stream-aquifer interactions. Results reveal that (1) the parent and child grids interact to shift the true head and flux solution to a different solution where the heads and fluxes of both grids are in equilibrium, (2) the locally refined model provided a solution for both heads and fluxes in the region of the refinement that was more accurate than a model without refinement only if iterations are performed so that both heads and fluxes are in equilibrium, and (3) the accuracy of the coupling is limited by the parent-grid size—a coarse parent grid limits correct representation of the hydraulics in the feedback from the child grid.     

Fitting a stochastic partial differential equation to aquifer data

The steady state two dimensional groundwater flow equation with constant transmissivities was studied by Whittle in 1954 as a stochastic Laplace equation. He showed that the correlation function consisted of a modified Bessel function of the second kind, order 1, multiplied by its argument. This paper uses this pioneering work of Whittle to fit an aquifer head field to unequally spaced observations by maximum likelihood. Observational error is also included in the model. Both the isotropic and anisotropic cases are considered. The fitted field is then calculated on a two dimensional grid together with its standard deviation. The method is closely related to the use of two-dimensional splines for fitting surfaces to irregularly spaced observations.     

Finite element numerical simulation of 2.5D direct current method based on mesh refinement and recoarsement

To deal with the problem of low computational precision at the nodes near the source and satisfy the requirements for computational efficiency in inversion imaging and finite-element numerical simulations of the direct current method, we propose a new mesh refinement and recoarsement method for a two-dimensional point source. We introduce the mesh refinement and mesh recoarsement into the traditional structured mesh subdivision. By refining the horizontal grids, the singularity owing to the point source is minimized and the topography is simulated. By recoarsening the horizontal grids, the number of grid cells is reduced significantly and computational efficiency is improved. Model tests show that the proposed method solves the singularity problem and reduces the number of grid cells by 80% compared to the uniform grid refinement.     

Grid cell distortion and MODFLOW's integrated finite-difference numerical solution

The ground water flow model MODFLOW inherently implements a nongeneralized integrated finite-difference (IFD) numerical scheme. The IFD numerical scheme allows for construction of finite-difference model grids with curvilinear (piecewise linear) rows. The resulting grid comprises model cells in the shape of trapezoids and is distorted in comparison to a traditional MODFLOW finite-difference grid. A version of MODFLOW-88 (herein referred to as MODFLOW IFD) with the code adapted to make the one-dimensional DELR and DELC arrays two dimensional, so that equivalent conductance between distorted grid cells can be calculated, is described. MODFLOW IFD is used to inspect the sensitivity of the numerical head and velocity solutions to the level of distortion in trapezoidal grid cells within a converging radial flow domain. A test problem designed for the analysis implements a grid oriented such that flow is parallel to columns with converging widths. The sensitivity analysis demonstrates MODFLOW IFD's capacity to numerically derive a head solution and resulting intercell volumetric flow when the internal calculation of equivalent conductance accounts for the distortion of the grid cells. The sensitivity of the velocity solution to grid cell distortion indicates criteria for distorted grid design. In the radial flow test problem described, the numerical head solution is not sensitive to grid cell distortion. The accuracy of the velocity solution is sensitive to cell distortion with error <1% if the angle between the nonparallel sides of trapezoidal cells is <12.5 degrees. The error of the velocity solution is related to the degree to which the spatial discretization of a curve is approximated with piecewise linear segments. Curvilinear finite-difference grid construction adds versatility to spatial discretization of the flow domain. MODFLOW-88's inherent IFD numerical scheme and the test problem results imply that more recent versions of MODFLOW 2000, with minor modifications, have the potential to make use of a curvilinear grid.     

An assembly model for simulation of large-scale ground water flow and transport

When managing large-scale ground water contamination problems, it is often necessary to model flow and transport using finely discretized domains--for instance (1) to simulate flow and transport near a contamination source area or in the area where a remediation technology is being implemented; (2) to account for small-scale heterogeneities; (3) to represent ground water-surface water interactions; or (4) some combination of these scenarios. A model with a large domain and fine-grid resolution will need extensive computing resources. In this work, a domain decomposition-based assembly model implemented in a parallel computing environment is developed, which will allow efficient simulation of large-scale ground water flow and transport problems using domain-wide grid refinement. The method employs common ground water flow (MODFLOW) and transport (RT3D) simulators, enabling the solution of almost all commonly encountered ground water flow and transport problems. The basic approach partitions a large model domain into any number of subdomains. Parallel processors are used to solve the model equations within each subdomain. Schwarz iteration is applied to match the flow solution at the subdomain boundaries. For the transport model, an extended numerical array is implemented to permit the exchange of dispersive and advective flux information across subdomain boundaries. The model is verified using a conventional single-domain model. Model simulations demonstrate that the proposed model operated in a parallel computing environment can result in considerable savings in computer run times (between 50% and 80%) compared with conventional modeling approaches and may be used to simulate grid discretizations that were formerly intractable.     

Simulating a lake as a high-conductivity variably saturated porous medium

One approach for simulating ground water–lake interactions is to incorporate the lake into the ground water solution domain as a high-conductivity region. Previous studies have developed this approach using fully saturated models. This study extends this approach to variably saturated models, so that ground water–lake interactions may be more easily simulated with commonly used or public domain variably saturated codes that do not explicitly support coupled lake–water balance modeling. General guidelines are developed for the choices of saturated hydraulic conductivity and moisture retention and relative permeability curves for the lake region. When applied to an example ground water–lake system, model results are very similar to those from a model in which the lake is represented as a specified head boundary continuously updated by a lake mass balance. The high-conductivity region approach is most suitable for relatively simple geometries and lakes with slower and smaller fluctuations when the overall flow pattern and system fluxes, rather than the detailed flow pattern around the intersection of the lake and land surfaces, are of interest.     

瞬变电磁三维FDTD正演多分辨网格方法

瞬变电磁三维时域有限差分（FDTD）正演的网格剖分受最小网格尺寸、时间步长、边界条件、目标尺寸、模型尺寸等的影响,结构化网格一直存在最小网格尺寸受限于异常目标尺寸的矛盾;尽管非均匀网格能够在保证模型尺寸的前提下尽可能的降低网格数量,但由于Yee网格结构的限制,非均匀网格不能无限制的扩大单一方向的尺寸,这是为了避免边界网格区域出现长宽比过大的畸形网格,影响计算精度甚至导致结果发散.在非均匀网格剖分的基础上,本文提出了瞬变电磁三维FDTD正演的多尺度网格方法,即首先使用较大尺寸的粗网格进行第一次剖分,然后在希望加密的区域进行二次剖分,使计算域中包含粗、细两套网格.尽管细网格包含在粗网格内部,但其具有Yee网格的全部属性,因而可以在网格中设置不同的电性参数模拟不同形状的目标.基于Maxwell方程组推导了细网格内电场和磁场的迭代公式,基于泰勒展开给出了设置粗、细网格后产生的内部边界条件,使电磁场的传播在粗、细网格和时间步进上得到统一.采用均匀半空间中包含三维低阻异常的经典模型和三维接触带复杂模型进行精度验证,发现多分辨网格方法计算结果满足精度要求.使用"L"型异常模型计算采用多分辨网格方法和不采用多分辨网格的传统FDTD方法对比计算效率,发现多分辨网格算法能够显著提高计算效率,并能够保证计算精度.     

Sensitivity and moment analyses of head in variably saturated regimes

A numerical approach for approximating statistical moments of hydraulic heads of variably saturated flows in multi-dimensional porous media is developed. The approximation relies on a first-order Taylor series expansion of a finite element flow model and an adjoint state numerical method for variably saturated flows to evaluate sensitivities. This approach can be employed to analyze uncertainties associated with predictions of head of steady-state or transient flows in variably saturated porous media, with any type of boundary and initial conditions. Limitations of stochastic analytical methods such as spectral/perturbation approaches and the time-consuming Monte Carlo simulation technique are thus alleviated. An example is given to demonstrate the utility of the approach and to investigate the temporal evolution of head variances in a variably saturated flow regime. Results show that the fluctuation of the water table can have significant impacts on the propagation of the head variance.     

Factors Influencing the Stream‐Aquifer Flow Exchange Coefficient

Knowledge of river gain from or loss to a hydraulically connected water table aquifer is crucial in issues of water rights and also when attempting to optimize conjunctive use of surface and ground waters. Typically in groundwater models this exchange flow is related to a difference in head between the river and some point in the aquifer, through a “coefficient.” This coefficient has been defined differently as well as the location for the head in the aquifer. This paper proposes a new coefficient, analytically derived, and a specific location for the point where the aquifer head is used in the difference. The dimensionless part of the coefficient is referred to as the SAFE (stream‐aquifer flow exchange) dimensionless conductance. The paper investigates the factors that influence the value of this new conductance. Among these factors are (1) the wetted perimeter of the cross‐section, (2) the degree of penetration of the cross‐section, and (3) the shape of the cross‐section. The study shows that these factors just listed are indeed ordered in their respective level of importance. In addition the study verifies that the analytical correct value of the coefficient is matched by finite difference simulation only if the grid system is sufficiently fine. Thus the use of the analytical value of the coefficient is an accurate and efficient alternative to ad hoc estimates for the coefficient typically used in finite difference and finite element methods.     

深度均匀采样梯形网格有限差分地震波场模拟方法

由于重力引起的岩石压实效应,一般来说,地震波传播速度由浅入深整体逐渐增大.梯形坐标系设计可耦合速度由浅入深逐渐增大的变化,该坐标系中均匀网格采样所对应的物理直角坐标系网格由浅入深逐渐增大,也即浅部低速区对应细网格,深部高速区对应粗网格.在梯形坐标系表征波动方程后利用有限差分求解,本文实现一种深度均匀采样、横向采样间隔随深度增加逐渐线性增大的有限差分地震波模拟方法.梯形坐标系波动方程离散后,仍采用常规均匀网格有限差分算法对其求解.由于横向网格大小由浅入深线性增加,本方法可避免不同大小网格区域过渡所产生的虚假反射.梯形坐标系波场模拟浅层精度高,深层横向响应范围广,可有效减少有限差分网格数量.本文提出的方法是在更广义的坐标系下利用有限差分求解波动方程,正交坐标系仅为该梯形坐标系之特例.本文旨在为大速度动态范围深地高效高精度地震波场模拟提供一种思路.     

Grid refinement for modeling multiphase flow in discretely fractured porous media

A study of the effects of grid discretization on the migration of DNAPL within a discrete-fracture network embedded in a porous rock matrix is presented. It is shown that an insufficiently fine discretization of the fracture elements can lead to an overprediction of the volume of DNAPL that continues to migrate vertically at the intersection of a vertical and horizontal fracture. Uniform discretization of elements at the scale of one centimetre (or less) accurately resolved the density and capillary pressure components of the head gradient in the DNAPL. An alternative, non-uniform method of discretization of elements within the discrete-fracture network is presented whereby only fracture elements immediately adjacent to fracture intersections are refined. To further limit the number of elements employed, the porous matrix elements adjacent to the fracture elements are not similarly refined. Results show this alternative method of discretization reduces the numerical error to an acceptable level, while allowing the simulation of field-scale DNAPL contamination problems. The results from two field-scale simulations of a DNAPL-contaminated carbonate bedrock site in Ontario, Canada are presented. These simulations compare different methods of grid discretization, and highlight the importance of grid refinement when simulating DNAPL migration problems in fractured porous media.     

Grid Refinement in Cartesian Coordinates for Groundwater Flow Models Using the Divergence Theorem and Taylor's Series

Grid refinement is introduced in a numerical groundwater model to increase the accuracy of the solution over local areas without compromising the run time of the model. Numerical methods developed for grid refinement suffered certain drawbacks, for example, deficiencies in the implemented interpolation technique; the non‐reciprocity in head calculations or flow calculations; lack of accuracy resulting from high truncation errors, and numerical problems resulting from the construction of elongated meshes. A refinement scheme based on the divergence theorem and Taylor's expansions is presented in this article. This scheme is based on the work of De Marsily (1986) but includes more terms of the Taylor's series to improve the numerical solution. In this scheme, flow reciprocity is maintained and high order of refinement was achievable. The new numerical method is applied to simulate groundwater flows in homogeneous and heterogeneous confined aquifers. It produced results with acceptable degrees of accuracy. This method shows the potential for its application to solving groundwater heads over nested meshes with irregular shapes.     

Comparison of finite difference and element models of internal tides on the Malin–Hebrides shelf

A three-dimensional baroclinic finite element model with a coarse and fine (i.e. local refinement along the shelf edge) grid is used to examine the influence of shelf edge grid refinement upon the internal tide generation and propagation off the west coast of Scotland. Comparisons are made with observations in the region and with a published solution using a finite difference model. The calculations show that provided that the finite element grid is refined in the internal tide generation area and the adjacent region through which the internal tide propagates, then a numerically accurate solution is obtained. In the regions of strong internal tide generation with a local grid refinement, internal wave energy can accumulate at small scales and must be removed by a scale-selective filter.     

An Analytical Method to Determine Ground Water Supply Well Network Designs

An analytical method is provided where the ground water practitioner can quickly determine the size (number of wells) and spacing of a well network capable of meeting a known ground water demand. In order to apply the method, two new parameters are derived that relate theoretical drawdown to the maximum drawdown that is achievable without mining the aquifer. The size of a well network is shown to be proportional to the ground water demand and inversely proportional to the transmissivity and available head. The spacing between wells in a supply well network is shown to be most sensitive to a derived parameter r _{HA/ 3}, which is related to the available head and the propagation of drawdown away from a theoretical well if the total ground water demand was applied to that well. The method can be used to quickly determine the required spacing between wells in well networks of various sizes that are completed in confined aquifers with no leakance.     

Simulation of Sampling and Hydraulic Tests to Assess a Hybrid Monitoring Well Design

Determination of the nature, extent, and rate of off-site chemical migration are common objectives of hazardous waste site investigations. Chemical analyses of water samples from monitoring wells and measurements of hydraulic head and hydraulic conductivity provide the basis for making these determinations. Accurate site assessment, therefore, depends upon the appropriate monitoring well design and sampling and testing procedures.
During the course of remedial investigations in Niagara Falls, New York, it has been necessary to evaluate the ground water quality and hydraulic characteristics of 5- to 30-feet thick overburden formations. Many of the monitoring wells completed to these formations consist of a partially penetrating screen (5 feet at the base of the formation) with a fully penetrating sandpack. Questions regarding how this well design influences the source of sampled ground water and hydraulic tests were examined using an extremely fine axisymmetric grid with SATURN, a two-dimensional, finite-element ground water model, and a particle tracking post-processor.
A discrete sensitivity analysis was made to determine how flow patterns induced by pumping at 1 gpm are affected by: different screen and sandpack configurations, the ratio of sandpack to formation hydraulic conductivities, heterogeneity, anisotropy, and sandpack thickness. The simulations show that the source (and chemistry given a non-uniform chemical distribution) of ground water sampled will vary considerably depending on a number of factors. Analysis of simulated drawdowns in the monitoring well during purging shows that calculated transmissivities for the range of well designs and conditions modeled will be accurate to within one-half order of magnitude.     

Delineating and quantifying ground water discharge zones using streambed temperatures

Streambed temperature mapping, hydraulic testing using minipiezometers, and geochemical analyses of interstitial water of the streambed were used to delineate the pattern of ground water discharge in a sandy streambed and to develop a flux-based conceptual model for ground water/surface water interactions. A new and simple empirical method was used to relate fluxes obtained from minipiezometer data to streambed temperatures. The relationship allowed flux to be calculated at locations where only streambed temperature measurements were made. Slug testing and potentiomanometer measurements at 34 piezometers indicated ground water discharge ranged from 0.03 to 446 L/m2/day (and possibly as high as 7060 L/m2/day) along a 60 m long by 11 to 14 m wide reach of river. Complex but similar plan-view patterns of flux were calculated for both summer and winter using hundreds of streambed temperatures measured on a 1 by 2 m grid. The reach was dominated by ground water discharge and 5% to 7% of the area accounted for approximately 20% to 24% of the total discharge. < 12% of the total area consisted of recharge zones or no-discharge zones. A conceptual model for ground water/surface water interactions consisting of five different behaviors was developed based on the magnitude and direction of flux across the surface of the streambed. The behaviors include short-circuit discharge (e.g., high-flow springs), high discharge (e.g., preferential flowpaths), low to moderate discharge, no discharge (e.g., horizontal hyporheic or ground water flow), and recharge. Geological variations at depth played a key role in determining which type of flow behavior occurred in the streambed.     

Effect of elevation resolution on evapotranspiration simulations using MODFLOW

Surface elevations represented in MODFLOW head-dependent packages are usually derived from digital elevation models (DEMs) that are available at much high resolution. Conventional grid refinement techniques to simulate the model at DEM resolution increases computational time, input file size, and in many cases are not feasible for regional applications. This research aims at utilizing the increasingly available high resolution DEMs for effective simulation of evapotranspiration (ET) in MODFLOW as an alternative to grid refinement techniques. The source code of the evapotranspiration package is modified by considering for a fixed MODFLOW grid resolution and for different DEM resolutions, the effect of variability in elevation data on ET estimates. Piezometric head at each DEM cell location is corrected by considering the gradient along row and column directions. Applicability of the research is tested for the lower Rio Grande (LRG) Basin in southern New Mexico. The DEM at 10 m resolution is aggregated to resampled DEM grid resolutions which are integer multiples of MODFLOW grid resolution. Cumulative outflows and ET rates are compared at different coarse resolution grids. Results of the analysis conclude that variability in depth-to-groundwater within the MODFLOW cell is a major contributing parameter to ET outflows in shallow groundwater regions. DEM aggregation methods for the LRG Basin have resulted in decreased volumetric outflow due to the formation of a smoothing error, which lowered the position of water table to a level below the extinction depth.     

Configuring Detection Wells Near Landfills

A graphical method was devised for designing contaminant detection monitoring networks in aquifers. The approach eliminates bias in detection efficiency among well pairs, thereby improving the overall efficiency of a ground water monitoring network. In the equidistant configurations derived by the graphical approach, all wells are located the same distance from a landfill, but the distance is measured parallel to ground water flow, Measured perpendicular to ground water flow, there is also an equal spacing between wells in an equidistant network. A simulation model was used to compare an equidistant network to a peripheral monitoring configuration, in which wells were spaced evenly along the downgradient boundaries of a landfill. The equidistant network yielded a 12.4% higher detection efficiency and also facilitated earlier release detection. In practice, the graphical approach that yields equidistant configurations can be used to identify candidate monitoring networks to detect potential releases from landfills.