We thank the reviewer for the kind comments and would like to take the opportunity to provide a short development history of Alquimia that both justifies the timing of this submission as well as the role of the perspective provided in the submission. In a sense, this contribution is a product of this long process. 

It was initially conceived and developed as a solution to provide proven geochemical capabilities to the code Amanzi during its development phase, specifically using PFLOTRAN geochemical subsroutine routines. The advantages of the generic stand-alone design, rather than being tied exclusively to PFLOTRAN, became apparent when CrunchFlow was added to the interface. As such, it became a representative example of code interoperability, at a time where lack of interoperability had been identified as an obstacle for code development in the Department of Energy (DOE) software ecosystem (from linear solvers to application codes). In developing Amanzi’s test suite, it also became apparent the possibilities that swapping drivers and engines opened up for code intercomparison efforts under way ( and references therein). Over the last few years, in addition to maintaining and refining aspects of the interface, most work has gone into implementing it into additional codes. In one case, it was done to replace a previous coupling to CrunchFlow; in another case, the availability of a well-defined and documented interface facilitated implementation into land models; in a third, the availability of a C-interface and implementation examples in the test drivers facilitated coupling to openFOAM, and lastly, ATS relied on the earlier Amanzi implementation due to their close relationship. While it was indeed a long-time coming, it would have not been possible to envision and present all these examples and applications back in 2013. 

The technical comments provided in the pdf will help us improve the manuscript by clarifying some terms that are not clear in the draft, such as call signature, and correct a number of errors. They are detailed in what follows, where the comment from the reviewer is provided first and then our discussion or way we are addressing it:

• Abstract: not clear what is a “call signature” 

Call signature refers to the function arguments and their types, whether input or output. We will replace it by a more explicit expression, avoiding the use of the word signature:

“The interface enforces the function arguments and their types for setting up,...”

• Introduction: “We present…” is not a good style for a scientific publication, rather “it is presented” should be used 

Yes, the sentence is unnecessarily long. We simplify it to 

“Alquimia v1.0 is a generic interface to geochemical solvers that facilitates development of multiphysics simulators by enabling code coupling, prototyping and benchmarking.”

• L 65: ci is the concentration of species i per unit volume of water? 

Yes:

“ci is the concentration of species i (mass per unit water)” 

• L 75: there is no Nc in Eq (2). Is Nc the number of primary species implicitly? 

Yes, we need to note more clearly that there are Nc equations in Eq. (2). 

• L 90: … does not stipulate that any specific mathematical form that is used. … 

To be removed. 

• L 114: Eqs (5-7) 

To be removed.

• L 117: Eqs (5-7) 

To be removed.

• L 133-135: What does it mean “enforcing a signature for geochemical subroutines using a single-cell model”? Why “signature” in this context? Is this a metadata structure? This needs more explanation… 

We replace the term and provide an explanation that clarifies what we mean: 

“Alquimia has two parts: (1) an engine-independent application programming interface (API) consisting of all relevant functions, data structures, constants, and their respective types (Table 1) and (2) an optional utility library.”

