SolveSAPHE (Solver Suite for Alkalinity-PH Equations)

SolveSAPHE (Solver Suite for Alkalinity-PH Equations) provides a collection of Fortran 90 routines to solve the total alkalinity–pH equation, which relates total alkalinity and pH for a given set of total concentrations of the acid–base systems that contribute to total alkalinity in a given water sample.

This version is a bug fix release of the previous versions. The original package (version 1.0.1) was published as a supplement to the paper:

Guy Munhoven (2013) Mathematics of the total alkalinity-pH equation – pathway to robust and universal solution algorithms:the SolveSAPHE package v1.0.1. Geosci. Model Dev., 6, 1367–1388, 2013 (URL: https://www.geosci-model-dev.net/6/1367/2013 or doi:10.5194/gmd-6-1367-2013).

During extension works on SolveSAPHE, limitations in the bracketing procedure of successive iterates for safeguarding solvers have been made out (in mod_acb_solvers.F90, mod_acbw_solvers.F90, mod_phsolvers.F90 and mod_phsolvers_logging.F90). With the Alk_T - DIC pair used in the published version of SolveSAPHE, these limitations did not have any effect; with DIC replaced by CO₂, the mathematical properties of the underlying fundamental equation change (the carbonate alkalinity term is not bounded any more), which causes problems with Newton iterates when the given CO₂ concentration is high. With this pair, the following shortcomings were diagnosed:

• an iterate obtained by Newton's method could differ from its predecessor by a factor of the exponential of several hundreds, leading to over- or underflow floating point exceptions, which prevented convergence;
• right after an update of the minimum and maximum bounds, the procedure adopted to keep Newton iterates within bounds could in some instances falsely trigger the convergence criterion.

Although we have not observed any similar effects during the stress-tests of SolveSAPHE with the Alk_T - DIC pair, the affected procedures have been amended to avoid these problems from showing up in the future:

• an iterate derived by Newton's method may now differ from its predecessor by a factor of 100 at most, and is rewritten to preclude floating point operation over- and underflows;
• if an iterate derived by Newton's method falls out of bounds, it is replaced by the geometric mean of the currently valid minimum and maximum bounds, and not considering the previous iterate any more for this purpose. The new iterate is then always different from both the minimum and the maximum bound, and can only trigger the convergence criterion if it is sufficiently close to its predecessor (even if this latter is equal to the minimum or the maximum bound).