Abstract. A simple numerical solution procedure – namely the method of lines combined with an off-the-shelf ordinary differential equation (ODE) solver – was shown in previous work to provide efficient, mass-conservative solutions to the pressure-head form of Richards' equation. We implement such a solution in our model openRE. We developed a novel method to quantify the boundary fluxes that reduce water balance errors without negative impacts on model runtimes – the solver flux output method (SFOM). We compare this solution with alternatives, including the classic modified Picard iteration method and the Hydrus 1D model. We reproduce a set of benchmark solutions with all models. We find that Celia's solution has the best water balance, but it can incur significant truncation errors in the simulated boundary fluxes, depending on the time steps used. Our solution has comparable runtimes to Hydrus and better water balance performance (though both models have excellent water balance closure for all the problems we considered). Our solution can be implemented in an interpreted language, such as MATLAB or Python, making use of off-the-shelf ODE solvers. We evaluated alternative SciPy ODE solvers that are available in Python and make practical recommendations about the best way to implement them for Richards' equation. There are two advantages of our approach: (i) the code is concise, making it ideal for teaching purposes; and (ii) the method can be easily extended to represent alternative properties (e.g., novel ways to parameterize the K(ψ) relationship) and processes (e.g., it is straightforward to couple heat or solute transport), making it ideal for testing alternative hypotheses.

Abstract. A simple numerical solution procedure – namely the method of lines combined with an off-the-shelf ordinary differential equation (ODE) solver – was shown in previous work to provide efficient, mass-conservative solutions to the pressure-head form of Richards' equation. We implement such a solution in our model openRE. We developed a novel method to quantify the boundary fluxes that reduce water balance errors without negative impacts on model runtimes – the solver flux output method (SFOM). We compare this solution with alternatives, including the classic modified Picard iteration method and the Hydrus 1D model. We reproduce a set of benchmark solutions with all models. We find that Celia's solution has the best water balance, but it can incur significant truncation errors in the simulated boundary fluxes, depending on the time steps used. Our solution has comparable runtimes to Hydrus and better water balance performance (though both models have excellent water balance closure for all the problems we considered). Our solution can be implemented in an interpreted language, such as MATLAB or Python, making use of off-the-shelf ODE solvers. We evaluated alternative SciPy ODE solvers that are available in Python and make practical recommendations about the best way to implement them for Richards' equation. There are two advantages of our approach: (i) the code is concise, making it ideal for teaching purposes; and (ii) the method can be easily extended to represent alternative properties (e.g., novel ways to parameterize the K(ψ) relationship) and processes (e.g., it is straightforward to couple heat or solute transport), making it ideal for testing alternative hypotheses.

Some systems of differential equations that model problems in science and engineering have natural splittings of the right-hand side into the sum of three parts, in particular, diffusion, reaction, and advection. Implicit-explicit (IMEX) methods treat these three terms with only two numerical methods, and this may not be desirable. Accordingly, this work gives a detailed study of 3-additive linear multi-step methods for the solution of diffusion-reaction-advection systems. Specifically, we construct new 3-additive linear multi-step methods that treat diffusion, reaction, and advection with separate methods. The stability of the new methods is investigated, and the order of convergence is tested numerically. A comparison of the new methods is made with some popular IMEX methods in terms of stability and performance. It is found that the new 3-additive methods have larger stability regions than the IMEX methods tested in some cases and generally outperform in terms of computational efficiency.

Some systems of differential equations that model problems in science and engineering have natural splittings of the right-hand side into the sum of three parts, in particular, diffusion, reaction, and advection. Implicit-explicit (IMEX) methods treat these three terms with only two numerical methods, and this may not be desirable. Accordingly, this work gives a detailed study of 3-additive linear multi-step methods for the solution of diffusion-reaction-advection systems. Specifically, we construct new 3-additive linear multi-step methods that treat diffusion, reaction, and advection with separate methods. The stability of the new methods is investigated, and the order of convergence is tested numerically. A comparison of the new methods is made with some popular IMEX methods in terms of stability and performance. It is found that the new 3-additive methods have larger stability regions than the IMEX methods tested in some cases and generally outperform in terms of computational efficiency.

