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.     

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.     

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.     

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.     

The role of hand calculations in ground water flow modeling

Most ground water modeling courses focus on the use of computer models and pay little or no attention to traditional analytic solutions to ground water flow problems. This shift in education seems logical. Why waste time to learn about the method of images, or why study analytic solutions to one-dimensional or radial flow problems? Computer models solve much more realistic problems and offer sophisticated graphical output, such as contour plots of potentiometric levels and ground water path lines. However, analytic solutions to elementary ground water flow problems do have something to offer over computer models: insight. For instance, an analytic one-dimensional or radial flow solution, in terms of a mathematical expression, may reveal which parameters affect the success of calibrating a computer model and what to expect when changing parameter values. Similarly, solutions for periodic forcing of one-dimensional or radial flow systems have resulted in a simple decision criterion to assess whether or not transient flow modeling is needed. Basic water balance calculations may offer a useful check on computer-generated capture zones for wellhead protection or aquifer remediation. An easily calculated "characteristic leakage length" provides critical insight into surface water and ground water interactions and flow in multi-aquifer systems. The list goes on. Familiarity with elementary analytic solutions and the capability of performing some simple hand calculations can promote appropriate (computer) modeling techniques, avoids unnecessary complexity, improves reliability, and is likely to save time and money. Training in basic hand calculations should be an important part of the curriculum of ground water modeling courses.     

Numerical resolution of well-balanced shallow water equations with complex source terms

This paper presents a well-balanced numerical scheme for simulating frictional shallow flows over complex domains involving wetting and drying. The proposed scheme solves, in a finite volume Godunov-type framework, a set of pre-balanced shallow water equations derived by considering pressure balancing. Non-negative reconstruction of Riemann states and compatible discretization of slope source term produce stable and well-balanced solutions to shallow flow hydrodynamics over complex topography. The friction source term is discretized using a splitting implicit scheme. Limiting value of the friction force is derived to ensure stability. This new numerical scheme is validated against four theoretical benchmark tests and then applied to reproduce a laboratory dam break over a domain with irregular bed profile.     

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.     

A unified formulation for the three-dimensional shallow water equations using orthogonal co-ordinates: theory and application

In this paper, the formulations of the primitive equations for shallow water flow in various horizontal co-ordinate systems and the associated finite difference grid options used in shallow water flow modelling are reviewed. It is observed that horizontal co-ordinate transformations do not affect the chosen co-ordinate system and representation in the vertical, and are the same for the three- and two-dimensional cases. A systematic derivation of the equations in tensor notation is presented, resulting in a unified formulation for the shallow water equations that covers all orthogonal horizontal grid types of practical interest. This includes spherical curvilinear orthogonal co-ordinate systems on the globe. Computational efficiency can be achieved in a single computer code. Furthermore, a single numerical algorithmic code implementation satisfies. All co-ordinate system specific metrics are determined as part of a computer-aided model grid design, which supports all four orthogonal grid types. Existing intuitive grid design and visual interpretation is conserved by appropriate conformal mappings, which conserve spherical orthogonality in planar representation. A spherical curvilinear co-ordinate solution of wind driven steady channel flow applying a strongly distorted grid is shown to give good agreement with a regular spherical co-ordinate model approach and the solution based on a β-plane approximation. Especially designed spherical curvilinear boundary fitted model grids are shown for typhoon surge propagation in the South China Sea and for ocean-driven flows through Malacca Straits. By using spherical curvilinear grids the number of grid points in these single model grid applications is reduced by a factor of 50–100 in comparison with regular spherical grids that have the same horizontal resolution in the area of interest. The spherical curvilinear approach combines the advantages of the various grid approaches, while the overall computational effort remains acceptable for very large model domains.     

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.     

Combined hillslope diffusion and sediment transport simulation applied to landscape dynamics modelling

We present a new numerical approach for simulating geomorphic and stratigraphic processes that combines open‐channel flow with non‐uniform sediment transport law and semi‐empirical diffusive mass wasting. It is designed to facilitate modelling of surface processes across multiple space‐ and time‐scales, and under a variety of environmental and tectonic conditions. The physics of open‐channel flow is primarily based on an adapted Lagrangian formulation of shallow‐water equations. The interaction between flow and surface geology is performed by a non‐uniform total‐load sediment transport law. Additional hillslope processes are simulated using a semi‐empirical method based on a diffusion approach. In the implementation, the resolution of flow dynamics is made on a triangulated grid automatically mapped and adaptively remeshed over a regular orthogonal stratigraphic mesh. These new methods reduce computational time while preserving stability and accuracy of the physical solutions. In order to illustrate the potential of this method for landscape and sedimentary system modelling, we present a set of three generic experiments focusing on assessing the influence of contrasting erodibilities on the evolution of an active bedrock landscape. The modelled ridges morphometrics satisfy established relationships for drainage network geometry and slope distribution, and provide quantitative information on the relative impact of hillslope and channel processes, sediment discharge and alluviation. Our results suggest that contrasting erodibility can stimulate autogenic changes in erosion rate and influence the landscape morphology and preservation. This approach offers new opportunities to investigate joint landscape and sedimentary systems response to external perturbations. The possibility to define and track a large number of materials makes the implementation highly suited to model source‐to‐sink problems where material dispersion is the key question that needs to be addressed, such as natural resources exploration and basin analysis. Copyright © 2014 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.