A nearly analytic symplectic partitioned Runge–Kutta method based on a locally one-dimensional technique for solving two-dimensional acoustic wave equations

In this paper, we develop a new nearly analytic symplectic partitioned Runge–Kutta method based on locally one-dimensional technique for numerically solving two-dimensional acoustic wave equations. We first split two-dimensional acoustic wave equation into the local one-dimensional equations and transform each of the split equations into a Hamiltonian system. Then, we use both a nearly analytic discrete operator and a central difference operator to approximate the high-order spatial differential operators, which implies the symmetry of the discretized spatial differential operators, and we employ the partitioned second-order symplectic Runge–Kutta method to numerically solve the resulted semi-discrete Hamiltonian ordinary differential equations, which results in fully discretized scheme is symplectic unlike conventional nearly analytic symplectic partitioned Runge–Kutta methods. Theoretical analyses show that the nearly analytic symplectic partitioned Runge–Kutta method based on locally one-dimensional technique exhibits great higher stability limits and less numerical dispersion than the nearly analytic symplectic partitioned Runge–Kutta method. Numerical experiments are conducted to verify advantages of the nearly analytic symplectic partitioned Runge–Kutta method based on locally one-dimensional technique, such as their computational efficiency, stability, numerical dispersion and long-term calculation capability.     

Second-order characteristic methods for advection–diffusion equations and comparison to other schemes

We develop two characteristic methods for the solution of the linear advection diffusion equations which use a second order Runge–Kutta approximation of the characteristics within the framework of the Eulerian–Lagrangian localized adjoint method. These methods naturally incorporate all three types of boundary conditions in their formulations, are fully mass conservative, and generate regularly structured systems which are symmetric and positive definite for most combinations of the boundary conditions. Extensive numerical experiments are presented which compare the performance of these two Runge–Kutta methods to many other well perceived and widely used methods which include many Galerkin methods and high resolution methods from fluid dynamics.     

An asynchronous solver for systems of ODEs linked by a directed tree structure

This paper documents our development and evaluation of a numerical solver for systems of sparsely linked ordinary differential equations in which the connectivity between equations is determined by a directed tree. These types of systems arise in distributed hydrological models. The numerical solver is based on dense output Runge–Kutta methods that allow for asynchronous integration. A partition of the system is used to distribute the workload among different processes, enabling a parallel implementation that capitalizes on a distributed memory system. Communication between processes is performed asynchronously. We illustrate the solver capabilities by integrating flow transport equations for a ∼17,000 km² river basin subdivided into 305,000 sub-watersheds that are interconnected by the river network. Numerical experiments for a few models are performed and the runtimes and scalability on our parallel computer are presented. Efficient numerical integrators such as the one demonstrated here bring closer to reality the goal of implementing fully distributed real-time flood forecasting systems supported by physics based hydrological models and high-quality/high-resolution rainfall products.     

Dynamic response of an elastic-viscoplastic system in modal co-ordinates

A dynamic analysis of elastic–viscoplastic systems, incorporating the modal co-ordinate transformation technique, is presented. The formulation results in uncoupled incremental equations of motion with respect to the modal co-ordinates. The elastic–viscoplastic model adopted allows the analysis not to involve yielding regions and loading/unloading processes. An implicit Runge–Kutta scheme together with the Newton–Raphson method are used to solve the non-linear constitutive equations. Stability and accuracy of the numerical solution are improved by utilizing a local time step sub-incrementing procedure. Applications of the analyses to multi-storey shear buildings show that good results can be obtained for the maximum displacement response by including only a few lower modes in the computation, but the prediction of the ductility factor response tends to underestimate the peak values when too few modes are used. In addition, stable and valid results can be obtained even with a sizable time step increment.     

High-order balanced CWENO scheme for movable bed shallow water equations