Fractional-step methods are a popular and powerful divide-and-conquer approach for the numerical solution of differential equations. When the integrators of the fractional steps are Runge--Kutta methods, such methods can be written as generalized additive Runge--Kutta (GARK) methods, and thus the representation and analysis of such methods can be done through the GARK framework. We show how the general Butcher tableau representation and linear stability of such methods are related to the coefficients of the splitting method, the individual sub-integrators, and the order in which they are applied. We use this framework to explain some observations in the literature about fractional-step methods such as the choice of sub-integrators, the order in which they are applied, and the role played by negative splitting coefficients in the stability of the method.

Fractional-step methods are a popular and powerful divide-and-conquer approach for the numerical solution of differential equations. When the integrators of the fractional steps are Runge--Kutta methods, such methods can be written as generalized additive Runge--Kutta (GARK) methods, and thus the representation and analysis of such methods can be done through the GARK framework. We show how the general Butcher tableau representation and linear stability of such methods are related to the coefficients of the splitting method, the individual sub-integrators, and the order in which they are applied. We use this framework to explain some observations in the literature about fractional-step methods such as the choice of sub-integrators, the order in which they are applied, and the role played by negative splitting coefficients in the stability of the method.

Abstract In this study, we manufacture an exact solution for a set of 2D thermochemical mantle convection problems. The derivation begins with the specification of a stream function corresponding to a non‐stationary velocity field. The method of characteristics is then applied to determine an expression for composition consistent with the velocity field. The stream function formulation of the Stokes equation is then applied to solve for temperature. The derivation concludes with the application of the advection‐diffusion equation for temperature to solve for the internal heating rate consistent with the velocity, composition, and temperature solutions. Due to the large number of terms, the internal heating rate is computed using Maple™, and code is also made available in Fortran and Python. Using the method of characteristics allows the compositional transport equation to be solved without the addition of diffusion or source terms. As a result, compositional interfaces remain sharp throughout time and space in the exact solution. The exact solution presented allows for precision testing of thermochemical convection codes for correctness and accuracy.

Abstract In this study, we manufacture an exact solution for a set of 2D thermochemical mantle convection problems. The derivation begins with the specification of a stream function corresponding to a non‐stationary velocity field. The method of characteristics is then applied to determine an expression for composition consistent with the velocity field. The stream function formulation of the Stokes equation is then applied to solve for temperature. The derivation concludes with the application of the advection‐diffusion equation for temperature to solve for the internal heating rate consistent with the velocity, composition, and temperature solutions. Due to the large number of terms, the internal heating rate is computed using Maple™, and code is also made available in Fortran and Python. Using the method of characteristics allows the compositional transport equation to be solved without the addition of diffusion or source terms. As a result, compositional interfaces remain sharp throughout time and space in the exact solution. The exact solution presented allows for precision testing of thermochemical convection codes for correctness and accuracy.

The susceptible-infected-recovered (SIR) model is perhaps the most basic epidemiological model for the evolution of disease spread within a population. Because of its direct representation of fundamental physical quantities, a true solution to an SIR model possesses a number of qualitative properties, such as conservation of the total population or positivity or monotonicity of its constituent populations, that may only be guaranteed to hold numerically under step-size restrictions on the solver. Operator-splitting methods with order greater than two require backward sub-steps in each operator, and the effects of these backward sub-steps on the step-size restrictions for guarantees of qualitative correctness of numerical solutions are not well studied. In this study, we analyze the impact of backward steps on step-size restrictions for guaranteed qualitative properties by applying third- and fourth-order operator-splitting methods to the SIR epidemic model. We find that it is possible to provide step-size restrictions that guarantee qualitative property preservation of the numerical solution despite the negative sub-steps, but care must be taken in the choice of the method. Results such as this open the door for the design and application of high-order operator-splitting methods to other mathematical models in general for which qualitative property preservation is important.

<abstract><p>This paper proposes a new very simply explicitly invertible function to approximate the standard normal cumulative distribution function (CDF). The new function was fit to the standard normal CDF using both MATLAB's Global Optimization Toolbox and the BARON software package. The results of three separate fits are presented in this paper. Each fit was performed across the range $ 0 \leq z \leq 7 $ and achieved a maximum absolute error (MAE) superior to the best MAE reported for previously published very simply explicitly invertible approximations of the standard normal CDF. The best MAE reported from this study is 2.73e–05, which is nearly a factor of five better than the best MAE reported for other published very simply explicitly invertible approximations.</p></abstract>

