Reviewer comments are in bold and our reply is in italic

This paper presents a very nicely executed groundwater flow model from a site with a remarkable set of field observations. This is an impressive model of an impressive field site, and I really enjoyed reading the paper. I made some notes in the pdf, but they are nothing more than minor word changes. Instead, my suggestions are ways to improve the focus of the paper, which should help improve its impact. This comes down to three suggestions:

Dear Dr. Wilson,

thank you very much for your taking the time reading through the manuscript and making these valuable suggestions.

1. Improve the focus of the introduction. As written, the justification of the work is pretty much that no one has ever made a realistic model of such a well-characterized field site. In the authors’ defense, it really is a remarkable model, and different people will take different things from it. As a coastal hydrogeologist, I am delighted to see such a detailed description of the model set-up, and I am looking forward to having my graduate students comb through the methods to make sure they understand them. As indicated more clearly in the abstract, there is a lot for biogeochemists in the mixing and variability at the fringes of the USP; for everyone, there is clear visualization of the flow paths that support observed temperatures, salinity and age dates; and the model provides a firm foundation for future modeling efforts that will include reactive transport. It needs to be published. So from some point of view “here is our model” is not a terrible justification for the paper. But still, it would be better (more citable, easier to get through review) if the research questions were clearer in the introduction.

To me, the obvious questions that this paper either answers or can answer are (1) where are the most active and transient zones of mixing surrounding the USP; (2) where and on what time scales should we sample in beach systems (possibly with added justification by posing this as an iterative, coupled investigation, where the initial observations were needed to create the model, and now the model (and its truly cool movie) can guide the location and timing of future field observations – with extra points for posing these locations in a way that is transferable to other sites, e.g. MLLW, MHHW, etc); (3) how much do major modifications of beach morphology change flow and mixing below the beach (the 2022 beach shortening event) and (4) how do the idealized (and typically steady-state) models that permeate the literature compare to a real field site (exactly how “bad” are they)?

 We fully understand your points here, and in response to that we will explicitly make clear that our major aims were to reconstruct the physical processes in the beach aquifer as a consequence of various superimposing hydro(geo)logical forces in a real high-energy beach setting and to identify potential hotspots of intense gradients and variability that may also be relevant for biogeochemical transformation processes. Second, we also would like to pick-up on your suggestion below regarding the fresh and saline submarine groundwater (SGD) fluxes and groundwater ages of this SGD, and how they seem to be dependent on the hydrological forcing and beach morphodynamics.

With respect to model-informed sampling as an objective: Given the rather complicated superposition of morphodynamics, spring-neap and daily tidal changes affecting the SGD fluxes, it will be rather hard, or even impossible to come up with an optimized, model-informed iteratively adapted strategy on sampling the different out-flowing SGD components, as conditions seem to quickly change in an unpredictable manner. In the beginning of the project this was not clear to us either. What can be said however, is that in such systems, dense monitoring at high resolution and over a few years are needed to come up with a reasonable conceptual model and parametrization that allows to visualize a more or less complete picture of the underground and which (i) can be used to identify specific locations for sampling where subsurface (reactive) hots-spots are likely prevailing on somewhat longer time-scales than just days, and (ii) serve as bases for further computing reactive transport. We will add this kind of elaboration to the Conclusion chapter.

Your suggestion on investigating the suitability of idealized models to capture the flow and transport patterns in the subsurface of a high-energy beach is an important point. Nevertheless, this would require more simulation runs, switching on and off the dynamics of individual boundary conditions, followed by a stringent analysis of the simulation results for the affected state-variables to be decided on, e.g., distribution, intensity and dynamics of hot-spot, fresh and saline SGD fluxes, residence time (distribution) etc. This requires careful work and planning, and is as we think, beyond the scope of this paper. It is certainly needed and will be worked on in the next phase of our project.

