Refined conceptualization of TOPMODEL for shallow subsurface flows

The TOPMODEL framework was used to derive expressions that account for saturated and unsaturated flow through shallow soil on a hillslope. The resulting equations were the basis for a shallow‐soil TOPMODEL (STOPMODEL). The common TOPMODEL theory implicitly assumes a water table below the entire watershed and this does not conceptually apply to systems hydrologically controlled by shallow interflow of perched groundwater. STOPMODEL provides an approach for extending TOPMODEL's conceptualization to apply to shallow, interflow‐driven watersheds by using soil moisture deficit instead of water table depth as the state variable. Deriving STOPMODEL by using a hydraulic conductivity function that changes exponentially with soil moisture content results in equations that look very similar to those commonly associated with TOPMODEL. This alternative way of conceptualizing TOPMODEL makes the modelling approach available to researchers, planners, and engineers who work in areas where TOPMODEL was previously believed to be unsuited, such as the New York City Watershed in the Catskills region of New York State. Copyright © 2002 John Wiley & Sons, Ltd.     

Generalized LSG models with spatially varying Coriolis parameter

In this paper, we derive and study approximate balance models for nearly geostrophic shallow water flow where the Coriolis parameter is permitted to vary across the domain as long as it remains nondegenerate. This situation includes, for example, the β-plane approximation to the shallow water equations at mid-latitudes. Our approach is based on changing configuration space coordinates in the underlying variational principle in such a way that consistent asymptotics in the transformed Lagrangian leads to a degenerate Lagrangian structure. In this paper, we restrict our attention to first-order models. We show that the resulting models can be formulated in terms of an advected potential vorticity with a nonlinear vorticity inversion relation. We study the associated solvability conditions and identify a subfamily of models for which these conditions are satisfied without additional restrictions on the data. Finally, we provide the link between our framework and the theory of constrained Hamiltonian systems.     

Two-dimensional numerical model of two-layer shallow water equations for confluence simulation

This study presents a finite-volume explicit method to solve 2D two-layer shallow water equations. This numerical model is intended to describe two-layer shallow flows in which the superposed layers differ in velocity, density and rheology in a two-dimensional domain. The rheological behavior of mudflow or debris flow is called the Bingham fluid. Thus, the shear stress on rigid bed can be derived from the constitutive equation. The computational approach adopts the HLL scheme, a novel approach for the purpose of computing a Godunov flux and solving the Riemann problem approximately proposed by Harten, Lax and van Leer, as a basic building block, treats the bottom slope by lateralizing the momentum flux, and refines the scheme using the Strang splitting to manage the frictional source term. This study successfully performed 2D two-layer shallow water computations on a rigid bed. The proposed numerical model can describe the variety of depths and velocities of substances including water and mud, when the hyperconcentrated tributary flows into the main river. The analytical results in this study will be valuable for further advanced research and for designing or planning hydraulic engineering structures.     

A mass conservative scheme for simulating shallow flows over variable topographies using unstructured grid

Most available numerical methods face problems, in the presence of variable topographies, due to the imbalance between the source and flux terms. Treatments for this problem generally work well for structured grids, but most of them are not directly applicable for unstructured grids. On the other hand, despite of their good performance for discontinuous flows, most available numerical schemes (such as HLL flux and ENO schemes) induce a high level of numerical diffusion in simulating recirculating flows. A numerical method for simulating shallow recirculating flows over a variable topography on unstructured grids is presented. This mass conservative approach can simulate different flow conditions including recirculating, transcritical and discontinuous flows over variable topographies without upwinding of source terms and with a low level of numerical diffusion. Different numerical tests cases are presented to show the performance of the scheme for some challenging problems.     

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.     

A robust well-balanced finite volume model for shallow water flows with wetting and drying over irregular terrain

