HETerogeneous vectorized or Parallel (HETPv1.0)HETP: An updated inorganic heterogeneous chemistry solver for the metastable state NH 4+ –Na + –Ca 2+ –K + –Mg 2+ –SO 42– –NO 3– –Cl – –H 2 O system based on ISORROPIA II

. We describe a new FORTRAN 90Fortran computer program to solve the the system of equations for the NH 4 + –Na + –Ca 2+ – 10 K + –Mg 2+ –SO 4 2– –NO 3 – –Cl – –H 2 O system, based on the algorithms of ISORROPIA II. Specifically, the code solves the system of equations describing the “forward” (gas + aerosol input) metastable state, , but with containing algorithm improvements and corrections. These algorithm changes allow the code to deliver more accurate solution results in formal evaluations of accuracy of the roots of the systems of equations, while reducing processing time in practical applications by about 50 %. The improved solution performance results from several implementation improvements relative to the original ISORROPIA 15 algorithms. These improvements include (i) the use of the ‘interpolate, truncate and project’ (ITP) root–finding approach rather than bisection, (ii) the allowance of search interval endpoints as valid roots at the onset of a search, (iii) the use of a more accurate method to solve polynomial subsystems of equations, (iv) the elimination of negative concentrations during iterative solutions, (v) corrections for mass conservation enforcement, and (vi) several code structure improvements. The new code may be run in either a “vectorization” mode wherein a global convergence criterion is used across multiple tests within 20 the same chemical subspace, or a “by gridpointcase-by-case” mode wherein individual test cases are solved with the same convergence criteria. The latter approach was found to be more efficient on the compiler tested here, but users of the code are recommended to test both options on their own systems. TWe also note that implementation of inorganic chemistry within chemical transport models should take care to retain residual or “free” mass of aerosol species remaining after partitioning, to ensure mass conservation – thehe new code has been constructed to explicitly conserve the input mass for all species considered 25 in the solver, and . The new code is provided as open–sourceopen source FORTRAN 90Fortran shareware.

The aerosols can reside in the crystalline solid phase or exist as an aqueous solution of ions, andions and may be in thermodynamic equilibrium with atmospheric gases.The partitioning of the inorganic species between the solid, gaseousgaseous, and aqueous phase is a complex computational problem, owing to the many nonlinearities involved.The equations describing high concentration (non--ideal) inorganic heterogeneous equilibrium between gases, ions and crystallized solid phases present a system of  equations in  unknowns, (where  is the number of chemical constituents).While these equations may be addressed through searching for roots of polynomials resulting from substitution of equations, the non-ideal nature of the problem manifests as corrections to the equilibrium constants in the equations, known as a (activity coefficients.The) which activity coefficientsin turn depend on concentrations in the condensed phase,increasing the nonlinearity of the system of equations, andequations and requiring the development of special techniques for their solution.
Several solvers have been developed to simulate the thermodynamic partitioning of inorganic species (see Zhang et al., 2000 andPye et al., 2020 for a detailed review of these solvers).AIM2 (Clegg and Pitzer, 1992;Wexler and Clegg, 2002), and GFEMN (Ansari and Pandis, 1999a, b) and UHAERO (Amundson et al., 2006) are considered the most rigorous solvers, in that they attempt to find a global minimum in the Gibbs free energy of the constituents.,However, thbut the downfall of this approach stems from the computational time and operator review required to discriminate between the true global minimum and (potentially many) local minima (Makar et al., 2003).This difficulty has prevented the use of these solvers in three dimensional (3D) chemical transport models (CTMs) to date.However, these models may be used to help determine subsystemssubsystems of equationslocal solution spaces where gas and aerosol partitioning will occur with a smaller number of constituentsand hence describe simplified systems that may be solved with more efficient methods.Inorganic heterogeneous chemistry implementations in CTMschemical transport models have relied on computationally efficient algorithms.These algorithms, which directly solve the system of inorganic heterogeneous chemistry equations by considering the species' chemical potentials within these predetermined subspaces of a smaller numbers of species, hence simplifying and reducing the number of equations and unknowns.The specific subspace to be solved is determined based on the input precursor species, and the ratio(s) of the total available cations to the total available sulfate (see Sect. 2).This approach effectively breaks the larger problem into several separate smaller problems.Solvers that apply this tactic include SCAPE (Kim et al., 1993a,b;Kim and Seinfeld, 1995;Meng et al., 1995), EQUILSOLV-II (Jacobson, 1999), ISORROPIA/ISORROPIA II/ISORROPIA-lite (Nenes et al., 1998;Fountoukis and Nenes, 2007;Kakavas et al., 2022), HETV (Makar et al., 2003) and HETP (presented herein).HETV (HETerogeneous Vectorized) was a vectorized solver (i.e., optimized for vectorized computer architecture) based on the original ISORROPIA algorithms (Nenes et al., 1998), but with numerical improvements related to more accurate evaluation of cubic and quadratic equations whose coefficients may vary by several orders of magnitude, coding structure changes to replace logical IF statements with mathematical equivalents, the elimination of redundant calculations, the replacement of intrinsic functions in activity coefficient calculations by high order Taylor series, and the gathering of similar problems within a single--subsystem for solutionto be solved using a global convergence criteria.
These modifications allowed HETV to perform calculations in 1/38 to 1/89 of the time required for ISORROPIA (v1.0), on a vector supercomputer (the fastest supercomputer architecture at the time the HETV code was created).-Mmore recent supercomputer architectures focus on parallel processing across multiple processors to reduce processing time).In 2007, an update to ISORROPIA was released that included 'crustal' species (Mg 2+ , K + , Ca 2+ ) and sea salt (Na + , Cl -), (referred to as ISORROPIA II (; Fountoukis and Nenes, 2007).More recently,, a simplified (and extended ) version of ISORROPIA II has been developed, (called ISORROPIA--lite.ISORROPIA-lite) that attempts to addresses the metastable state (i.e., it assumes (a supersaturated aqueous solution where crystalline states are ignored, except CaSO4)those subsystems in which liquid water is present), as well as effects of organic aerosols on the partitioning of the inorganic system.ISORROPIA--lite solves the same chemical subspaces as ISORROPIA II, but only for the metastable state option (i.e., efflorescence branch) and uses precalculated binary activity coefficients, resulting in a solver that executes about 35 % faster than ISORROPIA II (Kakavas et al., 2022).
The underlying issue driving the use of a metastable state assumption in regional air quality models for inorganic heterogeneous chemistry solvers is that the presence of water in the aerosol is not only controlled by the inorganic components, but also by other components within a mixed--phase aerosol.In the absence of these additional sources of aerosol water, the "pure" (i.e.only) inorganic aerosol thermodynamics can result in partitioning to the 'stable' aerosol phase as only crystalline solidssalts (no ions) or a mixture of crystalline salts and aqueous ions that are saturated with respect to the crystalline salts., The whereas the presence of the additional sources of aerosol water will ensure that some water is always presentand hence the subsystems of equations that have no water will not be encountered.It has been reported that metastable state aerosols may be 'ubiquitous' in the Earth's atmosphere, existing more than 50 % of the time when the relative humidity is between 45 and 75 % (Rood et al., 1989;Tang et al., 1995); this may be especially true in the case of dissolved impurities such as organic species (Fountoukis et al., 2009).Another issue driving the use of the metastable state assumption in regional air quality models is the need to track the RH history of aerosols to accurately predict their phase state, due to the hysteresis of salts.
Specifically, without knowing the RH history of the aerosol, it is not possible to determine whether the aerosol will exist as an aqueous solution of ions or as a crystalline salt between its efflorescence and deliquescence RH (Martin et al., 2004;Fountoukis et al., 2009).Given these reasons, aApplications of inorganic aerosol thermodynamics within 3D chemical transport modelsCTMs tthus tend to assume a metastable state as the most likely conditions in the troposphere, although absolutely stable aqueous aerosols are possible above the deliquescence RH.This assumption also reduces the number of chemical subspacessubspaces required to obtain a solution of the system of equations for inorganic heterogeneous chemistry, and additions such as formulae for the water activity associated with organic aerosols may be used to better simulate the aerosol water content (Kakavas et al, 2022).
In the different versions of ISORROPIA and HETV, the roots of sub-systemssubsystems of equilibrium equations are used to determine the thermodynamic equilibrium solution, the result being the concentrations of the inorganic ions and the partitioning gases.In ISORROPIA/ISORROPIA II/ISORROPIA-lite and HETV, convergence of these solutions to these systems of equations are obtained via a bisection search, while in SCAPE, Newton's method is employed.Newton's method is also used in ANISORROPIA where it is combined with the bisection method for chemical subspaces describing a neutral aerosol.ANISORROPIA performs a sensitivity analysis on each inorganic species considered in ISORROPIA (excluding Ca 2+ , Mg 2+ , K + ) with respect to the total input precursor species concentration (Capps et al., 2012).It is well known that Newton's method may fail to converge if the 'initial guess' of the root is too far away from the actual root (Burden and Faires, 2011).Unlike Newton's method, the bisection method is guaranteed to converge, al (though the convergence may be slow.) The bisection method, requirequiresring at most   = log 2 ( −

2𝜀
) function evaluations to locate the root () on the interval [, ] such that |  −  * | ≤ , where  is a set tolerance, and  * is the current estimate of the root, and   is the previous estimate of the root.In most cases, the bisection method will require all   function evaluations for convergence (Oliveira and Takahasi, 2021).Recently, Oliveira and Takahasi (2021) developed a modified bisection approach called "interpolate, truncate and project" (ITP), which may obtain superlinear convergence, therefore reducing the execution time required to obtain a solution with the same accuracy as the typical bisection method (note that the bisection method has linear convergence).To achieve an improved order of convergence, the ITP method incorporates a regula--falsi estimate into the bisection method.The 'typical' bisection method simply splits the original interval in half, with  * becoming the midpoint of ).-Tthis estimate defines the 'interpolation' aspect of the ITP method.By making use of these two estimates simultaneously (i.e.,  * and   ), ITP is able tocan outperform the typical bisection method for both convergence rate and accuracy.For well--behaved functions (i.e., there is onlexactlyy one root in the function's domain) ITP requires on average 24 to 37 % of the iterations required by bisection, and for ill--behaved functions (i.e., there are multiple roots in the function's domainn, or the function contains discontinuities) ITP requires on average 82 % of the iterations necessary forrequired by bisection.Oliveira and Takahasi (2021) also compared the performance of ITP against well-established alternative root-finding methods, such as Ridder's method, the Illinois method, Matlab's 'fzero' routine and the Secant method.For all mathematical functions evaluated for convergence, ITP required the least amount of function evaluations when compared against the other root-finding methods.For example, compared to Ridder's method, ITP requires an average of 20.2 function Theevaluations while Ridder's method requires an average of 26.1.The full mathematical details describing the ITP method (as well as pseudocode) are given in Oliveira and Takahasi (2021) and are not repeated herein.
In this work we present HETP (HETerogeneous vectorized or Parallel), a solver based on the "forward" (input precursor species as gas + aerosol) metastable state algorithms of ISORROPIA II.HETP , which hascan been optimized for vector processors where (i.e.similar problems for a subsystem are gathered and solved with a global convergence criterion,) or parallel processors, where (the latter employing local, by grid point case-by-case solutions to the system of equations are used to minimize processing time on parallel processors).HETP focuses exclusively on the metastable state (efflorescence branch) where some amount of liquid water is always assumed to be present in the aerosol, even at very low relative humidity.
; Tthe metastable state assumption is currently applied in various state--of--the--art global and regional chemical transport modelsCTMs, such as GEM--MACH, GEOS--ChemHEM and CMAQ.GEM-MACH uses HETV (Makar et al., 2018), while CMAQ (Wang et al., 2012) and GEOS-ChemHEM (Pye et al., 2009) use ISORROPIA II.HETP has been updated to improve its numerical stability and computational speed compared to ISORROPIA II, as will be discussed in detail below.Specifically, in addition to the numerical improvements associated with its predecessor, HETV, modifications have been made to incorporate base cations and chlorine, to ensure mass conservation, and to update the bisection method to ITP.In the following sections, we demonstrate that the implementation of ITP not only decreases the execution time of the solver, but it can also improve the final convergence of the chemical system by initializing the search with a species concentration (i.e., an initial guess) that is closer to the actual solution being sought (at thermodynamic equilibrium).Thus, we have developed a new solver (HETP) that has improved the accuracy and decreased the execution time compared to the original ISORROPIA II metastable state forward algorithms.Section 2 briefly outlines the background theory underpinning the solver, followed in Sect. 3 by a detailed list of modifications that are unique to HETP (relative to ISORROPIA II).The final sections provide a comprehensive comparison between ISORROPIA II and HETP, in terms of output results and computational speed, both of which are improved in the HETP algorithm.For brevity we will henceforth refer to ISORROPIA II as ISORROPIA in the remainder of this paper.

