Author's personal copy Operator-splitting errors in coupled reactive transport codes for transient variably saturated flow and contaminant transport in layered soil profiles D. Jacques a, ⁎ , J. Šimůnek b , D. Mallants a , M.Th. van Genuchten c a Waste and Disposal Department, SCK•CEN, Boeretang 200, B-2400 Mol, Belgium b Department of Environmental Sciences, A135 Bourns Hall, University of California Riverside, 900 University Avenue, Riverside, CA 92521, USA c George E. Brown, Jr. Salinity Laboratory, USDA, ARS, 450 W. Big Spring Rd, Riverside, CA 92507, USA Received 15 June 2005; received in revised form 20 June 2006; accepted 30 June 2006 Available online 17 August 2006 Abstract One possible way of integrating subsurface flow and transport processes with (bio)geochemical reactions is to couple by means of an operator-splitting approach two completely separate codes, one for variably-saturated flow and solute transport and one for equilibrium and kinetic biogeochemical reactions. This paper evaluates the accuracy of the operator-splitting approach for multicomponent systems for typical soil environmental problems involving transient atmospheric boundary conditions (precipitation, evapotranspiration) and layered soil profiles. The recently developed HP1 code was used to solve the coupled transport and chemical equations. For steady-state flow conditions, the accuracy was found to be mainly a function of the adopted spatial discretization and to a lesser extent of the temporal discretization. For transient flow situations, the accuracy depended in a complex manner on grid discretization, time stepping and the main flow conditions (infiltration versus evaporation). Whereas a finer grid size reduced the numerical errors during steady-state flow or the main infiltration periods, the errors sometimes slightly increased (generally less than 50%) when a finer grid size was used during periods with a high evapotranspiration demand (leading to high pressure head gradients near the soil surface). This indicates that operator-splitting errors are most significant during periods with high evaporative boundary conditions. The operator-splitting errors could be decreased by constraining the time step using the performance index (the product of the grid Peclet and Courant numbers) during infiltration, or the maximum time step during evapotranspiration. Several test problems were used to provide guidance for optimal spatial and temporal discretization. © 2006 Elsevier B.V. All rights reserved. Keywords: Coupled reactive transport; Unsaturated flow modelling; Geochemical reactions; Operator-splitting errors Journal of Contaminant Hydrology 88 (2006) 197 – 218 www.elsevier.com/locate/jconhyd ⁎ Corresponding author. E-mail address: djacques@sckcen.be (D. Jacques). 0169-7722/$ - see front matter © 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.jconhyd.2006.06.008