An unstructured Godunov-type finite volume model is developed for the numerical simulation of geometrically challenging two-dimensional shallow water flows with wetting and drying over convoluted topography. In the framework of sloping bottom model, a modified formulation of shallow water equations is used to preserve mass conservation during flooding and recession. The key ingredient of the model is the use of this combination of the sloping bottom model and the modified shallow water equations to provide a robust technique for wet/dry fronts tracking and, together with centered discretization of the bed slope source term, to exactly preserve the static flow on irregular topographies. The variable reconstruction technique ensures nonnegative reconstructed water depth and reasonable reconstructed velocity, and the friction terms are solved by semi-implicit scheme that does not invert the direction of velocity components. The robustness and accuracy of the proposed model are assessed by comparing numerical and reference results of extensive test cases. Moreover, the results of a dam-break flooding over real topography are presented to show the capability of the model on field-scale application.     

Sensitivity equations for hyperbolic conservation law-based flow models

The present paper focuses on the governing equations for the sensitivity of the variables to the parameters in flow models that can be described by one-dimensional scalar, hyperbolic conservation laws. The sensitivity is shown to obey a hyperbolic, scalar conservation law. The sensitivity is a conserved scalar except in the case of discontinuous flow solutions, where an extra, point source term must be added to the equations in order to enforce conservation. The propagation speed of the sensitivity waves being identical to that of the conserved variable in the original conservation law, the system of conservation laws formed by the original hyperbolic equation and the equation satisfied by the sensitivity is linearly degenerate. A consequence on the solution of the Riemann problem is that rarefaction waves for the variable of the original equation result in vacuum regions for the sensitivity. The numerical solution of the hyperbolic conservation law for the sensitivity by finite volume methods requires the implementation of a specific shock detection procedure. A set of necessary conditions is defined for the discretisation of the source term in the sensitivity equation. An application to the one-dimensional kinematic wave equation shows that the proposed numerical technique allows analytical solutions to be reproduced correctly. The computational examples show that first-order numerical schemes do not yield satisfactory numerical solutions in the neighbourhood of moving shocks and that higher-order schemes, such as the MUSCL scheme, should be used for sharp transients.     

Inversion of ocean-bottom seismometer (OBS) waveforms for oceanic crust structure: a synthetic study

The waveform inversion method is applied—using synthetic ocean-bottom seismometer(OBS) data—to study oceanic crust structure. A niching genetic algorithm(NGA) is used to implement the inversion for the thickness and P-wave velocity of each layer, and to update the model by minimizing the objective function, which consists of the misfit and cross-correlation of observed and synthetic waveforms. The influence of specific NGA method parameters is discussed, and suitable values are presented.The NGA method works well for various observation systems, such as those with irregular and sparse distribution of receivers as well as single receiver systems. A strategy is proposed to accelerate the convergence rate by a factor of five with no increase in computational complexity; this is achieved using a first inversion with several generations to impose a restriction on the preset range of each parameter and then conducting a second inversion with the new range. Despite the successes of this method,its usage is limited. A shallow water layer is not favored because the direct wave in water will suppress the useful reflection signals from the crust. A more precise calculation of the air-gun source signal should be considered in order to better simulate waveforms generated in realistic situations; further studies are required to investigate this issue.     

Multiple objective evaluation of a simple phosphorus transfer model

A simple phosphorus (P) transfer model of the Welland catchment, UK, is evaluated against multiple objective functions using a Monte Carlo approach that combines calibration, identifiability, sensitivity and uncertainty analysis. The model is based on simple conceptual rainfall‐runoff and river routing components, combined with estimates of the daily non‐point source load derived from annual landuse‐based export coefficients, disaggregated as a function of the runoff. The model has limited data requirements, consistent with data availability, and is parsimoneous with respect to the number of parameters identified through inverse modelling. The best performing parameter sets capture the main aspects of the observed flow and total P (TP) concentrations and provide a suitable basis for a decision‐support tool. However, a trade‐off is evident between matching the observed flow peaks, flow recessions and TP concentrations simultaneously, highlighting some limitations of the model structure and/or calibration data. Model analysis indicates that daily non‐point source load cannot be described as a function of near‐surface runoff and land use alone, but that other influences, including seasonality, are important. However, further model development to improve performance is likely to introduce additional complexity (in terms of parameter numbers), and hence additional problems of parameter identifiability and output uncertainty, which in turn raises issues of the information content of the available data. Copyright © 2004 John Wiley & Sons, Ltd.     