Background theory
HETP is based on the algorithms of ISORROPIA, which are in turn based on Gibbs free energy minimizations to define subspaces of systems of equations for inorganic heterogeneous chemistry.ISORROPIA solves two types of problems, referred to as the 'forward' or 'reverse' problem.The forward problem requires known input precursor concentrations (total gas + aerosol), along with a relative humidity and air temperature, to predict the equilibrium state.HETP does not consider the reverse problem where the relative humidity, air temperature and aqueous aerosol species concentrations are known (i.e., no gaseous species are included in the input precursor concentrations), and a solution is sought to determine the resulting equilibrium and gas concentrations.For measured data, the reverse problem is typically not recommended since it lacks the inclusion of gas phase speciation in the input, making its predictions highly sensitive to measurement errors.For example, Hennigan et al., (2015) show that a ±10 % measurement error in NH4 + can alter the pH predicted by the reverse mode by more than 1 pH unit.Furthermore, Song et al., (2018) found that the aerosol pH predicted by the reverse mode may result in a bimodal pH distribution; in their study a negative ion balance gave highly acidic conditions while a positive ion balance gave near-neutral conditions.We note that the reverse mode is used in CMAQ to perform mass transfer with the coarse mode (Pye et al., 2022), but other CTMs that employ ISORROPIA use only the forward mode.
The ISORROPIA solvers have been used in a large number of chemical transport modelCTM applications (i.e., ISORROPIA: 1250 citations; ISORROPIA II: 1245 citations), and have been a key component in these models, allowing inorganic heterogeneous chemistry calculations to be carried out in a timely fashion.Here, we build on those solvers, and would like to acknowledge their important contribution to air--quality modelling science.As stated in Sect. 1, HETP assumes a metastable state, (where some liquid water is always present even at low relative humidity).The required input precursor species are the total sulfate (TS, expressed as molar equivalent H2SO4), total ammonium (TA, expressed as molar equivalent NH3), total nitrate (TN, expressed as molar equivalent HNO3), total sodium (TNa, expressed as molar equivalent Na + ), total chloride (TCl, expressed as molar equivalent HCl), total magnesium (TMg, expressed as molar equivalent Mg 2+ ), total potassium (TK, expressed as molar equivalent K + ) and total calcium (TCa, expressed as molar equivalent Ca 2+ ).Units of these net precursor species are mol m -3 air upon input into both ISORROPIA and HETP.For some input conditions ISORROPIA will adjust the input precursor concentrations prior to determining the subroutine that should be entered.Specifically, ISORROPIA will adjust TA and TCl so that they are no less than 1×10 -10 mol m -3 , and if (TNa + TS + TN) < 1×10 -10 mol m - 3 , then ISORROPIA will adjust TNa and TN so that they are no less than 1×10 -10 mol m -3 (note these are applicable only to Branch 3 and 4; see Fig. 1).These adjustments performed within a chemical transport modelCTM result in output speciation that violates mass conservation, since mass is created for TA, TN, TCl and TNa.For example, for 50,000 unique sets of input conditions executing Branch 4 subroutines (i.e., winter input from Sect.4.2), performing these adjustments results in a median of 1.09×10 -3 ug m -3 of TCl being created by the solver.On a relative scale ( output mass input mass × 100 %) this represents a median increase in TCl mass by 42.7 %; for 25 % of these input conditions the relative increase in TCl mass ≥ 4414 %.In a CTM these mass violations would occur at a single timestep, therefore the impact would increase as the simulation progresses.As a resultConsidering these mass violations, ISORROPIA currently used in GEOS-ChemHEM v14.0.0 (GEOS-CHEMChem, 2022) does not perform these mass adjustments.; Iit should be noted that GEOS-CHEM Chem v14.0.0 uses ISORROPIA v2.2 which contains minor bug fixes compared to ISORROPIA II (v2.0).CMAQv5.4 which also uses ISORROPIA v.2.2 (CMAS, 2016; USEPA, 2022), does perform these initial mass adjustments;, however, any output that results from input data that are mass adjusted are flagged.HETP adopts the approach of GEOS-ChemHEM and likewise does not perform these initial mass adjustments.Therefore, ISORROPIA v2.2 used herein (obtained from CMAQv5.4;USEPA, 2022) has been modified so that it also does not perform the aforementioned mass adjustments.Other than this modification, the branches and chemical subspaces (shown in Fig. 1) are identical to ISORROPIA.
Table 1 lists the entire set of equilibrium reactions (ER1 to ER7) that are solved in various chemical subspaces of the metastable state 'forward' option of both ISORROPIA and HETP.The decision tree (outlined at the end of this section) used to select the appropriate chemical subspace, as well as the equilibrium reactions shown in Table 1, are identical to ISORROPIA (Fountoukis and Nenes, 2007).ER1 to ER7 are solved by introducing additional relationships for mass conservation, electroneutrality (i.e., a charge balance equation), aerosol water activity, and mean molality-based activity coefficients () to represent ion--ion interactions in non--ideal aqueous solutions ( → 1 as the solution becomes more dilute, i.e. more "ideal").
Given in Table S1 are the equilibrium reactions that form the basis of dry salt partitioning (ER8 to ER25) that is completed during the initialization of several metastable state subspaces.It should be noted that ER8 to ER25 are not solved directlyinstead the input precursor species are partitioned into various salts based on these equilibrium reactions.
The exact salts that form (i.e., which anions are matched by which cations) depends on the specific chemical subspace that is entered and whether the subspace is 'sulfate rich', 'sulfate super--rich' or 'sulfate poor'; these classifications are determined by the relative amounts of the input cations to the total available sulfate.For example, in CALCP13 (the algorithm branch describing a sulfate poor case with base cations present) calcium, potassium and magnesium first react with the sulfates to produce CaSO4, K2SO4 and MgSO4 respectively, and sodium and chloride react to form NaCl. Any free calcium will then react with nitrate and free chloride to form Ca(NO3)2 and CaCl2 respectively.Next, free magnesium will react with free nitrate and free chloride to form Mg(NO3) and CaCl2, respectively, and then free sodium will then react with free nitrate to form NaNO3. Finally, free potassium will react with free chloride and free nitrate to form KCl and KNO3, respectively.The order of dry salt partitioning in the remaining chemical subspaces, (where applicable,) are provided in Table S2 of the Supplemental Information, and are identical to ISORROPIA (except for CALCL9, discussed in Sect.3).Depending on the amount of anions and cations present for this initial partitioning stage, some of these input components may be in excess of the amount which can be partitioned into salts.This excess mass, beyond that required to create a set of salts, is referred to as the "free" amount of the given component.The salts created in this initial stage of partitioning are then assumed to undergo deliquescence in each of the problems to be solved, resulting in an aqueous phase speciation that is then used as the initial conditions for which a thermodynamic solution is required.In addition to the free amounts generated within a chemical subspace during dry salt partitioning, free amounts may also be generated during the initialization of HETP and ISORROPIA, prior to entering a chemical subspace.Specifically, automatic adjustments are applied if the input precursor species are nonelectroneutral.In this case, any excess cations are ignored, and free amounts of Na, Ca, K and Mg may be created.These automatic adjustments help constrain the particle alkalinity of the equilibrium solution, ensuring that it does not exceed the pH of dissolved particulate calcium carbonate (Pye et al., 2020).The "free" mass must therefore be treated carefully in the context of the application of thermodynamic solvers within chemical transport modelsCTMs.A key requirement for chemical transport modelsCTMs is that they conserve the mass of transported species, within process representation such as inorganic thermodynamics.Solvers such as ISORROPIA conserve mass for the "captured" or "non--free" portion of the input chemical speciation, but not the "free" portion. .However, the "free" mass must be retained by the program accessing the solver, to prevent loss of mass of species such as Na, Mg, K, and Ca; the free mass must be added back to the captured mass partitioned by the solver prior to returning to the program accessing the inorganic heterogeneous chemistry solver.Currently, ISORROPIA only outputs the aqueous, solid, or gaseous species that result after partitioning at thermodynamic equilibrium, and not 'free' amounts.If the non-volatile species (Ca 2+ , Mg 2+ , K + , Na + ) output by the solver are used by the CTM, and the 'free' amounts are not retained and used to conserve mass, then inputs to the solver which result in 'free' species will be lost in the solver call.Some of the chemical subsystem solvers in ISORROPIA retain the free amounts, while others do not; Wwe note, however, thatt the free amounts are not being tracked in someCTMs community regional chemical transport models such as CMAQ v5.4 and GEOS-Chem v14.0.2 avoid this potential problem by only allowing the semi-volatile species employing ISORROPIA (i.e., CMAQv5.4,GEOS-CHEM v14.0.2)(i.e., Cl --HCl, NO3 --HNO3, NH4 + -NH3) to be modified on output from the solver.The semi-volatile species are then saved and transferred back to the model.The non-volatile species are not used after chemical partitioning and are not transferred back to the model calling ISORROPIA.Therefore, any non-volatile free mass that was created in ISORORPIA is not lost in the solver call in these CTMs (aerosol mass is conserved).; these implementations may be inadvertently losing aerosol mass due to this issue, where the free amounts were not being retained and hence inorganic aerosol mass may sometimes inadvertently be lost in these regional models.IIn HETP , the free amounts have been retained in all cases and are returned to the calling code for completeness.The manner in which the initial salt concentrations are determined, including the "free" amounts, is provided in detail in Table S2 (Supplemental information).; HETP tracks all free amounts explicitly, otherwise, the initial dry salt concentrations outlined in Table S2 are determined identical to ISORROPIA (except CALCL9 which is discussed in Sect.3).
Table 1: Equilibrium reactions (ER) considered in the metastable state chemical subspaces, identical to ISORROPIA (Fountoukis and Nenes, 2007).These reactions are solved directly within the appropriate major system.∆  0 , ∆  0 and ∆  0 are the standard molar Gibbs free energy, enthalpy of formation and heat capacity at standard pressure,  = 8.314 J mol -1 K -1 is the universal gas constant, and  0 = 298.15K is the reference temperature.Each species concentration The ions denoted in square brackets […] (i.e., [H + ], [SO4 2-], [HSO4 -], etc.) with unitsrefer to of molalities with units of mol kg -1 mol m -3 is converted to a molality using the aerosol liquid water content in kg m -3 .Here,  is a multicomponent activity coefficient and  is a gas partial pressure.
Theoretically, equilibrium constants are unitless since each pressure or concentration should be normalized by a standard state; here standard states are neglected.

