the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
HETerogeneous vectorized or Parallel (HETPv1.0): an updated inorganic heterogeneous chemistry solver for the metastablestate NH_{4}^{+}–Na^{+}–Ca^{2+}–K^{+}–Mg^{2+}–SO_{4}^{2−}–NO_{3}^{−}–Cl^{−}–H_{2}O system based on ISORROPIA II
Paul A. Makar
Colin J. Lee
We describe a new Fortran computer program to solve the system of equations for the NH${}_{\mathrm{4}}^{+}$–Na^{+}–Ca^{2+}–K^{+}–Mg^{2+}–SO${}_{\mathrm{4}}^{\mathrm{2}}$–NO${}_{\mathrm{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 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 algorithms. These improvements include (i) the use of the “interpolate, truncate and project” (ITP) rootfinding 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 the same chemical subspace or a “by casebycase” 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. The new code has been constructed to explicitly conserve the input mass for all species considered in the solver and is provided as opensource Fortran shareware.
 Article
(3608 KB)  Fulltext XML

Supplement
(840 KB)  BibTeX
 EndNote
Anthropogenic atmospheric particulate matter (aerosols) can negatively impact the Earth's climate and biosphere – aerosols can alter the atmosphere's radiative forcing (Jacobson, 2001; Schmale et al., 2021), contribute to acid rain (Irwin and Williams, 1988), reduce atmospheric visibility (Quan et al., 2015), and cause morbidity in humans (Atkinson et al., 2014) and other plant and animal species (Lovett et al., 2009). Atmospheric particulate matter is comprised of organic and inorganic species, with 25 % to 60 % of particulate matter being inorganic by mass (Harrison and Pio, 1983; Heintzenberg, 1989). The inorganic portion of atmospheric particulate matter consists primarily of sulfate (SO${}_{\mathrm{4}}^{\mathrm{2}}$), nitrate (NO${}_{\mathrm{3}}^{}$), ammonium (NH${}_{\mathrm{4}}^{+}$), chloride (Cl^{−}), calcium (Ca^{2+}), potassium (K^{+}), magnesium (Mg^{2+}), sodium (Na^{+}), and water (H_{2}O) (Harrison and Pio, 1983; Wang et al., 2003). Along coastlines and within marine air masses, inorganic bromide (Br^{−}) (Sander et al., 2003) and iodide (I^{−}) (SaizLopez et al., 2011) may also be common. Ca^{2+}, K^{+}, Mg^{2+}, Na^{+}, and Cl^{−} exist principally in coarsemode aerosols (particle diameter >2.5 µm), and these species are particularly important to the partitioning of ammonium and nitrate (Metzger et al., 2006). As an example, coarsemode particle nitrate may form via adsorption of nitric acid (HNO_{3}) onto sea salt (Savoie and Prospero, 1982). It should be noted that a considerable amount of K^{+} may also be present in finemode aerosols (particle diameter <2.5 µm) when it is generated during biomass burning events, termed “pyrogenic potassium” (Metzger et al., 2006). The transfer of cation and anion mass between the gas and particulate phase is crucially dependent on inorganic thermodynamic partitioning. For example, observations have indicated that base cations (Ca^{2+}, K^{+}, Mg^{2+}, Na^{+}) and NH${}_{\mathrm{4}}^{+}$ can compete for uptake of HNO_{3} (the former residing in coarsemode and the latter in finemode particle nitrate formation) (Makar et al., 1998; Anlauf et al., 2006).
The aerosols can reside in the crystalline solid phase or exist as an aqueous solution of ions and may be in thermodynamic equilibrium with atmospheric gases. The partitioning of the inorganic species between the solid, gaseous, and aqueous phase is a complex computational problem, owing to the many nonlinearities involved. The equations describing highconcentration (nonideal) inorganic heterogeneous equilibrium between gases, ions, and crystallized solid phases present a system of N equations in N unknowns, where N is the number of chemical constituents. While these equations may be addressed through searching for roots of polynomials resulting from substitution of equations, the nonideal nature of the problem manifests as corrections to the equilibrium constants in the equations, known as activity coefficients. The activity coefficients depend on concentrations in the condensed phase, increasing the nonlinearity of the system of equations 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; Pye et al., 2020, for a detailed review of these solvers). AIM2 (Clegg and Pitzer, 1992; Wexler and Clegg, 2002), 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, 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 threedimensional (3D) chemical transport models (CTMs) to date. However, these models may be used to help determine subsystems of equations – local solution spaces where gas and aerosol partitioning will occur with a smaller number of constituents – and hence describe simplified systems that may be solved with more efficient methods. Inorganic heterogeneous chemistry implementations in CTMs have relied on computationally efficient algorithms. These algorithms 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), EQUILSOLVII (Jacobson, 1999), ISORROPIA/ISORROPIA II/ISORROPIAlite (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 highorder Taylor series, and the gathering of similar problems within a single subsystem to be solved using global convergence criteria. These modifications allowed HETV to perform calculations in $\mathrm{1}/\mathrm{38}$ to $\mathrm{1}/\mathrm{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). More 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 ISORROPIAlite. ISORROPIAlite addresses the metastable state (i.e., it assumes a supersaturated aqueous solution where crystalline states are ignored, except CaSO_{4}), as well as effects of organic aerosols on the partitioning of the inorganic system. ISORROPIAlite solves the same chemical subspaces as ISORROPIA II but only for the metastablestate option 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 metastablestate assumption in regional airquality 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 mixedphase 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 salts (no ions) or a mixture of crystalline salts and aqueous ions that are saturated with respect to the crystalline salts. The presence of the additional sources of aerosol water will ensure that some water is always present – and hence the subsystems of equations that have no water will not be encountered. It has been reported that metastablestate 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 metastablestate assumption in regional airquality 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, applications of inorganic aerosol thermodynamics within CTMs 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 subspaces 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 subsystems 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, and ISORROPIAlite and HETV, convergence of these solutions to these systems of equations is 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, although the convergence may be slow. The bisection method requires at most ${n}_{\mathrm{eval}}={\mathrm{log}}_{\mathrm{2}}\left(\frac{ba}{\mathrm{2}\mathit{\epsilon}}\right)$ function evaluations to locate the root (x) on the interval [a, b] such that $\left{x}_{i}{x}^{\ast}\right\le \mathit{\epsilon}$, where ε is a set tolerance, x^{∗} is the current estimate of the root, and x_{i} is the previous estimate of the root. In most cases, the bisection method will require all n_{eval} 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 x^{∗} becoming the midpoint of this interval (${x}^{\ast}={x}_{\mathrm{1}/\mathrm{2}}=\mathrm{0.5}(a+b)$); a new interval is then chosen (i.e., $[a,{x}^{\ast}]$ or $[{x}^{\ast},b]$) based on the sign change. The regula falsi estimate, however, is determined by fitting a straight line through the identified interval using the function values at each endpoint (i.e., ${x}_{\mathrm{f}}=\left[bf\left(a\right)af\left(b\right)\right]/\left[f\left(a\right)f\left(b\right)\right])$. This estimate defines the “interpolation” aspect of the ITP method. By making use of these two estimates simultaneously (i.e., x^{∗} and x_{f}), ITP can outperform the typical bisection method for both convergence rate and accuracy. For wellbehaved functions (i.e., there is exactly one root in the function's domain) ITP requires on average 24 % to 37 % of the iterations required by bisection, and for illbehaved functions (i.e., there are multiple roots in the function's domain or the function contains discontinuities) ITP requires on average 82 % of the iterations required by bisection. Oliveira and Takahasi (2021) also compared the performance of ITP against wellestablished alternative rootfinding 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 lowest number of function evaluations when compared against the other rootfinding methods. For example, compared to Ridder's method, ITP requires an average of 20.2 function evaluations, 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) metastablestate algorithms of ISORROPIA II. HETP has been optimized for vector processors where similar problems for a subsystem are gathered and solved with a global convergence criterion or parallel processors, where local casebycase solutions to the system of equations are used to minimize processing time. HETP focuses exclusively on the metastable state where some amount of liquid water is always assumed to be present in the aerosol, even at very low relative humidity. The metastablestate assumption is currently applied in various stateoftheart global and regional CTMs, such as GEMMACH, GEOSChem, and CMAQ. GEMMACH uses HETV (Makar et al., 2018), while CMAQ (Wang et al., 2012) and GEOSChem (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 metastablestate 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.
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 gasphase 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 NH${}_{\mathrm{4}}^{+}$ 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 nearneutral conditions. We note that the reverse mode is used in CMAQ to perform mass transfer with the coarse mode (Pye et al., 2020), but other CTMs that employ ISORROPIA use only the forward mode.
The ISORROPIA solvers have been used in a large number of CTM 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 airquality modeling 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 total sulfate (TS; expressed as molar equivalent H_{2}SO_{4}), total ammonium (TA; expressed as molar equivalent NH_{3}), total nitrate (TN; expressed as molar equivalent HNO_{3}), 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 $\mathrm{1}\times {\mathrm{10}}^{\mathrm{10}}$ mol m^{−3}, and if (TNa + TS + TN) $<\mathrm{1}\times {\mathrm{10}}^{\mathrm{10}}$ mol m^{−3}, then ISORROPIA will adjust TNa and TN so that they are no less than $\mathrm{1}\times {\mathrm{10}}^{\mathrm{10}}$ mol m^{−3} (note these are applicable only to Branch 3 and 4; see Fig. 1). These adjustments performed within a CTM 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 $\mathrm{1.09}\times {\mathrm{10}}^{\mathrm{3}}$ µg m^{−3} of TCl being created by the solver. On a relative scale $\left(\frac{\text{output mass}}{\text{input mass}}\times \mathrm{100}\phantom{\rule{0.125em}{0ex}}\mathit{\%}\right)$ 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 time step; therefore the impact would increase as the simulation progresses. Considering these mass violations, ISORROPIA currently used in GEOSChem v14.0.0 (GEOSChem, 2022) does not perform these mass adjustments. It should be noted that GEOSChem 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 is flagged. HETP adopts the approach of GEOSChem 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 metastablestate “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, is 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 molalitybased activity coefficients (γ) to represent ion–ion interactions in nonideal 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 metastablestate subspaces. It should be noted that ER8 to ER25 are not solved directly – instead 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) depend on the specific chemical subspace that is entered and whether the subspace is “sulfate rich”, “sulfate superrich”, 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 CaSO_{4}, K_{2}SO_{4}, and MgSO_{4} respectively, and sodium and chloride react to form NaCl. Any free calcium will then react with nitrate and free chloride to form Ca(NO_{3})_{2} and CaCl_{2} respectively. Next, free magnesium will react with free nitrate and free chloride to form Mg(NO_{3}) and CaCl_{2}, respectively, and then free sodium will react with free nitrate to form NaNO_{3}. Finally, free potassium will react with free chloride and free nitrate to form KCl and KNO_{3}, respectively. The order of dry salt partitioning in the remaining chemical subspaces, where applicable, is provided in Table S2 of the Supplement and is 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 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 CTMs. A key requirement for CTMs 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 “nonfree” portion of the input chemical speciation but not the free portion. Currently, ISORROPIA only outputs the aqueous, solid, or gaseous species that result after partitioning at thermodynamic equilibrium and not free amounts. If the nonvolatile 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. We note, however, that CTMs such as CMAQ v5.4 and GEOS–Chem v14.0.2 avoid this potential problem by only allowing the semivolatile species (i.e., Cl^{−}–HCl, NO${}_{\mathrm{3}}^{}$–HNO_{3}, NH${}_{\mathrm{4}}^{+}$–NH_{3}) to be modified on output from the solver. The semivolatile species are then saved and transferred back to the model. The nonvolatile species are not used after chemical partitioning and are not transferred back to the model calling ISORROPIA. Therefore, any nonvolatile free mass that was created in ISORROPIA is not lost in the solver call in these CTMs (aerosol mass is conserved). In 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. HETP tracks all free amounts explicitly, otherwise the initial dry salt concentrations outlined in Table S2 are determined to be identical to those in ISORROPIA (except CALCL9 which is discussed in Sect. 3).
Note that $\frac{{\mathit{\gamma}}_{{\mathrm{H}}^{+}}{\mathit{\gamma}}_{{\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}}}}{{\mathit{\gamma}}_{{\mathrm{HSO}}_{\mathrm{4}}^{}}}=\frac{{\mathit{\gamma}}_{{\mathrm{H}}^{+}}^{\mathrm{2}}{\mathit{\gamma}}_{{\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}}}}{{\mathit{\gamma}}_{{\mathrm{H}}^{+}}{\mathit{\gamma}}_{{\mathrm{HSO}}_{\mathrm{4}}^{}}}=\frac{{\mathit{\gamma}}_{{\mathrm{H}}_{\mathrm{2}}{\mathrm{SO}}_{\mathrm{4}}}^{\mathrm{3}}}{{\mathit{\gamma}}_{\mathrm{H}{\mathrm{HSO}}_{\mathrm{4}}}^{\mathrm{2}}}$ (Kim and Seinfeld, 1993b).
The equilibrium constants are calculated from the Van't Hoff equation, where ΔH^{0}(T_{0}) is approximated for a small temperature range (Denbigh, 1981) as
where K_{0} is the equilibrium constant at a reference temperature of T_{0}=298.15 K, R=8.314 J mol^{−1} K^{−1} is the universal gas constant, $\mathrm{\Delta}{C}_{p}^{\mathrm{0}}$ (J mol^{−1} K^{−1}) is the change of molar heat capacity of products minus reactants, and ΔH^{0} (J mol^{−1}) is the enthalpy change of the reaction at temperature T_{0} (K). K_{0} is determined as
where $\mathrm{\Delta}{G}_{\mathrm{f}}^{\mathrm{0}}$ (J mol^{−1}) is the standard molar Gibbs free energy of formation at T_{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 is calculated following Meissner and Peppas (1973). HETP (as in ISORROPIA) assumes that OH${}_{\left(\mathrm{aq}\right)}^{}$ is small compared to other species, and hence it is not used in the calculation of ionic strength. HETP only allows online calculation of activity coefficients and does not use precalculated lookup tables.
Aerosol liquid water content in kg m^{−3} air is calculated according to the Zdanovskii–Stokes–Robinson (ZSR) relation (Robinson and Stokes, 1965), as
where M_{i} is the concentration of species i in mol m^{−3} air, and m_{i} is the molality (mol kg^{−1}) of an aqueous solution of i at the same water activity (a_{w}) as the mixture. 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 H_{2}O does not affect the ambient water vapor pressure (i.e., no effect on the ambient RH). Therefore, equilibrium between the vapor (gas) and liquid (aerosol) phase is assumed with a_{w}= RH (Seinfeld and Pandis, 2016).
There are other simplifications and assumptions applied to the metastable state in HETP and ISORROPIA, including the following: (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 (CaSO_{4}) 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 NH${}_{\mathrm{4}}^{+}$, NO${}_{\mathrm{3}}^{}$, 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 sulfaterich 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 NO${}_{\mathrm{3}}^{}$ 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 are 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., 2003) and the original release of ISORROPIA (Nenes et al., 1998). Both codes follow the same procedure, first creating three sulfate ratios: the “total sulfate ratio” (R_{1}), “crustal species and sodium ratio” (R_{2}), and “crustal species ratio” (R_{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 R_{1}, R_{2}, and R_{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 floatingpoint arithmetic limitations. This limitation occurs since HETP does not apply the mass modification that resets TCl to a floor value of $\mathrm{1}\times {\mathrm{10}}^{\mathrm{10}}$ mol m^{−3}, as discussed near the start of the section.
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 the following:

An updated root finding algorithm, referred to as “interpolate, truncate and project (ITP)” (Oliveira and Takahashi, 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 fewer iterations. The increased rate of convergence can affect the activity coefficients. In 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).

All bisection subroutines in ISORROPIA employ a rootbracketing approach to obtain an initial interval [x_{a},x_{b}], where f(x_{a})f(x_{b})<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 f(x_{a})=0 or f(x_{b})=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 rootbracketing stage to identify cases when x_{a} or x_{b} 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.

In some cases that require ITP (or bisection in ISORROPIA) to obtain an equilibrium solution, the independent variable (i.e., x) converges, but the function being evaluated at x (i.e., y=f(x)) oscillates between a negative and positive value, and thus y does not converge to zero as expected if x is a root (despite convergence of x). This oscillating behavior of y may indicate (i) that x is a discontinuity; (ii) that there is significant nonlinearity in the partitioning solution; or (iii) that the accepted tolerance on x is too loose for convergence, and hence x is not an accurate solution to the system of equations at the targeted tolerance level for x. For all subroutines requiring ITP, HETP will track the species concentrations, activity coefficients, and value of x that are found to minimize y during the iterative process. If after convergence of x it is determined that y is not minimized compared to all earlier iterations, then HETP will “reset” and instead use the x value, species concentrations, and activity coefficients that were found to minimize y – 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.

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 $f\left(x\right)=a{x}^{\mathrm{2}}+bx+c$, where the solution corresponding to f(x)=0 is usually expressed as the standard quadratic formula $x=\frac{b\pm \sqrt{{b}^{\mathrm{2}}\mathrm{4}ac}}{\mathrm{2}a}$. x has two possible solutions, x_{1} and x_{2}, determined by the sign in front of the radical. As identified in Makar et al. (2003) in the original version of HETV, when the coefficient “b” differs by several orders of magnitude from coefficients “a” or “c”, floatingpoint arithmetic can fail to give an accurate answer for x when using the standard root formula. For example, if $\sqrt{{b}^{\mathrm{2}}\mathrm{4}ac}\approx b$, then addition in the quadratic formula may be problematic since we are subtracting two nearly equal numbers (i.e., $\approx b+b)$. To avoid this issue, HETP uses the analytic formula given in Press et al. (2007) to solve the quadratic equation: $q=\frac{\mathrm{1}}{\mathrm{2}}\left(b+\text{sign}\left(b\right)\sqrt{{b}^{\mathrm{2}}\mathrm{4}ac}\right)$, with roots ${x}_{{p}_{\mathrm{1}}}=\frac{c}{q}$ and ${x}_{{p}_{\mathrm{2}}}=\frac{q}{a}$. Care must be taken when applying this formula since the appropriate choice of ${x}_{{p}_{\mathrm{1}}}$ and ${x}_{{p}_{\mathrm{2}}}$ depends simultaneously on the chosen solution (i.e., x_{1} or x_{2}) and the sign of the b coefficient, as described in Table S3. 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 “b” and “c” differ by orders of magnitude and hence when the numerical precision issues as described above are likely to occur (note that a=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 are from Spiegel et al. (2009) and are 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., ISORROPIAlite) did not address these outstanding numerical issues.

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 NH${}_{\mathrm{4}}^{+}$ can occur when solving the minor system NH${}_{\mathrm{3}\left(\mathrm{g}\right)}+{\mathrm{H}}_{\left(\mathrm{aq}\right)}^{+}\leftrightarrow {\mathrm{NH}}_{{\mathrm{4}}_{\left(\mathrm{aq}\right)}}^{+}$ for thermodynamic equilibrium. In this case, HETP and ISORROPIA will solve a quadratic equation to determine the concentration of ammonia gas (NH_{3}). From the concentration of NH_{3}, the ammonium cation is determined as NH${}_{\mathrm{4}}^{+}={\mathrm{NH}}_{{\mathrm{4}}_{i}}^{+}{\mathrm{NH}}_{\mathrm{3}}$, where NH${}_{{\mathrm{4}}_{i}}^{+}$ is the ammonium concentration determined from the major system (see Table S2). If partitioning (after solving the quadratic equation) at this stage gives NH_{3}≈ NH${}_{{\mathrm{4}}_{i}}^{+}$, then subtraction of two nearly identical numbers may lead to a floatingpoint arithmetic error and a final concentration of NH${}_{\mathrm{4}}^{+}<\mathrm{0}$ in the original ISORROPIA equations. In HETP, negative output is strictly prohibited. To accomplish this, 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).

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, potassium, and sodium, in some cases. In HETP we have slightly modified the initial dry salt partitioning of CALCL9 (see Table S2) to ensure mass conservation holds for all cases. Any free TA that may result in L9 is assumed to be in the gas phase as NH_{3} 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 SO_{4}, 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.

Mass conservation may not hold in ISORROPIA when the input precursor concentrations are near the lower limit for species concentrations, “tiny” ($\mathrm{1}\times {\mathrm{10}}^{\mathrm{20}}$ mol m^{−3}), used in the solver. The same lower limit used to bound the input precursor concentrations is also 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.

The subroutine “adjust” performs a postconvergence mass balance adjustment for ammonium, sulfate, nitrate, 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, 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 CaSO_{4} 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 CaSO_{4} in the mass balance adjustment of sulfate.

Improvements to the overall code structure and efficiency include the following:
 (a)
use of modern Fortran compared to FORTRAN 77 in ISORROPIA;
 (b)
use of explicit declarations only (all subroutines now start with an “implicit none” statement, and all common blocks have been removed);
 (c)
removing all GOTO statements and instead using modern 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 postconvergence mass balance adjustment (adjust), with the merging of functions and some short subroutines allowing several variables to be calculated once and reused throughout the iterative process, reducing computational time;
 (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)
precalculating 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 vectorizationbygrid point approach (Makar, 1995), which may reduce the call factor overhead on some compilers.
 (a)
4.1 Casebycase 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 solutions – adjacent 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 the two solvers. For activity coefficients, ${\mathit{\u03f5}}_{\mathrm{act}}=\mathrm{1}\times {\mathrm{10}}^{\mathrm{6}}$ and maxit_{act}=4, where ϵ_{act} is the relative error limit between successive iterations of activity coefficient calculations, and maxit_{act} is the maximum number of allowed iterations. For bisection or ITP, $\mathit{\epsilon}=\mathrm{1}\times {\mathrm{10}}^{\mathrm{9}}$, maxit_{bsec}=100, and ndiv=5, where ε is defined in Sect. 1, maxit_{bsec} is the maximum number of allowed iterations, and “ndiv” 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. The ISORROPIA code used in this comparison is the base version (ISORROPIA v2.2) used in the CMAQ airquality model (USEPA, 2022). ISORROPIA throughout this paper has been compiled using the “r8” flag that converts all real variables 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 “doubleprecision” variables, but some singleprecision 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 nontrivial 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 randomaccess memory. The compiler used was an Intel compiler (IFORT) version 2021.5.0.2021109.
Figure 2 displays the output from ISORROPIA and HETP for two example chemical subspaces: panels (a) and (b) display CALCO7, and panels (c) and (d) show CALCM8. These chemical subspaces involve the presence of at least one of Ca^{2+}, K^{+}, and Mg^{2+}, and so they were not included in the original HETV package, which was designed for the SO_{4}–NO_{3}–NH_{4}–H_{2}O system. Furthermore, these two subspaces are frequently called in practical CTM 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 humidity (RH) was set to 35 % and the air temperature (T) to 306 K, conditions typical of a hot summer's 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 $\mathrm{2.1}\times {\mathrm{10}}^{\mathrm{5}}$ and $\mathrm{2.4}\times {\mathrm{10}}^{\mathrm{5}}$ mol m^{−3}, where visual differences begin to appear, particularly for H^{+}, HSO${}_{\mathrm{4}}^{}$, and NH_{3}. 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, the ISORROPIA solution shows the effects of numerical instability in the bisection rootfinding procedure. The ISORROPIA algorithm used in CALCM8 is designed so that the variable being bisected is proportional to Cl^{−} (see Table S2). 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. 2c and d, where the differences between ISORROPIA and HETP are related only to the choice of rootfinding 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 behavior that is demonstrated in the ISORROPIA simulation shown in Fig. 2d. It should be noted that these differences are due to the choice of rootfinding 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 logarithmic difference between the “calculated” equilibrium constant (K_{calc}) and the “true” equilibrium constant (K_{true}); that is, $\mathit{\xi}=\mathrm{log}\left({K}_{\mathrm{calc}}\right)\mathrm{log}\left({K}_{\mathrm{true}}\right)$. K_{calc} 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 K_{true} 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 ${\mathit{\xi}}^{\prime}=\left\mathit{\xi}\right$. A logarithmic difference is used herein (instead of a percent difference, for example) since the difference between K_{calc} and K_{true} can span several orders of magnitude. In this way, a difference on the order of 1 implies that K_{calc} and K_{true} differ by an order of magnitude, while a difference on the order of $\mathrm{1}\times {\mathrm{10}}^{\mathrm{2}}$ implies K_{calc} and K_{true} differ starting at the second or third digit, when written in scientific notation. The error analysis has been completed using the casebycase implementation of HETP (see Sect. 4.3). Ideally, ${\mathit{\xi}}^{\prime}=\mathrm{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 K_{calc} and K_{true}. The accuracy of K_{true} calculated from Eq. (1) (used in both solvers) is limited to three significant digits due to the variable $\mathrm{\Delta}{H}_{\mathrm{f}}^{\mathrm{0}}/\left(R{T}_{\mathrm{0}}\right)$. Therefore when ${\mathit{\xi}}^{\prime}<\mathrm{1}\times {\mathrm{10}}^{\mathrm{3}}$ in either solver, we can conclude that K_{calc} after convergence is identical to K_{true} within its known accuracy. However, in practical applications (i.e., within a CTM), the value of K_{true} 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 $\ll \mathrm{1}\times {\mathrm{10}}^{\mathrm{3}}$. Hence, we seek a solver that obtains ξ^{′} as close to zero as possible. Table 2 gives the median, the maximum, and the 25th and 75th 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 convergence criteria. For example, for K_{HCl} HETP has a median ${\mathit{\xi}}^{\prime}\approx \mathrm{1.77}\times {\mathrm{10}}^{\mathrm{8}}$, while ISORROPIA has a median ${\mathit{\xi}}^{\prime}\approx \mathrm{0.39}$, with similar results for ${K}_{{\mathrm{HNO}}_{\mathrm{3}}}$. The superior performance of HETP for this set of initial conditions can also be confirmed visually by comparing Fig. 2c to 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 demonstrate that the ξ^{′} values are linked to the poor convergence performance of ISORROPIA 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 25th and 75th percentiles are considered (i.e., for K_{HCl} the 75th percentile of ξ^{′} is $\mathrm{4.88}\times {\mathrm{10}}^{\mathrm{7}}$ and $\mathrm{4.40}\times {\mathrm{10}}^{\mathrm{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 =65 % and T=263 K. The main difference here is that CALCO7 performs slightly worse in HETP than ISORROPIA, as determined from the median and 75th percentile of ξ^{′}. 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 $\mathrm{1.0}\times {\mathrm{10}}^{\mathrm{4}}$, implying that the difference between K_{calc} and K_{true} occurs in the fourth digit when written in scientific notation. As a result, the differences between HETP and ISORROPIA 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 T=263 K), HETP has an unstable behavior in the output speciation for TS between $\mathrm{1.6}\times {\mathrm{10}}^{\mathrm{7}}$ and $\mathrm{2.3}\times {\mathrm{10}}^{\mathrm{7}}$ mol m^{−3}, while ISORROPIA has an unstable behavior for all TS $>\mathrm{0.7}\times {\mathrm{10}}^{\mathrm{7}}$ mol m^{−3} (see Fig. S1). This poor performance in CALCM8 for these meteorological conditions is demonstrated in the statistics of ξ^{′} shown in Table 2.
Figure 3 displays a comparison of HETP and ISORROPIA, where now TS and TA are varied simultaneously while all other input precursor species are held constant. Figure 3 displays output generated from n=2 000 000 unique test cases. These test cases are divided into two tests, 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). 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 $<\mathrm{1}\times {\mathrm{10}}^{\mathrm{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 $\mathrm{1}\times {\mathrm{10}}^{\mathrm{10}}$ mol m^{−3}, thereby creating mass. This adjustment has not been applied here.
The colors in Fig. 3a, b, d, and e represent the amount of gaseous NH_{3} after partitioning between the gas and aerosol phase. The test input spans across all of Branch 4 (O7, M8, P13, L9, and K4), using the same convergence criteria as Figs. 1 and 2. The colors shown in Fig. 3c and f give the absolute percent difference between Fig. 3a and b and Fig. 3d and e, respectively, calculated relative to HETP as HETPISO/HETP × 100 %. The color contour intervals in Fig. 3 are on a logarithmic scale. In each figure panel, dashed black lines separate the different chemical subspaces, with the particular subspace label superimposed. In Fig. 3a, b, d, and 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. 3a, b), 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 fewer call to calculate activity coefficients in HETP for some test cases, specifically those test cases that have no convergence of activity coefficients after completing the maximum number of allowed iterations. The largest absolute differences of 100 % to 600 % are in L9 and are predominantly due to (ii), where for some input conditions ISORROPIA creates dry salt mass for TA, TS, and TK. Specifically in ISORROPIA, 6.02 %, 0.05 %, and 5.97 % of the test input conditions shown in Fig. 3a and 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., ${\mathrm{species}}_{\mathrm{out}}{\mathrm{species}}_{\mathrm{in}}>\mathrm{9.999}\times {\mathrm{10}}^{\mathrm{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. 3a–c, there is a large amount of “noise” in K4 for $\mathrm{TS}>\mathrm{1}\times {\mathrm{10}}^{\mathrm{5}}$ mol m^{−3} and $\mathrm{TA}<\mathrm{12}\times {\mathrm{10}}^{\mathrm{6}}$ mol m^{−3} in ISORROPIA that is not present in HETP. This 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 %.
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 T=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., HSO${}_{\mathrm{4}}^{}$). 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. Without the modifications applied in Fig. 4c and d, the output from HETP and ISORROPIA is 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^{+}–HSO${}_{\mathrm{4}}^{}$–SO${}_{\mathrm{4}}^{\mathrm{2}}$, requiring a quadratic root with a large variation in coefficient magnitudes to be derived – and therefore an error in H^{+} will propagate through to the minor systems that are solved thereafter (see Table S2). It should be noted that the y 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 HSO${}_{\mathrm{4}}^{}$ 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, updated methodologies are used to solve the quadratic equation that avoid numerical inaccuracies due to catastrophic cancellation. 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 NO${}_{\mathrm{3}}^{}$. 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 floatingpoint 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., $\mathrm{1}\times {\mathrm{10}}^{\mathrm{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 K_{HCl} and ${K}_{{\mathrm{HNO}}_{\mathrm{3}}}$ in the chemical subspaces I6, J3, L9 and K4 for some sets of initial conditions. Statistics 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 K_{HCl} from 13.0 to $\mathrm{1.46}\times {\mathrm{10}}^{\mathrm{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^{−}, NO${}_{\mathrm{3}}^{}$, H^{+}, and NH_{3} near $\mathrm{1.4}\times {\mathrm{10}}^{\mathrm{12}}$ mol m^{−3} – and hence these species now show a smooth transition over the entire range of TS. HETP has a limiting precision of $\mathrm{1}\times {\mathrm{10}}^{\mathrm{28}}$ mol m^{−3}, which is the likely cause of the HSO${}_{\mathrm{4}}^{}$ concentration becoming zero in Fig. 4c and d when TS is between about $\mathrm{2.15}\times {\mathrm{10}}^{\mathrm{12}}$ and $\mathrm{2.4}\times {\mathrm{10}}^{\mathrm{12}}$ mol m^{−3}.
4.2 Comparison using input from the GEMMACH airquality 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 GEMMACH airquality model (Makar et al., 2018). In this section, 20 000 unique sets of input data (“test cases”) from GEMMACH are investigated for each chemical subspace, with 10 000 test cases obtained from summer days and 10 000 test cases obtained from winter days. These test cases were selected from input conditions generated from a 10 km resolution simulation with a domain covering North America. The test cases were chosen randomly so that the selected set of test cases spans across a broad range of temperatures and relative humidity, typical of actual tropospheric conditions. Table 4 gives the relative frequency of calls to each chemical subspace as a percentage of the total calls in GEMMACH determined from 4 d (2 in the winter and 2 in the summer). It should be noted that subspaces A2, B4, and C2 all require that TN be formally zero. A low number limit in the GEMMACH model prevents true zero conditions from occurring; hence the given subroutines are not called in this practical application test. The majority of calls are to the subspaces O7, M8, and L9 which comprise more than 75 % of the total calls on these 4 d. Therefore, most situations encountered in GEMMACH over North America have a nonzero amount of base cation species present (K^{+}, Mg^{2+} and Ca^{2+}).
Figure 5 displays a scatter plot of Cl^{−} and HCl (left panels) and NO${}_{\mathrm{3}}^{}$ and HNO_{3} (right panels) output from ISORROPIA (y axis) and HETP (x axis). Figure 5a and b display CALCM8:summer (hereafter M8:S), and Fig. 5c and d show CALCG5:winter (hereafter G5:W). The dashed black lines give a 1 : 1 relationship, denoting where HETP and ISORROPIA agree exactly. There is relatively good agreement between the two solvers for M8:S, despite the differences noted for this subspace in Sect. 4.1. However, for G5:W a large amount of scatter exists, demonstrating disagreement between the two solvers for some test cases. This disagreement is likely related to the choice of rootfinding method and/or other numerical updates that have been made to the HETP code, as described in Sect. 3. The differences between the two solvers noted for Cl^{−} and HCl in Fig. 5c are only for very low concentrations, which likely would not be impactful in practical airquality applications.
As in Sect. 4.1, statistics of ξ^{′} are calculated from the output of each solver to judge the accuracy of the equilibrium solution. This is especially important since the test cases in this section cannot be plotted in a regular fashion (as in Sect. 4.1), to graphically reveal obvious numerical instabilities. Figure 6 displays a boxandwhisker plot of ξ^{′} for the chemical subspaces G5, H6, O7, M8, and P13. These subspaces all require bisection or ITP and must have chloride present, with K_{HCl} providing the “final convergence check” (except for H6). The statistics shown in Fig. 6 include the data shown in Fig. 5 for subspaces M8:S and G5:W. Fig. 6a and b show ξ^{′} for K_{HCl} and ${K}_{{\mathrm{HNO}}_{\mathrm{3}}}$ respectively. Each panel shows ξ^{′} for both seasons with summer having a “:S” label and winter having a “:W” label. In the box plot, the 25th percentile, the median, and the 75th 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 give 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 ${K}_{{\mathrm{HNO}}_{\mathrm{3}}}$ is smaller in HETP than ISORROPIA for all subspaces shown in Fig. 6b. For K_{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 ${\mathit{\xi}}^{\prime}>\mathrm{0.5}$ for all equilibrium constants (but with marginally worse performance in HETP than in ISORROPIA). For example, in H6:S for ${K}_{{\mathrm{HNO}}_{\mathrm{3}}}$ the 75th percentile in HETP is 31.8, and in ISORROPIA it is 13.4. H6 is unique relative to the other subspaces requiring a rootfinding 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 K_{HCl} and ${K}_{{\mathrm{HNO}}_{\mathrm{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 HNO_{3} and NO${}_{\mathrm{3}}^{}$ between the two solvers in G5:W (Fig. 5d), it is clear from the statistics of ξ^{′} for ${K}_{{\mathrm{HNO}}_{\mathrm{3}}}$ and K_{HCl} shown in Fig. 6 that both solvers are producing output that spans a broad range of accuracy. The 75th percentiles of ${K}_{{\mathrm{HNO}}_{\mathrm{3}}}$ and K_{HCl} are 2 orders of magnitude lower in HETP than ISORROPIA. For K_{HCl} the 75th percentile of ξ^{′} is $\mathrm{6.93}\times {\mathrm{10}}^{\mathrm{2}}$ and 4.35 in HETP and ISORROPIA respectively. However, the maximum ξ^{′} values 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 75th percentile of ξ^{′} (for ${K}_{{\mathrm{HNO}}_{\mathrm{3}}}$ and K_{HCl}) by an order of magnitude in both solvers. For HETP the median ξ^{′} for ${K}_{{\mathrm{HNO}}_{\mathrm{3}}}$ reduces to $\mathrm{4.89}\times {\mathrm{10}}^{\mathrm{8}}$ (from $\mathrm{4.60}\times {\mathrm{10}}^{\mathrm{7}}$), and for ISORROPIA the median ξ^{′} reduces to $\mathrm{2.72}\times {\mathrm{10}}^{\mathrm{6}}$ (from $\mathrm{5.59}\times {\mathrm{10}}^{\mathrm{5}}$). The modification to account for oscillatory behavior has the effect of reducing ξ^{′} for the flagged test cases in HETP compared to ISORROPIA – for the 14.02 % of test cases affected, the median ξ^{′} for ${K}_{{\mathrm{HNO}}_{\mathrm{3}}}$ is 0.28 for HETP, but for ISORROPIA it is 2.65. Furthermore, 75.3 % of the flagged test cases are times when Cl^{−} is predicted to be $<\mathrm{1}\times {\mathrm{10}}^{\mathrm{16}}$ mol m^{−3} (Cl^{−} is the bisected variable in G5), and all flagged test cases have TCl $<\mathrm{1}\times {\mathrm{10}}^{\mathrm{10}}$ mol m^{−3}. For test cases where the output from each solver agrees well and falls along the 1 : 1 line in Fig. 5c and d, ξ^{′} for ${K}_{{\mathrm{HNO}}_{\mathrm{3}}}$ and K_{HCl} is minimized in each solver. The statistics of ξ^{′} for other subspaces not discussed here are summarized in Table S4 (summer) and Table S5 (winter) of the Supplement.
4.3 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 is given in Table 5. The timing tests have an estimated uncertainty of ±1 %. For HETP, two sets of timing tests are reported. Test 1, labeled “T_{HETV}”, refers to timing using a global convergence criterion for all tests within a given chemical subspace, representing a “vectorized” test where all n 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, labeled “T_{HETP}”, refers to a casebycase test where the solver is called individually for each test case (i.e., the solver is called n times). In the latter test, the time associated with subroutine calls is offset by the number of iterations becoming test 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. 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 (T_{ISO}) requires a casebycase implementation and cannot solve n 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) and removing 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 rootfinding 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 unique input data derived from GEMMACH 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 and likewise in the summer for subspaces A2, B4, and C2. 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 GEMMACH 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, B4, and C2 respectively, except with TN = 0.
The CPU timing results demonstrate that all subspaces (except H6, winter and summer) execute faster in HETP's vectorized T_{HETV} implementation than ISORROPIA. In some cases the speedup is significant: for CALCO7 the speedup is about a factor of 1.74 to 1.83 when using T_{HETV}. An even greater speedup can be achieved by using the casebycase T_{HETP} implementation for some subspaces, specifically those that require bisection (A2, D3, G5, H6, O7, M8, and P13). Unlike T_{HETV}, all chemical subspaces execute faster in T_{HETP} than ISORROPIA. For the sets of test cases investigated in this work, the bestcase performance is found in P13:S, where T_{HETV} executes in about 0.38 s, but T_{HETP} executes in about 0.18 s (the latter being ∼4.3 times faster than ISORROPIA). The speedup afforded by HETP for this subcase is largely the result of HETP's rootfinding 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 rootfinding methodology of each solver to achieve convergence of the major systems are given in Table 6, for the same input data used to generate the timing tests shown in Table 5. For P13:S which has the bestcase performance, ITP in HETP requires on average 8.0 iterations for convergence, while bisection in ISORROPIA requires on average 42.5 iterations. Thus, HETP's rootfinding method requires about 19 % of the iterations required by ISORROPIA for this set of input conditions, while executing in about 23 % of the time using the casebycase mode. The overall performance for the tests in GEMMACH (bottom row of Table 5) shows the average performance of HETP operating in the casebycase mode results in a speedup relative to ISORROPIA of a factor of 2.24 times for the summer tests and 2.12 times 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.
The difference between T_{HETP} and T_{HETV} becomes even more apparent and in favor of T_{HETP} if a significant amount of test cases do not require bisection. While T_{HETV} includes a return statement to reorder the problem, removing those test cases that have converged or have no solution prior to entering ITP, the rootbracketing stage in T_{HETV} will nonetheless need to be repeated a second time for all test cases that do require ITP. Note that the rootbracketing stage identifies an interval where the objective function has a sign change. Assuming a continuous function, this sign change signifies that a root exists within the interval. Furthermore, in T_{HETV} some test cases may iterate in the rootbracketing 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 criterion do not), thereby introducing excess computations into T_{HETV} that do not exist in T_{HETP}. This is especially true as the variable ndiv, 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 casebycase T_{HETP} 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.3 times faster than ISORROPIA, with an average performance increase in a practical application between 2.12 times and 2.24 times (using the casebycase mode). The subspace H6 which executes slower in T_{HETV} than T_{HETP} and is also less accurate than ISORROPIA for most input test cases accounts for <1 % of the test cases on the days sampled (see Table 4).
In this work we have presented HETP, an updated solver to perform thermodynamic equilibrium calculations of the H^{+}–SO${}_{\mathrm{4}}^{\mathrm{2}}$–NH${}_{\mathrm{4}}^{+}$–NO${}_{\mathrm{3}}^{}$–Cl^{−}–Na^{+}–Ca^{2+}–K^{+}–Mg^{2+}–H_{2}O chemical system, based on the algorithms of ISORROPIA in the forward metastable state. HETP has been updated in several ways to improve both the computational speed and accuracy compared to ISORROPIA. For most input conditions HETP produces equivalent results to ISORROPIA, but for some input conditions the output from the solvers can diverge. Analysis of the output from each solver suggests that HETP's use of ITP, instead of bisection, improves the accuracy of its equilibrium solution for some input conditions by obtaining a more accurate initial estimate of the root prior to the commencement of the ITP search. At the same time, ITP can reduce the number of iterations required for convergence. The differences may be formally linked to reduced accuracy of the ISORROPIA solver's output due to several numerical issues as described in the sections above. In addition to providing more accurate output for most test cases, HETP, when implemented to solve n test cases simultaneously, may execute 1.3 to 2.1 times faster than ISORROPIA (except for CALCH6), based on input from the CTM GEMMACH. Alternatively, when HETP is implemented as a casebycase solver (the solver is called n times), then HETP is 1.3 to 4.3 times faster than ISORROPIA for individual chemical subspaces and 2.1 to 2.2 times faster than ISORROPIA on average, with the speedup being most significant in subspaces that require the application of a rootfinding method for convergence.
The data used in the analysis presented herein, and the HETP code, are available online at https://doi.org/10.5281/zenodo.8164704 (Miller, 2024).
The supplement related to this article is available online at: https://doi.org/10.5194/gmd1721972024supplement.
SJM developed the HETP code, performed the analysis presented herein, and wrote the manuscript, with significant contributions from PAM. PAM also assisted with the development of the HETP code. CJL assisted with inorganic heterogeneous chemistry module interfacing with GEMMACH, GEMMACH tests, and incorporation of the new module into the regional transport model.
The contact author has declared that none of the authors has any competing interests.
Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.
This work was funded under the Oil Sands Monitoring (OSM) Program. It is independent of any position of the OSM Program. We would like to thank Wanmin Gong for reviewing the manuscript prior to submission for peer review and Verica SavicJovcic for providing advice during the peerreview process.
This research has been supported by the Environment and Climate Change Canada (grant no. Oil Sands Monitoring (OSM) Program, subproject “Integrated Atmospheric Deposition”, APD62324).
This paper was edited by Andrea Stenke and reviewed by three anonymous referees.
Amundson, N. R., Caboussat, A., He, J. W., Martynenko, A. V., Savarin, V. B., Seinfeld, J. H., and Yoo, K. Y.: A new inorganic atmospheric aerosol phase equilibrium model (UHAERO), Atmos. Chem. Phys., 6, 975–992, https://doi.org/10.5194/acp69752006, 2006.
Anlauf, K., Li, S.M., Leaitch, R., Brook, J., Hayden, K., Toom–Sauntry, D., and Wiebe, A.: Ionic composition and size characteristics of particles in the lower fraser valley: Pacific 2001 field study, Atmos. Environ., 40, 2662–2675, https://doi.org/10.1016/j.atmosenv.2005.12.027, 2006.
Ansari, A. S. and Pandis, S. N.: Prediction of multicomponent inorganic atmospheric aerosol behavior, Atmos. Environ., 33, 745–757, https://doi.org/10.1016/s13522310(98)002210, 1999a.
Ansari, A. S. and Pandis, S. N.: An analysis of four models predicting the partitioning of semivolatile inorganic aerosol components, Aerosol Sci. Technol., 31, 129–153, https://doi.org/10.1080/027868299304200, 1999b.
Atkinson, R. W., Mills, I. C., Walton, H. A., and Anderson, H. R.: Fine particle components and health – a systematic review and metaanalysis of epidemiological time series studies of Daily Mortality and hospital admissions, J. Expo. Sci. Env. Epid., 25, 208–214, https://doi.org/10.1038/jes.2014.63, 2014.
Bromley, L. A.: Thermodynamic properties of strong electrolytesin aqueous solutions, AIChE J., 19, 313–320, 1973.
Burden, R. L. and Faires, J. D.: Numerical analysis (9th ed.), Cengage Learing, Boston, MA, USA, 861 pp., ISBN 9780538733519, 2011.
Capps, S. L., Henze, D. K., Hakami, A., Russell, A. G., and Nenes, A.: ANISORROPIA: the adjoint of the aerosol thermodynamic model ISORROPIA, Atmos. Chem. Phys., 12, 527–543, https://doi.org/10.5194/acp125272012, 2012.
Clegg, S. L. and Pitzer, K. S.: Thermodynamics of multicomponent, miscible, ionic solutions: Generalized equations for symmetrical electrolytes, J. Phys. Chem., 96, 3513–3520, https://doi.org/10.1021/j100187a061, 1992.
Community Modeling and Analysis System (CMAS, 2016): https://www.airqualitymodeling.org/index.php/CMAQv5.1_Isorropia, last access: 17 July 2023.
Denbigh, K.: The principles of chemical equilibrium, 4th edn., Cambridge University Press, Cambridge, 520 pp., ISBN 9780521281508, 1981.
Fountoukis, C. and Nenes, A.: ISORROPIA II: a computationally efficient thermodynamic equilibrium model for K^{+}–Ca^{2+}–Mg^{2+}–NH${}_{\mathrm{4}}^{+}$–Na^{+}–SO${}_{\mathrm{4}}^{\mathrm{2}\mathrm{}}$–NO${}_{\mathrm{3}}^{\mathrm{}}$–Cl^{−}–H_{2}O aerosols, Atmos. Chem. Phys., 7, 4639–4659, https://doi.org/10.5194/acp746392007, 2007.
Fountoukis, C., Nenes, A., Sullivan, A., Weber, R., Van Reken, T., Fischer, M., Matías, E., Moya, M., Farmer, D., and Cohen, R. C.: Thermodynamic characterization of Mexico City aerosol during MILAGRO 2006, Atmos. Chem. Phys., 9, 2141–2156, https://doi.org/10.5194/acp921412009, 2009.
GEOSChem 14.0.0, Zenodo [code], https://doi.org/10.5281/zenodo.7254288, 2022.
Harrison, R. M. and Pio, C. A.: Major ion composition and chemical associations of Inorganic Atmospheric Aerosols, Environ. Sci. Technol., 17, 169–174, https://doi.org/10.1021/es00109a009, 1983.
Heintzenberg, J.: Fine particles in the global troposphere a review, Tellus B, 41B, 149–160, https://doi.org/10.1111/j.16000889.1989.tb00132.x, 1989.
Hennigan, C. J., Izumi, J., Sullivan, A. P., Weber, R. J., and Nenes, A.: A critical evaluation of proxy methods used to estimate the acidity of atmospheric particles, Atmos. Chem. Phys., 15, 2775–2790, https://doi.org/10.5194/acp1527752015, 2015.
Irwin, J. G. and Williams, M. L.: Acid rain: Chemistry and transport, Environ. Pollut., 50, 29–59, https://doi.org/10.1016/02697491(88)901844, 1988.
Jacobson, M. Z.: Studying the effects of calcium and magnesium on size–distributed nitrate and ammonium with EQUISOLV II, Atmos. Environ., 33, 3635–3649, https://doi.org/10.1016/s13522310(99)001053, 1999.
Jacobson, M. Z.: Global direct radiative forcing due to multicomponent anthropogenic and natural aerosols. J. Geophys. Res.Atmos., 106, 1551–1568, https://doi.org/10.1029/2000jd900514, 2001.
Kakavas, S., Pandis, S. N., and Nenes, A.: ISORROPIA–Lite: A comprehensive atmospheric aerosol thermodynamics module for Earth System Models, Tellus B, 74, 1, https://doi.org/10.16993/tellusb.33, 2022.
Kim, Y. P. and Seinfeld, J. H.: Atmospheric gas–aerosol equilibrium: III. Thermodynamics of crustal elements Ca^{2+}, K^{+}, and Mg^{2+}, Aerosol Sci. Technol., 22, 93–110, https://doi.org/10.1080/02786829408959730, 1995.
Kim, Y. P., Seinfeld, J. H., and Saxena, P.: Atmospheric gas–aerosol equilibrium I. Thermodynamic model, Aerosol Sci. Technol., 19, 157–181, https://doi.org/10.1080/02786829308959628, 1993a.
Kim, Y. P., Seinfeld, J. H., and Saxena, P.: Atmospheric gas–aerosol equilibrium II. Analysis of common approximations and activity coefficient calculation methods, Aerosol Sci. Technol., 19, 182–198, https://doi.org/10.1080/02786829308959629, 1993b.
Kusik, C. L. and Meissner, H. P.: Electrolyte activity coefficients in inorganic processing, AIChE Symp. Series, 173, 14–20, 1978.
Lovett, G. M., Tear, T. H., Evers, D. C., Findlay, S. E. G., Cosby, B. J., Dunscomb, J. K., Driscoll, C. T., and Weathers, K. C: Effects of air pollution on ecosystems and biological diversity in the Eastern United States, Annals of the New York Academy of Sciences, 1162, 99–135, https://doi.org/10.1111/j.17496632.2009.04153.x, 2009.
Makar, P. A.: Fast use chemical numerics methods: the use of “vectorization by gridpoint”, in: Air Pollution III, Vol. 1, edited by: Moussiopoulos, H. N. and Brebbia, C. A., Computational Mechanics Publications, Southampton, 434 pp., 1995.
Makar, P. A., Wiebe, H. A., Staebler, R. M., Li, S. M., and Anlauf, K: Measurement and modeling of particle nitrate formation, J. Geophys. Res.Atmos., 103, 13095–13110, https://doi.org/10.1029/98jd00978, 1998.
Makar, P. A., Bouchet, V. S., and Nenes, A.: Inorganic Chemistry calculations using HETV – a vectorized solver for the SO${}_{\mathrm{4}}^{\mathrm{2}}$–NO${}_{\mathrm{3}}^{}$–NH${}_{\mathrm{4}}^{+}$ system based on the ISORROPIA algorithms, Atmos. Environ., 37, 2279–2294, https://doi.org/10.1016/s13522310(03)000748, 2003.
Makar, P. A., Akingunola, A., Aherne, J., Cole, A. S., Aklilu, Y.A., Zhang, J., Wong, I., Hayden, K., Li, S.M., Kirk, J., Scott, K., Moran, M. D., Robichaud, A., Cathcart, H., Baratzedah, P., Pabla, B., Cheung, P., Zheng, Q., and Jeffries, D. S.: Estimates of exceedances of critical loads for acidifying deposition in Alberta and Saskatchewan, Atmos. Chem. Phys., 18, 9897–9927, https://doi.org/10.5194/acp1898972018, 2018.
Martin, S. T., Hung, H.M., Park, R. J., Jacob, D. J., Spurr, R. J. D., Chance, K. V., and Chin, M.: Effects of the physical state of tropospheric ammoniumsulfatenitrate particles on global aerosol direct radiative forcing, Atmos. Chem. Phys., 4, 183–214, https://doi.org/10.5194/acp41832004, 2004.
Meissner, H. P. and Peppas, N. A.: Activity coefficients – aqueous Solutions of polybasic acids and their salts, AIChE Journal, 19, 806–809, 1973.
Meng, Z., Seinfeld, J. H., Saxena, P., and Kim, Y. P.: Atmospheric gas–aerosol equilibrium: IV. Thermodynamics of carbonates, Aerosol Sci. Technol., 23, 131–154, https://doi.org/10.1080/02786829508965300, 1995.
Metzger, S., Mihalopoulos, N., and Lelieveld, J.: Importance of mineral cations and organics in gasaerosol partitioning of reactive nitrogen compounds: case study based on MINOS results, Atmos. Chem. Phys., 6, 2549–2567, https://doi.org/10.5194/acp625492006, 2006.
Miller, S.: HETP: An updated inorganic heterogeneous chemistry solver for metastable state based on ISORROPIA II, Zenodo [code], https://doi.org/10.5281/zenodo.8164704, 2024.
Nenes, A., Pandis, S. N., and Pilinis, C.: ISORROPIA: A New Thermodynamic Equilibrium Model for Multiphase Multicomponent Inorganic Aerosols, Aquat. Geochem., 4, 123–152, https://doi.org/10.1023/a:1009604003981, 1998.
Oliveira, I. F. and Takahashi, R. H.: An enhancement of the bisection method average performance preserving Minmax optimality, ACM Transactions on Mathematical Software, 47, 1–24, https://doi.org/10.1145/3423597, 2021.
Press, W. H., Teukolsky, S. A., Vetterling, W. T., Flannery B. P.,: Numerical Recipes The Art of Scientific Computing, 3rd edn., Cambridge University Press, Cambridge, UK, 1235 pp., ISBN 9780511335556, 2007.
Pye, H. O., Liao, H., Wu, S., Mickley, L. J., Jacob, D. J., Henze, D. K., and Seinfeld, J. H.: Effect of changes in climate and emissions on future sulfatenitrateammonium aerosol levels in the United States, J. Geophys. Res.Atmos., 114, D01205, https://doi.org/10.1029/2008jd010701, 2009.
Pye, H. O. T., Nenes, A., Alexander, B., Ault, A. P., Barth, M. C., Clegg, S. L., Collett Jr., J. L., Fahey, K. M., Hennigan, C. J., Herrmann, H., Kanakidou, M., Kelly, J. T., Ku, I.T., McNeill, V. F., Riemer, N., Schaefer, T., Shi, G., Tilgner, A., Walker, J. T., Wang, T., Weber, R., Xing, J., Zaveri, R. A., and Zuend, A.: The acidity of atmospheric particles and clouds, Atmos. Chem. Phys., 20, 4809–4888, https://doi.org/10.5194/acp2048092020, 2020.
Quan, J., Liu, Q., Li, X., Gao, Y., Jia, X., Sheng, J., and Liu, Y.: Effect of heterogeneous aqueous reactions on the secondary formation of inorganic aerosols during haze events, Atmos. Environ., 122, 306–312, https://doi.org/10.1016/j.atmosenv.2015.09.068, 2015.
Robinson R. A. and Stokes R. H.: Electrolyte solutions (revised 2nd ed.), Butterworths, London, 608 pp., ISBN 0486422259, 1965.
Rood, M. J., Shaw, M. A., Larson, T. V., and Covert, D. S.: Ubiquitous nature of ambient metastable aerosol, Nature, 337, 537–539, https://doi.org/10.1038/337537a0, 1989.
SaizLopez, A., Plane, J. M., Baker, A. R., Carpenter, L. J., von Glasow, R., Gómez Martín, J. C., McFiggans, G., and Saunders, R. W.: Atmospheric chemistry of iodine, Chem. Rev., 112, 1773–1804, https://doi.org/10.1021/cr200029u, 2011.
Sander, R., Keene, W. C., Pszenny, A. A. P., Arimoto, R., Ayers, G. P., Baboukas, E., Cainey, J. M., Crutzen, P. J., Duce, R. A., Hönninger, G., Huebert, B. J., Maenhaut, W., Mihalopoulos, N., Turekian, V. C., and Van Dingenen, R.: Inorganic bromine in the marine boundary layer: a critical review, Atmos. Chem. Phys., 3, 1301–1336, https://doi.org/10.5194/acp313012003, 2003.
Savoie, D. L. and Prospero, J. M.: Particle size distribution of nitrate and sulfate in the marine atmosphere, Geophys. Res. Lett., 9, 1207–1210, https://doi.org/10.1029/gl009i010p01207, 1982.
Schmale, J., Zieger, P., and Ekman, A. M.: Aerosols in current and future arctic climate, Nat. Clim. Change, 11, 95–105, https://doi.org/10.1038/s41558020009695, 2021.
Seinfeld, J. H. and Pandis, S. N.: Atmospheric Chemistry and physics: From air pollution to climate change, Wiley & Sons, 1152 pp., ISBN 9781118947401, 2016.
Shaw, M. A. and Rood, M. J.: Measurement of the crystallization humidities of ambient aerosol particles, Atmos. Environ. A, 24, 1837–1841, https://doi.org/10.1016/09601686(90)90516p, 1990.
Song, S., Gao, M., Xu, W., Shao, J., Shi, G., Wang, S., Wang, Y., Sun, Y., and McElroy, M. B.: Fineparticle pH for Beijing winter haze as inferred from different thermodynamic equilibrium models, Atmos. Chem. Phys., 18, 7423–7438, https://doi.org/10.5194/acp1874232018, 2018.
Spiegel, M. R., Lipschutz, S., and Liu, J.: Schaum's outlines – Mathematical Handbook of Formulas and Tables, 3rd Edn., The McGrawHill Companies, USA, 289 pp., ISBN 0071548564, 2009.
Tang, I. N., Fung, K. H., Imre, D. G., and Munkelwitz, H. R.: Phase transformation and metastability of hygroscopic microparticles, Aerosol Sci. Technol., 23, 443–453, https://doi.org/10.1080/02786829508965327, 1995.
United States Environmental Protection Agency (USEPA): CMAQ (Version 5.4), Zenodo [code], https://doi.org/10.5281/zenodo.7218076, 2022.
Wang, G., Wang, H., Yu, Y., Gao, S., Feng, J., Gao, S., and Wang, L.: Chemical characterization of water–soluble components of PM_{10} and PM_{2.5} atmospheric aerosols in five locations of Nanjing, China, Atmos. Environ., 37, 2893–2902, https://doi.org/10.1016/S13522310(03)002711, 2003.
Wang, K., Zhang, Y., Nenes, A., and Fountoukis, C.: Implementation of dust emission and chemistry into the Community Multiscale Air Quality modeling system and initial application to an Asian dust storm episode, Atmos. Chem. Phys., 12, 10209–10237, https://doi.org/10.5194/acp12102092012, 2012.
Wexler, A. S. and Clegg, S. L.: Atmospheric aerosol models for systems including the ions H^{+}, NH${}_{\mathrm{4}}^{+}$, Na^{+}, SO${}_{\mathrm{4}}^{\mathrm{2}}$, NO${}_{\mathrm{3}}^{}$, Cl^{−}, Br^{−}, and H_{2}O, J. Geophys. Res., 107, ACH14, https://doi.org/10.1029/2001jd000451, 2002.
Zhang, Y., Seigneur, C., Seinfeld, J. H., Jacobson, M., Clegg, S. L., and Binkowski, F. S.: A comparative review of inorganic aerosol thermodynamic equilibrium modules: Similarities, differences, and their likely causes, Atmos. Environ., 34, 117–137, https://doi.org/10.1016/s13522310(99)002368, 2000.