169,831 research outputs found

    Picard and Newton Linearization For the Coupled Model of Saltwater Intrusion In Aquifers

    No full text
    Difficulties in the numerical solution of the partial differential equations governing seawater intrusion in aquifers arise from the coupling between the flow and transport equations and from the nonlinear aspects of this coupling. Several linearization approaches are discussed for the solution of the nonlinear system which results from a finite element discretization of the coupled equations. It is first shown that the most commonly used solution method can be viewed as a Picard linearization applied to the transport equation, with the coupling resolved by iteration over the two governing equations. The full Newton scheme for solving the coupled problem produces a Jacobian of size 2N x 2N, where N is the number of nodes in the discretization of both the flow and transport equations. To reduce the size and complexity of the full Newton scheme, a partial Newton method is proposed, which, like the Picard approach, produces matrix systems of size N x N. This scheme applies Newton linearization to the transport equation, and conventional iteration to resolve the coupling. Results from two- and three-dimensional test simulations show that the partial Newton scheme gives improved convergence and robustness compared to Picard linearization, especially for highly advective problems or large density ratios. Both approaches involve the solution of a symmetric (flow) and a nonsymmetric (transport) system of equations, and it is shown that the per iteration CPU cost for the partial Newton method is not significantly greater than that of the Picard scheme

    A Comparison of Picard and Newton Iteration In the Numerical-solution of Multidimensional Variably Saturated Flow Problems

    No full text
    Picard iteration is a widely used procedure for solving the nonlinear equation governing flow in variably saturated porous media. The method is simple to code and computationally cheap, but has been known to fail or converge slowly. The Newton method is more complex and expensive (on a per-iteration basis) than Picard, and as such has not received very much attention. Its robustness and higher rate of convergence, however, make it an attractive alternative to the Picard method, particularly for strongly nonlinear problems. In this paper the Picard and Newton schemes are implemented and compared in one-, two-, and three-dimensional finite element simulations involving both steady state and transient flow. The eight test cases presented highlight different aspects of the performance of the two iterative methods and the different factors that can affect their convergence and efficiency, including problem size, spatial and temporal discretization, initial solution estimates, convergence error norm, mass lumping, time weighting, conductivity and moisture content characteristics, boundary;conditions, seepage fades, and the extent of fully saturated zones in the soil. Previous strategies for enhancing the performance of the Picard and Newton schemes are revisited, and new ones are suggested. The strategies include chord slope approximations for the derivatives of the characteristic equations, relaxing convergence requirements along seepage faces, dynamic time step control, nonlinear relaxation, and a mixed Picard-Newton approach. The tests show that the Picard or relaxed Picard schemes are often adequate for solving Richards' equation, but that in cases where these fail to converge or converge slowly, the Newton method should be used. The mixed Picard-Newton approach can effectively overcome the Newton scheme's sensitivity to initial solution estimates, while comparatively poor performance is reported for the various chord slope approximations. Finally,given the reliability and efficiency of current conjugate gradient-like methods for solving linear nonsymmetric systems; the only real drawback of using Newton rather than Picard iteration is the algebraic complexity and computational cost of assembling the derivative terms of the Jacobian matrix, and it is suggested that both methods can be effectively implemented and used in numerical models of Richards' equation

    Assessment of adaptive vs. heuristic time stepping for variably saturated flow

    No full text
    The performance of improved initial estimates and 'heuristic' and 'adaptive' techniques for time step control in the iterative solution of Richards equation is evaluated. The so-called heuristic technique uses the convergence behaviour of the iterative scheme to estimate the next time step whereas the adaptive technique regulates the time step on the basis of an approximation of the local time truncation error. The sample problems used to assess these various schemes are characterized by nonuniform (in time) boundary conditions, sharp gradients in the infiltration fronts, and discontinuous derivatives in the soil hydraulic properties. It is found that higher order initial solution estimates improve the convergence of the iterative scheme for both the heuristic and adaptive techniques, with greater overall performance gains for the heuristic scheme, as could be expected. It is also found that the heuristic technique outperforrns the adaptive method under strongly nonlinear conditions. Previously reported observations suggesting that adaptive techniques perform best when accuracy requirements on the numerical solution are very stringent are confirmed. Overall both heuristic and adaptive techniques have their limitations, and a more general or mixed time stepping strategy combining truncation error and convergence criteria is recommended for complex problem

    Mass-conservative reconstruction of Galerkin velocity fields for transport simulations

    No full text
    Accurate calculation of mass-conservative velocity fields from numerical solutions of Richards’ equation is central to reliable surface–subsurface flow and transport modeling, for example in long-term tracer simulations to determine catchment residence time distributions. In this study we assess the performance of a local Larson-Niklasson (LN) post-processing procedure for reconstructing mass-conservative velocities from a linear (P₁) Galerkin finite element solution of Richards’ equation. This approach, originally proposed for a-posteriori error estimation, modifies the standard finite element velocities by imposing local conservation on element patches. The resulting reconstructed flow field is characterized by continuous fluxes on element edges that can be efficiently used to drive a second order finite volume advective transport model. Through a series of tests of increasing complexity that compare results from the LN scheme to those using velocity fields derived directly from the P₁ Galerkin solution, we show that a locally mass-conservative velocity field is necessary to obtain accurate transport results. We also show that the accuracy of the LN reconstruction procedure is comparable to that of the inherently conservative mixed finite element approach, taken as a reference solution, but that the LN scheme has much lower computational costs. The numerical tests examine steady and unsteady, saturated and variably saturated, and homogeneous and heterogeneous cases along with initial and boundary conditions that include dry soil infiltration, alternating solute and water injection, and seepage face outflow. Typical problems that arise with velocities derived from P₁ Galerkin solutions include outgoing solute flux from no-flow boundaries, solute entrapment in zones of low hydraulic conductivity, and occurrences of anomalous sources and sinks. In addition to inducing significant mass balance errors, such manifestations often lead to oscillations in concentration values that can moreover cause the numerical solution to explode. These problems do not occur when using LN post-processed velocities.</p

    Impact of grid resolution on the integrated and distributed response of a coupled surface-subsurface hydrological model for the des Anglais catchment, Quebec

    No full text
    Digital elevation models (DEMs) at different resolutions (180, 360, and 720 m) are used to examine the impact of different levels of landscape representation on the hydrological response of a 690-km2 catchment in southern Quebec. Frequency distributions of local slope, plan curvature, and drainage area are calculated for each grid size resolution. This landscape analysis reveals that DEM grid size significantly affects computed topographic attributes, which in turn explains some of the differences in the hydrological simulations. The simulations that are then carried out, using a coupled, process-based model of surface and subsurface flow, examine the effects of grid size on both the integrated response of the catchment (discharge at the main outlet and at two internal points) and the distributed response (water table depth, surface saturation, and soil water storage). The results indicate that discharge volumes increase as the DEM is coarsened, and that coarser DEMs are also wetter overall in terms of water table depth and soil water storage. The reasons for these trends include an increase in the total drainage area of the catchment for larger DEM cell sizes, due to aggregation effects at the boundary cells of the catchment, and to a decrease in local slope and plan curvature variations, which in turn limits the capacity of the watershed to transmit water downslope and laterally. The results obtained also show that grid resolution effects are less pronounced during dry periods when soil moisture dynamics are mostly controlled by vertical fluxes of evaporation and percolation

    Examination of the seepage face boundary condition in subsurface and coupled surface/subsurface hydrological models

    No full text
    A seepage face is a nonlinear dynamic boundary that strongly affects pressure head distributions, water table fluctuations, and flow patterns. Its handling in hydrological models, especially under complex conditions such as heterogeneity and coupled surface/subsurface flow, has not been extensively studied. In this paper, we compare the treatment of the seepage face as a static (Dirichlet) versus dynamic boundary condition, we assess its resolution under conditions of layered heterogeneity, we examine its interaction with a catchment outlet boundary, and we investigate the effects of surface/subsurface exchanges on seepage faces forming at the land surface. The analyses are carried out with an integrated catchment hydrological model. Numerical simulations are performed for a synthetic rectangular sloping aquifer and for an experimental hillslope from the Landscape Evolution Observatory. The results show that the static boundary condition is not always an adequate stand-in for a dynamic seepage face boundary condition, especially under conditions of high rainfall, steep slope, or heterogeneity; that hillslopes with layered heterogeneity give rise to multiple seepage faces that can be highly dynamic; that seepage face and outlet boundaries can coexist in an integrated hydrological model and both play an important role; and that seepage faces at the land surface are not always controlled by subsurface flow. The paper also presents a generalized algorithm for resolving seepage face outflow that handles heterogeneity in a simple way, is applicable to unstructured grids, and is shown experimentally to be equivalent to the treatment of atmospheric boundary conditions in subsurface flow models.ECH

    Comparison of data assimilation techniquesfor a coupled model of surface and subsurface flow

    No full text
    Data assimilation in the geophysical sciences refers to methodologies to optimally merge model predictions and observations. The ensemble Kalman filter (EnKF) is a statistical sequential data assimilation technique explicitly developed for nonlinear filtering problems. It is based on a Monte Carlo approach that approximates the conditional probability densities of the variables of interest by a finite number of randomly generated model trajectories. In Newtonian relaxation or nudging (NN), which can be viewed as a special case of the classic Kalman filter, model variables are driven toward observations by adding to the model equations a forcing term, or relaxation component, that is proportional to the difference between simulation and observation. The forcing term contains four-dimensional weighting functions that can, ideally, incorporate prior knowledge about the characteristic scales of spatial and temporal variability of the state variable(s) being assimilated. In this study, we examined the EnKF and NN algorithms as implemented for a complex hydrologic model that simulates catchment dynamics, coupling a three-dimensional finite element Richards' equation solver for variably saturated porous media and a finite difference diffusion wave approximation for surface water flow. We report on the retrieval performance of the two assimilation schemes for a small catchment in Belgium. The results of the comparison show that nudging, while more straightforward and less expensive computationally, is not as effective as the ensemble Kalman filter in retrieving the true system state. We discuss some of the strengths and weaknesses, both physical and numerical, of the NN and EnKF schemes
    corecore