Step-on versus step-off signals in time-domain controlled source electromagnetic methods using a grounded electric dipole

The time-domain controlled source electromagnetic method is a geophysical prospecting tool applied to image the subsurface resistivity distribution on land and in the marine environment. In its most general set-up, a square-wave current is fed into a grounded horizontal electric dipole, and several electric and magnetic field receivers at defined offsets to the imposed current measure the electromagnetic response of the Earth. In the marine environment, the application often uses only inline electric field receivers that, for a 50% duty-cycle current waveform, include both step-on and step-off signals. Here, forward and inverse 1D modelling is used to demonstrate limited sensitivity towards shallow resistive layers in the step-off electric field when transmitter and receivers are surrounded by conductive seawater. This observation is explained by a masking effect of the direct current signal that flows through the seawater and primarily affects step-off data. During a step-off measurement, this direct current is orders of magnitude larger than the inductive response at early and intermediate times, limiting the step-off sensitivity towards shallow resistive layers in the seafloor. Step-on data measure the resistive layer at times preceding the arrival of the direct current signal leading to higher sensitivity compared to step-off data. Such dichotomous behaviour between step-on and step-off data is less obvious in onshore experiments due to the lack of a strong overlying conductive zone and corresponding masking effect from direct current flow. Supported by synthetic 1D inversion studies, we conclude that time-domain controlled source electromagnetic measurements on land should apply both step-on and step-off data in a combined inversion approach to maximize signal-to-noise ratios and utilize the sensitivity characteristics of each signal. In an isotropic marine environment, step-off electric fields have inferior sensitivity towards shallow resistive layers compared to step-on data, resulting in an increase of non-uniqueness when interpreting step-off data in a single or combined inversion.     

An automated routing methodology to enable direct rainfall in high resolution shallow water models

Recent high profile flood events have highlighted the need for hydraulic models capable of simulating pluvial flooding in urban areas. This paper presents a constant velocity rainfall routing scheme that provides this ability within the LISFLOOD‐FP hydraulic modelling code. The scheme operates in place of the shallow water equations within cells where the water depth is below a user‐defined threshold, enabling rainfall‐derived water to be moved from elevated features such as buildings or curbstones without causing instabilities in the solution whilst also yielding a reduction in the overall computational cost of the simulation. Benchmarking against commercial modelling packages using a pluvial and point‐source test case demonstrates that the scheme does not impede the ability of LISFLOOD‐FP to match both predicted depths and velocities of full shallow water models. The stability of the scheme in conditions unsuitable for traditional two‐dimensional hydraulic models is then demonstrated using a pluvial test case over a complex urban digital elevation model containing buildings. Deterministic single‐parameter sensitivity analyses undertaken using this test case show limited sensitivity of predicted water depths to both the chosen routing speed within a physically plausible range and values of the depth threshold parameter below 10 mm. Local instabilities can occur in the solution if the depth threshold is >10 mm, but such values are not required even when simulating extreme rainfall rates. The scheme yields a reduction in model runtime of ~25% due to the reduced number of cells for which the hydrodynamic equations have to be solved. Copyright © 2012 John Wiley & Sons, Ltd.     

Inversion of dynamic production data for permeability: fast streamline-based computation of sensitivity coefficients of fractional flow rate

