Improved numerical solvers for coupled chemistry-climate model simulations

A key and expensive part of coupled atmospheric chemistry-climate model simulations is the integration of gas phase chemistry, which involves dozens of species and hundreds of reactions. These species and reactions form a highly-coupled network of Differential Equations (DEs). There exists orders of magnitude variability in the lifetimes of the different species present in the atmosphere and so solving these DEs to obtain robust numerical solutions poses a “stiff problem”. With newer models having more species and increased complexity it is now becoming increasingly important to have chemistry solving schemes that reduce time but maintain accuracy.

A sound way to handle stiff systems is by using implicit DE solvers but the computational costs for such solvers are high due to internal iterative algorithms (e.g., Newton-Raphson (NR) methods). Here we propose an approach for implicit DE solvers that improves their convergence speed and robustness with relatively small modification in the code. We achieve this by using Quasi-Newton (QN) methods. We test our approach with numerical experiments on the UK Chemistry and Aerosol (UKCA) model, part of the UK Met Office Unified Model suite, run in both an idealized box-model environment and under realistic 3D atmospheric conditions. The box model tests reveal that the proposed method reduces the time spent in the solver routines significantly, with each QN call costing 27% of a call to the full NR routine. A series of experiments over a range of chemical environments was conducted with the box-model to find the optimal iteration steps to call the QN routine which result in the greatest reduction in the total number of NR iterations whilst minimising the chance of causing instabilities and maintaining solver accuracy. The 3D simulations show that our method for the chemistry solver, speeds up the chemistry routines by around 13%, resulting in a net improvement in overall run-time of the full model by approximately 3% with negligible loss in the accuracy (relative error of order 10^-7) . The QN method also improves the robustness of the solver by significantly reducing (40% ) the number of grid cells which fail to converge hence avoiding unnecessary timestep adjustments.