This manuscript aims to study the effects of beach morphological changes, hydrodynamic boundary conditions, and hydrogeologic properties on flow regimes, salt distribution, and the potential for mixing controlled chemical reactions in high-energy beach aquifers. The authors achieve this objective by building 2-D transient density-driven groundwater flow and transport numerical models that systematically explore boundary conditions and hydrogeological parameters in a changing beach morphology setting. The authors then conclude that changing beach morphology causes the migration of infiltration and exfiltration locations along the beach transect, leading to moving flow and salt patterns in the subsurface and enhancing mixing-controlled reactions.

I have read the manuscript with great interest. My overall opinion is that the manuscript is well-written but needs some moderate revisions and clarifications. Below, I have listed comments, hoping they may help improve the manuscript’s quality.

Specific Comments

1. Some clarification is needed regarding the 2-D numerical model approach.
   1. The beach boundary, where freshwater inflow occurs, is close to where flow, salinity, and mixing variations occur. Have the authors considered moving this boundary further inland to remove possible boundary effects?
   2. The description of the freshwater inflow boundary is missing from the methods section. Also, why is the value of 0.5 m³/d/m decided for most cases? Is this value related to the groundwater discharge in the area where the model is based?
   3. The authors mention that meteoric groundwater recharge was applied at the upper beach slope above the mean high water line (MHWL) (lines 123 to 124). What is the value of this recharge (350 or 400 mm/a)? These values are missing in Table 1.
   4. What is the initial concentration distribution in the model? Did the authors run a warm-up period before calculating the variations?
   5. In my mind, when the beach morphology changes, part of the model geometry also changes. How do the authors estimate the variations in salinity and mixing patterns with varying geometry?

1. Throughout the manuscript, the authors highlight the importance of changes in beach morphology in the dynamics of upper salinity plumes (USPs). However, the modeling choice to account for these changes needs further clarification.
   1. It is unclear how the surface interpolation was made, and given that the authors express the importance of this feature in the modeling choice, more information should be provided on this matter. Thus, I recommend adding more information on the methods of how the interpolation was made and refining the plot in Figure 1 so that the reader can picture how beach morphology changes over time.
   2. This reviewer understands the study's limitations. Between lines 351 and 355, the authors reference the limitations of linearly interpolating the beach morphology. However, how do the authors reconcile the idea of cycling beach morphology when changes in the hydrologic forcings and hydraulic properties occur? Have the authors considered changing the order of the interpolation in the beach profiles to see if different patterns emerge? Also, the authors only considered a stable case with and without storm floods. Have they considered exploring the variations in hydrogeologic parameters and boundary conditions with stable cases to compare them with the changes in beach morphology?

1. Regarding mixing-controlled reactions, the authors used the initial salinity patterns for the Rf and Rs concentration patterns. This modeling choice could create misleading mixing patterns because these chemical constituents were assumed conservative, and the model needs to stabilize first. So, did the authors consider warming up the model before setting the initial concentration patterns?

Technical Corrections

Besides the comments described above, I have a few technical recommendations for the manuscript.

1. In line 17 (in the abstract), I recommend being more specific in the sentence: “The objective was to investigate their individual effects…,” whose individual effects?
2. I recommend rewording the sentence between lines 55 and 58. It is difficult to understand what the other studies found and their limitations.
3. A comma is missing between “…model Greskowiak…” in line 62.
4. Review the sentence between lines 67 and 69. It sounds redundant.
5. There is a double comma in line 86.
6. The term variability of total dissolved solutes (TDS) (in line 88) is ambiguous. Are there variations in the concentration values of TDS, the spatial distribution or which variability?
7. The statement generic modelling approach (in line 93) is also ambiguous. What do the authors mean by “generic modelling”?
8. Line 95 is missing a space in “…conditions(Hayes, 1979)….”
9. Consider rewording the sentence and citation between lines 100 and 101.
10. In Figure 1,
    1. There is a small box in (a) that is not described in the caption. Is this the area where changes in beach morphology occur? If so, I recommend stating it in the caption.
    2. Is the hydraulic conductivity presented in (c) through (f) the horizontal (K_h) or the vertical (K_v) component?
    3. There is a space missing between “…scans(Grünenbaum et al., 2020a)…” in the caption.
11. The concentration of saltwater is missing in line 134.
12. This might conflict with the journal's requirements, but I suggest changing the equations to a math format as it is easier to read.
13. In Table 1,
    1. What do the bold values mean? Are they to highlight the changing parameters? If so, I recommend stating this in the caption and keeping the bold values consistent throughout the rows.
    2. What is the `-“-` symbol in α_L in cases 4 and 5?
    3. What is the acronym SF? Is it Storm Flood? If so, please state it somewhere in the text.
    4. I recommend capitalizing the “based case” in the description of case 3 to be consistent with the rest of the table.
    5. I recommend correcting the spacing in the description of case 15 to be consistent.
14. I consider the clustering analysis a neat exercise for examining similar models. However, switching between Figure 6 and Figures 3 through 5 can be cumbersome. I recommend adding boxes around the subplots in Figures 3 to 5 that represent the colors of the clusters; that way, it is easier for the reader to visualize which models are clustered together.

The main contribution of this work is that the authors showed that geomorphological changes drive changes in salinization and reaction potential in beach aquifers using general 2-D numerical models. They showed that there is greater mixing of fresh and saltwater due to geomorphic change and this leads to greater reaction potential. They also demonstrated that hydraulic conductivity, dispersivity, and storm floods control the extent of the mixing zone and salinization under changing geomorphic conditions.

            The authors provide an important contribution by incorporating geomorphic change into coastal groundwater models. The paper is technically sound and well written, but it could benefit from a few clarifications in the methods section. Therefore, I recommend minor revisions in accordance with the comments below.

Minor Comments

Line 96 What is the datum that MHWL and MLWL are referenced to?

Line 122 What is the specified flux that was used?

Line 123 What is the flux for the meteoric recharge boundary?

Lines 124 Why was a general head boundary used?

Lines 128-130 How was the linear interpolation performed from the lidar scans (i.e. was each point in the domain sampled from the lidar scan or was it sampled at a different resolution), and how was the morphological change represented in the model (i.e. how was the interpolation done between each morphological realization, if erosion occurred was salt mass removed from the model and if deposition occurred was it assumed to be saline or fresh)?

Line 133 Why was the simulated salinity assigned to the water discharging across the ocean boundary?

Line 134 Was the simulation run to a steady-state salinity distribution before running the transient model?

153-154 It seems like the initial distribution of the reactants is unrealistic given that there was no spin up to a steady-state salinity distribution. This should be mentioned as a limitation.

Line 180 I think it would help the reader to describe the difference between the base case and the stable case before discussing the results.

Line 193-195 I am confused by Figure 5 case 3. How is it normalized? If it is normalized to the base case then there should be no variation in the base case figure.  A description of the normalization should be included in the methods.

Technical Corrections

Line 60 led should be lead

Lines 67-69 This sentence does not make sense. Move citation to end of sentence.

Table 1 Case 2 Description What does SF stand for? I think storm flood but this is not defined anywhere.

Line 261 reached should be reach and mowing should be moving

Line 289 Where is the black cross?

Figure 6 Some of the points are hard to see. It would be helpful if the scenarios were labeled so that each number could be read.

Line 308 lead should be led

Line 322 Remove “related to”

Line 331 Rephrase “impact significantly on the shape”

Line 339 Rephrase “made that simulations ended”