A terrain-following grid transform and preconditioner for parallel,large-scale,integrated hydrologic modeling

A terrain-following grid formulation (TFG) is presented for simulation of coupled variably-saturated subsurface and surface water flow. The TFG is introduced into the integrated hydrologic model, ParFlow, which uses an implicit, Newton Krylov solution technique. The analytical Jacobian is also formulated and presented and both the diagonal and non-symmetric terms are used to precondition the Krylov linear system. The new formulation is verified against an orthogonal stencil and is shown to provide increased accuracy at lower lateral spatial discretization for hillslope simulations. Using TFG, efficient scaling to a large number of processors (16,384) and a large domain size (8.1 Billion unknowns) is shown. This demonstrates the applicability of this formulation to high-resolution, large-spatial extent hydrology applications where topographic effects are important. Furthermore, cases where the analytical Jacobian is used for the Newton iteration and as a non-symmetric preconditioner for the linear system are shown to have faster computation times and better scaling. This demonstrates the importance of solver efficiency in parallel scaling through the use of an appropriate preconditioner.     

New preconditioning strategy for Jacobian-free solvers for variably saturated flows with Richards’ equation

We develop a new approach for solving the nonlinear Richards’ equation arising in variably saturated flow modeling. The growing complexity of geometric models for simulation of subsurface flows leads to the necessity of using unstructured meshes and advanced discretization methods. Typically, a numerical solution is obtained by first discretizing PDEs and then solving the resulting system of nonlinear discrete equations with a Newton-Raphson-type method. Efficiency and robustness of the existing solvers rely on many factors, including an empiric quality control of intermediate iterates, complexity of the employed discretization method and a customized preconditioner. We propose and analyze a new preconditioning strategy that is based on a stable discretization of the continuum Jacobian. We will show with numerical experiments for challenging problems in subsurface hydrology that this new preconditioner improves convergence of the existing Jacobian-free solvers 3-20 times. We also show that the Picard method with this preconditioner becomes a more efficient nonlinear solver than a few widely used Jacobian-free solvers.     

Improvement of surface run‐off in the hydrological model ParFlow by a scale‐consistent river parameterization

We propose an improvement of the overland‐flow parameterization in a distributed hydrological model, which uses a constant horizontal grid resolution and employs the kinematic wave approximation for both hillslope and river channel flow. The standard parameterization lacks any channel flow characteristics for rivers, which results in reduced river flow velocities for streams narrower than the horizontal grid resolution. Moreover, the surface areas, through which these wider model rivers may exchange water with the subsurface, are larger than the real river channels potentially leading to unrealistic vertical flows. We propose an approximation of the subscale channel flow by scaling Manning's roughness in the kinematic wave formulation via a relationship between river width and grid cell size, following a simplified version of the Barré de Saint‐Venant equations (Manning–Strickler equations). The too large exchange areas between model rivers and the subsurface are compensated by a grid resolution‐dependent scaling of the infiltration/exfiltration rate across river beds. We test both scaling approaches in the integrated hydrological model ParFlow. An empirical relation is used for estimating the true river width from the mean annual discharge. Our simulations show that the scaling of the roughness coefficient and the hydraulic conductivity effectively corrects overland flow velocities calculated on the coarse grid leading to a better representation of flood waves in the river channels.     

Application of Jacobian-free Newton–Krylov with physics-based preconditioning to biogeochemical transport

Modern multicomponent geochemical transport models require the use of parallel computation for carrying out three-dimensional, field-scale simulations due to extreme memory and processing demands. However, to fully exploit the advanced computational power provided by today’s supercomputers, innovative parallel algorithms are needed. We demonstrate the use of Jacobian-free Newton–Krylov (JFNK) within the Newton–Raphson method to reduce memory and processing requirements on high-performance computers. We also demonstrate the use of physics-based preconditioners, which are often necessary when using JFNK since no explicit Jacobian matrix is ever formed. We apply JFNK to simulate enhanced in situ bioremediation of a NAPL source zone, which entails highly coupled geochemical and biodegradation reactions. The algorithm’s performance is evaluated and compared with conventional solvers and preconditioners. We found that JFNK provided substantial saving in memory (i.e. 30–60%) on problems utilizing up to 512 processors on LANL’s ASCI Q. However, the performance based on wallclock time was less advantageous, coming out on par with conventional techniques. In addition, we illustrate deficiencies in physics-based preconditioner performance for biogeochemical transport problems with components that undergo significant sorption or form a local quasi-stationary state.     

基于高斯牛顿法的频率域可控源电磁三维反演研究