Generation of permeability field in a reservoir model that matchs historical dynamic production data requires an inverse calculation. A gradient method is typically used to solve the inverse minimization problem and requires sensitivity coefficients of reservoir responses, e.g. fractional flow rate or pressure, with respect to the change in the permeability. This paper presents a novel semi-analytical streamline-based method for computing such sensitivity coefficients under the framework of two-phase (oil-water) flow conditions. This method is shown to be significantly faster and generate permeability fields with lower objective function than the traditional perturbation method. The method decomposes the multiple-dimensional full flow problem into multiple 1D problems along streamlines. The sensitivity of fractional flow rate at the production well is directly related to the sensitivity of time-of-flight (TOF) along each individual streamline and the sensitivity of pressure at grid cells along the streamline. The sensitivity of TOF of a streamline can be obtained analytically. The sensitivity of pressure is obtained as part of a fast single phase flow simulation. The proposed method is implemented in a geostatistically based inverse technique, called the sequential self-calibration (SSC) method. Results for fractional flow rate sensitivities are presented and compared with the traditional perturbation method. This new method can be easily extended to compute sensitivity coefficients of saturation (concentration) data.     

Analytic Element Modeling of Steady Interface Flow in Multilayer Aquifers Using AnAqSim

This paper presents the analytic element modeling approach implemented in the software AnAqSim for simulating steady groundwater flow with a sharp fresh‐salt interface in multilayer (three‐dimensional) aquifer systems. Compared with numerical methods for variable‐density interface modeling, this approach allows quick model construction and can yield useful guidance about the three‐dimensional configuration of an interface even at a large scale. The approach employs subdomains and multiple layers as outlined by Fitts (2010) with the addition of discharge potentials for shallow interface flow (Strack 1989). The following simplifying assumptions are made: steady flow, a sharp interface between fresh‐ and salt water, static salt water, and no resistance to vertical flow and hydrostatic heads within each fresh water layer. A key component of this approach is a transition to a thin fixed minimum fresh water thickness mode when the fresh water thickness approaches zero. This allows the solution to converge and determine the steady interface position without a long transient simulation. The approach is checked against the widely used numerical codes SEAWAT and SWI/MODFLOW and a hypothetical application of the method to a coastal wellfield is presented.     

A weakly nonlinear stability of centrally symmetric magnetohydrodynamic systems to perturbations involving large scales

The problem of weakly nonlinear stability of 3-D centrally symmetric magnetohydrodynamic systems to perturbations involving large scales is considered. It is assumed that large space-time scales are absent in the magnetohydrodynamic state under study, which is stable with respect to perturbations whose scales are as small as those of the state itself. Equations derived by asymptotic methods for average fields of perturbations generalize the Navier-Stokes and magnetic induction equations. They include a combined eddy diffusion operator, generally anisotropic and not necessarily negative definite, and additional quadratic terms. An effective method is proposed for the calculation of coefficients of eddy diffusion and advection in equations governing average fields.     

Theoretical estimation of free and entrapped nonwetting–wetting fluid interfacial areas in porous media

Fluid–fluid interfacial areas play important roles in numerous subsurface processes such as dissolution, volatilization, and adsorption. Integral expressions have been derived to estimate both entrapped (discontinuous) and free (continuous) nonwetting fluid–wetting fluid specific interfacial areas in porous media. The expressions, compatible with widely used capillary head-saturation and entrapment models, require information on capillary head-saturation relation parameters, porosity, and fluid-pair interfacial tension. In addition, information on the maximum entrapped nonwetting fluid saturation as well as the main drainage branch reversal point for water and total liquid saturations is necessary to estimate entrapped fluid interfacial areas. Implementation of the interfacial area equations in continuum-based multifluid flow simulators is straightforward since no additional parameters are needed than those required by the simulators to complete the multifluid flow computations. A limited sensitivity analysis, based on experimentally obtained parameter values, showed that imposed variations resulted in logical and consistent changes in predicted specific interfacial areas for both entrapped and free nonwetting fluid–wetting fluid systems. A direct comparison with published experimental work to test the derived expressions was limited to free air–water systems and yielded reasonable results. Such comparisons are often not possible because of the lack of information given on retention parameters, and variables used to determine nonwetting fluid entrapment. This contribution is dedicated to John W. Cary.     

Improved bounds on linear instability of barotropic zonal flow within the shallow water equations