Equation
No.

Equilibrium reactions and values of
with  H + = 1 and  OH − = 1

ER5
HNO 3 : HNO 3 (g) ⇌ H (aq) Note: (Kim and Seinfeld, 1993b) The equilibrium constants are calculated from the Van't Hoff equation, where ∆ 0 ( 0 ) is approximated for a small temperature range (Denbigh, 1981) as where  0 is the equilibrium constant at a reference temperature of  0 = 298.15K,  = 8.314 J mol -1 K -1 is the universal gas constant, ∆  0 (J mol -1 K -1 ) is the change of molar heat capacity of products minus reactants and ∆ 0 (kJ mol -1 ) is the enthalpy change of the reaction at temperature  0 (K). 0 is determined as where ∆  0 (kJ mol -1 ) is the standard molar Gibbs free energy of formation at  0 .
The mean activity coefficients are calculated following the same methodology as in ISORROPIA: multicomponent activity coefficients are calculated according to Bromley's formula (Bromley, 1973), binary activity coefficients are determined from the Kusik--Meissner relationship (Kusik and Meissner, 1978), and the temperature dependence of the multicomponent activity coefficients areis calculated following Meissner and Peppas (1973).HETP (as in ISORROPIA) assumes that OH - (aq) is small compared to other species, and hence it is not used in the calculation of ionic strength.HETP only allows on--line calculation of activity coefficients and does not use precalculated look--up tables.
Aerosol liquid water content in kg m -3 air is calculated according to the Zdanovskii-Stokes-Robinson (ZSR) correlation relation (Robinson and Stokes, 1965), as where   is the concentration of species  in mol m -3 air and   is the molality (mol kg -1 ) of an aqueous solution of  at the same water activity (  ) as the mixturethe water activity (  ) is equal to the fractional relative humidity (0 to 1 scale).It is assumed that there are negligible effects from droplet curvature (i.e., Kelvin effect), and that the growth of an aerosol by uptake of H2O does not affect the ambient water vapor pressure (i.e., no effect on the ambient relative humidityRH).Therefore, equilibrium between the vapor (gas) and liquid (aerosol) phase is assumed with   = RH (Seinfeld and Pandis, 2016).
There are other simplifications and assumptions applied to the metastable state in HETP and ISORROPIA including: (i) sulfuric acid, sodium, magnesium, calcium and potassium are assumed to only exist in the aerosol phase (i.e., no sulfuric acid gas), (ii) calcium sulfate (CaSO4) never dissolves and will only be present as a solid species, (iii) in cases that are sulfate rich (B4, C2, E4, F2, I6, J3, L9, K4), the ions NH4 + , NO3 -and Cl -are "assumed to be minor species that do not significantly perturb the [thermodynamic] equilibrium" (Fountoukis and Nenes, 2007) the partitioning problem to be solved for these ions in sulfate--rich cases is referred to as a "minor system".All minor systems are solved after convergence of the major system has been achieved.Practically, for point (iii) above, this implies that NO3 -and Cl -within the minor system will not affect the charge balance or the activity coefficients of the major system.The concentration of H + determined from the major system is used as the basis to perform the partitioning between the aerosol and gas phase in the minor system(s), (using the equilibrium reaction(s) in Table 1 which describe the minor system(s) to be solved).
The system of equations and order of the operations to create a solution is identical between ISORROPIA and HETP using the same chemical subspaces.The subspace that will be entered, (and therefore the speciation that will be present,) is determined based on the input precursor species.If crustal species (TK, TMg and TCa), TNa and TCl are all near zero, then the set of chemical subspaces reduces to those used in HETV (Makar et al., 2002) and the original release of ISORROPIA (Nenes et al., 1998).Both codes follow the same procedure, first creating three sulfate ratios used to determine the chemical subspace for solution: the "total sulfate ratio" ( 1 ), "crustal species and sodium ratio" ( 2 ) and "crustal species ratio" ( 3 ), These ratios are used as the basis to determine the appropriate chemical subspace that is entered and solved, with 15 possible metastable subspaces in total.The possible subspaces given the input ratios 1, 2 and 3 are summarized in Fig. 1, along with the resulting speciation (aqueous, gaseous, and solid).The bold font species are solved in the major system while regular font species are solved in the minor system.Four unique 'branches' exist: in Branch 1 only TS and TA are present, in Branch  These ratios are used as the basis to determine the appropriate chemical subspace that is entered (15 possible metastable subspaces in total).The possible subspaces (given the input ratios 1, 2 and 3) are summarized in Fig. 1, along with the resulting speciation (aqueous, gaseous and solid).The bold font species are solved in the major system, while regular font species are solved in the minor system.Four unique 'branches' exist: in Branch 1 only TS and TA are present, in Branch 2 only TS, TA and TN are present, in Branch 3 TS, TA and TN are present, and at least one of TNa or TCl, and in Branch 4 TS, TA and TN are present, and at least one of TCa, TK or TMg.The branches are further subdivided into subcases depending on input concentrations.It should be noted that the subcases G5, H6, O7, M8 and P13 require that TCl be present (along with the aforementioned requisite species), otherwise a solution is not possible due to small numbers and floating point arithmetic limitations; this limitation occurs since HETP does not apply the mass modification that sets TCl = max(TCl, 1×10 -10 ), as discussed near the start of the section.

