We would like to thank the anonymous referee for his thorough review of our manuscript and for his constructive comments. We have carefully considered all remarks and revised the manuscript accordingly. Below, we provide a point-by-point response to each comment:

1. Section 2.1 has been reduced accordingly. The description of the governing equations has been condensed and only modified or newly implemented equations are now shown in the manuscript, while the remaining formulations are referenced to previous publications. Additionally, small modifications were implemented in Section 2.2 (coupling modes) to further reduce its length. All changes are highlighted in the revised manuscript version.

2. All eleven issues identified in the manuscript have been fully addressed in the updated solver. At the same time, we acknowledge that the specific test case selected for this study does not explicitly exercise every one of these improvements. As you correctly identified, points 3, 10, and 11 are not directly demonstrated because they are not relevant to the configuration of the chosen benchmark:

   • Point 3 (periodic boundary conditions): This functionality is implemented in the new solver, but it applies only to closed-domain simulations. The test case used in the manuscript is an open system, so periodic boundaries are not invoked.
   • Point 10 (interface boundary conditions for one‑way or single‑domain simulations): The presented test case focuses on fully two‑way coupled surface–subsurface flow. As a result, the one‑way or single‑domain configurations covered by this improvement are not part of the scenario.
   • Point 11 (runtime reinitialization of data fields): Although this capability is included in the new solver, the selected test case does not require any field reinitialization during the simulation.

   To avoid any ambiguity, we have added a table in Section 2.3 which explicitly clarifies that these three points are implemented in the solver but are not triggered by the particular test case used in this study. This ensures that readers understand both the completeness of the improvements and the scope of the demonstrated results.

3. The sentence referred to was misleading at the original location and has, therefore, been relocated to line 143 in the revised manuscript, where it now more explicitly describes the pseudo-transient coupling procedure and the assumptions behind it. Specifically, the sentence clarifies that, owing to the slower response of groundwater compared to surface water, the subsurface solution is assumed to adjust quasi-instantaneously to the surface solution within each time step, thereby avoiding the need for multiple coupling iterations to ensure consistency between surface and subsurface exchanges at the same time level.

4. The volume fraction (alpha) was missing in the diffusion term of Eq. 15 and has now been added in the revised manuscript to ensure consistency with the code and the correct computation of solute diffusion restricted to the water phase.

5. Analogously to the solver’s code, the data associated with the test case was provided as part of the assets for review (https://editor.copernicus.org/EGUsphere/ms-records/egusphere-2025-3521). The link to the data repository is also included in the reference list of the manuscript as “Pardo-Álvarez, Á.: Code and data associated with Pardo-Álvarez et al., 2025, Zenodo [data set], https://doi.org/10.5281/zenodo.16283559, 2025c.”

6. Quantitative validation of variables other than temperature (e.g., flow velocities, water levels or solute concentrations) could not be performed for the selected test case due to the lack of corresponding in-situ observations. Note that the acquisition of velocity fields and the exact, spatially distributed water table constitutes a major observational challenge. Consequently, additional validation plots for these variables cannot be provided. Nevertheless, the simulated flow, temperature and concentration fields are qualitatively consistent with results reported in previous studies conducted at the same site. In addition, the simulated thermal plume compares well with that from drone-based thermal imagery acquired on-site. Together, these qualitative comparisons support the correct functioning of the newly implemented solver.

We would like to thank the anonymous referee for the time and effort he devoted to reviewing our manuscript. We have carefully considered all his valuable comments and suggestions and revised the manuscript accordingly. Below, we address each comment point by point:

1. A quantitative validation metric has been added in Section 4.2 by calculating the root mean square error (RMSE) of the simulated temperature field. This metric complements the qualitative comparison with drone-based thermal imagery and further supports the correct functioning of the newly implemented heat transfer module, based on the spatial heat pattern of the plume generated by the discharged groundwater. Quantitative comparisons with additional variables (e.g., flow velocities, water levels or solute concentrations) could not be performed due to the lack of corresponding field measurements for the selected test case, as already indicated in the manuscript.

2. A quantitative criterion to assess when steady-state conditions can approximate transient system response to periodic stresses has been added in Section 5, following the guidelines proposed by Haitjema (2006). This approach calculates a dimensionless aquifer response time, τ, given by:

   τ = (1/4*P)*(S*L^2/K*b)

   where L is the characteristic system length defined as the distance between major surface water features (m), P is the period of the forcing function (=365 days for seasonal fluctuations), K is the hydraulic conductivity (m/s), b is the saturated thickness (m) and S the is aquifer storativity (specific yield, S_y, for an unconfined aquifer).

   According to the resulting value of , a steady-state approximation will be appropriate for τ > 1.0, while a transient model will be required for 0.1 < τ < 1.0 and a bounding or successive steady-state solution may be used for τ < 0.1.

   An expanded discussion of this timescale classification is provided in the manuscript to clarify the conditions under which the steady-state assumption is valid.

3. A new table has been added in Section 2.3 including the suggested information, namely: (1) the limitations of the original solver, (2) which of these limitations are addressed in the new code, (3) the type of update (new or extended functionality) and (4) whether or not each new feature is applied in the test case of this study. The added table clearly distinguishes between extensions and newly introduced features, thereby clarifying the novel contributions of the upgraded solver.

4. A comprehensive and systematic computational cost analysis, including direct performance comparisons with previous versions of the code and with other solvers of the same class, is beyond the scope of the present study. The primary focus of this paper is on the description and verification of the newly introduced modelling capabilities. Nevertheless, to address practical considerations related to computational cost, Section 5 provides guidance on the computational resources required to run efficient simulations with the new solver, thereby assisting users in selecting appropriate hardware for their applications.

5. Figure 1 was designed as a flow chart to illustrate the workflow of the newly developed solver. We have reviewed the figure carefully and believe that its current form already fulfils this purpose. Could you, please, clarify this point?