三维反演解释是电磁法勘探发展的重要趋势,而如何提高三维反演的可靠性、稳定性和计算效率是算法开发者们目前的研究重点.本文实现了一种频率域可控源电磁(CSEM)三维反演算法.其中正演基于拟态有限体积法离散化,利用直接矩阵分解技术来求解大型线性系统方程,不仅准确、稳定,而且特别有利于含有大量发射场源位置的CSEM勘探情况;对目标函数的最优化采用高斯牛顿法(GN),具有近似二次的收敛性;使用预条件共轭梯度法(PCG)求解每次GN迭代所得到的法方程,避免了显式求解和存储灵敏度矩阵,减小了计算量.以上这些方法的结合应用,使得本文的三维反演算法准确、稳定且高效.通过陆地和海洋CSEM勘探场景中的典型理论模型的反演测试,验证了本文算法的有效性.     

On Solving Groundwater Flow and Transport Models with Algebraic Multigrid Preconditioning

Iterative solvers preconditioned with algebraic multigrid have been devised as an optimal technology to speed up the response of large sparse linear systems. In this work, this technique was implemented in the framework of the dual delineation approach. This involves a single groundwater flow linear solution and a pure advective transport solution with different right-hand sides. The new solver was compared with other preconditioned iterative methods, the MODFLOW's GMG solver, and direct sparse solvers. Test problems include two- and three-dimensional benchmarks spanning homogeneous and highly heterogeneous and anisotropic formations. For the groundwater flow problems, using the algebraic multigrid preconditioning speeds up the numerical solution by one to two orders of magnitude. The algebraic multigrid preconditioner efficiency was preserved for the three dimensional heterogeneous and anisotropic problem unlike for the MODFLOW's GMG solver. Contrarily, a sparse direct solver was the most efficient for the pure advective transport processes such as the forward travel time simulations. Hence, the best sparse solver for the more general advection-dispersion transport equation is likely to be Péclet number dependent. When equipped with the best solvers, processing multimillion grid blocks by the dual delineation approach is a matter of seconds. This paves the way for its routine application to large geological models. The paper gives practical hints on the strategies and conditions under which algebraic multigrid preconditioning would remain competitive for the class of nonlinear and/or transient problems.     

Efficient scheme for the shallow water equations on unstructured grids with application to the Continental Shelf

In this paper, a shallow-water flow solver is presented, based on the finite-volume method on unstructured grids The method is suitable for flows that occur in rivers, channels, sewer systems (1D), shallow seas, rivers, overland flow (2D), and estuaries, lakes and shelf breaks (3D). We present an outline of the numerical approach and show three 2D test cases and an application of tidal propagation on the Continental Shelf. The benefits of applying an unstructured grid were explored by creating an efficient model network that aims at keeping the number of grid cells per wavelength constant. The computational speed of our method was compared with that of WAQUA/TRIWAQ and Delft3D (the commonly used structured shallow-flow solvers in The Netherlands), and comparable performance was found.     

A parallel,fully coupled,fully implicit solution to reactive transport in porous media using the preconditioned Jacobian-Free Newton-Krylov Method

Modeling large multicomponent reactive transport systems in porous media is particularly challenging when the governing partial differential algebraic equations (PDAEs) are highly nonlinear and tightly coupled due to complex nonlinear reactions and strong solution-media interactions. Here we present a preconditioned Jacobian-Free Newton-Krylov (JFNK) solution approach to solve the governing PDAEs in a fully coupled and fully implicit manner. A well-known advantage of the JFNK method is that it does not require explicitly computing and storing the Jacobian matrix during Newton nonlinear iterations. Our approach further enhances the JFNK method by utilizing physics-based, block preconditioning and a multigrid algorithm for efficient inversion of the preconditioner. This preconditioning strategy accounts for self- and optionally, cross-coupling between primary variables using diagonal and off-diagonal blocks of an approximate Jacobian, respectively. Numerical results are presented demonstrating the efficiency and massive scalability of the solution strategy for reactive transport problems involving strong solution-mineral interactions and fast kinetics. We found that the physics-based, block preconditioner significantly decreases the number of linear iterations, directly reducing computational cost; and the strongly scalable algebraic multigrid algorithm for approximate inversion of the preconditioner leads to excellent parallel scaling performance.     

Parametric study of a physically-based,plot-scale hillslope hydrological model through virtual experiments

Abstract