Algorithm design and improvements
During the development of HETP, several improvements related to the mathematical techniques were incorporated relative to ISORROPIA (and HETV), as well as additional modifications related to mass balance.These modifications and improvements include: (1) An updated root finding algorithm, referred to as 'interpolate, truncate and project (ITP)' (Oliveria et al., 2021), has been used instead of the bisection method in HETP.ITP has the advantage of 'superlinear convergence', and hence may obtain a root with the same accuracy as bisection, but in less iterations.The increased rate of convergence can affect the activity coefficients.;Iin some cases, the faster convergence of ITP can alter the ionic strength, resulting in different activity coefficients being calculated early on in the iterative process than would be determined from the bisection algorithm used in ISORROPIA.The new approach may also contribute to an improved formal accuracy performance for estimating the roots, for the same convergence criteria level (see Sect. 4.1).
(2) All bisection subroutines in ISORROPIA employ a root bracketing approach to obtain an initial interval [  ,   ] where (  )(  ) < 0, signifying that a root exists within the interval according to the intermediate value theorem, (assuming a continuous function).We have found that ISORROPIA does not check to determine if either endpoint is a valid root, that is, if (  ) = 0 or (  ) = 0. Instead, ISORROPIA will proceed to the next interval, continuing its search for a root and, potentially locating a different root than expected (the code seeks the smallest positive real root in the case of multiple roots in the search domain), or a slower convergence towards the start or end of the root interval than might otherwise be the case.In HETP we have included a check during the root bracketing stage to identify cases when   or   is a valid root.If an endpoint is a root, then HETP will return since an equilibrium solution has been found.It should be noted that the occurrence of an endpoint as a valid root is extremely rare and hence neglecting this modification will have no effect on most output from the solver, but nonetheless we have included this possibility in HETP for completeness and accuracy.
(3) In some cases that require ITP (or bisection in ISORROPIA) to obtain an equilibrium solution, the independent variable (i.e., ) converges, but the function being evaluated at  (i.e.,  = ()) oscillates between a negative and positive value, and thus || does not converge to zero as expected if  is a root (despite convergence of ).This oscillating behavior of  may indicate (i) that  is a discontinuity, (ii) or that there is significant non-linearity in the partitioning solution, or (iii) that the accepted tolerance on  is too loose for convergence, and hence  is not an accurate solution to the system of equations at the targeted tolerance level for .For all subroutines requiring ITP, HETP will track the species concentrations, activity coefficients and the value of  that are found to minimize || during the iterative process.If after convergence of  it is determined that || is not minimized compared to all earlier iterations, then HETP will 'reset', and instead use the  value, species concentrations and activity coefficients that were found to minimize ||this is chosen as the solution of the system.The effect of this modification on the output from HETP is discussed in Sect 4.2.
(4) In all chemical subspaces, a quadratic equation must be solved for a subsystem of the equations, while in some cases a cubic equation will be solved.Quadratic equations have the form () =  2 +  + , where the solution (i.e., corresponding to ( ) = 0) i, is usually expressed as the standard quadratic formula  = .
has two possible solutions,  1 and  2 , determined by the sign in front of the radical.As identified in Makar et al., (2003) and implemented in the original version of HETV, when the coefficient '' differs by several orders of magnitude from coefficients '' or '', floating--point arithmetic can fail to give an accurate answer for  when using the standard root formula.For example, if √ 2 − 4 ≈ , then addition in the quadratic formula may be problematic since we are subtracting two nearly equal numbers (i.e., ≈ − + ).. To avoid this issue, HETP uses the In HETP (and HETV), a Taylor series expansion of the quadratic formula is used instead, to approximate the root for times when the coefficients '' and '' differ by orders of magnitude (note that  = 1 in all subroutines; formulae were normalized)analytic formula given in Press et al., (2007)  .Care must be taken when applying this formula since the appropriate choice of   1 and   2 depends simultaneously on the chosen solution (i.e.,  1 or  2 ) and the sign of the  coefficient, as described in Table S3 of the supplement.In addition to the analytic formula from Press et al., (2007), HETP also includes code (which is commented out) to solve the quadratic equation using a Taylor series expansion of the quadratic formula.In this code, the Taylor series expansion is only applied when the coefficients '' and '' differ by orders of magnitude, and hence when the numerical precision issues as described above are likely to occur (note that  = 1 in all subroutines; formulae were normalized). .Both methods produce very similar results, but the analytic formula provided by Press et al., (2007) is superior to the Taylor expansion since it provides an exact solution, giving lower error metrics (i.e., Sect.4).For cases where a cubic equation must be solved, HETP will employ an ITP search to obtain an estimate of the smallest positive real root if an exact analytic solution is not possible.The generic formulae describing the exact analytic solution of a cubic polynomial is from Spiegel et al., (2009) and is used in ISORROPIA.It should be noted that the requirement to solve a cubic equation occurs only during the solution procedure of the minor systems of I6, J3, L9 and K4.For example, the call to solve a cubic equation occurs on line 130 of subroutine 'mach_hetp_calchclhno3'.The most recent version of ISORROPIA (i.e., ISORROPIA--lite) did not address these outstanding numerical issues.
(5) During the development of HETP we have identified several cases where a negative ion or gas concentration can be output from ISORROPIA.For example, a negative concentration of NH4 + can occur when solving the minor system NH3(g) + H + (aq) ↔ NH4 + (aq) for thermodynamic equilibrium.In this case, HETP and ISORROPIA will solve a quadratic equation to determine the concentration of ammonia gas (NH3).From the concentration of NH3, the ammonium cation is determined as NH4 + = NH4 + i -NH3, where NH4 + i is the ammonium concentration determined from the major system (see Table S2, Supplemental Information).If partitioning (after solving the quadratic equation) at this stage gives NH3 ≈ NH4 + i, then subtraction of two nearly identical numbers may lead to a floating-point arithmetic error and a final concentration of NH4 + < 0 (in the original ISORROPIA equations).In HETP, negative output is strictly prohibited.To accomplish thisthis, we have utilized max statements that force any negative concentrations to zero, in conjunction with the more accurate evaluation of the quadratic formula (i.e., point 4 above).
(6) In ISORROPIA, the initial dry salt partitioning that is completed at the commencement of chemical subspace L9 may fail to conserve mass for sulfate, ammonium, potassiumpotassium, and sodium, in some cases.In HETP we have slightly modified the initial dry salt partitioning of CALCL9 (see Table S2, Supplemental Information) to ensure mass conservation holds for all cases.;Aany free TA that may result in L9 is assumed to be in the gas phase as NH3, andNH3 and is added back to the final equilibrium solution after convergence of both the major and minor systems.
As discussed in Sect.2, the free amounts of SO4, Na, Mg, K and Ca are explicitly tracked in HETP for all chemical subspaces and returned to the calling code to prevent a loss of mass in the output speciation.
(7) Mass conservation may not hold in ISORROPIA when the input precursor concentrations are near the lower limit for species concentrations, "tiny" (1×10 -20 mol m -3 ), used in the solver.The same lower limit used to bound the input precursor concentrations is also is used throughout ISORROPIA to bound the species concentrations during and after chemical partitioning.In HETP we use the same lower limit as ISORROPIA to bound the input precursor species (i.e., tiny), but during and after partitioning the lower limit for gaseous speciation is reduced to tiny2 = 1×10 -28 mol m -3 .This reduction of the lower limit for gaseous speciation during the iterative process improves mass conservation for the limiting case when the input precursor concentrations are near the lower limit of tiny.
(8) The subroutine 'adjust' performs a post--convergence mass balance adjustment for ammonium, sulfate, nitratenitrate, and chloride, with the goal of ensuring mass conservation holds to machine precision.Specifically, this subroutine checks only for excess mass relative to the input totals.If identified, the excess mass is removed first from the aqueous phase, and then from the solid phase, and finally from the gaseous phase, until no excess remains.However, the mass adjustment of sulfate in ISORROPIA does not include CaSO4 in the mass balance calculations, and therefore in some cases, ISORROPIA will fail to properly conserve mass to machine precision.In HETP we have included CaSO4 in the mass balance adjustment of sulfate.
(9) Improvements to the overall code structure and efficiency include: (a) Use of modern Fortran 90 (compared to Fortran FORTRAN 77 in ISORROPIA), (b) Use of explicit declarations onlyall subroutines now start with an 'implicit none' statement and all common blocks have been removed, (c) Removing all GOTO statements, and instead using Fortran 90modern Fortran constructs such as 'do while' loops, (d) Removing function and short subroutine calls, except for process calls to calculate activity coefficients (calcact), to solve a cubic equation (poly3), to solve minor systems, and to perform a post--convergence mass balance adjustment (adjust).The merging of functions and some short subroutines allowed several variables to be calculated once and reused throughout the iterative process, reducing computational time,i.e.reducing the call factor overhead for individual subroutine calls to the largest extent possible, (e) Moving expressions being recalculated unnecessarily within loops to take place prior to the loop, and removing calculations that serve no purpose to the actual solution being sought, (f) Pre--calculating constant values which are then stored as variables to be used later in the subroutine and, (g) Designing the code to include an optional use of a vectorization--by--grid point approach (Makar, 1995), which may reduce the call factor overhead on some compilers.