In this work the numerical integration of 1D shallow water equations (SWE) over movable bed is performed using a well-balanced central weighted essentially non-oscillatory (CWENO) scheme, fourth-order accurate in space and in time. Time accuracy is obtained following a Runge–Kutta (RK) procedure, coupled with its natural continuous extension (NCE). Spatial accuracy is obtained using WENO reconstructions of conservative variables and of flux and bed derivatives. An original treatment for bed slope source term, which maintains the established order of accuracy and satisfies the property of exactly preserving the quiescent flow (C-property), is introduced in the scheme. This treatment consists of two procedures. The former involves the evaluation of the point-values of the flux derivative, considered as a whole with the bed slope source term. The latter involves the spatial integration of the source term, analytically manipulated to take advantage from the expected regularity of the free surface elevation. The high accuracy of the scheme allows to obtain good results using coarse grids, with consequent gain in terms of computational effort. The well-balancing of the scheme allows to reproduce small perturbations of the free surface and of the bottom otherwise of the same order of magnitude of the numerical errors induced by the non-balancing. The accuracy, the well-balancing and the good resolution of the model in reproducing free surface flow over movable bed are tested over analytical solutions and over numerical results available in literature.     

Modelling overland flow based on Saint‐Venant equations for a discretized hillslope system

Mathematical modelling of overland flow is a critical task in simulating transport of water, sediment and other pollutants from land surfaces to receiving waters. In this paper, an overland flow routing method is developed based on the Saint‐Venant equations using a discretized hillslope system for areas with high roughness and steep slope. Under these conditions, the momentum equation reduces to a unique relationship between the flow depth and discharge. A hillslope is treated as a system divided into several subplanes. A set of first‐order non‐linear differential equations for subsequent subplanes are solved analytically using Chezy's formula in lieu of the momentum equation. Comparison of the analytical solution of the first‐order non‐linear ordinary differential equations and a numerical solution using the Runge‐Kutta method shows a relative error of 0·3%. Using runoff data reported in the literature, comparison between the new approach and a numerical solution of the full Saint‐Venant equations showed a close agreement. Copyright © 2002 John Wiley & Sons, Ltd.     

Novel coupling Rosenbrock‐based algorithms for real‐time dynamic substructure testing

Real‐time testing with dynamic substructuring is a novel experimental technique capable of assessing the behaviour of structures subjected to dynamic loadings including earthquakes. The technique involves recreating the dynamics of the entire structure by combining an experimental test piece consisting of part of the structure with a numerical model simulating the remainder of the structure. These substructures interact in real time to emulate the behaviour of the entire structure. Time integration is the most versatile method for analysing the general case of linear and non‐linear semi‐discretized equations of motion. In this paper we propose for substructure testing, L‐stable real‐time (LSRT) compatible integrators with two and three stages derived from the Rosenbrock methods. These algorithms are unconditionally stable for uncoupled problems and entail a moderate computational cost for real‐time performance. They can also effectively deal with stiff problems, i.e. complex emulated structures for which solutions can change on a time scale that is very short compared with the interval of time integration, but where the solution of interest changes on a much longer time scale. Stability conditions of the coupled substructures are analysed by means of the zero‐stability approach, and the accuracy of the novel algorithms in the coupled case is assessed in both the unforced and forced conditions. LSRT algorithms are shown to be more competitive than popular Runge–Kutta methods in terms of stability, accuracy and ease of implementation. Numerical simulations and real‐time substructure tests are used to demonstrate the favourable properties of the proposed algorithms. Copyright © 2007 John Wiley & Sons, Ltd.     

Numerical modeling of 3D fully nonlinear potential periodic waves

A simple and exact numerical scheme for long-term simulations of 3D potential fully nonlinear periodic gravity waves is suggested. The scheme is based on the surface-following nonorthogonal curvilinear coordinate system. Velocity potential is represented as a sum of analytical and nonlinear components. The Poisson equation for the nonlinear component of velocity potential is solved iteratively. Fourier transform method, the second-order accuracy approximation of vertical derivatives on a stretched vertical grid and the fourth-order Runge–Kutta time stepping are used. The scheme is validated by simulation of steep Stokes waves. A one-processor version of the model for PC allows us to simulate evolution of a wave field with thousands degrees of freedom for hundreds of wave periods. The scheme is designed for investigation of nonlinear 2D surface waves, generation of extreme waves, and direct calculations of nonlinear interactions.     

A limiter-based well-balanced discontinuous Galerkin method for shallow-water flows with wetting and drying: One-dimensional case