Many mathematical models of natural phenomena are described by partial differential equations (PDEs) that consist of additive contributions from different physical processes. Classical methods for the numerical solution to such equations are monolithic in that all processes are treated with a single method. Additive numerical methods, in contrast, apply distinct methods to each additive term. There are, however, different ways mathematically to specify the additive terms, and it is not always clear which ways (if any) offer advantages over monolithic methods. This study compares the performance of two different additive splitting techniques (physics-based splitting and dynamic linearization) on a suite of eight test problems that involve advection, reaction, and diffusion with various 2-additive Runge–Kutta methods and (monolithic) Runge–Kutta–Chebyshev (RKC) methods. Results show that dynamic linearization generally outperforms physics-based splitting and so should be preferred as the splitting technique when splitting is required or otherwise desirable. RKC methods are the best performers on three of the eight problems, especially at coarse tolerances, but they can also be prone to severe underperformance.

Despite the proliferation of computer-based research on hydrology and water resources, such research is typically poorly reproducible. Published studies have low reproducibility due to incomplete availability of data and computer code, and a lack of documentation of workflow processes. This leads to a lack of transparency and efficiency because existing code can neither be quality controlled nor reused. Given the commonalities between existing process-based hydrologic models in terms of their required input data and preprocessing steps, open sharing of code can lead to large efficiency gains for the modeling community. Here, we present a model configuration workflow that provides full reproducibility of the resulting model instantiations in a way that separates the model-agnostic preprocessing of specific data sets from the model-specific requirements that models impose on their input files. We use this workflow to create large-domain (global and continental) and local configurations of the Structure for Unifying Multiple Modeling Alternatives (SUMMA) hydrologic model connected to the mizuRoute routing model. These examples show how a relatively complex model setup over a large domain can be organized in a reproducible and structured way that has the potential to accelerate advances in hydrologic modeling for the community as a whole. We provide a tentative blueprint of how community modeling initiatives can be built on top of workflows such as this. We term our workflow the “Community Workflows to Advance Reproducibility in Hydrologic Modeling” (CWARHM; pronounced “swarm”).

Distributed hydrological models predict the spatial variability in processes that govern observed mass and energy fluxes. A challenge associated with the use of these models is the computational burden associated with representing the Earth's (sub)surface via millions of computational elements. This burden is exacerbated as more complex process representations are included because their parameterizations involve computationally intensive mathematical functions. Lookup tables (LUTs) approximate a mathematical function by interpolating precomputed values of the function. Highly accurate approximations are possible for substantially reduced computational costs. In this work, a general methodology using the C++ LUT library FunC is applied to identify and replace computationally intensive mathematical function evaluations in the Canadian Hydrological Model (CHM). The use of LUTs introduces a pointwise relative error below 10 − 8 and provides a reduction in run time of almost 20%. This work shows how LUTs can be implemented with relatively little pain and yield significant computational savings for distributed hydrological models. • The Canadian Hydrological Model (CHM) is profiled and expensive mathematical functions identified. • FunC was used to replace the expensive mathematical functions in CHM with lookup tables. • The run-time performance of CHM was improved by approximately 20% on two realistic simulations. • A general methodology for using FunC to replace expensive mathematical functions with lookup tables is given.

Distributed hydrological models predict the spatial variability in processes that govern observed mass and energy fluxes. A challenge associated with the use of these models is the computational burden associated with representing the Earth's (sub)surface via millions of computational elements. This burden is exacerbated as more complex process representations are included because their parameterizations involve computationally intensive mathematical functions. Lookup tables (LUTs) approximate a mathematical function by interpolating precomputed values of the function. Highly accurate approximations are possible for substantially reduced computational costs. In this work, a general methodology using the C++ LUT library FunC is applied to identify and replace computationally intensive mathematical function evaluations in the Canadian Hydrological Model (CHM). The use of LUTs introduces a pointwise relative error below 10 − 8 and provides a reduction in run time of almost 20%. This work shows how LUTs can be implemented with relatively little pain and yield significant computational savings for distributed hydrological models. • The Canadian Hydrological Model (CHM) is profiled and expensive mathematical functions identified. • FunC was used to replace the expensive mathematical functions in CHM with lookup tables. • The run-time performance of CHM was improved by approximately 20% on two realistic simulations. • A general methodology for using FunC to replace expensive mathematical functions with lookup tables is given.