Case--by--case comparison
In this section the output from HETP is compared to ISORROPIA for a set of 10,000 artificially generated input 'test cases' that span the domain of each chemical subspace.The test cases have all precursor species held constant except the total sulfate (TS) which is slowly varied (linearly) over the range of the chemical subspace.Tests of this nature demonstrate the stability of numerical solutionsadjacent tests along the same axis of variation in general are expected to be smoothly varying (Makar et al., 2003).The convergence criteria are consistent between both solvers.For activity coefficients,   = 1×10 -6 ^(-6) and   = 4 = 4,, where   is the relative error limit between successive iterations of activity coefficient calculations, and   is the maximum number of allowed iterations.For bisection or ITP,  = 1×10 -9 ^(-9),   = 100 = 100 and  = 5 = 5,, where  is defined in Sect. 1,   is the maximum number of allowed iterations, and  is the number of subdivisions searched for an interval containing a root (i.e., sign change) prior to the start of bisection or ITP.All output from HETP, (in this section and those presented hereafter,) includes the modifications outlined in Sect. 3 unless stated otherwise.,T while the ISORROPIA code used in this comparison is the base version (ISORROPIA v2.2) used in the CMAQ air--quality model (USEPA, 2022).ISORROPIA throughout this paper has been compliled using the '--r8' flag that (converts allall real variables converted to double precision,) to ensure the precision of both solvers is consistent (HETP uses double precision throughout).It should be noted that ISORROPIA is coded to use mostly 'double precision' variables, but some single precision variables exist (i.e., declared as 'real', either explicitly or by default under Fortran variable naming conventions).While compiling ISORROPIA with the intel Fortran compiler flag '-r8' does not have a large impact on the execution time, it may in some cases produce non-trivial differences in the output, compared to output produced without the '-r8' flag.Aside from the '-r8' flag, no other compilation flags were used in this work.All numerical tests herein were executed on a Lenovo ThinkSystem SD650v2 DWC computer, which uses an Intel ® Xeon ® Platinum 8380 CPU running at a clock speed of 2.30 GHz, with 512 GB of available random accessrandom-access memory.The compiler used was an intel compiler (IFORT) version 2021.5.0.2021109.and so they were not included in the original HETV package, which was designed for the SO4-NO3-NH4-H2O system.Furthermore, these two subspaces are frequently called in practical chemical transport modelCTM applications (see Sect. 4.2) and hence are used to compare HETP against ISORROPIA in this section.For the test cases shown in Fig. 2, the relative humidithumidity (RH) y (RH) was set to 35 % and the air temperature () to 306 K, conditions typical of a hot summer day in central North America.The output for CALCO7 is nearly identical between the two solvers, with a difference of < 1 % between HETP (Fig. 2a) and ISORROPIA (Fig. 2b), except for TS between 2.1×10 -5 and 2.4×10 -5 mol m -3 , where visual differences begin to appear,, particularly for H + , HSO4 -and NH3.In the case of CALCM8, the output from HETP (Fig. 2c) is vastly different from ISORROPIA (Fig. 2d), for the same initial conditions and convergence criteria.For these initial conditions, , and the ISORROPIA solution shows the effects of numerical instability in the bisection root--finding procedure.
The ISORROPIA algorithm used in CALCM8 is designed so that the variable being bisected is proportional to Cl -(see Table S2, Supplemental Information).At the same time, the multicomponent activity coefficients are dependent on the ionic strength of the aqueous aerosol, determined from the molar concentration of all ions present, including Cl -.Both of these iterative procedures are completed simultaneously, and impact each other in a nonlinear fashion.The choice of Cl -during the first iteration of bisection (or ITP) may considerably impact the final equilibrium solution, by altering the initial ionic strength, and as a result, the convergence of the multicomponent activity coefficients.This effect is demonstrated for CALCM8 in Fig. 2(c--d), where the differences between ISORROPIA and HETP are related only to the choice of root--finding methodology.In fact, if the ITP approach within HETP is reverted to the same bisection algorithm used in ISORROPIA, then the output from HETP begins to show the same unstable behaviour that is demonstrated in the ISORROPIA simulation shown in Fig. 2d.It should be noted that these differences are due to the choice of root--finding methodology and are not the result of allowing the ends of the interval to be potentially valid roots (i.e., point 2 in Sect.3).
The accuracy of each solver can be assessed directly by introducing an error term (), determined as the absolute logarithmic difference between the 'calculated' equilibrium constant (  ) and the 'true' equilibrium constant (  ), that is,  = log(  ) − log(  ).  is determined from the species concentrations (converted to molalities using the aerosol liquid water content in kg m -3 ) and activity coefficients after convergence of the major or minor system (i.e., from the equations in Table 1), while   is calculated from the Van't Hoff equation (Eq.1).The parameter  thus provides a direct measure of each solver's proximity to the actual root of the system of equations, for a given level of convergence criteria employed in both solvers.For statistical characterization of , the absolute value of the difference is used, so that  ′ = ||.A logarithmic difference is used herein (instead of a percent difference, for example) since the difference between   and   can span several orders of magnitude.In this way, a difference on the order of 1 implies that   and   differ by an order of magnitude different, while a difference on the order of 1×10 -2 implies   and   differ starting at the second or third digit, (when written in scientific notation).The error analysis has been completed using the case--by--case implementation of HETP (see Sect. 4.3).Ideally, ′ = 0, signifying that the problem has converged to a solution whose concentrations and activity coefficients satisfy the equilibrium equations of the major (and minor systems) precisely.In reality, however, there may be some magnitude of difference between   and   .The accuracy of   calculated from Eq. 1 (used in both solvers) is limited to 3 significant digits due to the variable −∆  0 ( 0 ) ⁄ .Therefore when ′ < 1×10 -3 in either solver, we can conclude that   after convergence is identical to   within its known accuracy.However, in practical applications (i.e., within a chemical-transport modelCTM), the value of   calculated from Eq. 1 will retain all digits as determined by the precision of the code (i.e., double precision in HETP) and therefore ′ may be ≪ 1×10 -3 .HenceHence, we seek a solver that obtains ′ as close to zero as possible.Table 2 gives the median, the maximum, and the 25 th and 75 th percentiles of  ′ for HETP and ISORROPIA, corresponding to the data presented in Fig. 2. For CALCM8, the median  ′ is lower in HETP than ISORROPIA for all equilibrium constants, which suggests that HETP is obtaining a more accurate solution for this set of input conditions.The difference in median  ′ between the two solvers is large and indicates that HETP values are more accurate than ISORROPIA by many orders of magnitude, for the same level of the convergence criteria.For example,, i.e., for  HCl , HETP has a median  ′ ≈ 1.77×10 -8 , while ISORROPIA has a median ′ ≈ 0.39, with (similar results are found for  HNO 3 for CALCM8).The superior performance of HETP for this set of initial conditions can also be confirmed visually by comparing Fig. 2(c) to Fig. 2(d).For all species present in this subspace, HETP shows a smooth transition with incremental change in TS, but this is not the case for ISORROPIA.In CALCM8, the very large differences in median ′ between the two codes demonstrates that the ′ values are linked to the poor convergence performance of ISORROPIA, andISORROPIA and are associated with the high degree of sensitivity of that algorithm's use of bisection towards initial conditions.
In CALCO7 (Fig. 2a--b), the median ′ for all equilibrium constants is lower in HETP than ISORROPIA, but the difference between the two solvers is marginal, especially when the 25 th and 75 th percentiles are considered (i.e., for  HCl the 75 th percentile of ′ is 5.174.88×10 - and 4.40×10 -7 for HETP and ISORROPIA respectively).Table 2 also gives statistics of ′ for the same set of input precursor concentrations, but now with a RH a  = 65 % and  = 263 K.The main difference here is that CALCO7 performs slightly worse in HETP than ISORROPIA, (as determined from the median and 75 th percentile of ′).D However, despite this worse statistical performance in HETP, there are no visual differences between them when the output from each solver is plotted (see Fig. S1).In this case, the median ′ of both solvers is on the order of 1.0×10 -4 , implying that the difference between   and   occurs in the 4 th digit (when written in scientific notation).As a result, the differences between HETP and ISORROPIARORPIA do not become apparent unless the graph is zoomed in very close to the data points.For CALCM8 at these new meteorological conditions ( =(RH = 65 % and  = 263 K), HETP has an unstable behavior in the output speciation for TS between 1.6×10 -7 mol m -3 and 2.3×10 -7 mol m -3 , while ISORROPIA has an unstable behavior for all TS > 0.7×10 -7 mol m -3 (see Fig. S1).; Tthis poor performance in CALCM8 for these meteorological conditions is demonstrated in the statistics of ′ shown in Table 2.
These test cases are divided into tTwo tests, were conducted denoted as a high Mg 2+ -Ca 2+ -K + -Na + case (Fig. 3a-c) and a low Mg 2+ -Ca 2+ -K + -Na + case (Fig. 3d-f).th at span the same range of TS and TA, The input conditions used to generate Fig. 3 are summarized in Table 3.It should be noted that in the unaltered version of ISORROPIA, TCl < 1×10 -14 mol m -3 would have necessitated a mass adjustment at the commencement of the solver.In this case, TCl would be reset to a floor value of 1×10 - 10 mol m -3 , thereby creating mass.This adjustment has not been applied here.
In Fig. 3(a-b) and Fig. 3(d-e) the output compares well between HETP and ISORROPIA for the subspaces O7, M8 and P13, with absolute differences typically < 0.1 % and no obvious visual differences between the two solvers.However, in the high Mg 2+ -Ca 2+ -K + -Na + case (Fig. 3(a-b), for the subspaces K4, and particularly L9, there are some noticeable visual differences between the two solvers for the subspaces K4 and particularly L9..The differences in L9 between the two solvers result from (i) the updated methodology within HETP to calculate polynomial roots, (ii) a correction within HETP to the initial dry salt partitioning to ensure mass conservation, and (iii) one less call to calculate activity coefficients in HETP for some test cases, (specifically thosethose test cases that have no convergence of activity coefficients after completing the maximum number of allowed iterations).The largest absolute differences of (i.e., 100 %to -600 %) are in L9, and are predominantly due to (ii), where for some input conditions ISORROPIA creates dry salt mass for TA, TSTS, and TK.Specifically in ISORROPIA, 6.02%, 0.05% and 5.97% of the test input conditions shown in Fig. 3(a-b) create mass for TS, TA, and TK respectively that cannot be attributed to machine precision near the lower limit used in the solver (i.e., Species  − Species  > 9.999×10 -19 mol m -3 ).The median relative mass created for these input conditions is 22.6% for TS, 0.24% for TA and 2.93×10 10 % for TK.
In K4, (ii) is not applicable, so the differences are thus due to (i) and (iii).As demonstrated in Fig 3(a-c), there is a large amount of 'noise' in K4 for TS > 1×10 -5 mol m -3 and TA < 12×10 -6 mol m -3 in ISORROPIA that is not present in HETP.-Tthis 'noise' shows up as speckling in the percent difference plots and is due mainly to (i).If the noise in ISORROPIA is neglected for K4, then the output from ISORROPIA is quite similar to HETP, with differences < 1%.
Figure 3: Regular variation (linear) of the total available sulfate (TS) and the total available ammonium (TA), while holding all other input precursor species constant (see the main text for a description of precursor species that are held constantTable 3 for a summary of the input conditions).The colors in figure panels (a-b) and (d-e) give the amount of gaseous NH3 after chemical partitioning at thermodynamic equilibrium (the color scale is logarithmic and identical in these panels).Panel (c)(c) and (f) show the percent difference of (a-b) and (d-e) respectively.Each panel set (i.e., a-c and d-f) includes 1,000,000 unique input test cases, with the same convergence criteria as used to generate Fig. 2 and 3. Superimposed on each panel are dashed black lines denoting the boundariesy between different chemical subspaces.; Tthe actual subspace contained within a set of dashed lines is given as a text label.
An additional concern identified in Makar et al., (2003) is the potential impact of the inaccurate evaluation of the quadratic and cubic formula (i.e., analytic formulae to obtain an 'exact' solution), which remains present in subsequent iterations of ISORROPIA since the development of HETV (see Sect 3, point 4).An example showing the incremental improvement of the quadratic and cubic solution procedure on the output speciation is displayed in Fig. 4, which depicts the output of CALCI6 from ISORROPIA (Fig. 4b) and HETP (Fig. 4a, c--d).This case illustrates differences that would occur at rather low temperatures and relative humidity, in this case  = 243 K and RH = 5 %.While such a combination of air temperature and relative humidity is likely to be rare in the lower troposphere, it is not uncommon for surface air temperatures to reach 243 K or lower in the winter in Canada, and at similar higher latitudes in other parts of the world.The choice of RH = 5 % here is used to highlight the numerical issues present in ISORROPIA, which occur more frequently and are more pronounced at low RH.However, the numerical issues highlighted in Fig. 4 continue to be present in the output from ISORROPIA, but to a lesser extent, even as the RH is increased to 35 % for the same set of initial conditions.At an ambient RH of 5 % the assumption of a supersaturated aqueous phase may be less justified and is more likely to be representative of a very hypothetical case.Nonetheless, observations from southern California have indicated (although in warmer air temperatures than investigated here) that crystallization of some ambient aerosols may not occur until a RH as low as 4 % (Shaw and Rood, 1990), suggesting that in some atmospheric conditions metastable aerosols are possible even at a very low RH of 5 %.
In Fig. 4a, HETP has been executed without any modifications to improve the accuracy of polynomial root calculations, so that the only improvement over ISORROPIA is that HETP will not allow negative species concentrations (i.e., HSO4 -).In Fig. 4c, HETP now includes an improved methodology to calculate roots of quadratic polynomials ('analytic quad'), in addition to the improvement related to negative species concentrations of Fig. 4a.Lastly, In Fig. 4d, HETP now includes an ITP search to determine the roots of cubic polynomials, in addition to the improvements of Fig. 4a and c. Figure 4 follows the same procedure as Fig. 2 that is, an incremental variation of the input TS while holding all other precursor species constant.This case illustrates differences that would occur at rather low temperatures and relative humidity, in this case  = 243 K and RH = 5%.Without the modifications applied in Fig. 4c and d, the output from HETP and ISORROPIA are quite similar.However, as numerical improvements are incrementally applied to HETP, clear visual differences between HETP and ISORROPIA become apparent for most chemical species in this subspace.In CALCI6, the major system being solved is H + -HSO4 --SO4 2-, requiring a quadratic root with a large variation in coefficient magnitudes to be derivedand therefore an error in H + will propagate through to the minor systems that are solved thereafter (see Table S2, Supplemental Information).It should be noted that the --axis in Fig. 4 is logarithmic, so negative values are not shown in the figure panels.
Nonetheless, there are many instances when ISORROPIA outputs a negative concentration of HSO4 -for this subspace (Fig. 4b), as a result of the use of the standard (and under these circumstances inaccurate) formula for the roots of a quadratic equation for H + in this subspace.In HETP we have included a updated methodologies are used to solve the quadratic equation, that avoid numerical inaccuracies due to catastrophic cancellation.Taylor expansion of the quadratic formula, which is applied when numerical precision is likely to cause erroneous output.The result of this modification (as demonstrated in Fig. 4c) is the removal of the numerical instability present in the output of HETP for this set of initial conditions shown in Fig. 4a.
Numerical instability caused by the erroneous evaluation of the quadratic formula appears to be most prevalent at a low relative humidity (low aerosol water mass).
Following convergence of the major system in CALCI6, the minor systems are solved, one of which requires the roots of a cubic polynomial to be identified; the smallest positive real root determines the concentration of Cl -and NO3 -.In HETP, an ITP search is employed to determine the smallest positive real root of the cubic polynomial when an exact analytic solution from the cubic root formulae is not possible, (due to a large range in the magnitude of the coefficients of the cubic polynomial, which may lead to floating-point arithmetic errors).For the set of input conditions shown in Fig. 4, including an ITP search to solve cubic polynomials results in about 72 % more roots being identified in HETP than in ISORROPIA.If ISORROPIA is unable to determine a valid root from the cubic formula, it will assume that the root is a tiny value (i.e., 1×10 - 20 mol m -3 )this is the procedure that was applied to generate the output shown in Fig. 4a-c.The effect of including an ITP search to solve cubic polynomials is a very large reduction in ′ for  HCl and  HNO 3 in the chemical subspaces I6, J3, L9 and K4 for some sets of initial conditions.(Sstatistics of ′ corresponding to CALCI6 shown in Fig. 4 are given at the bottom of Table 2).For example, in Fig. 4d, HETP has been implemented with an ITP search to solve cubic polynomials, and as shown in Table 2, this implementation leads to a large reduction in the median ′ for  HCl from 13.01 to 1.46×10 -9 .The difference here is a solution that is accurate versus one that is not.The output shown in Fig. 4d demonstrates that including an ITP search to solve cubic polynomials removes discontinuities that occur in Cl -, NO3 -, H + and NH3 near 1.4 mol m -3and hence these species now show a smooth transition over the entire range of TS.HETP has a limiting precision of 1×10 -28 mol m -3 , which is the likely cause of the HSO4 -concentration becoming zero in Fig. 4(c-d) when TS is between about 2.15×10 -12 and 2.4×10 - solves quadratic equations using the 'standard' quadratic formula, and attempts to find an exact analytic solution of cubic equations.All input precursor species are held constant, except the total available sulfate (TS) which is varied over 10,000 sets of initial conditions.The air temperature and relative humidity are 243 K and 5% respectively, for all test cases in the figure.The convergence criteria are consistent between both solvers (see text).

Comparison using input from the GEM--MACH air--quality model
Aside from generating artificial sets of input data to evaluate HETP (Sect.4.1), the value of which is to demonstrate relative solution stability across small increments in input conditions, a comparison between HETP and ISORROPIA can be completed using more realistic input conditions obtained from the GEM--MACH air--quality model (Makar et al., 2018).In this section, percentile of ′ correspond to the bottom of the box, center line in the box, and top of the box respectively.The bottom and top whisker of each box gives the minimum and maximum of ′ respectively; if the bottom whisker extends off the graph, then the minimum ′ is zero.Except for G5:-S, H6:-S, H6:-W, M8:-S and M8:-W, the median ′ of  HNO 3 is smaller in HETP than ISORROPIA for all subspaces shown in Fig. 6b.For  HCl , all subspaces except H6:-S and H6:-W have a smaller median ′ in HETP than ISORROPIA.We note that despite HETP having lower median ′ than ISORROPIA for some subspaces, the magnitude of ′ suggests that ISORROPIA is nevertheless providing sufficiently accurate output for most test cases.For the input data investigated here, the subspace H6 is performing poorly in both solvers, with a (median  ′ > 0.5 for all equilibrium constants) (but with marginally worse, but the performance is marginally worse in HETP than in ISORROPIA).FA (for example, in H6:-S for  HNO 3 the 75 th percentile in HETP is 31.8,andbut in ISORROPIA it is 13.4).H6 is unique relative to the other subspaces requiring a root-finding method (i.e., G5, O7, M8 and P13), since the objective function used to determine the root of the system of equations does not include H + explicitly.The expressions for ′ used in Fig. 6, however, explicitly evaluate the convergence of H + relative to  HCl and  HNO 3 equilibria.The relatively poor performance of the H6 algorithm when evaluated using ′ thus tells us that although the other ions and gases in the H6 chemical subspace have converged with the existing solution procedure, convergence with respect to H + remains poor.
Returning to the scatter noted in HNO3 and /NO3 -between the two solvers in G5:-W (Fig. 5d), it is clear from the statistics of  ′ for  HNO 3 and  HCl shown in Fig. 6 that both solvers are producing output that spans a broad range of accuracy.
The 75 th percentile of  HNO 3 and  HCl are 2 orders of magnitude lower in HETP than ISORROPIA.(Ffor  HCl the 75 th percentile of ′ is 6.93×10 -2 and 4.35 in HETP and ISORROPIA respectively.),However, but the maximum ′ are a similar magnitude in each solver.This suggests that both solvers are struggling with partitioning between the aqueous and gaseous phase for some test cases investigated here.Of the 10,000 test cases analyzed in G5:-W, 14.02 % are identified in HETP as having 'oscillatory behavior' (see Sect. 3, point 3).These flagged test cases generally have large  ′ for all equilibrium constants (in both solvers), which is related to poor convergence during the iterative process.Removing these flagged test cases reduces the median and 75 th percentile of  ′ (for  HNO 3 and  HCl ) by an order of magnitude in both solvers.; Ffor HETP the median  ′ for  HNO 3 reduces to 4.9089×10 -8 (from 4.60×10 -7 ) and for ISORROPIA the median  ′ reduces to 2.72×10 -6 (from 5.59×10 -5 ).The modification to account for 'oscillatory behavior' has the effect of reducing  ′ for the flagged test cases in HETP compared to ISORROPIA -(i.e., for the 14.02 % of test cases affected, the median  ′ for  HNO 3 is 0.28 for HETP, but for ISORROPIA it is 2.65).Furthermore, 98.875.3 % of the flagged test cases are times when [Cl − ]Cl -is predicted to be < 1×10 -16 mol m -3 (Cl -note that [Cl − ] is the bisected variable in G5), and all flagged test cases have TCl < 1×10 -10 mol m -3 .
For test cases where the output from each solver agrees well and (i.e., falls along the one-to-one line in Fig. 5c-d),  ′ for  HNO 3 and  HCl are minimized in each solver.The statistics of  ′ for other subspaces not discussed here are summarized in Table S43 (summer) and Table S54 (winter) of the supplemental information.

Computational time
The mean time (determined from 10 repeated samples) required for the central processing unit (CPU) of a Lenovo SV650v2 DWC computer to solve the test cases from Sect.4.2 (for each season and subspace) are given in Table 54.; Tthe timing tests have an estimated uncertainty of ± 1 %.For HETP, two sets of timing tests are reported.Test 1, labelled 'THETV', refers to timing using a global convergence criteriaglobal convergence criterion for all tests within a given chemical subspace, representing ; a "vectorized" test where all  test cases for a given subspace are solved simultaneously.This is the methodology used in Makar et al., (2003), where the great reduction in processing time associated with vectorization on a vector compiler was used to offset the fact that the number of iterations was determined by the single test case with the worst convergence behavior.Test 2,Test 2, labelled 'THETP', refers to a case--by--case test where the solver is called individually for each test case (i.e., the solver is called  times).In the latter test, the time associated with subroutine calls is offset by the number of iterations becoming test-specifictest specific.The first strategy may be more efficient, (aside from vectorization architecture gains,) when the convergence criteria are relatively similar across grid--cells, that is, all input problems converge with the same number of iterations.Twhile the second strategy may be more efficient when the distribution of convergence is more heterogeneous, with some test cases requiring many more iterations than others.ISORROPIA (TISO) requires a case--by-case implementation, andimplementation and cannot solve  cases simultaneously.The convergence criteria are identical to those used in the previous sections (Sect.4.1 and 4.2).In the case of ISORROPIA, it is important to reaffirm that the '--r8' flag was used during compilation, forcing all calculations to be performed in double precision ((as in the default implementation of HETP) andthat is, the precision of the solver has been removedremoving precision as a possible cause for differences in performance.).For the subspaces D3, G5, H6, O7, M8 and P13 all test cases investigated were chosen so that they require the application of a root--finding method for convergence, since these are the most computationally intensive cases encountered by the solver.As noted above, not all chemical subspaces have 10,000 a sufficient amount of unique input data derived from GEM--MACH simulations for the days sampled from each season.Specifically, in the winter the subspaces A2, B4, C2 and F2 do not have enough suitable input data from which to draw 10,000 unique samples, while iand likewise in the summer for , the subspaces A2, B4 and, C2 .haveinsufficient input data.For winter, input data from J3 are used for F2, except with TNa = 0 and TCl = 0, the aim here being to provide timing tests across a realistic range of initial conditions.It should be noted again that the subspaces A2, B4 and C2 were not executed by GEM--MACH on either day for the reasons noted in Sect.4.2.Therefore, like F2 (winter), the input data used to analyze D3, E4 and F2 are used to analyze A2, B42 and C2 respectively, except with TN = 0.

Table 54:
The average computational time () (calculated from 10 samples) required to solve 10,000 unique sets of input conditions (from summer and winter), using ISORROPIA (TISO), the vectorized solver of HETP (THETV) and the case--by--case solver of HETP (THETP).Input conditions were obtained from the GEM--MACH air--quality model, and the convergence criteria are consistent between both solvers (see text).The speed up is a dimensionless quantity, with the non--bracketed value representing TISO/THETV and the bracketed value representing TISO/THETP.; Aa value > 1 implies that HETP (or HETV) is computationally faster, while a value < 1 implies that ISORROPIA is computationally faster.In the first three columns of each season, the bolded value denotes the fastest execution time between each of the 820 solvers.The bolded value in the speed up column shows which solver style is computationally faster (i.e., HETP or HETV); an underlined value in this column signifies that HETV is computationally slower than ISORROPIA for that subcase (row).The CPU timing results demonstrate that all subspaces (except H6, winter and summer) execute faster in HETP's vectorized THETV implementation than ISORROPIA.-Iin some cases the speed -up is significant: f (i.e., for CALCO7 the speed up is about a factor of 1.745 to 1.835 when using THETV).An even more significantgreater speed up can be achieved by using the case--by--case THETP implementation for some subspaces, specifically those that require bisection (A2, D3, G5, H6, O7, M8 and P13).Unlike THETV, all chemical subspaces execute faster in in THETP than ISORROPIA.For the sets of test cases investigated in this work, the best--case performance is found in P13:-S, where THETV executes in about 0.38 s, but THETP executes in about 0.189 s (the latter being about ~4.32x faster than ISORROPIA).The speed -up afforded by HETP for this subcase is largely the result of HETP's updated root--finding methodology (ITP), which requires fewer iterations on average to obtain a solution with an equivalent (or better) level of accuracy as ISORROPIA.The statistics related to the number of iterations required by the root--finding methodology of each solver to achieve convergence (of the major systems) are given in Table 56, for the same input data used to generate the timing tests shown in Table 45.For P13:-S which has the best--case performance, ITP in HETP requires on average 8.02 iterations for convergence, while bisection in ISORROPIA requires on average 42.5 iterations.Thus, HETP's root--finding method requires about 19 % of the iterations required by ISORROPIA for this set of input conditions, while executing in about 23 4% of the time (using the case--by--case mode).The overall performance for the tests in GEM--MACH (bottom row of Table 45) show the average performance of HETP operating in the case--by--case mode results in a speed up relative to ISORROPIA of a factor of 2.2407x for the summer tests, and 2.1218x for the winter tests.The inclusion of an ITP search for the smallest positive real root of cubic equations in I6, J3, L9 and K4 substantially increases the execution time of the solver for these chemical subspaces relative to no ITP search, but despite this, HETP still executes in less time than ISORROPIA for these subcases.

Subroutine
The difference between THETP and THETV becomes even more apparent, and in favor of THETP, if a significant amount of test cases do not require bisection.While THETV includes a return statement to reorder the problem, (removing those test cases that have converged or have no solution prior to entering ITP), the root bracketing stage in THETV will nonetheless need to be repeated a second time for all test cases that do require ITP.Note that the root bracketing stage identifies an interval where the objective function has a sign change.; Aassuming a continuous function, this sign change signifies that a root exists within the interval.Furthermore, in THETV some test cases may iterate in the root--bracketing stage more times than necessary (i.e., one test case has an identified interval, but other test cases within the same chemical subspace being solved by a global convergence criteriaglobal convergence criterion do not), thereby introducing excess computations into in THETV that do not exist in THETP.This is especially true as the variable , (which controls the number of subdivisions searched for a sign change, ) is increased.Thus, in most applications, and for the computer architecture tested here, the case--by--case THETP implementation will be preferred.Both options are available as separate versions of code, and we recommend users test both options of the code on their own system to determine the best performance.
The results presented herein have demonstrated that HETP is able to provide output for these subspaces that is more accurate overall, while executing up to 4..32x faster than ISORROPIA, with an average performance increase in a practical application between 2.1207x and 2.2418x (using the case--by--case mode).The subspace H6 which executes slower in THETV this interval ( * =  1 2 ⁄ = 0.5( + ));a new interval is then chosen (i.e., [,  * ] or [ * , ]) based on the sign change.The regula--falsi estimate, however, is determined by fitting a straight line through the identified interval by using the function values at each endpoint (i.e.,   = [() − ()] [() − ()] ⁄ 2 only TS, TA and TN are present, in Branch 3 TS, TA and TN are present, and at least one of TNa or TCl, and in Branch 4 TS, TA and TN are present, and at least one of TCa, TK or TMg.The branches are further subdivided into subcases depending on input concentrations.It should be noted that the subcases G5, H6, O7, M8 and P13 require that TCl be present, along with the aforementioned requisite species, otherwise a solution is not possible due to small numbers and floating-point arithmetic limitations.This limitation occurs since HETP does not apply the mass modification that resets TCl to a floor value of 1×10 - 10 mol m -3 , as discussed near the start of the section.

Figure 11 :
Figure 11: Domains of the systems of equations, based on ISORROPIA.For Branch 3, each of TS, TA and TN > tiny, as well as one (or both) of TNa and TCl.For Branch 4, each of TS, TA and TN > tiny, as well as one (or all) of TMg, TK and TMgthus Branch 4 does not necessarily require TNa or TCl > tiny.However, it should be noted that for a solution to be possible, subcases CALCH6, CALCG5, CALCM8, CALCO7 and CALCP13 do require TCl > tiny.The dashed lines in the figure implies that the domain extends infinintely in the direction of increasing R1 or R2. ; Ffor example, in Branch 1, 0 ≤  1 < ∞, but in the figure  1 only extends to 4, and subcase CALCA2 extends for all TA/TS >2.
to solve the quadratic equation:

Figure 22 :
Figure 22: A side--by--side comparison of the output from HETP (left) and ISORROPIA (right), for the chemical subspace CALCO7 (ab) and CALCM8 (c-d).All input species are held constant, except the total available sulfate (TS) which is varied over 10,000 sets of initial conditions.The air temperature and relative humidity are 306 K and 35 % respectively, for all test cases in the figure.The convergence criteria are consistent between both solvers (see text).

Figure 2
Figure 2 displays the output from ISORROPIA and HETP for two example chemical subspaces: (a--b) displays CALCO7 and (c--d) shows CALCM8.These chemical subspaces involve the presence of at least one of Ca 2+ , K + and Mg 2+

Figure 4 :
Figure 4: A side--by--side comparison of the output from HETP (a, c, d) and ISORROPIA (b) for CALCI6.In (a), HETP does not include any methodological improvements to polynomial root calculations.In (c) HETP may apply a Taylor series expansion touses an updated methodology to calculate polynomial roots.In (d), HETP uses an updated methodology to calculate polynomial rootsmay apply a Taylor series expansion to calculate polynomial roots, as well as an ITP search to determine cubic polynomial roots.ISORROPIA shown in (b)

Figure 5 :
Figure 5: A scatter plot of the output concentrations (mol m -3 ) from ISORROPIA (y--axis) compared agaagainstinist HETP (x--axis) for M8:-Summer (M8:S) (a,b) and G5:-Winter (G5:W) (c,d), calculated from 10,000 input test cases obtained from the GEM-MACH air-quality model.The solid black line gives a one--to--one relationship.Speciation is given in the legend shown in panel a.

Figure 6 :
Figure 6: A box and whisker plot of the absolute error  ′ = |log(  ) − log(  )| for (a)  HCl and (b)  HNO 3 .The summer season is denoted by ':-S' and the winter season is denoted by ':-W' in the x-axis labels. ′ is calculated from a set of 10,000 test cases in each season (obtained from the GEM--MACH air--quality model). ′ shown in the figure for M8 and G5 correspond to the scatter plots shown in Fig. 5.The median  ′ is represented by the solid black line in the center of each box, and the 25 th and 75 th percentiles correspond to the bottom 775

Table 2 :
Theoretical error (′) for  = 10,000 generated input conditions corresponding to the chemical subspaces O7, M8 and I6.Statistics of ′ for two sets of atmospheric conditions are presented (temperature,  and relative humidity, RH)RH).The bolded values denotes the smallest median error for that equilibrium constant (i.e., row) between HETP and ISORROPIA.