Stefan Miller, Paul Makar and Colin Lee describe a new and improved solver method termed HETP for solving the thermodynamic inorganic aerosol gas–liquid partitioning problem based on the ISORROPIA II model and its approach. This comprehensive study provides an overview of inorganic aerosol thermodynamics and the typically considered set of species, a description of how the ISORROPIA algorithm and its division of the thermodynamic equilibrium problem into subspaces is designed, and where, how, and why HETP differs in some of the details. It implies that the authors have acquired a detailed understanding of the ISORROPIA algorithm and potential numerical limitations, which they address to a large extent with their HETP model.

Overall, the paper is well written, the tables and figures are clear, useful, and well described, and the evaluation of HETP (and HETV) against ISORROPIA II is done in a meaningful way in terms of numerical accuracy, reliability, and computing time. This reviewer appreciates the distinct sets of comparisons discussed in Sections 4.1 and 4.2, with the latter drawing realistic examples of species inputs from an air quality model (GEM-MACH) run for the North America domain. I find this work to be a valuable effort both in terms of this article and associated evaluations and documentation, as well as in terms of the accompanying Fortran source code, which should enable relatively easy adoption into various air quality models.

My comments are minor. They should help the authors to further clarify a few details of this work prior to final publication.

General comments:

1. Abstract and Section 3, a wording detail: The HETP code is described as being written in “FORTRAN 90”. This is a bit of an oddity, not because it is coded in Fortran, which is perfectly adequate for this computer program, but because of the spelling and implied unnecessarily narrow meaning of the language. If anything, the correct spelling is “Fortran 90” (not capitalized, unlike FORTRAN 77 and prior). However, why use/refer to the code as being in Fortran 90 and why not 95, 2003, 2008 or a more recent version of the standard of the language, such as Fortran 2018 or the upcoming Fortran 2023 standard? By that I do not mean that the authors should make use of all the latest language features, but that by writing “Fortran 90” they may simply, yet incorrectly, mean what is often referred to as “modern Fortran”. That is, a Fortran version that uses free-format source code, a modular procedural (or object-oriented) programming paradigm or a combination of both, portable variable kinds, dynamic allocatable arrays, etc (see, e.g. https://fortran-lang.org). Note that the “.f90” file extension (or .ftn90 in the linked HETP package) does not specifically refer to nor restrict code to Fortran 90; rather, it is a file extension that is (to date) used for any modern free-form Fortran code (since Fortran 90). Also note that most (nearly all) of Fortran 90 is a subset of the most recent standard (since Fortran is backwards-compatible, except for a few very rare, depreciated features).
   
   Therefore, I strongly suggest revising the sentences in this paper when Fortran 90 is mentioned, since it points to an outdated standard of the language.

Specific comments:

1. Line 29: “The inorganic portion ….” Mentioning water as a key inorganic aerosol component would be adequate.
2. Line 32: Aside from bromide, some marine regions also show notable iodine species amounts (such as I- and IO3- ions in aerosols), e.g. Saiz-Lopez et al. (2011, Chem. Reviews).
3. Line 33: “in the coarse mode”; clarify for non-experts: coarse mode of what?
4. Line 50: Another, more recent overview of inorganic aerosol equilibrium models, including E-AIM, MOSAIC, AIOMFAC-GLE, ISORROPIA II and EQUISOLV II, is provided in the review paper by Pye et al. (2020), https://acp.copernicus.org/articles/20/4809/2020/ . Perhaps worth mentioning here.
5. Line 51: Regarding rigorous solvers: another rigorous and sophisticated solver approach is the UHAERO model/method by Amundson et al. (2006), atmos-chem-phys.net/6/975/2006/, suggested to be cited here too.
6. Line 76: “the metastable state (those subsystems in which liquid water is present)”. The text in parenthesis describes a side effect of aqueous inorganic aerosols in metastable state, not the state itself. It is a potentially misleading description of what metastable state means and should be improved. The metastability refers to the fact that certain cation-anion combinations (and potentially hydrates) are present in a supersaturated aqueous solution relative to a more stable state involving presence of some crystalline salts aside from a remaining aqueous solution that is saturated with respect to those salts. In metastable state mode, the thermodynamic solution approach ignores potential formation of crystalline states (perhaps with some exceptions like CaSO4).
7. Line 89: “assume a metastable state as the most likely”; yes often, but note that an aqueous inorganic phase could also be absolutely stable, not just metastable, e.g. at high RH > 80%. Except that there may be minerals present or forming in aerosols that remain crystalline (e.g. gypsum) due to extremely low solubility in water. This was also assumed in this work for CaSO4.
8. Line 96: Perhaps of interest: in ANISORROPIA and related changes to ISORROPIA Capps et al. (2012), doi:10.5194/acp-12-527-2012, have introduced a combination of bisection and Newton's method (primarily needed for the adjoint model, not for better speed).
9. Line 117: “by grid point solutions…”; It should be explained here what is meant by grid point in this context. Is this about a grid point in a 3D air quality model?
10. Line 140: “HETP does not consider the reverse problem…” Please provide an explanation why not? Is that case irrelevant in the applications for which HETP is designed – or for another reason?
11. Line 149: “expressed as molar equivalent Na, …” For those non-volatile cations, I guess this should be stated as Na^+, K^+, Mg^2+ etc. – or is the charge omitted for a reason?
12. Line 156: “As a result, …” This seems not to be a “result” but rather a choice or consideration, right?
13. Line 168: “mean activity coefficients”; Are these molality-based mean activity coefficients or molarity-based or mole-fraction-based? It should be mentioned somewhere. In the following line, “non-ideal solutions” should be non-ideal aqueous solutions because water is used as reference solvent for the ions.
14. Lines 185–188: In this context, do the free amounts of present components mean that charge balance is violated. If charge balance is violated, is additional H^+ or OH^- considered to balance the charges and allow for all "free" amounts to be accounted for in aqueous solution?
15. Line 191: “… to prevent loss of mass of species such as Na, Mg, K, and Ca; the free mass must be added back to the captured mass partitioned by the solver prior to returning to the program accessing the inorganic heterogeneous chemistry solver.” When such cases occur, does the presence of free mass not indicate a problem in the inputs? In other words, if the inputs are not charge balanced, and presumably ions like H+ and OH- not directly measured or traced, would ignoring such free mass really be an important issue or rather reflect minor problems with the accuracy and consistency of input species amounts? Please discuss briefly.
16. Line 198: related to previous comment; “these implementations may be inadvertently losing aerosol mass due to this issue…”; but that loss is presumably only happening at initialization, since afterwards all species should be part of the “captured” portion, or not? Please comment.
17. Table 1: ER1 row on right, please check the notation for the activity coefficients. It is unclear what gamma_H2SO4 means compared to gamma_H-HSO4.
    
    Furthermore, presumably the [H+] notation refers to molalities of H+, not molarity, please clarify in the caption? Also, actual molalities should be normalized by unit molality in such equilibrium relationships, which means that any equilibrium constant is a dimensionless quantity, but scale-dependent; see e.g. Pye et al. (2020). The often-stated units of the equilibrium constants arise from a less rigorous (but common) use of dropping any unit molalities and unit pressure values. I think this should be mentioned in a footnote to this table.
18. Lines 218 and 223: (kJ mol^-1); Shouldn't these stated units be (J mol^-1) to remain consistent with the notation and values used in Table 1? While such values are typically tabulated in kJ/mol, in the application here they are properly converted to J/mol and that is the SI unit that I'd recommend being stated.
19. Line 230–231: “correlation” should be “relation”. And, clarify “where the water activity (𝑎𝑤) is equal to the fractional relative humidity (0 to 1 scale).” This sentence could be improved. The ZSR relation is the assumption that the water uptake by individual components of a mixture equilibrated at a constant RH is assumed to be additive. It does not state that a_w is equal to RH, which is a statement of vapor-liquid equilibrium (modified Raoult's law) for bulk solutions.
20. Figure 1, caption text: “subcases H6, G5, M8, …”. Please point out that subcase H6 refers to CALCH6 in the figure; this wasn't obvious on a first read.
21. Line 278: “ITP has the advantage of ‘superlinear convergence’, and hence …”. So do Brent's method or Ridder's method, e.g., outlined in Numerical Recipes (https://search.worldcat.org/search?q=bn:0521880688). Have the authors considered such well-established alternatives?
22. Line 299: “and hence 𝑥 is not an accurate solution to the system of equations.” … at the targeted tolerance level for x (I assume). Could it be in some cases that the accepted tolerance on x is simply too loose to determine an adequate |f(x)| ~= 0.0 value (a problem with numerical computations) even if the function is continuous?
23. Lines 308–313: Regarding solving quadratic equations with the standard formula, this is a known problem; another option is to use the method outlined in Section 5.6 of the Numerical Recipes book.
24. Line 355: “Removing function and subroutine calls, except for process calls …”. Unclear; do you mean that a lot of the subroutine contents have been merged into a single subroutine? That would seem to be a disadvantage in code readability and unnecessary given that modern Fortran allows for assumed-shape array association via subroutine and function interfaces that avoid the generation of temporary local copies... Please clarify this a bit more.
25. Line 376: “… complied using the ‘-r8’ flag…” Spelling of compiled; also clarify which compiler is used with this flag, since it may differ amongst various options. Also, if ISORROPIA is not used with double precision floating point numbers in its usual implementation within models like GEOS-Chem, would you expect the move from single to double precision to slow down ISORROPIA substantially?
26. Related to above, other compiler flags, such as for fast math or a certain optimization level may also matter a lot in terms of run time performance; were those consistent between ISORROPIA and HETP compilations? It would be good to state somewhere in the SI the compiler flags used for the model comparison runs.
27. Figure 4, Line 560: “The air temperature and relative humidity are 243 K and 5% respectively,”.
    
    Was this combination of temperature and RH used because the demonstrated numerical issues only surface under such thermodynamic conditions? At least in the lower to mid troposphere, such RH and temperature conditions seem rare and when present, how likely would it be that crystallization of certain salts could safely be ignored at 5 % RH? In other words, while the demonstration is good, one might wonder whether these issues within ISORROPIA matter in applications within air quality models – and whether the solution shown with HETP in panel (d) would represent a realistic case or a very hypothetical strongly supersaturated metastable state.
28. Figure 5: The case shown in panel (d) seems to matter most in terms of numerical inaccuracies, while panels (a)–(c) show excellent to sufficient agreement. Slight disagreements in panel (c) only become notable for very low molar amounts, which probably do not matter from an air quality point of view, would you agree?
29. Figure 6: Looking at the relatively poor performance of both models for the H6 category: are these just indicative of an attempt to solve equations for which double precision arithmetic is insufficient? If one were to set a threshold for accepted convergence, say of 10^-3 in terms of the xi’ metric, how many out of all computations would violate this? Such an implications-oriented quantification (thinking of getting the inorganic aerosol mass concentration and composition roughly right) may be of interest for flagging the fraction of insufficiently solved HETP or ISORROPIA predictions in 3D models.
    
    
    
    

Comments on supplement and code:

 30.  HETP package code: I had a look at the Fortran source code and have some related minor comments:

• Instead of using fixed kinds for reals and integers (like kind=8), it is better for clarity and portability to use a module that sets kind parameters, e.g. often something like a “dp” or “wp” kind for double precision (64-bit floating point numbers); a simple example would be: integer, parameter, public :: wp = kind(0.0D0) and use that parameter to declare all real(wp) variables. This will work with any compiler and is thus better than using kind=8 (which does not exist in all Fortran compilers, some using a different kind numbering system). Related, literal values would then be written e.g. as tstd = 298.15_wp instead of 298.15d0 and 82.0567d-6 would be 82.0567E-6_wp. Similar for integers and logicals (one could simply use default integer and logical (no specific kind stated). See https://fortran-lang.org/learn/best_practices/floating_point/ .
• In module mach_hetp_mod, instead of using the “save” attribute module-wide, it would seem unnecessary if the various water activity parameters, such as awsc and awss were defined as parameter arrays; e.g. real(wp), parameter :: awsc(100) = real( [28.16, 28.16, 28.16, …, 0.1], kind=wp)
• It seems that most subroutines are written as quasi-external, standalone procedures provided in the same file. Ideally, with modern Fortran, one would place them inside of one or several modules, so that their interfaces are implicitly known within the module’s use scope and assumed-shape arrays could be used. Many of the subroutines may also qualify as “pure subroutines”.
• Subroutine mach_hetp_calcd3 includes logical variable statement evaluations like
  
  “if (soln == .false. .and. frst == .true.) then”. This is non-standard Fortran; the correct, standard way of writing the same would be: “if ((.not. soln) .and. frst) then”.

In this work the authors present their HETP inorganic thermodynamic partitioning solver, an evolution of the ISORROPIA algorithms commonly embedded in many current atmospheric chemistry models. Based on the stability, accuracy, and efficiency metrics presented here, this solver is an important development that should see rapid adoption within the atmospheric modeling community. On the whole I find this manuscript to be very well planned and composed, with clear figures and text. I have just a few suggestions on how to strengthen an otherwise excellent paper, and I recommend publication following just a few minor changes.

• Formatting of "GEOS-Chem" should be fixed to remove all-caps from the final "hem".
• Lines 155-156 and other similar references mention mass conservation issues in ISORROPIA. Obviously mass conservation is an important modeling goal for a number of reasons, but it would be helpful to have some context on the relative and absolute scale of these violations. How much mass are we talking about here, and how impactful might it be in common practice?
• Figure 3: This figure could use some help to make it more legible and comprehensible. Label text sizes (especially axis labels) are unreasonably small relative to manuscript text. Also, it should be much easier to find and interpret the conditions that differentiate the top row (panels a-c) from the bottom row (panels d-f). I suggest that the values in lines 479-480 be relocated to a small table for easy reference, perhaps with some form of qualitative description for each (High/Low NaCaKMg?) for easier description.
• Figure 4: Panel titles may lead to confusion as written, as "HETP - base" could be interpreted as subtraction (i.e. a difference plot). I suggest switching to a colon ("HETP: base") for the HETP panels.

Please see the comments in the attached pdf.