Thank you very much for your constructive comments and suggestions, which have been extremely helpful to our research. After received your comments, we have undergone two rounds of model improvements by including nitrite and dissolved oxygen as state variables. Here, we present the responses to each of your comments. The related revisions will be incorporated into the upcoming new version of the manuscript.

Comment 1:

Study area. The study area should be described before the model. Key contextual information is missing, including hydrography, sediment characteristics, primary production, and especially dissolved oxygen (O₂) dynamics. Given the central role of O₂ in the biogeochemical reaction network, available O₂ data should be presented (e.g., in Fig. 3).

Response 1:

For the convenience of introducing the model input conditions, we placed the study area after the model description. Of course, we can position the study area before the model, which does not change the contents of this paper. It is true that the study area section lacks an introduction to hydrological conditions, primary production, and other water environment factors, which will be added in revision. Sediment characteristics were mentioned in other places but will be moved to the study area section. In this paper, since we initially avoid modeling dissolved oxygen, we did not provide the dissolved oxygen conditions. In response to your suggestions, we will add a systematic description in the revised version.

Comment 2:

Inconsistent units. Units are inconsistent throughout the text, tables, figures, and equations. For example, diffusion coefficients in Eqs. 1–2 are given in m² s⁻¹, while sedimentation rates are reported in mm yr⁻¹, without explicit unit conversions. Table 2 mixes mm, cm, and m, and concentrations are variously expressed as µmol L⁻¹, mmol m⁻³, or mass-based units. Units should be standardized throughout.

Response 2:

The use of units in this article is not sufficiently rigorous; we will standardize the units in the revised version.

Comment 3:

Conceptual model and state variables. The conceptual diagram in Fig. 1 is incomplete. The fluff layer is not shown, and NO₂⁻ is absent as a state variable despite appearing in the anammox rate formulation (Table 2). The Supplement suggests that NO₂⁻ is prescribed as a constant, but the rationale and implementation are unclear. If NO₂⁻ is important for anammox, it should be treated as a dynamic state variable. This may partly explain the unusually large modelled contribution of anammox to fixed N loss (Fig. 5d), which exceeds typical values reported for coastal sediments (e.g., Dalsgaard et al., 2005; doi:10.1016/j.resmic.2005.01.011). Although porewater profiles in Fig. 4 are reasonably reproduced, this may be a result of unrealistic internal N cycling if key intermediates are not modelled explicitly.

Response 3:

The conceptual diagram mainly distinguishes three parts: the oxygenated zone, the hypoxic zone, and the area outside the model domain. The fluff layer is only the surface layer of the model and is part of the oxygenated zone. Therefore, the fluff layer is not listed alongside the other parts and was not labeled. In the revised version, we will add a label for the fluff layer.

Since many sediment models do not describe nitrification in separate steps and do not model nitrite, this study treated nitrite as a constant. In fact, we have recently introduced nitrite as a state variable in the model. The new results show significant changes in the relationship between anammox and denitrification. Clearly, in the previous version, fixing nitrite concentration as a constant to drive the anammox led to an overestimation of its rate. In the new results, the denitrification rate is significantly higher than the anammox rate (Figure 1c, compared with Figure 5c in the manuscript). Therefore, we greatly appreciate this suggestion and its contribution to our research. We will replace the results of old model using those from our new model with nitrite in the revised manuscript.

Comment 4:

Treatment of oxygen. Several reaction terms depend on O₂, yet O₂ is not included as a state variable. Instead, an oxygen penetration depth (OPD) appears to be imposed. It is unclear whether O₂ concentrations are assumed constant above the OPD. Given the strong control of O₂ on N transformations, this approach is difficult to justify. Including O₂ as a dynamic, seasonally varying state variable would substantially improve the model. At minimum, a fixed O₂ upper boundary condition would be preferable to a static OPD. Again, without doubt this will have an important influence on the anammox rate. Seasonal variability in O2 concentrations over depth versus time would be a key plot to show in the main manuscript.

Response 4:

In this study, we focus on the seasonal variations of nitrogen cycling in sediments, primarily relying on model results validated against observational data. Since there is no measurement data of metal elements or sulfide concentrations in the sediment, we gave up the simulation of DO and instead specified the dissolved oxygen (DO) profile as an input condition (handled similarly to nitrite) in our model. We are also well aware of the critical role that dissolved oxygen plays in sediment nitrogen cycling. Therefore, after model validation, we quantitatively described the DO profile using two indicators—maximum concentration and oxygen penetration depth (OPD)—and conducted two groups of numerical experiments to explore the impact of changes in the DO profile on the model results. As reported in Section 4.2, when OPD remains unchanged, variations in the maximum DO concentration have no effect on the results; however, when the maximum concentration is constant, changes in OPD below 3 mm significantly affect the model results. Overall, variations in the DO profile have a limited impact on the model results, mainly affecting nitrate and its related processes.