An important part in the numerical simulation of tsunami and storm surge events is the accurate modeling of flooding and the appearance of dry areas when the water recedes. This paper proposes a new algorithm to model inundation events with piecewise linear Runge–Kutta discontinuous Galerkin approximations applied to the shallow water equations. This study is restricted to the one-dimensional case and shows a detailed analysis and the corresponding numerical treatment of the inundation problem.The main feature is a velocity based “limiting” of the momentum distribution in each cell, which prevents instabilities in case of wetting or drying situations. Additional limiting of the fluid depth ensures its positivity while preserving local mass conservation. A special flux modification in cells located at the wet/dry interface leads to a well-balanced method, which maintains the steady state at rest. The discontinuous Galerkin scheme is formulated in a nodal form using a Lagrange basis. The proposed wetting and drying treatment is verified with several numerical simulations. These test cases demonstrate the well-balancing property of the method and its stability in case of rapid transition of the wet/dry interface. We also verify the conservation of mass and investigate the convergence characteristics of the scheme.     

Robust numerical solution of the reservoir routing equation

The robustness of numerical methods for the solution of the reservoir routing equation is evaluated. The methods considered in this study are: (1) the Laurenson–Pilgrim method, (2) the fourth-order Runge–Kutta method, and (3) the fixed order Cash–Karp method. Method (1) is unable to handle nonmonotonic outflow rating curves. Method (2) is found to fail under critical conditions occurring, especially at the end of inflow recession limbs, when large time steps (greater than 12 min in this application) are used. Method (3) is computationally intensive and it does not solve the limitations of method (2). The limitations of method (2) can be efficiently overcome by reducing the time step in the critical phases of the simulation so as to ensure that water level remains inside the domains of the storage function and the outflow rating curve. The incorporation of a simple backstepping procedure implementing this control into the method (2) yields a robust and accurate reservoir routing method that can be safely used in distributed time-continuous catchment models.     

A discontinuous Galerkin method for two-dimensional flow and transport in shallow water

A discontinuous Galerkin (DG) finite element method is described for the two-dimensional, depth-integrated shallow water equations (SWEs). This method is based on formulating the SWEs as a system of conservation laws, or advection–diffusion equations. A weak formulation is obtained by integrating the equations over a single element, and approximating the unknowns by piecewise, possibly discontinuous, polynomials. Because of its local nature, the DG method easily allows for varying the polynomial order of approximation. It is also “locally conservative”, and incorporates upwinded numerical fluxes for modeling problems with high flow gradients. Numerical results are presented for several test cases, including supercritical flow, river inflow and standard tidal flow in complex domains, and a contaminant transport scenario where we have coupled the shallow water flow equations with a transport equation for a chemical species.     

A novel semi-analytical approach for non-uniform vegetated flows

In this paper a semi-analytical approach is proposed to understand the mechanism by which a non-uniform vegetated flowpasses over a finite thick soil layer covered with grass. The flow region is divided into three layers: a homogenous water layer, a mixed water-grass layer, and a finite thick soil layer (hereafter referred to as the water layer, the grass layer, and the soil layer). The flow of the water layer is governed by the Navier–Stokes equations. Both the grass and soil layers are regarded as porous media and the Biot’s theory of poroelasticity is applied to the porous medium flow. The semi-closed solutions are then obtained by the Runge–Kutta method. The drag force induced by the flow through the grass layer and the flow profiles of three patterns: submerged grass, emergent grass and mixed type are also discussed.     

Improved numerical solvers for implicit coupling of subsurface and overland flow