A physically-based hillslope hydrological model with shallow overland flow and rapid subsurface stormflow components was developed and calibrated using field experiments conducted on a preferential path nested hillslope in northeast India. Virtual experiments were carried out to perform sensitivity analysis of the model using the automated parameter estimation (PEST) algorithm. Different physical parameters of the model were varied to study the resulting effects on overland flow and subsurface stormflow responses from the theoretical hillslopes. It was observed that topographical shapes had significant effects on overland flow hydrographs. The slope profiles, surface storage, relief, rainfall intensity and infiltration rates primarily controlled the overland flow response of the hillslopes. Prompt subsurface stormflow responses were mainly dominated by lateral preferential flow, as soil matrix flow rates were very slow. Rainfall intensity and soil macropore structures were the most influential parameters on subsurface stormflow. The number of connected soil macropores was a more sensitive parameter than the size of macropores. In hillslopes with highly active vertical and lateral preferential pathways, saturation excess overland flow was not evident. However, saturation excess overland flow was generated if the lateral macropores were disconnected. Under such conditions, rainfall intensity, duration and preferential flow rate governed the process of saturation excess overland flow generation from hillslopes.

Editor D. Koutsoyiannis; Associate editor C. Perrin     

A finite element model for simulating surface run‐off and unsaturated seepage flow in the shallow subsurface

This work introduces water–air two‐phase flow into integrated surface–subsurface flow by simulating rainfall infiltration and run‐off production on a soil slope with the finite element method. The numerical model is formulated by partial differential equations for hydrostatic shallow flow and water–air two‐phase flow in the shallow subsurface. Finite element computing formats and solution strategies are presented to obtain a numerical solution for the coupled model. An unsaturated seepage flow process is first simulated by water–air two‐phase flow under the atmospheric pressure boundary condition to obtain the rainfall infiltration rate. Then, the rainfall infiltration rate is used as an input parameter to solve the surface run‐off equations and determine the value of the surface run‐off depth. In the next iteration, the pressure boundary condition of unsaturated seepage flow is adjusted by the surface run‐off depth. The coupling process is achieved by updating the rainfall infiltration rate and surface run‐off depth sequentially until the convergence criteria are reached in a time step. A well‐conducted surface run‐off experiment and traditional surface–subsurface model are used to validate the new model. Comparisons with the traditional surface–subsurface model show that the initiation time of surface run‐off calculated by the proposed model is earlier and that the water depth is larger, thus providing values that are closer to the experimental results.     

Storm runoff processes and subcatchment characteristics in a new zealand hill country catchment

The quickflow responses of six subcatchment areas in a small hill country catchment in the Craigieburn Range, South Island, New Zealand, were compared for a range of storm sizes, rainfall intensities and antecedent wetness conditions. Topography and soil characteristics suggested that all subcatchments would receive subsurface stormflow input, but that some would receive larger saturation overland flow inputs than others. Quickflow yields and response ratios were positively correlated with storm size and antecedent wetness conditions in the subcatchment most suited to producing saturation overland flow. In subcatchments more likely to be dominated by subsurface flow, quickflow yields and response ratios were positively correlated with storm size, but were either not correlated, or negatively correlated, with antecedent wetness. Quickflow responses were either not significantly or negatively correlated with rainfall intensity variables. Quickflow from the subcatchment most suited to produce saturation overland flow providing an increasing proportion of total catchment quickflow in larger storms and as antecedent conditions became wetter. Subcatchment responses varied greatly in space and time and there was less pattern to the variation than had been expected. Where topographic and pedologic conditions permit substantial responses to storm rainfall by both saturation overland flow and subsurface stormflow, simple topographic and soil indicators may not be useful guides to the relative importance of runoff mechanisms, or to the identification of runoff-source areas.     

HYDROLOGICALLY DRIVEN MECHANISMS OF HEADCUT DEVELOPMENT

An analytic investigation of the effect of surface seal mechanical properties, overland flow, and subsurface hydrology was performed on headcut development. Headcut growth rates on upland areas have been observed to be quite small （less than 0.00015 meter per second） and that they occur in increments in which chips break off at points where cracks have developed in surface seals. The substrate soil under the seal collapses and is removed by the flow. This mode of headcut development is the result of a strong interaction between the surface and the subsurface processes. The surface process is energetically controlled by the mechanical features of the seal whereas the subsurface process is hydrologically controlled. The analysis yields estimates of the temporal scale of headcut velocities. In cases of infiltration from the vertical gully wall into the substrate, the flexural wave velocity （seismic sound velocity） was found to inversely affect headcut velocity.     

Efficient simulation of geothermal processes in heterogeneous porous media based on the exponential Rosenbrock–Euler and Rosenbrock-type methods

