On the comparison of numerical methods for the integration of kinetic equations in atmospheric chemistry and transport models

The integration of systems of ordinary differential equations (ODEs) that arise in atmospheric photochemistry is of significant concern to tropospheric and stratospheric chemistry modelers. As a consequence of the stiff nature of these ODE systems, their solution requires a large fraction of the total computational effort in three-dimensional chemical model simulations. Several integration techniques have been proposed and utilized over the years in an attempt to provide computationally efficient, yet accurate, solutions to chemical kinetics ODES. In this work, we present a comparison of some of these techniques and argue that valid comparisons of ODE solvers must take into account the trade-off between solution accuracy and computational efficiency. Misleading comparison results can be obtained by neglecting the fact that any ODE solution method can be made faster or slower by manipulation of the appropriate error tolerances or time steps. Comparisons among ODE solution techniques should therefore attempt to identify which technique can provide the most accurate solution with the least computational effort over the entire range of behavior of each technique. We present here a procedure by which ODE solver comparisons can achieve this goal. Using this methodology, we compare a variety of integration techniques, including methods proposed by Hesstvedt et al. (1978, Int. J. Chem. Kinet. 10, 971–994), Gong and Cho (1993, Atmospheric Environment 27A, 2147–2160), Young and Boris (1977, J. phys. Chem. 81, 2424–2427) and Hindmarsh (1983, In Scientific Computing (edited by Stepleman R. S. et al.), pp. 55–64. North-Holland, Amsterdam). We find that Gear-type solvers such as the Livermore Solver for ordinary differential equations (LSODE) and the sparse-matrix version of LSODE (LSODES) provide the most accurate solution of our test problems with the least computational effort.