Due to complex dynamics inherent in the physical models, numerical formulation of subsurface and overland flow coupling can be challenging to solve. ParFlow is a subsurface flow code that utilizes a structured grid discretization in order to benefit from fast and efficient structured solvers. Implicit coupling between subsurface and overland flow modes in ParFlow is obtained by prescribing an overland boundary condition at the top surface of the computational domain. This form of implicit coupling leads to the activation and deactivation of the overland boundary condition during simulations where ponding or drying events occur. This results in a discontinuity in the discrete system that can be challenging to resolve. Furthermore, the coupling relies on unstructured connectivities between the subsurface and surface components of the discrete system, which makes it challenging to use structured solvers to effectively capture the dynamics of the coupled flow. We present a formulation of the discretized algebraic system that enables the use of an analytic form of the Jacobian for the Newton–Krylov solver, while preserving the structured properties of the discretization. An effective multigrid preconditioner is extracted from the analytic Jacobian and used to precondition the Jacobian linear system solver. We compare the performance of the new solver against one that uses a finite difference approximation to the Jacobian within the Newton–Krylov approach, previously used in the literature. Numerical results explores the effectiveness of using the analytic Jacobian for the Newton–Krylov solver, and highlights the performance of the new preconditioner and its cost. The results indicate that the new solver is robust and generally outperforms the solver that is based on the finite difference approximation to the Jacobian, for problems where the overland boundary condition is activated and deactivated during the simulation. A parallel weak scaling study highlights the efficiency of the new solver.     

A model based on Hirano-Exner equations for two-dimensional transient flows over heterogeneous erodible beds

In order to study the morphological evolution of river beds composed of heterogeneous material, the interaction among the different grain sizes must be taken into account. In this paper, these equations are combined with the two-dimensional shallow water equations to describe the flow field. The resulting system of equations can be solved in two ways: (i) in a coupled way, solving flow and sediment equations simultaneously at a given time-step or (ii) in an uncoupled manner by first solving the flow field and using the magnitudes obtained at each time-step to update the channel morphology (bed and surface composition). The coupled strategy is preferable when dealing with strong and quick interactions between the flow field, the bed evolution and the different particle sizes present on the bed surface. A number of numerical difficulties arise from solving the fully coupled system of equations. These problems are reduced by means of a weakly-coupled strategy to numerically estimate the wave celerities containing the information of the bed and the grain sizes present on the bed. Hence, a two-dimensional numerical scheme able to simulate in a self-stable way the unsteady morphological evolution of channels formed by cohesionless grain size mixtures is presented. The coupling technique is simplified without decreasing the number of waves involved in the numerical scheme but by simplifying their definitions. The numerical results are satisfactorily tested with synthetic cases and against experimental data.     

From Satellite Gradiometry Data to Subcrustal Stress Due to Mantle Convection

Subcrustal stress induced by mantle convection can be determined by the Earth’s gravitational potential. In this study, the spherical harmonic expansion of the simplified Navier–Stokes equation is developed further so satellite gradiometry data (SGD) can be used to determine the subcrustal stress. To do so, we present two methods for producing the stress components or an equivalent function thereof, the so-called S function, from which the stress components can be computed numerically. First, some integral estimators are presented to integrate the SGD and deliver the stress components and/or the S function. Second, integral equations are constructed for inversion of the SGD to the aforementioned quantities. The kernel functions of the integrals of both approaches are plotted and interpreted. The behaviour of the integral kernels is dependent on the signal and noise spectra in the first approach whilst it does not depend on extra information in the second method. It is shown that recovering the stress from the vertical–vertical gradients, using the integral estimators presented, is suitable, but when using the integral equations the vertical–vertical gradients are recommended for recovering the S function and the vertical–horizontal gradients for the stress components. This study is theoretical and numerical results using synthetic or real data are not given.     

Numerical prediction of ground vibrations induced by high-speed trains including wheel–rail–soil coupled effects

A simplified analytical model including the coupled effects of the wheel–rail–soil system and geometric irregularities of the track is proposed for evaluation of the moving train load. The wheel–rail–soil system is simulated as a series of moving point loads on an Euler–Bernoulli beam resting on a visco-elastic half-space, and the wave-number transform is adopted to derive the 2.5D finite element formulation. The numerical model is validated by published data in the literature. Numerical predictions of ground vibrations by using the proposed method are conducted at a site on the Qin-Shen Line in China.     

Investigating the use of a rational Runge Kutta method for transport modelling