2. Having just suggested the beach shortening event as a possible focus of the paper, it would be helpful to illustrate the changes that surrounded that event in a more focused way. As it is, the shortening event comes up a lot, but it is hard to see it in the figures, and it is dispersed among different sections. It is nicely visible in the movie, but I wonder if some kind of summary figure could be developed, with a focused discussion of the event all in one place. Part of this may also be that I thought I saw some interesting things presented as part of the calibration (e.g. illustrating the ability of the model to reproduce transients), when some of the findings could actually be interesting results in addition to calibration targets. (Perhaps sub-headings within the results section would help focus readers’ attention on the most salient points?)

Although the beach shortening event was mentioned several times, it was not meant to be a focus of the paper, and we would like to keep it that way in the revised manuscript. One of the major aims was to reconstruct the physical processes in the beach aquifer as a consequence of the superposition of various impacting hydro(geo)logical forces, where the sudden beach erosion event is only one of out of many effects. We will make this aim clearer in the introduction (see our reply to point #1). Therefore, we would like to refrain from adding a summary figure on the beach erosion event.

3. Some of this may be intended for a later paper, but it would be very helpful to pull in the SGD crowd with some fluxes. The model is the perfect way to reconstruct volumetric fluxes through different parts of the model (the USP, the FDT, the total variation in fresh groundwater discharge from the dunes). This would actually be a great way to illustrate the potential importance of the 2022 beach shortening, where line graphs could show the width of the beach and the changes in the various fluxes over time. Perhaps with the addition of some kind of average age of water discharging from different flow systems (the USP, FDT), or the average area of each flow system (admittedly harder to automate than fluxes or ages)? Presenting information that is hard to collect in the field (fluxes), expensive to analyze (age dates), or that happened before extensive monitoring began (beach shortening) is a fantastic justification for creating a model. Again, some of this might be for the next paper, but all of it would be interesting.

We agree, and will pick up on your suggestion with respect to the SGD fluxes. Interestingly, our pre-analyses already reveal that the sudden beach shortening and corresponding abrupt change in the mean beach slope, appear to have little effect on the SGD fluxes. Rather, on the longer time-scales, the spatial variability of the beach slope (in terms of standard deviation) in the intertidal zone seem to control the SGD fluxes. Also, the spring-neap cycle play a role (as also previously shown in Heiss & Michael, 2014, WRR), and even on a daily time-scale, the sub-daily changes in semi-diurnal tide range seem to have a big impact. We will add a figure on SGD fluxes and mean groundwater ages, respectively, and a corresponding discussion in the Results and Discussion Chapter to elaborate on this in detail.

In summary, I think this work is interesting for many reasons, and I hope that an improved focus will help a broader range of readers appreciate its contributions.

Thank you for your kind evaluation.

The reviewer comments are given in bold and our reply in italic font.

The authors investigate the groundwater flow system, the heat and mass transport processes along a two-dimensional model of a high-energy subterranean estuary. During a year and a half of field observations — with sampling varying in time and space — they recorded the hydraulic head, temperature, salinity, and water age. A significant part of the model parameters (49 in total!) was determined by calibration, performing approximately 10,000 simulations and "selecting the best solution" using the particle swarm optimization (PSO) method. I consider the work performed to be very enormous, the techniques used (field measurements, numerical simulation, automated calibration) to be very comprehensive, and the results achieved and presented to be realistic. In light of all this, I consider the manuscript — after some clarification and completion — suitable for publication in the HESS journal.

Dear reviewer,

thank you very much for your careful evaluation of our manuscript and the kind suggestions for improvement. This is very much appreciated.

Major comments/concerns

The description of the numerical model is rather sporadic, which does not facilitate the work of modelers, a deeper understanding of the article, or its potential reproducibility.

We agree and will add information in the manuscript/supporting information with respect to your comments, as specified below.

Therefore, I suggest:

a. the equations used in the simulations should be presented in the article if possible, or at least in the supplementary material, where it should be clarified that water density depends purely on concentration, while viscosity is constant.

We will include the governing equations in the supporting information.

b. Due to the complexity of the processes, the boundary conditions are also very complicated (Section 2.2) and difficult to follow. It would be useful for the reader if they were summarized in a figure, emphasizing their time dependence. This would also be useful because the total length of the model is 1450 m (Page5Line134), but the entire model area is never shown in the manuscript, only between 500–1300 m and 800–1300 m. This would also eliminate the shortcoming that I could not find any information on the salinity, temperature, and water age boundary conditions for the vertical and lower boundaries.

We agree, we will include a figure of the entire model domain, including the definition of all model boundaries and their time-dependence in the supporting information.

c. Fixed model parameters (e.g. R, Ss, D) should be tabled too.

Yes, we will include them in Table A1.

2. The PSO algorithm is indeed suitable for finding the minimum of highly non-linear inversion problems, although this is a major challenge in the case of a model involving almost 50 unknowns, time-dependent, free topohaline convection. Could you show a figure that illustrates the convergence of the method? What object function did PSO use to quantify the individual simulations at the pilot points?

We will include figures in the supporting information that show how the parameters converge to their final values over the course of the optimization iterations. The objective function was the weighted sum of squared residuals, of all observations, i.e., one single value. We will elaborate on the weighting scheme in our reply to the next comment.

Regarding pilot points: From your comment we get the impression that our explanations on ‘pilot points’ were perhaps misleading. Pilot points are not used for evaluating the difference between observed and simulates state-variables, e.g., salinity, and thus they do not appear in the objective function. Instead, they are predefined locations where hydraulic conductivities are estimated and updated during each optimization iteration and which serve then as reference points to fill the grid-cells in between them via interpolation. We will make this clearer in the manuscript where we describe the pilot point approach.

3. The cross plots in Figure 3 show the relation between the measured data and the simulated results for different calibration quantities. The figure is not discussed in the manuscript, which could answer some of the arising questions. Are the parameters of the best-fitting model shown on the y-axes? Were these diagrams used by PSO to qualify each simulation? Are the quantities at different pilot points and at different times shown here? Are all measurement points given equal weight? Could this explain why the root mean square error of the smaller number of water age data is much greater than that of the large number of e.g. head data? I would appreciate an explaining paragraph regarding Figure 3.

Yes, this figure was just to show how the model performs with respect to the different state-variables individually. All the displayed numerical performance measures were calculated by using the unweighted residuals between simulation and observation. We feel that the figure is not as important as it first appeared, as the model’s performance is also very well visible in the follow-up figures in the manuscript. Figure 3 is merely meant as an add-on. Thus, we have decided to shift this figure into the supporting information and add some more explanations in the figure caption.

Weights: In the objective function, the weights are usually set simply as the inverse of the measurements standard deviation. However, in a multi-objective optimization framework, where different types of state-variables and with an uneven number of data points are used for calibration, weighting is not that straight forward any more. The weights need to balance the uneven amount of data points in the individual datasets of different state-variables and their often vastly different numerical values – And with that, their uneven contribution to the objective function. During the optimization process, these contributions likely change, as the goodness of fit for one or the other dataset may become superior. Thus, to compensate for that and keep the right balance, an iterative approach would be ideal to evaluate the objective function after each iteration in the optimization process and to adapt the weights ‘on-the-fly’. However, there is still a debate on the weight problem in the scientific community when it comes to multi-objective optimization. It is difficult and we would refer to in-depth elaborations on the topic by e.g. Dai and Samper (2004) and to the PEST++ documentation (https://github.com/usgs/pestpp).

In the PSO framework of PEST++, however, there is no automatic on-the-fly update of weights. Instead, one has to stick with fixed weights. This leads to the problem that we do not end up with only one single optimization run, but with many trial runs, in which suitable weights have to be manually found via trial-and-error – A daunting task. While 10000 simulations were carried out for the ‘final’ optimization run, more than 10 of those runs had to be carried out prior to that in order to find suitable weights which could balance the different contributions of the individual datasets to the objective function in a way that finally a reasonable fit could be achieved.

We will provide more details on the weights used in the objective function for each evaluated state-variable. This information will go into the supporting information.

Dai, Z., and J. Samper (2004), Inverse problem of multicomponent reactive chemical transport in porous media: Formulation and applications, Water Resour. Res., 40, W07407, doi:10.1029/2004WR003248.

4. The thermal retardation factor was fixed during the simulations, R=2. This value seems very high, although there is some ambiguity in the literature concerning the definition of the R value. What porosity, density and specific heat values give this factor?

Yes, in the literature different values of the thermal retardation factor R were reported. We have not calculated it from the porosity, density and specific heat capacity but simply used the average value of values reported in previous studies (referenced in the manuscript) that employed heat transport simulations in sandy aquifers. These ranged between 1.8 and 2.24.

With specific heat capacities of quartz and water of 710 J/kg/K and 4183 J/kg/K, respectively, and specific densities of 2650 kg/m³ and 1025 kg/m³, respectively, the porosity would need to be 0.32 when assuming and R of 2. The calibration results suggested effective porosities between 0.25 and 0.28. Given the uncertainty of those model-based estimates at field-scale, we think that the assumption of R=2 is acceptable.

Minor comments

Fig. 2 It would be useful to indicate the north-south direction.

We will include this.

P7L181 Dez → Dec

Will be corrected

P15L326 3He → 3H

Will be corrected

P15L332 3He → 3H

Will be corrected

P16L362 permeable material → permeable medium/sediment…

We will change this.

Fig. 7 caption I guess the figure illustrates the standard deviation of simulated salinity, groundwater age and temperature time series along…

Yes. We will add the term ‘time series’.

P17L382 Figure 5 → Figure 5.b

Will be corrected.

P20L457 Reference is missing

We will include appropriate references. The sentence will then read: “Here, mixing-controlled reactions, e.g., the oxidation of dissolved ferrous iron and precipitation ofiron(hydr)oxides as so-called iron-curtain (Charette and Sholkovitz, 2002), or nitrification due to ammonium oxidation (Spiteri et al., 2008) may be promoted.”

Charette, M.A., Sholkovitz, E.R., 2002. Oxidative precipitation of groundwater-derived ferrous iron in the subterranean estuary of a coastal bay. Geophysical Research Letters 29, 85-1-85–4. https://doi.org/10.1029/2001GL014512

Spiteri, C., Slomp, C.P., Tuncay, K., Meile, C., 2008. Modeling biogeochemical processes in subterranean estuaries: Effect of flow dynamics and redox conditions on submarine groundwater discharge of nutrients. Water Resources Research 44. https://doi.org/10.1029/2007WR006071

P20L460 Reference is missing

We will include appropriate reference. The sentence will then read: “In these zones, biogeochemical turnover and secondary reactions may be controlled by the varying degradability of dissolved organic matter as a function of groundwater age (e.g., Waska et al., 2021).”

Waska, H., Simon, H., Ahmerkamp, S., Greskowiak, J., Ahrens, J., Seibert, S.L., Schwalfenberg, K., Zielinski, O., Dittmar, T., 2021. Molecular Traits of Dissolved Organic Matter in the Subterranean Estuary of a High-Energy Beach: Indications of Sources and Sinks. Front. Mar. Sci. 8. https://doi.org/10.3389/fmars.2021.607083

It is somewhat surprising that the terms “biogeochemical” and “reactive” appear 16 and 11 times in the manuscript, resp., even though the simulation does not involve biogeochemical or reactive transport modeling…

Yes, there is a large scientific community that is interested in unravelling the hydrobiogeochemical transformation processes in subterranean estuaries (or beach aquifers) aiming at estimating element fluxes to the ocean via submarine groundwater discharge (SGD). Flow and transport dynamics critically affect biogeochemical turnover and thus reactive transport in the subsurface. As elaborated in our conclusion, the here identified hotspots of high temporal variability with respect to salinity and groundwater age conditions will likely be also hot spots for reactive processes. Thus, we think that our results are also important for hydro- and biogeochemists in this research field.

Nevertheless, maybe somewhat related to this comment, we will take up the suggestion of reviewer #1 and will extract and present fresh and saline SGD fluxes from the model in order to make this paper also potentially interesting for scientists that focus on the SGD water fluxes, rather than solely biogeochemistry.

Overall, the work presented in the manuscript is very substantial, multifaceted, and fits well with the subject matter of the HESS journal, so, after taking into account the above clarifications, I support its publication.

Thank you for your kind evaluation.