We describe a method for analyzing the structure of a system of nonlinear integro-differential–algebraic equations (IDAEs) that generalizes the Σ -method for the structural analysis of differential–algebraic equations. The method is based on the sparsity pattern of the IDAE and the ν -smoothing property of a Volterra integral operator. It determines which equations and how many times they need to be differentiated to determine the index, and it reveals the hidden constraints and compatibility conditions in order to prove the existence of a solution. The success of the Σ -method is indicated by the non-singularity of a certain Jacobian matrix. Although it is likely the Σ -method can be directly applied with success to many problems of practical interest, it can fail on some solvable IDAEs. Accordingly, we also present two techniques for addressing these failures.

We describe a method for analyzing the structure of a system of nonlinear integro-differential–algebraic equations (IDAEs) that generalizes the Σ -method for the structural analysis of differential–algebraic equations. The method is based on the sparsity pattern of the IDAE and the ν -smoothing property of a Volterra integral operator. It determines which equations and how many times they need to be differentiated to determine the index, and it reveals the hidden constraints and compatibility conditions in order to prove the existence of a solution. The success of the Σ -method is indicated by the non-singularity of a certain Jacobian matrix. Although it is likely the Σ -method can be directly applied with success to many problems of practical interest, it can fail on some solvable IDAEs. Accordingly, we also present two techniques for addressing these failures.

Abstract The intent of this paper is to encourage improved numerical implementation of land models. Our contributions in this paper are two-fold. First, we present a unified framework to formulate and implement land model equations. We separate the representation of physical processes from their numerical solution, enabling the use of established robust numerical methods to solve the model equations. Second, we introduce a set of synthetic test cases (the laugh tests) to evaluate the numerical implementation of land models. The test cases include storage and transmission of water in soils, lateral sub-surface flow, coupled hydrological and thermodynamic processes in snow, and cryosuction processes in soil. We consider synthetic test cases as “laugh tests” for land models because they provide the most rudimentary test of model capabilities. The laugh tests presented in this paper are all solved with the Structure for Unifying Multiple Modeling Alternatives model (SUMMA) implemented using the SUite of Nonlinear and DIfferential/Algebraic equation Solvers (SUNDIALS). The numerical simulations from SUMMA/SUNDIALS are compared against (1) solutions to the synthetic test cases from other models documented in the peer-reviewed literature; (2) analytical solutions; and (3) observations made in laboratory experiments. In all cases, the numerical simulations are similar to the benchmarks, building confidence in the numerical model implementation. We posit that some land models may have difficulty in solving these benchmark problems. Dedicating more effort to solving synthetic test cases is critical in order to build confidence in the numerical implementation of land models.

Abstract The intent of this paper is to encourage improved numerical implementation of land models. Our contributions in this paper are two-fold. First, we present a unified framework to formulate and implement land model equations. We separate the representation of physical processes from their numerical solution, enabling the use of established robust numerical methods to solve the model equations. Second, we introduce a set of synthetic test cases (the laugh tests) to evaluate the numerical implementation of land models. The test cases include storage and transmission of water in soils, lateral sub-surface flow, coupled hydrological and thermodynamic processes in snow, and cryosuction processes in soil. We consider synthetic test cases as “laugh tests” for land models because they provide the most rudimentary test of model capabilities. The laugh tests presented in this paper are all solved with the Structure for Unifying Multiple Modeling Alternatives model (SUMMA) implemented using the SUite of Nonlinear and DIfferential/Algebraic equation Solvers (SUNDIALS). The numerical simulations from SUMMA/SUNDIALS are compared against (1) solutions to the synthetic test cases from other models documented in the peer-reviewed literature; (2) analytical solutions; and (3) observations made in laboratory experiments. In all cases, the numerical simulations are similar to the benchmarks, building confidence in the numerical model implementation. We posit that some land models may have difficulty in solving these benchmark problems. Dedicating more effort to solving synthetic test cases is critical in order to build confidence in the numerical implementation of land models.