Simulation of geothermal systems is challenging due to coupled physical processes in highly heterogeneous media. Combining the exponential Rosenbrock–Euler method and Rosenbrock-type methods with control-volume (two-point flux approximation) space discretizations leads to efficient numerical techniques for simulating geothermal systems. In terms of efficiency and accuracy, the exponential Rosenbrock–Euler time integrator has advantages over standard time-discretization schemes, which suffer from time-step restrictions or excessive numerical diffusion when advection processes are dominating. Based on linearization of the equation at each time step, we make use of matrix exponentials of the Jacobian from the spatial discretization, which provide the exact solution in time for the linearized equations. This is at the expense of computing the matrix exponentials of the stiff Jacobian matrix, together with propagating a linearized system. However, using a Krylov subspace or Léja points techniques make these computations efficient.The Rosenbrock-type methods use the appropriate rational functions of the Jacobian of the ODEs resulting from the spatial discretization. The parameters in these schemes are found in consistency with the required order of convergence in time. As a result, these schemes are A-stable and only a few linear systems are solved at each time step. The efficiency of the methods compared to standard time-discretization techniques are demonstrated in numerical examples.     

A sensitivity analysis of Hortonian flow

We present a sensitivity analysis for infiltration excess (Hortonian) overland flow based on a classic laboratory experiment by Smith and Woolhiser [Smith RE, Woolhiser DA. Overland flow on an infiltrating surface. Water Resour Res 1971;7(4):899–913]. The model components of the compartment approach are comprised of a diffusive wave approximation to the Saint–Venant equations for overland flow, a Richards model for flow in the variably saturated zone, and an interface coupling concept that combines the two components. In the coupling scheme a hydraulic interface is introduced to allow the definition of an exchange flux between the surface and the unsaturated zone. The effects of friction processes, soil capillarity, hydraulic interface, and vertical soil discretization on both infiltration and runoff prediction are investigated in detail. The corresponding sensitivity analysis is conducted using a small-perturbation method. As a result the importance of the hydraulic processes and related parameters are evaluated for the coupled hydrosystem.     

A distributed TOPMODEL for modelling impacts of land‐cover change on river flow in upland peatland catchments

There is global concern about headwater management and associated impacts on river flow. In many wet temperate zones peatlands can be found covering headwater catchments. In the UK there is major concern about how environmental change, driven by human interventions, has altered the surface cover of headwater blanket peatlands. However, the impact of such land‐cover changes on river flow is poorly understood. In particular, there is poor understanding of the impacts of different spatial configurations of bare peat or well‐vegetated, restored peat on river flow peaks in upland catchments. In this paper, a physically based, distributed and continuous catchment hydrological model was developed to explore such impacts. The original TOPMODEL, with its process representation being suitable for blanket peat catchments, was utilized as a prototype acting as the basis for the new model. The equations were downscaled from the catchment level to the cell level. The runoff produced by each cell is divided into subsurface flow and saturation‐excess overland flow before an overland flow calculation takes place. A new overland flow module with a set of detailed stochastic algorithms representing overland flow routing and re‐infiltration mechanisms was created to simulate saturation‐excess overland flow movement. The new model was tested in the Trout Beck catchment of the North Pennines of England and found to work well in this catchment. The influence of land cover on surface roughness could be explicitly represented in the model and the model was found to be sensitive to land cover. Copyright © 2014 John Wiley & Sons, Ltd.     

大地电磁有限差分数值解对比

基于网格精度、形成系统方程的方式、边界条件以及预条件线性算子等,文中对大地电磁（MT）有限差分数值解作了对比.对不同网格剖分方式下的三个均匀半空间模型的一维MT响应对比显示,在降低首层厚度的同时保持层间厚度变化在合理范围可以同时提高主场和辅助场的精度.在利用正常中心网格法（主场和辅助场都定义在单元顶面的中心）计算二维（2-D）TM模式响应时,应该从Maxwell一次差分方程开始组建二次差分方程,这样可以更充分考虑模型电阻率的变化.在对边界值如何影响数值解的测试表明,仅仅提高一维（1-D）边界值的精度对提高2-D MT有限差分数值解的精度是有限的.线性算子对提高MT解的效率十分重要,简单的对比进一步表明合适的预条件再配合好的线性算子（如文中求解2-D MT时所采用的DILU-BICGSTAB方法）不仅可以加速收敛,而且可以降低迭代次数.     

Comparison of an algebraic multigrid algorithm to two iterative solvers used for modeling ground water flow and transport