Here we develop mathematical results to describe the location of linear instability of a parallel mean flow within the framework of the shallow water equations; growth estimates of near neutral modes (for disturbances subcritical with respect to gravity wave speed) in the cases of non-rotating and rotating shallow water. The bottom topography is taken to be one-dimensional and the isobaths are parallel to the mean flow. In the case of a rotating fluid, the isobaths and the mean flow are assumed to be zonal. The flow is front-like: there is a monotonic increase of mean flow velocity. Our results show that for barotropic flows the location of instabilities will be a semi-ellipse region in the complex wave velocity plane, that is based on the wave-number, Froude number, and depth of the fluid layer. We also explore the instability region for the case of spatially unbounded mean velocity profiles for non-rotating shallow water.     

The GeoClaw software for depth-averaged flows with adaptive refinement

Many geophysical flow or wave propagation problems can be modeled with two-dimensional depth-averaged equations, of which the shallow water equations are the simplest example. We describe the GeoClaw software that has been designed to solve problems of this nature, consisting of open source Fortran programs together with Python tools for the user interface and flow visualization. This software uses high-resolution shock-capturing finite volume methods on logically rectangular grids, including latitude-longitude grids on the sphere. Dry states are handled automatically to model inundation. The code incorporates adaptive mesh refinement to allow the efficient solution of large-scale geophysical problems. Examples are given illustrating its use for modeling tsunamis and dam-break flooding problems. Documentation and download information is available at www.clawpack.org/geoclaw.     

Estimating flow parameter distributions using ground-penetrating radar and hydrological measurements during transient flow in the vadose zone

Methods for estimating the parameter distributions necessary for modeling fluid flow and contaminant transport in the shallow subsurface are in great demand. Soil properties such as permeability, porosity, and water retention are typically estimated through the inversion of hydrological data (e.g., measurements of capillary pressure and water saturation). However, ill-posedness and non-uniqueness commonly arise in such non-linear inverse problems making their solutions elusive. Incorporating additional types of data, such as from geophysical methods, may greatly improve the success of inverse modeling. In particular, ground-penetrating radar (GPR) methods have proven sensitive to subsurface fluid flow processes and appear promising for such applications. In the present work, an inverse technique is presented which allows for the estimation of flow parameter distributions and the prediction of flow phenomena using GPR and hydrological measurements collected during a transient flow experiment. Specifically, concepts from the pilot point method were implemented in a maximum a posteriori (MAP) framework to allow for the generation of permeability distributions that are conditional to permeability point measurements, that maintain specified patterns of spatial correlation, and that are consistent with geophysical and hydrological data. The current implementation of the approach allows for additional flow parameters to be estimated concurrently if they are assumed uniform and uncorrelated with the permeability distribution. (The method itself allows for heterogeneity in these parameters to be considered, and it allows for parameters of the petrophysical and semivariogram models to be estimated as well.) Through a synthetic example, performance of the method is evaluated under various conditions, and some conclusions are made regarding the joint use of transient GPR and hydrological measurements in estimating fluid flow parameters in the vadose zone.     

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

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

Adaptive volume penalization for ocean modeling

The development of various volume penalization techniques for use in modeling topographical features in the ocean is the focus of this paper. Due to the complicated geometry inherent in ocean boundaries, the stair-step representation used in the majority of current global ocean circulation models causes accuracy and numerical stability problems. Brinkman penalization is the basis for the methods developed here and is a numerical technique used to enforce no-slip boundary conditions through the addition of a term to the governing equations. The second aspect to this proposed approach is that all governing equations are solved on a nonuniform, adaptive grid through the use of the adaptive wavelet collocation method. This method solves the governing equations on temporally and spatially varying meshes, which allows higher effective resolution to be obtained with less computational cost. When penalization methods are coupled with the adaptive wavelet collocation method, the flow near the boundary can be well-resolved. It is especially useful for simulations of boundary currents and tsunamis, where flow near the boundary is important. This paper will give a thorough analysis of these methods applied to the shallow water equations, as well as some preliminary work applying these methods to volume penalization for bathymetry representation for use in either the nonhydrostatic or hydrostatic primitive equations.