During revision, however, after introducing nitrite as a state variable, numerical experiments revealed that the updated model is highly sensitive to changes in DO profiles. This sensitivity likely results from the decisive role of DO on nitrite, which is an intermediate variable. Consequently, we have recently taken a further step by including DO as a state variable in the simulation. Due to the lack of any data on metal elements and sulfide concentrations, the treatment of processes related to DO follows Soetaert’s scheme (Soetaert et al., 1996). The model with DO is currently under debugging, and related updates will be presented in the revised manuscript.

Soetaert, K., Herman, P. M. J., and Middelburg, J. J.: A model of early diagenetic processes from the shelf to abyssal depths, Geochim. Cosmochim. Acta, 60, 1019-1040, 10.1016/0016-7037(96)00013-0, 1996.

Comment 5:

N₂O is reported as a model output (Line 277), but the corresponding governing equations are not included in Table 1. These must be provided.

Response 5:

Perhaps we did not clearly explain the calculation of nitrous oxide. Its flux is calculated only through the denitrification. The related calculation equation is described in the denitrification formula in Table 1, and the proportional parameters are listed in Table 2.

Comment 6:

Table 1 lists a diffusive boundary layer thickness of 3 m, which is likely a typo (∼3 mm would be more realistic). This should be clarified. If a diffusive boundary layer is included, Fig. 1 and Eq. 9 should be revised to reflect flux continuity at the sediment–water interface.

Response 6:

Thanks for pointing this typo and its unit should be mm. We will correct it in the revised version.

Comment 7:

It is unclear how accumulation and erosion of the fluff layer are treated, particularly in relation to the advection and diffusion terms in Eqs. 1–2. This needs explicit explanation.

Response 7:

In this study, the accumulation and erosion of the fluff layer are reflected through the exchange processes at the sediment-water interface, including the sinking-resuspension of particulate matter and the diffusion exchange of dissolved matter. These processes are described in detail in the boundary conditions section. The exchange between the fluff layer and the interior of the sediment is expressed through the convective-diffusion terms in the governing equations 1-2, consistent with the processes between the other sediment layers. Additionally, the mineralization rate in the fluff layer is set to be several tens of times higher than that in the sediment.

Comment 8:

Table 1. PO4 does not need to be included in the model description since it is not simulated.

Response 8:

Exactly, this study did not focus on phosphorus. We will remove the related expression of phosphorus in Table 1 in the revised version.

Comment 9:

Reaction equations in Table 1 are not balanced with respect to H, O, or charge. Each POC degradation term includes an additional limitation factor (1/lim), even though limitation terms are already specified. The stoichiometric coefficients (x, y) could be used directly in the mass-balance equations (Eqs. 5–7), potentially removing the need for the rCN parameter. The basis for the chosen x and y values should be explained. The assumed oxidation state of organic carbon (apparently zero) should also be stated explicitly.

Response 9:

Due to our oversight in checking the chemical equations in the references, there was a problem with unbalanced equations. We will correct it in the revised manuscript.

To make the expressions clearer, it is better to retain the carbon-to-nitrogen ratio as a parameter. We exactly omitted the introduction of the values for this parameter. In fact, we also measured the total carbon concentration profile in the sediments. Therefore, we calculated the average carbon-to-nitrogen concentration ratio to use as the value for this parameter.

We will also include the assumptions regarding the chemical states of the organic matter in the revised version.

Comment 10:

Anaerobic solutes are represented by a lumped oxygen demand unit (ODU), but ODU is not treated as a state variable. Consequently, O₂ and NO₃⁻ consumption during ODU oxidation is not represented. Explicit inclusion of ODU would improve internal consistency and confidence in the model output.

Response 10:

Since we initially tried to avoid simulating DO, we did not set ODU as a state variable. Recently we started simulating changes in DO concentration in the model and have also simulated ODU. Thanks for your suggestion.

Comment 11:

Why was bioirrigation not included in the model? I would assume that this would be a major solute transport term in coastal sediments, even if hypoxic (Dale et al., 2013; doi:10.5194/bg-10-629-2013). This needs careful justification.

Response 11:

Through previous researches, we found that the process of bioirrigation can be represented in two ways: as a diffusion term (Berner, 1980) or as a source-sink term (Dale et al., 2011). Following your suggestions, we now include it in the model using the source-sink form. This change will be added to the revised version.

Berner, R.: Early diagenesis: a theoretical approach: Princeton, 1980.

Dale, A. W., Sommer, S., Bohlen, L., Treude, T., Bertics, V. J., Bange, H. W., Pfannkuche, O., Schorp, T., Mattsdotter, M., and Wallmann, K.: Rates and regulation of nitrogen cycling in seasonally hypoxic sediments during winter (Boknis Eck, SW Baltic Sea): Sensitivity to environmental variables, Estuar. Coast. Shelf Sci., 95, 14-28, 10.1016/j.ecss.2011.05.016, 2011.

Comment 12:

The denominator in Eq. 4 should read 2 ln(porosity) instead of 2.02 ln(porosity).

Response 12:

We may not have adequately cited the reference of this expression. In fact, we referenced earlier research as the origin of the expression, but the value of 2.02 ln(porosity) was taken from another source (Radtke et al., 2019). We believe the difference between the two is little. We will add a citation note at the expression accordingly in the revision.