The hydrography of the Canadian Prairies and adjacent northern US Great Plains is unusual in that the landscape is flat and recently formed due to the effects of pleistocene glaciation and a semi-arid climate since holocene deglaciation. Therefore, there has not been sufficient energy, time, or runoff water to carve typical dendritic surface water drainage networks in many locations. In these regions, runoff is often detented and sometimes stored by the millions of depressions (known locally as “potholes” or “sloughs”) that cover the landscape.

The hydrography of the Canadian Prairies and adjacent northern US Great Plains is unusual in that the landscape is flat and recently formed due to the effects of pleistocene glaciation and a semi-arid climate since holocene deglaciation. Therefore, there has not been sufficient energy, time, or runoff water to carve typical dendritic surface water drainage networks in many locations. In these regions, runoff is often detented and sometimes stored by the millions of depressions (known locally as “potholes” or “sloughs”) that cover the landscape.

Blowing snow is ubiquitous in cold, windswept environments. In some regions, blowing snow sublimation losses can ablate a notable fraction of the seasonal snowfall. It is advantageous to predict alpine snow regimes at the spatial scale of snowdrifts (≈1 to 100 m) because of the role of snow redistribution in governing the duration and volume of snowmelt. However, blowing snow processes are often neglected due to computational costs. Here, a three‐dimensional blowing snow model is presented that is spatially discretized using a variable resolution unstructured mesh. This represents the heterogeneity of the surface explicitly yet, for the case study reported, gained a 62% reduction in computational elements versus a fixed‐resolution mesh and resulted in a 44% reduction in total runtime. The model was evaluated for a subarctic mountain basin using transects of measured snow water equivalent (SWE) in a tundra valley. Including blowing snow processes improved the prediction of SWE by capturing inner‐annual snowdrift formation, more than halved the total mean bias error, and increased the coefficient of variation of SWE from 0.04 to 0.31 better matching the observed CV (0.41). The use of a variable resolution mesh did not dramatically degrade the model performance. Comparison with a constant resolution mesh showed a similar CV and RMSE as the variable resolution mesh. The constant resolution mesh had a smaller mean bias error. A sensitivity analysis showed that snowdrift locations and immediate up‐wind sources of blowing snow are the most sensitive areas of the landscape to wind speed variations.

BACOLI is a Fortran software package for solving one-dimensional parabolic partial differential equations (PDEs) with separated boundary conditions by B-spline adaptive collocation methods. A distinguishing feature of BACOLI is its ability to estimate and control error and correspondingly adapt meshes in both space and time. Many models of scientific interest, however, can be formulated as multiscale parabolic PDE systems, that is, models that couple a system of parabolic PDEs describing dynamics on a global scale with a system of ordinary differential equations describing dynamics on a local scale. This article describes the Fortran software eBACOLI, the extension of BACOLI to solve such multiscale models. The performance of the extended software is demonstrated to be statistically equivalent to the original for purely parabolic PDE systems. Results from eBACOLI are given for various multiscale models from the extended problem class considered.

The speed of mathematical function evaluations can significantly contribute to the overall performance of numerical simulations. Two common approaches to evaluate a mathematical function are by dir...

Abstract Unstructured triangular meshes are an efficient and effective landscape representation that are suitable for use in distributed hydrological and land surface models. Their variable spatial resolution provides similar spatial performance to high-resolution structured grids while using only a fraction of the number of elements. Many existing triangulation methods either sacrifice triangle quality to introduce variable resolution or maintain well-formed uniform meshes at the expense of variable triangle resolution. They are also generally constructed to only fulfil topographic constraints. However, distributed hydrological and land surface models require triangles of varying resolution to provide landscape representations that accurately represent the spatial heterogeneity of driving meteorology, physical parameters and process operation in the simulation domain. As such, mesh generators need to constrain the unstructured mesh to not only topography but to other important surface and sub-surface features. This work presents novel multi-objective unstructured mesh generation software that allows mesh generation to be constrained to an arbitrary number of important features while maintaining a variable spatial resolution. Triangle quality is supported as well as a smooth gradation from small to large triangles. Including these additional constraints results in a better representation of spatial heterogeneity than from classic topography-only constraints.