Efficient operator splitting for modelling transport and transformations of multiple nitrogen species in a variably-saturated subsurface environment

In order to study the potential impacts of nitrogen leaching from agricultural land into the vadose zone and shallow groundwater systems that drain into New Zealand’s oligotrophic lakes, we coupled a multi-species reaction model with an existing 3D variably-saturated finite element flow and transport model, FEMWATER (Lin et al., 1997). The resulting code is called FEMWATER-N (Wang et al., 2003). FEMWATER solves the transient advection-dispersion partial differential equation (PDE) using a Lagrangian-Eulerian method. In FEMWATER-N this transport solver was coupled with a reaction ordinary differential equation (ODE) system solver using an operator splitting approach. Initially the reaction ODE was solved using an iterative exponential Euler method, and the coupled transport-reaction PDE-ODE system was solved using standard sequential iterative operator splitting. Improvements to FEMWATER-N described in this paper include implementing the popular LSODE (Livermore solver for ODEs) code to solve the non-linear reaction ODE system at every mesh node, and a symmetric iterative operator splitting method to solve the coupled advection-dispersion-reaction PDE-ODE system. Solution time was dominated by calls to the transport solver. The symmetric operator split allowed more rapid convergence of the coupled solvers, and, by preconditioning the operator split with an additional reaction call, a single iteration was sufficient to obtain an accurate solution to the coupled equations after the first few time steps. This approach achieved a 47% reduction in CPU time in our example. The symmetric operator splitting method is also free from mass balance error.