Numerical solution of large-scale ground water flow and transport problems is often constrained by the convergence behavior of the iterative solvers used to solve the resulting systems of equations. We demonstrate the ability of an algebraic multigrid algorithm (AMG) to efficiently solve the large, sparse systems of equations that result from computational models of ground water flow and transport in large and complex domains. Unlike geometric multigrid methods, this algorithm is applicable to problems in complex flow geometries, such as those encountered in pore-scale modeling of two-phase flow and transport. We integrated AMG into MODFLOW 2000 to compare two- and three-dimensional flow simulations using AMG to simulations using PCG2, a preconditioned conjugate gradient solver that uses the modified incomplete Cholesky preconditioner and is included with MODFLOW 2000. CPU times required for convergence with AMG were up to 140 times faster than those for PCG2. The cost of this increased speed was up to a nine-fold increase in required random access memory (RAM) for the three-dimensional problems and up to a four-fold increase in required RAM for the two-dimensional problems. We also compared two-dimensional numerical simulations of steady-state transport using AMG and the generalized minimum residual method with an incomplete LU-decomposition preconditioner. For these transport simulations, AMG yielded increased speeds of up to 17 times with only a 20% increase in required RAM. The ability of AMG to solve flow and transport problems in large, complex flow systems and its ready availability make it an ideal solver for use in both field-scale and pore-scale modeling.     

A hydrochemical approach to quantify the role of return flow in a surface flow‐dominated catchment

Stormflow generation in headwater catchments dominated by subsurface flow has been studied extensively, yet catchments dominated by surface flow have received less attention. We addressed this by testing whether stormflow chemistry is controlled by either (a) the event‐water signature of overland flow, or (b) the pre‐event water signature of return flow. We used a high‐resolution hydrochemical data set of stormflow and end‐members of multiple storms in an end‐member mixing analysis to determine the number of end‐members needed to explain stormflow, characterize and identify potential end‐members, calculate their contributions to stormflow, and develop a conceptual model of stormflow. The arrangement and relative positioning of end‐members in stormflow mixing space suggest that saturation excess overland flow (26–48%) and return flow from two different subsurface storage pools (17–53%) are both similarly important for stormflow. These results suggest that pipes and fractures are important flow paths to rapidly release stored water and highlight the value of within‐event resolution hydrochemical data to assess the full range and dynamics of flow paths.     

A particle-tracking scheme for simulating pathlines in coupled surface-subsurface flows

A Lagrangian particle tracking scheme has been extended to simulate advective transport through coupled surface and subsurface flows. This extended scheme assumes a two-dimensional flow field for the overland domain and a three-dimensional flow field for the subsurface domain. Moreover it is assumed that the flow fields are simulated by a cell centered finite difference method. Pathlines in both the subsurface and the overland domain are simulated by classical particle tracking methods. Exchange of particles between the two domains is simulated by newly-developed algorithms presented in this study. Different algorithms are used depending on the direction of the exchange across the interface between the two domains. In the subsurface domain knowledge about a particle’s pathline is enough to detect a transfer to the surface domain and the solution is straightforward. However, in the two-dimensional overland domain pathlines are parallel to the land surface. Therefore the velocity field in the overland domain alone cannot be used to detect a transfer to the subsurface. We propose a relatively simple algorithm to estimate the probability of transfer to the subsurface domain. It is shown that this algorithm can also be used to handle the transfer from the overland domain to the atmosphere domain. The algorithm to estimate the transfer probabilities is based on the mass balance of water on a streamtube aligned with the particle’s pathline. This newly developed technique ensures that transit time distributions can be simulated accurately. These new relationships are implemented in an existing particle tracking code and are verified using analytical solutions for transit times.     

A New Equation Solver for Modeling Turbulent Flow in Coupled Matrix‐Conduit Flow Models

Karst aquifers represent dual flow systems consisting of a highly conductive conduit system embedded in a less permeable rock matrix. Hybrid models iteratively coupling both flow systems generally consume much time, especially because of the nonlinearity of turbulent conduit flow. To reduce calculation times compared to those of existing approaches, a new iterative equation solver for the conduit system is developed based on an approximated Newton–Raphson expression and a Gauß–Seidel or successive over‐relaxation scheme with a single iteration step at the innermost level. It is implemented and tested in the research code CAVE but should be easily adaptable to similar models such as the Conduit Flow Process for MODFLOW‐2005. It substantially reduces the computational effort as demonstrated by steady‐state benchmark scenarios as well as by transient karst genesis simulations. Water balance errors are found to be acceptable in most of the test cases. However, the performance and accuracy may deteriorate under unfavorable conditions such as sudden, strong changes of the flow field at some stages of the karst genesis simulations.