• L 183: Would it be nice to provide a link to GitHub repository (https://github.com/LBLEESA/alquimia-dev) here? 

To be added

• L 215 on: Would it be nice to point at “call signature” or “signature for geochemical subroutines” in code block 1 or 2, to give the reader some idea of what it really is? 

As noted earlier, we will be using a different wording. When referring to these code blocks at line 277, we write:

“While Alquimia’s approach to engine data sometimes introduces more detail in the code, its flexibility allows Alquimia to accommodate the needs of very different engines. Specifically, the subroutine ReactionStepOperatorSplit shares the same arguments for any engine (lines 1 and 2, code blocks 1 and 2).”

• P 19, Figure 2 caption: in the last 3 lines, the text regarding Langmuir and Freundlich isotherms seems to miss something, e.g. “The Langmuir and Freundlich sorption isotherms presented for CrunchFlow as they are not directly available in CrunchFlow, …” 

To be corrected:

“The Langmuir and Freundlich sorption isotherms are not presented for CrunchFlow. While they can be simulated via a surface complexation model and a single sorbing species, no specific keyword in the input deck is available and this was not pursued further here.”

• L 653: what can be non-zero? (R_i_n from eq 12?) 

Yes, the rate can be non-zero.

• L 779: The statement “Alquimia can be used even when there is no transport” sounds like a paradox in the context of Multiphysics and needs more explanation.

The only thing that this means is that there is no requirement when using Alquimia to include any and all processes in a multiphysics simulation or be solving a spatially distributed problem. In the same way, that e.g. one can use PHREEQC, PFLOTRAN or CrunchFlow as batch simulators. It was stated in the wrong way. This will be rephrased to make it clearer:

“For example, Alquimia can be used to perform simple batch chemistry calculations of the kind routinely carried out by geochemical models such as PHREEQc, without consideration of other processes that involve fluxes over a spatially discretized domain. This may be necessary for batch-scale laboratory experiments or, as in the land surface model example, to expand the range of reactive processes considered (Sulman et al., 2022, 2020).”

We thank the reviewer for the kind comments and we are glad to see that the messages that we wanted to convey came through in his/her reading as well. We address the comments below to improve the article, where the comments from the reviewer are presented first and then our discussion and our approach to address it

• It is correctly stated that chemistry is embarassingly parallel and  that alquimia is designed with a "single cell" approach which makes it possible to flexibly deal with hardware resources, however no further comment nor code example is given in the paper about it. I would like the authors to elaborate on this aspect since it is in my  opinion crucial in the scope of RTM. Is the parallelization strategy completely determined by the driver (e.g., reusing domain partitions as the AMR example of code listing n. 4 seems to implicitly suggest) or does alquimia explicitly allow to introduce ad-hoc parallelization schemes for the geochemical engine in the coupling interface? As an exemplary use case, the driver (transport code) would use N domain partitions on N CPUs to solve flow and transport, but chemistry could be solved on M > N CPUs since linear scaling is expected for geochemistry. Could for example a round robin parallelization independent on the one used in the transport code be implemented at the alquimia level, as e.g. implemented by https://doi.org/10.5194/gmd-14-7391-2021 ?

Yes, we need to clarify this point. 

By single-cell model we mean that structures hold data for one cell of the domain and functions apply to these data. Alquimia is a 0-dimensional model that has no notion of the spatial problem. It also implies well-mixed geochemical conditions in this cell, or batch reactor (e.g. beaker). The geochemical problem is solved for this cell, and this solution is independent from all other grid cells. 

The driver determines the parallelization strategy for the solution of the multiphysics problem. Because the geochemical (ODE) equations are uncoupled across the domain, there is considerable freedom in staging their integration and load balancing. All code examples (except the land model) are multithreaded CPU-based implementations and involve MPI parallelization by the driver. Code block 4 shows how Amanzi-S is also explicitly exploiting multi-threading within a Fortran Array Box.  The work to integrate all the cells in the domain can be distributed across the available threads with no race or synchronization concerns. As suggested by the reviewer, load balancing of chemical calculations may be at odds with load balancing strategies that incorporate stencil operations (such as advection or diffusion). While none of the examples presented does it, it is possible to use different load balancing strategies for transport and chemistry when using Alquimia. The driver could choose to redistribute the state data between each stage of the time-split integration differently, including one that seeks better load balancing when sharp geochemical fronts are present in areas of the domain (and thus computations are more expensive) such as round-robin approaches (De Lucia et al., 2021). The situation is more complex in GPU implementations, but because this is out of scope for this work (we do not have the ability to address this specifically with the current engines and drivers), we recommend the discussion in (Balos et al., 2025) about the topic, including on how to implement parallelization of time-split operators with tools such as SUNDIALS (Balos et al., 2025). 

Although it is not the case for PFLOTRAN or CrunchFlow, the geochemical engine could also implement a form of task parallelism for the geochemical solution within a grid cell. It has been suggested that such a form of parallelism could speed up geochemical calculations in systems with a very large number of species and reactions.  Anecdotal reports indicate that this has been effective for the codes Geochemist’s Workbench (GWB) ChemPlugin (Bethke, 2024) and TOUGHREACT (Sonnenthal et al., 2021). 

We will add this discussion in a revised version of the paper, with an emphasis on the flexibility provided by Alquimia, when it comes to the parallelization of the solution of the time-split problems. 

• Formal references to eqs solved by Alquimia and geochemical engine respectively are to be checked: sometimes (5-7), e.g., lines 114,  117, but often (8-5), e.g. lines 140, 142. Please carefully review these references to avoid confusion.

Agreed. We need to ensure correct reference to equations. 

Additional references provided in this reply are given here

Balos, C. J., Day, M., Esclapez, L., Felden, A. M., Gardner, D. J., Hassanaly, M., et al. (2025). SUNDIALS time integrators for exascale applications with many independent systems of ordinary differential equations. The International Journal of High Performance Computing Applications, 39(1), 123–146. https://doi.org/10.1177/10943420241280060

Bethke, C. M. (2024). The Geochemist’s Workbench® Release 17 ChemPluginTM User’s Guide. Champaign, Illinois: Aqueous Solutions, LLC. Retrieved from https://www.gwb.com/pdf/GWB/ChemPluginUsersGuide.pdf

De Lucia, M., Kühn, M., Lindemann, A., Lübke, M., & Schnor, B. (2021). POET (v0.1): speedup of many-core parallel reactive transport simulations with fast DHT lookups. Geoscientific Model Development, 14(12), 7391–7409. https://doi.org/10.5194/gmd-14-7391-2021

Sonnenthal, E., Spycher, N., Xu, T., & Zheng, L. (2021). TOUGHREACT V4.12-OMP and TReactMech V1.0 Geochemical and Reactive-Transport User Guide. Retrieved from https://escholarship.org/uc/item/8945d2c1

Steefel, C. I., Appelo, C. A. J., Arora, B., Jacques, D., Kalbacher, T., Kolditz, O., et al. (2015). Reactive transport codes for subsurface environmental simulation. Computational Geosciences, 19(3), 445–478. https://doi.org/10.1007/s10596-014-9443-x