An unconditionally stable explicit time integrator has recently been developed for parabolic systems of equations. This rational Runge Kutta (RRK) method, proposed by Wambecq¹ and Hairer², has been applied by Liu et al.³ to linear heat conduction problems in a time-partitioned solution context. An important practical question is whether the method has application for the solution of (nearly) hyperbolic equations as well.In this paper the RRK method is applied to a nonlinear heat conduction problem, the advection-diffusion equation, and the hyperbolic Buckley-Leverett problem. The method is, indeed, found to be unconditionally stable for the linear heat conduction problem and performs satisfactorily for the nonlinear heat flow case. A heuristic limitation on the utility of RRK for the advection-diffusion equation arises in the Courant number; for the second-order accurate one-step two-stage RRK method, a limiting Courant number of 2 applies. First order upwinding is not as effective when used with RRK as with Euler one-step methods. The method is found to perform poorly for the Buckley-Leverett problem.     

Eulerian-Lagrangian formulation of the equations for groundwater denitrification using bacterial activity

A mathematical model for groundwater denitrification using bacterial activity is presented. The model includes the momentum and mass balance equations for water and nitrogen, substrate and bacteria, and chemical reactions between them. The resulting multiphase, multicomponent, flow and transport governing equations, are coupled and nonlinear. A Eulerian-Lagrangian formulation of the equations is developed. The water and gas flow and transport equations are split into forward advection along characteristics, and a residual at a fixed frame of reference. Discontinuities, sharp fronts and steep gradients of the dependent variables are imposed on the advection mode and solved exactly. It is believed that this novel method will avoid numerical artifacts for the solution of the multiphase flow equations (e.g., upstream permeability) and numerical dispersion for the transport equation.     

An integral-balance nonlinear model to simulate changes in soil moisture,groundwater and surface runoff dynamics at the hillslope scale

We present a system of ordinary differential equations (ODEs) capable of reproducing simultaneously the aggregated behavior of changes in water storage in the hillslope surface, the unsaturated and the saturated soil layers and the channel that drains the hillslope. The system of equations can be viewed as a two-state integral-balance model for soil moisture and groundwater dynamics. Development of the model was motivated by the need for landscape representation through hillslopes and channels organized following stream drainage network topology. Such a representation, with the basic discretization unit of a hillslope, allows ODEs-based simulation of the water transport in a basin. This, in turn, admits the use of highly efficient numerical solvers that enable space–time scaling studies. The goal of this paper is to investigate whether a nonlinear ODE system can effectively replicate observations of water storage in the unsaturated and saturated layers of the soil. Our first finding is that a previously proposed ODE hillslope model, based on readily available data, is capable of reproducing streamflow fluctuations but fails to reproduce the interactions between the surface and subsurface components at the hillslope scale. However, the more complex ODE model that we present in this paper achieves this goal. In our model, fluxes in the soil are described using a Taylor expansion of the underlying storage flux relationship. We tested the model using data collected in the Shale Hills watershed, a 7.9-ha forested site in central Pennsylvania, during an artificial drainage experiment in August 1974 where soil moisture in the unsaturated zone, groundwater dynamics and surface runoff were monitored. The ODE model can be used as an alternative to spatially explicit hillslope models, based on systems of partial differential equations, which require more computational power to resolve fluxes at the hillslope scale. Therefore, it is appropriate to be coupled to runoff routing models to investigate the effect of runoff and its uncertainty propagation across scales. However, this improved performance comes at the expense of introducing two additional parameters that have no obvious physical interpretation. We discuss the implications of this for hydrologic studies across scales.     

Mode coupling analysis and differential rotation in a flow driven by a precessing cylindrical container

We present a theoretical weakly nonlinear analysis of the dynamics of an inviscid flow submitted to both rotation and precession of an unbounded cylindrical container, by considering the coupling of two Kelvin (inertial) waves. The parametric centrifugal instability known for this system is shown to saturate when one expands the Navier–Stokes equation to higher order in the assumed small precession parameter (ratio of precession to rotation frequencies) with the derivation of two coupled Landau equations suitable to describe the dynamics of the modes. It is shown that an azimuthal mean flow with differential rotation is generated by this modes coupling. The time evolution of the associated dynamical system is studied. These theoretical results can be compared with water experiments and also to some numerical simulations where viscosity and finite length effects cannot be neglected.