Radtke, H., Lipka, M., Bunke, D., Morys, C., Woelfel, J., Cahill, B., Böttcher, M. E., Forster, S., Leipe, T., Rehder, G., and Neumann, T.: Ecological ReGional Ocean Model with vertically resolved sediments (ERGOM SED 1.0): coupling benthic and pelagic biogeochemistry of the south-western Baltic Sea, Geoscientific Model Development, 12, 275-320, 10.5194/gmd-12-275-2019, 2019.

Comment 13:

According to Fig. 4, NH4 fluxes are out of the sediment, opposite to NO3, yet in Fig. 5b the NO3 and NH4 fluxes have the same sign.

Response 13:

The misunderstanding here may be due to the key information in Fig. 5f being too small. In fact, when the image is enlarged, it can be seen that the nitrate concentration first increases and then decreases from the surface layer to the subsurface layer. Therefore, the nitrate flux is also outward into the bottom water.

Comment 14:

The numerical code used for the model should be specified, and model mass-balance performance should be reported. Analytical methods are insufficiently described. Finally, both the model code and the empirical data should be made publicly available in an online repository for scrutiny by the reviewers.

Response 14:

We may not have clearly introduced information related to the model code. The model is written in Fortran. The mass conservation issue can be verified through the nitrogen budget in section 3.3. Additionally, we discuss the sensitivity to various parameters in section S1 of the supplementary materials to describe the characteristic and reliability of the model. 

We will upload the updated model code and experimental data in subsequent revisions.

We sincerely thank you for your constructive comments, which have helped us improve the quality of the manuscript. All comments have been carefully addressed.

Comment 1：

According to my view, the drawbacks encompass the underestimation of bioturbation (in the model) and the lack of verification of some important parameters: PON (PN) accumulation rate and benthic fluxes (possibly benthic chambers).

Response 1:

We are not very certain about the magnitude of the bioturbation intensity. In fact, we can only verify and confirm this intensity through the concentration profile of particulate organic nitrogen (PON). The diffusion of the PON concentration profile includes only the bioturbation diffusion coefficient. By adjusting this parameter, we can easily bring the vertical distribution characteristics of the simulation results closer to the observed data, which is the method we used to determine the diffusion coefficient value. Moreover, this value falls within the range used in some previous studies (Dale et al., 2011; Soetaert et al., 1996). Therefore, we consider the bioturbation intensity to be reasonable.

Regarding the validation of PON accumulation rates and fluxes, we indeed lack data. We do not have the capability to perform in situ measurements using benthic chambers. We hope to validate the model’s reliability by examining the seasonal variations of various nitrogen species concentrations and dissolved inorganic nitrogen (DIN) fluxes. For the nitrogen system in the sediments studied here, if the internal concentration simulations are accurate and the DIN fluxes, as the input-output fluxes, are also simulated correctly, we believe the model can be able to reasonably represent the actual conditions of the system.

Dale, A. W., Sommer, S., Bohlen, L., Treude, T., Bertics, V. J., Bange, H. W., Pfannkuche, O., Schorp, T., Mattsdotter, M., and Wallmann, K.: Rates and regulation of nitrogen cycling in seasonally hypoxic sediments during winter (Boknis Eck, SW Baltic Sea): Sensitivity to environmental variables, Estuar. Coast. Shelf Sci., 95, 14-28, 10.1016/j.ecss.2011.05.016, 2011.

Soetaert, K., Herman, P. M. J., and Middelburg, J. J.: A model of early diagenetic processes from the shelf to abyssal depths, Geochim. Cosmochim. Acta, 60, 1019-1040, 10.1016/0016-7037(96)00013-0, 1996.

Comment 2：

Please, use PON (PN) for water column particulate N and Norg. (Ntot.) for sediment solid phase.

Response 2:

This may be a matter of usage habits. In many sediment studies, particulate organic matter (POM) is commonly used in chemical equations to represent organic substances or detritus (Paraska et al., 2014), such as particulate organic carbon (POC). In this study, since we focus on the nitrogen cycle, we use “PON” to represent particulate organic nitrogen. In fact, some studies also use this abbreviation to denote particulate organic nitrogen in sediments (Lefebvre et al., 2001; Wang et al., 2015). To distinguish particulate organic nitrogen in bottom water, we use a superscript notation to the right of the term, with "bw" indicating particulate organic nitrogen in bottom water (PON^bw in section 4.1). We believe this approach may also achieve clear expression.

Paraska D W, Hipsey M R, Salmon S U. Sediment diagenesis models: Review of approaches, challenges and opportunities[J]. Environmental modelling & software, 61: 297-325, 2014.

Lefebvre S, Bacher C, Meuret A, et al. Modeling approach of nitrogen and phosphorus exchanges at the sediment–water interface of an intensive fishpond system[J]. Aquaculture, 195(3-4): 279-297, 2001.

Wang C, Jiang R, Mao X, et al. Estimating sediment and particulate organic nitrogen and particulate organic phosphorous yields from a volcanic watershed characterized by forest and agriculture using SWAT model[C]. Annales de Limnologie-International Journal of Limnology. EDP Sciences, 51(1): 23-35, 2015.