Seasonal variations and controlling factors of nitrogen fluxes at the sediment-water interface in a semi-enclosed inland sea

Nitrogen fluxes across the sediment-water interface and nitrogen removal from sediments are essential components of nitrogen cycle in semi-enclosed inland seas. However, the difficulty in observational sampling hinders continuous data availability that is necessary to understand their seasonal variations and underlying mechanisms. To address this issue, we developed a one-dimensional sediment nitrogen-cycle model and used sensitivity experiments to quantify the relative roles of environmental and biogeochemical drivers. Model results indicate that 48 % of particulate organic nitrogen (PON) settling into sediments is returned to the bottom water as dissolved inorganic nitrogen (DIN), while 6 % is removed via N-loss flux (dinitrogen gas and nitrous oxide). The seasonal variations of PON and DIN fluxes are controlled by fundamentally different mechanisms, with PON flux primarily regulated by bottom-water PON concentration and bottom stress, while DIN flux is mainly governed by temperature-dependent biogeochemical transformations and nitrogen availability. These contrasting responses reveal a decoupling between particulate and dissolved nitrogen fluxes, reflecting the buffering capacity of sediments. Denitrification controls nitrogen removal but its amount is not large because the oligotrophic conditions of the study site limits the nitrate availability.

1 Introduction

Sediments play a crucial role in regulating the storage and removal of nitrogen within semi-enclosed inland seas, which significantly influences the balance and stability of the ecosystem (Huettel et al., 2014; Zhang et al., 2020; Kalvelage et al., 2013; Thamdrup, 2012). On one hand, the sediment can bury substantial amounts of nitrogen, potentially on par with river inputs (Liu et al., 2021; Han et al., 2021; Lønborg and Markager, 2021). On the other hand, because of remineralization, the sediment continuously releases dissolved inorganic nitrogen (DIN) into the overlying water (Devol, 2015; Boynton et al., 2017; Liu et al., 2012). In semi-enclosed inland seas, vertical mixing facilitates the contribution of DIN from sediments as a vital nutrient source for primary production because of shallow water depth (Huettel et al., 2014; Yi et al., 2023; Leng et al., 2023; Liu and Yin, 2007; Berelson et al., 2013; Leynaert et al., 2011; Mei Liu et al., 2011). Notably, the sediment serves as the primary anoxic environment in these regions, where key nitrogen removal processes, such as denitrification and anammox, occur. These processes are essential for regulating nitrogen inventories and influencing nutrient structure in these areas. Some research indicates that continental shelves, which constitute merely 7.5 % of the global seafloor, are responsible for 44 % of fixed nitrogen losses (Devol, 2015; Mctigue et al., 2016; Huang et al., 2021; Jäntti and Hietanen, 2012; Khalil and Rasmussen, 2012). Furthermore, nitrous oxide, one of the gaseous forms of nitrogen removal, is a significant greenhouse gas that contributes to global warming (Quick et al., 2019; Wilson et al., 2020; Yang et al., 2022).

To clarify the function of the sediment in ecological and environmental changes, it is necessary to conduct a quantitative analysis of nitrogen fluxes at the sediment-water interface in semi-enclosed inland seas, along with their seasonal variations. Direct observation of particulate organic nitrogen (PON) flux at this interface poses challenges due to the influence of resuspension processes. PON flux measurement techniques are confined to in situ observations, with relatively limited studies. In contrast, DIN flux can be obtained not only through in situ methods but also via laboratory experiments. Consequently, numerous studies focused on DIN flux and its seasonal variations primarily through in situ observations, incubation experiments, or estimations based on Fick's First law (De Vittor et al., 2012; Mu et al., 2017; Chen et al., 2023). For example, studies in the northern Adriatic Sea and the Baltic Sea have demonstrated the same conclusion that DIN flux exhibits seasonal variations, with higher levels in summer and lower levels in winter, likely influenced by water temperature (De Vittor et al., 2012; Niemistö et al., 2018). There are also studies identifying organic matter deposition as an important factor influencing sedimentary nitrogen cycling and DIN fluxes (Albert et al., 2021; Ratmaya et al., 2022). However, the relative contributions of organic matter supply and temperature to seasonal variations in nitrogen fluxes remain poorly quantified. Furthermore, the directions and controlling mechanisms of ammonium (${\mathrm{NH}}_{\mathrm{4}}^{+}$) and nitrate (${\mathrm{NO}}_{\mathrm{3}}^{-}$) fluxes vary considerably among coastal environments, with sediments acting either as sources or sinks of these nitrogen species. In addition to DIN exchange, nitrogen removal processes such as denitrification and anammox have also been reported to exhibit pronounced seasonal variability, generally with higher rates in summer and lower rates in winter (Zhang et al., 2018; Rich et al., 2020; Teng and Lin, 2024). Previous studies suggest that dissolved oxygen (DO), temperature, organic matter availability, bioturbation, and hydrodynamic conditions can all influence these processes (Mu et al., 2017; De Vittor et al., 2015; Dale et al., 2022; Bohlen et al., 2011; Ratmaya et al., 2022; Zhou et al., 2022).

However, sedimentary nitrogen cycling involves a cascade of interconnected transformations linking particulate organic nitrogen deposition, mineralization, nitrification, denitrification, and nitrogen removal. As a result, changes in organic matter input do not necessarily lead to proportional changes in DIN release or nitrogen loss. Different nitrogen pools and fluxes may respond to environmental forcing through distinct pathways and timescales, potentially causing decoupled responses between particulate and dissolved nitrogen fluxes. Consequently, despite extensive observational studies, a mechanistic framework capable of explaining the seasonal variability of sedimentary nitrogen cycling and the contrasting responses of different nitrogen fluxes remains lacking. Quantifying nitrogen fluxes at the sediment–water interface in semi-enclosed inland seas and elucidating their seasonal dynamics and controlling mechanisms are crucial for understanding benthic nutrient recycling, evaluating the contribution of sediments to ecosystem functioning, and improving predictions of biogeochemical responses to environmental change.

Numerical modelling is a valuable tool for understanding the factors that govern nitrogen fluxes at the sediment-water interface. Typically, nitrogen fluxes at this interface derived from empirical formulas serve as essential boundary conditions for modelling the aquatic nitrogen cycle in semi-enclosed inland seas (Lønborg and Markager, 2021). However, this approach largely overlooks the changes occurring within the sediment and the corresponding flux responses to these alterations. An alternative approach is the implementation of a box model, which investigates nitrogen processes within the sediment. Nonetheless, the box model's treatment of all sediment as a single domain, coupled with its failure to account for the vertical distribution of chemical substances – particularly oxidants such as DO – limits its capacity to provide a comprehensive understanding of the nitrogen cycling processes within the sediment (Yang et al., 2022). A one-dimensional vertical sediment model addresses this limitation, albeit at the cost of increased complexity and computational demands (Radtke et al., 2019; Umlauf et al., 2023). It not only reproduces seasonal variations in nitrogen fluxes but also provides a framework for disentangling the relative roles of environmental forcing and internal biogeochemical processes in regulating sedimentary nitrogen cycling. Despite the clear advantages of this modelling approach over empirical formulas and box models, well-validated applications remain relatively scarce.

This study focuses on the unclear seasonal variations in the magnitude and direction of nitrogen fluxes at the sediment–water interface in semi-enclosed inland seas, alongside the insufficiently understood underlying mechanisms with controlling factors driving these changes. We developed a one-dimensional vertical sediment model to simulate nitrogen cycling and accurately quantify nitrogen fluxes at the interface. Supported by monthly continuous observations conducted in the bottom water and sediment of Harima Nada in Japan, a typical semi-enclosed inland sea characterized by restricted water exchange, seasonal stratification, and active sediment-water nutrient cycling, we obtained a robust dataset for model validation. The model successfully reproduced the seasonal variations in the concentrations of three nitrogen species within the sediment and nitrogen fluxes at the sediment-water interface. Additionally, numerical experiments were further conducted to disentangle the relative roles of organic matter supply, temperature, DO, and benthic processes in controlling seasonal nitrogen cycling. Particular attention was given to the contrasting responses of particulate and dissolved nitrogen fluxes and to the mechanisms that regulate the transfer of deposited organic nitrogen into DIN and nitrogen loss pathways. By identifying the dominant controls on different nitrogen pools and fluxes, this study provides new insights into the mechanisms governing sedimentary nitrogen cycling in semi-enclosed coastal systems and improves our ability to predict benthic responses to environmental change and eutrophication management.

2 Materials and Methods

2.1 Study area

Harima Nada is situated within Seto Inland Sea, the largest semi-enclosed sea in Japan (Fig. 1). It has an average water depth of 26 m and is connected to Osaka Bay and Hiuchi Nada on the eastern and western sides, respectively. To the south, it connects to the Pacific Ocean through the Naruto Strait and Kii Channel (Zhu et al., 2019; Chang et al., 2009). Harima Nada is strongly controlled by seasonal variations in wind forcing and density structure (Tong-u-dom et al., 2023). Its hydrographic circulation exhibits clear seasonal patterns: in winter, a wind-driven downwind circulation dominates under the prevailing north-westerly winds, whereas in summer, a density-driven cyclonic circulation develops due to the horizontal change in stratification. The water column structure also shows pronounced seasonality. Stratification intensifies in summer, promoting the formation of bottom cold water, while in winter, the water column becomes vertically well mixed. These seasonal changes in circulation and stratification jointly regulate the spatial distribution of temperature, salinity, and currents within the basin.

Figure 1 Study area of the Harima Nada, Japan. The red star represents the site where the observational data used in this study were collected.

Based on the dataset from the Japanese Ministry of the Environment (https://water-pub.env.go.jp/water-pub/mizu-site/mizu/kouiki/dataMap.asp, last access: 2 September 2026), DO concentrations in Harima Nada range from 180 to 300 mmol m⁻³ (Fig. S1 in the Supplement). Surface DO exhibits clear seasonal variability, with the lowest values observed in autumn and the highest in winter. Controlled by hydrodynamic conditions, bottom DO reach its minimum in summer due to the intensification of stratification, while vertical mixing in winter leads to the highest bottom DO concentrations.

Chlorophyll concentrations and primary production in the water column are primarily regulated by nutrient availability (Leng et al., 2023; Tong-u-dom et al., 2025). During summer, elevated chlorophyll levels are observed in the estuarine region due to increased riverine nutrient input (Fig. S2 in the Supplement). In autumn and winter, enhanced vertical mixing redistributes DO, promoting higher chlorophyll concentrations in bottom waters. A previous study estimated the annual primary production in Harima Nada to be approximately 320 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ (Tada, 2021).

Most of the sediments in the Harima Nada are muddy, with a mud content exceeding 70 %. Such fine-grained sediments support abundant macrobenthic communities throughout Harima Nada (Tsujino, 2018; Umehara et al., 2019). Previous studies have shown that sediment properties, particularly mud content, organic matter content, and redox conditions, strongly influence the abundance and community structure of benthic organisms in this region (Nishijima et al., 2015). Through bioturbation, bioirrigation, and organic matter processing, benthic organisms can affect sediment biogeochemical cycling and nutrient transformations in this region.

Over the past several decades, the environment of the Seto Inland Sea has undergone significant changes (Ishii et al., 2014; Yamamoto et al., 2021). With the continuous observations of sediments in Harima Nada by our collaborators, we have realized that the DIN input from sediments in this area is comparable to that from rivers, exhibiting notable seasonal variations (Leng et al., 2023; Nakakuni et al., 2024). However, the driving factors and mechanisms behind these seasonal changes remain unclear. Therefore, Harima Nada serves as an excellent study area, and these observational data can provide a sufficient validation to ensure the reliability of model results (Nakakuni et al., 2024). This allows us to examine the seasonal variations in sediment nitrogen cycling and interfacial nitrogen fluxes, as well as their controlling mechanisms. Additionally, this study lays the groundwork for future exploration of how sediment nitrogen cycling responds to long-term changes in the aquatic environment.

2.2 Model descriptions

Based on the theory of early diagenesis (Berner, 1980), we have developed a biogeochemical model for sediment (Fig. 2). This model encompasses two components: solid matter and pore water. The state variables consist of the concentrations of two types of PON (fast-decayed and slow-decayed ones) for the solid matter and the concentrations of ${\mathrm{NH}}_{\mathrm{4}}^{+}$, ${\mathrm{NO}}_{\mathrm{2}}^{-}$, ${\mathrm{NO}}_{\mathrm{3}}^{-}$, DO, and ODU (the oxygen demand units) in the pore water. Throughout this study, PON refers to particulate organic nitrogen in the sediment solid phase and excludes dissolved nitrogen species in porewater. ODU is a lumped variable representing reduced metabolites generated by anoxic mineralization (e.g. reduced Fe, Mn, and sulfur species), expressed in oxygen-equivalent units to quantify their potential oxygen demand upon re-oxidation (Soetaert et al., 1996). Additionally, we set a fluff layer at the top of sediment, where the active benthos and material exchange processes with the above water column occur (Lee et al., 2002; Kuhrts et al., 2006; Laima et al., 2002).

Figure 2 Physical and biogeochemical processes for nitrogen cycle in the sediment model. Six state variables in the model are PON, ${\mathrm{NH}}_{\mathrm{4}}^{+}$, ${\mathrm{NO}}_{\mathrm{2}}^{-}$, ${\mathrm{NO}}_{\mathrm{3}}^{-}$, DO, and ODU (the oxygen demand units). PON is a solid substance (enclosed by solid lines), while ${\mathrm{NH}}_{\mathrm{4}}^{+}$, ${\mathrm{NO}}_{\mathrm{2}}^{-}$, ${\mathrm{NO}}_{\mathrm{3}}^{-}$, DO, and ODU are solutes in the pore water (enclosed by dash lines). The physical processes and biogeochemical processes are denoted by black arrows and coloured arrows, respectively. The numbers with circle denote the transport processes included in the model. The panel on the right represents the sediment separated into two distinct phases: solid and dissolved components, with porosity diminishing with increasing depth. The surface layer of the model is identified as the fluff layer.

In the solid part, the fast-decayed and slow-decayed PON are treated with a large and small mineralization rate, respectively. Both of them undergo identical vertical transport processes, including downward sedimentation driven by gravity and vertical diffusion caused by bioturbation. Their mineralization occurring within the sediment provides a flux of ${\mathrm{NH}}_{\mathrm{4}}^{+}$ into the pore water.

In the pore water, the dissolved ${\mathrm{NH}}_{\mathrm{4}}^{+}$, ${\mathrm{NO}}_{\mathrm{2}}^{-}$, ${\mathrm{NO}}_{\mathrm{3}}^{-}$, DO, and ODU are transported and diffused vertically, as corresponding to both bioturbation and molecular diffusion. The biogeochemical processes occurring within the pore water include nitrification, denitrification, anammox, dissimilatory nitrate reduction to ammonium (DNRA), ODU oxidation, and ODU burial. ODU oxidation denotes the oxygen-mediated re-oxidation of reduced compounds represented by ODU, resulting in oxygen consumption. ODU burial represents the permanent removal of these reduced compounds from the active sediment layer through downward sediment transport into deeper sediments (Soetaert et al., 1996).

At the sediment-water interface, the exchange of solid matter is controlled by the bottom shear stress. The solid matter in the water column sinks into the fluff layer when the bottom stress is below a critical value, while resuspension occurs when the bottom stress exceeds this threshold. The exchange fluxes of ${\mathrm{NH}}_{\mathrm{4}}^{+}$, ${\mathrm{NO}}_{\mathrm{2}}^{-}$, ${\mathrm{NO}}_{\mathrm{3}}^{-}$, DO, and ODU at the interface are determined by the product of the diffusion coefficient, including both bioturbation and molecular diffusion, and their concentration gradient between the fluff layer and the overlying water.

In the deepest layer of sediment, only the slow-decayed PON is allowed to move downward and leave the model domain via burial. The exchange fluxes of ${\mathrm{NH}}_{\mathrm{4}}^{+}$, ${\mathrm{NO}}_{\mathrm{2}}^{-}$, and ODU between the model domain and the underlying pore water are calculated by the product of the diffusion coefficient and the concentration gradient derived from the model results at the nearest grid point and the observations in the underlying pore water. Because ${\mathrm{NO}}_{\mathrm{3}}^{-}$ and DO are assumed absent in the underlying deep anoxic sediments, no diffusive exchange across the lower boundary is allowed, and their bottom-boundary fluxes are therefore set to zero.

2.3 Equations and parameters

The equations for the concentrations of solid matter ($C_{\mathrm{s}},\mathrm{mmol}{\mathrm{m}}^{-\mathrm{3}}$) and dissolved matter ($C_{\mathrm{d}},\mathrm{mmol}{\mathrm{m}}^{-\mathrm{3}}$) are as follows:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle \frac{\partial \left[\left(\mathrm{1}-\mathrm{\Phi }\right)\cdot C_{\mathrm{s}}\right]}{\partial t}}={\displaystyle \frac{\partial \left[\left(\mathrm{1}-\mathrm{\Phi }\right)\cdot D_{\mathrm{B}}\cdot \frac{\partial C_{\mathrm{s}}}{\partial z}\right]}{\partial z}}\\ \text{(1)}& {\displaystyle }& {\displaystyle }-{\displaystyle \frac{\partial \left[\mathit{\omega }\cdot \left(\mathrm{1}-\mathrm{\Phi }\right)\cdot C_{\mathrm{s}}\right]}{\partial z}}+\left(\mathrm{1}-\mathrm{\Phi }\right)\cdot \sum R_{\mathrm{s}}{\displaystyle }& {\displaystyle \frac{\partial \left[\mathrm{\Phi }\cdot C_{\mathrm{d}}\right]}{\partial t}}={\displaystyle \frac{\partial \left[\mathrm{\Phi }\cdot \left(D_{\mathrm{B}}+D_{\mathrm{M}}\right)\cdot \frac{\partial C_{\mathrm{d}}}{\partial z}\right]}{\partial z}}\\ \text{(2)}& {\displaystyle }& {\displaystyle }-{\displaystyle \frac{\partial \left[\mathit{\omega }\cdot \mathrm{\Phi }\cdot C_{\mathrm{d}}\right]}{\partial z}}+\mathrm{\Phi }\cdot \mathit{\alpha }\cdot \left(C_{\mathrm{d}}^{\mathrm{w}}-C_{\mathrm{d}}\right)+\mathrm{\Phi }\cdot \sum R_{\mathrm{d}}\end{array}$$

Where z (m) is vertical axis for sediment depth, t (s) is time, Φ is the sediment porosity, D_B (m² s⁻¹) is the bioturbation coefficient, D_M (m² s⁻¹) is the molecular diffusion coefficient, ω (m s⁻¹) is the sedimentation velocity, α (s⁻¹) is the bioirrigation coefficient, $C_{\mathrm{d}}^{\mathrm{w}}$ (mmol m⁻³) is the concentration of dissolved matter in the bottom water, R_s is the biochemical reaction rate of PON, and R_d is the biochemical reaction rate of dissolved matter, including ${\mathrm{NH}}_{\mathrm{4}}^{+}$, ${\mathrm{NO}}_{\mathrm{2}}^{-}$, ${\mathrm{NO}}_{\mathrm{3}}^{-}$, DO, and ODU. Since the sediment type in the study area is muddy, the advection of pore water in sandy sediments controlled by pressure is neglected.

The vertical profile of D_B is expressed as Eq. (3) (Radtke et al., 2019).

$$\begin{array}{}\text{(3)}& {\displaystyle }{\displaystyle }{D}_{\mathrm{B}}\left(z\right)=\left\{\begin{array}{ll}{D_{\mathrm{B}}}_{max},& \mathrm{0}<z<z_{max}\\ {D_{\mathrm{B}}}_{max}\cdot \exp \left(-\frac{z-z_{max}}{z_{\mathrm{d}}}\right),& z\ge z_{max}\end{array}\right.\end{array}$$

Where z_max (m) is the depth down to which the maximum bioturbation coefficient (${D_{\mathrm{B}}}_{max}$) is applied, z_d (m) is the decaying length scale of ${D_{\mathrm{B}}}_{max}$ (m² s⁻¹).

D_M of dissolved matter in the pore water is a function of porosity, and is expressed as follows (Boudreau et al., 1998; Radtke et al., 2019):

$$\begin{array}{}\text{(4)}& {\displaystyle }{\displaystyle }{D}_{\mathrm{M}}\left(z\right)={\displaystyle \frac{D^T}{\mathrm{1}-\mathrm{2.02}\ln \left(\mathrm{\Phi }\right)}}\end{array}$$

where D^T (m² s⁻¹) is the diffusion coefficient of the dissolved matter in the particle-free liquid at T °C. The relationship between D^T and temperature (T) is expressed as (Soetaert et al., 1996):

$$\begin{array}{}\text{(5)}& {\displaystyle }{\displaystyle }{D}^T=D_{\mathrm{0}}+a\cdot T\end{array}$$

where D₀ (m² s⁻¹) is the diffusion coefficient of the dissolved matter in the particle-free liquid at 0 °C, a (${\mathrm{m}}^{\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}(\mathrm{°}\mathrm{C}{)}^{-\mathrm{1}}$) is an ion-specific coefficient.

Bioirrigation characterizes the dissolved matter exchange caused by biological burrowing activities, and the coefficient decreases with depth (Dale et al., 2011):

$$\begin{array}{}\text{(6)}& {\displaystyle }{\displaystyle }\mathit{\alpha }\left(z\right)={\mathit{\alpha }}_{\mathrm{1}}\cdot {\displaystyle \frac{\exp ({\mathit{\alpha }}_{\mathrm{2}}-z)}{\mathrm{1}+\exp ({\mathit{\alpha }}_{\mathrm{2}}-z)}}\end{array}$$

The parameters α₁ (s⁻¹) and α₂ (m) represent the magnitude of bioirrigation and the depth at which it decreases to half of α₁, respectively.

The biochemical reactions are expressed in Eqs. (7)–(12), and their formulations are summarized in Table 1 (Soetaert et al., 1996; Capet et al., 2016; Akbarzadeh et al., 2018). The reaction rates are regulated by DO availability, electron acceptor competition, and temperature. Following Soetaert et al. (1996), the limitation functions associated with different mineralization pathways are normalized by their sum ($\sum lim$) to partition organic matter mineralization among oxic mineralization, denitrification, and anoxic mineralization, ensuring that the combined mineralization rate does not exceed the prescribed first-order degradation rate. In the model, temperature correction is applied only to organic matter mineralization, which serves as the primary source of reactive nitrogen and reduced substances in the sediment. The influence of temperature on other nitrogen transformation processes is represented indirectly through coupled biogeochemical interactions and temperature-dependent molecular diffusion. Following Akbarzadeh et al. (2018), several pathway-partitioning and leakage coefficients are introduced to account for the incomplete conversion of nitrogen intermediates and the distribution of products among nitrate reduction pathways. All related coefficients and parameters are listed in Table 2. The sensitivity analysis of model parameters is presented in Fig. S3 and Sect. S1 in the Supplement.

$$\begin{array}{ll}\text{(7)}& {\displaystyle }& {\displaystyle }\sum R_{\text{PON}}=-R_{\text{Min}}^{\text{OM}}-R_{\text{Den}}-R_{\text{Min}}^{\text{AM}}-R_{\text{DNRA}}{\displaystyle }& {\displaystyle }\sum R_{{\mathrm{NH}}_{\mathrm{4}}^{+}}=R_{\text{Min}}^{\text{OM}}+R_{\text{Den}}+R_{\text{Min}}^{\text{AM}}-R_{\text{Nit}\mathrm{1}}\\ \text{(8)}& {\displaystyle }& {\displaystyle }+(\mathrm{1}+\mathrm{0.5}\cdot \left(\mathrm{1}-\mathit{\delta }\right)\cdot r_{\mathrm{C}:\mathrm{N}})R_{\text{DNRA}}-R_{\text{Ana}}{\displaystyle }& {\displaystyle }\sum R_{{\mathrm{NO}}_{\mathrm{2}}^{-}}=\mathrm{0.8}\cdot \mathit{\alpha }\cdot r_{\mathrm{C}:\mathrm{N}}\cdot R_{\text{den}}+R_{\text{Nit}\mathrm{1}}-R_{\text{Nit}\mathrm{2}}\\ \text{(9)}& {\displaystyle }& {\displaystyle }+\mathrm{0.5}\cdot \mathit{\delta }\cdot r_{\mathrm{C}:\mathrm{N}}\cdot R_{\text{DNRA}}-R_{\text{Ana}}{\displaystyle }& {\displaystyle }\sum R_{{\mathrm{NO}}_{\mathrm{3}}^{-}}=R_{\text{Nit}\mathrm{2}}-\mathrm{0.8}\cdot r_{\mathrm{C}:\mathrm{N}}\cdot R_{\text{den}}\\ \text{(10)}& {\displaystyle }& {\displaystyle }-\mathrm{0.5}\cdot r_{\mathrm{C}:\mathrm{N}}\cdot R_{\text{DNRA}}{\displaystyle }& {\displaystyle }\sum R_{{\mathrm{O}}_{\mathrm{2}}}=-r_{\mathrm{C}:\mathrm{N}}\cdot R_{\text{Min}}^{\text{OM}}-\mathrm{1.5}\cdot R_{\text{Nit}\mathrm{1}}\\ \text{(11)}& {\displaystyle }& {\displaystyle }-\mathrm{0.5}\cdot R_{\text{Nit}\mathrm{2}}-R_{\text{ODUox}}\text{(12)}& {\displaystyle }& {\displaystyle }\sum R_{\text{ODU}}=r_{\mathrm{C}:\mathrm{N}}\cdot (\mathrm{1}-r_{\text{ODUsolid}})\cdot R_{\text{Min}}^{\text{AM}}-R_{\text{ODUox}}\end{array}$$

Table 1 Reaction formulas and rate expressions.

$\sum lim=\frac{{\mathrm{O}}_{\mathrm{2}}}{{\mathrm{O}}_{\mathrm{2}}+{\text{Ks}}_{{\mathrm{O}}_{\mathrm{2}}}^{\text{OM}}}+\frac{{\mathrm{NO}}_{\mathrm{3}}^{-}}{{\mathrm{NO}}_{\mathrm{3}}^{-}+{\text{Ks}}_{{\mathrm{NO}}_{\mathrm{3}}^{-}}^{\text{Den}}}\cdot \frac{{\text{Kin}}_{{\mathrm{O}}_{\mathrm{2}}}^{\text{Den}}}{{\mathrm{O}}_{\mathrm{2}}+{\text{Kin}}_{{\mathrm{O}}_{\mathrm{2}}}^{\text{Den}}}+\frac{{\text{Kin}}_{{\mathrm{O}}_{\mathrm{2}}}^{\text{AM}}}{{\mathrm{O}}_{\mathrm{2}}+{\text{Kin}}_{{\mathrm{O}}_{\mathrm{2}}}^{\text{AM}}}\cdot \frac{{\text{Kin}}_{{\mathrm{NO}}_{\mathrm{3}}^{-}}^{\text{AM}}}{{\mathrm{NO}}_{\mathrm{3}}^{-}+{\text{Kin}}_{{\mathrm{NO}}_{\mathrm{3}}^{-}}^{\text{AM}}}$. $\sum lim$ is used to rescale competing mineralization pathway, ensuring that the combined mineralization rate does not exceed the maximum first-order degradation rate (Soetaert et al., 1996).
x:y denotes the molar C:N ratio of organic matter.

Table 2 Model parameters. (L: based on literature; M: constrained with the model, I: independently determined from observations).

^a Ichimi et al. (2005); ^b Soetaert et al. (1996); ^c Akbarzadeh et al. (2018); ^d Dale et al. (2011); ^e Capet et al. (2016); ^f Laurent et al. (2016); ^g Pastor et al. (2011); ^h Ding et al. (2020); ⁱ Wang et al. (2016)

Due to the high benthic biomass and rapid organic matter turnover within the fluff layer, mineralization in this layer was assumed to be substantially enhanced relative to the underlying sediment (Lee et al., 2002; Kuhrts et al., 2006; Laima et al., 2002). The enhancement factor was treated as a tuning parameter and was finally set to 50 based on model calibration against observed sediment biogeochemical characteristics.

At the sediment-water interface, the sinking and resuspension of PON depend on the bottom stress. The exchange flux (F_s) is expressed as (Wang, 2002):

$$\begin{array}{}\text{(13)}& {\displaystyle }{\displaystyle }{F}_{\mathrm{s}}=\left\{\begin{array}{ll}w\cdot C_{\mathrm{s}}^{\mathrm{w}}\cdot \left(\frac{\mathit{\tau }}{{\mathit{\tau }}_{\mathrm{c}}}-\mathrm{1}\right),& \mathit{\tau }<{\mathit{\tau }}_{\mathrm{c}}\\ E_{\mathrm{s}}\cdot \left(\frac{\mathit{\tau }}{{\mathit{\tau }}_{\mathrm{c}}}-\mathrm{1}\right),& \mathit{\tau }\ge {\mathit{\tau }}_{\mathrm{c}}\end{array}\right.\end{array}$$

where w (m s⁻¹) is the sinking velocity of particles in the bottom water, $C_{\mathrm{s}}^{\mathrm{w}}$ (mmol m⁻³) is concentration of PON in the bottom water, τ (N m⁻²) is the bottom stress, τ_c (N m⁻²) is the critical bottom stress, E_s ($\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$) is the resuspension coefficient.

The exchange of solutes (${\mathrm{NH}}_{\mathrm{4}}^{+}$, ${\mathrm{NO}}_{\mathrm{2}}^{-}$, ${\mathrm{NO}}_{\mathrm{3}}^{-}$, DO, and ODU) are mainly through bioturbation, bioirrigation and molecular diffusion, and the flux (F_d) is expressed as (Boudreau et al., 1998):

$$\begin{array}{}\text{(14)}& {\displaystyle }{\displaystyle }{F}_{\mathrm{d}}=-\left(D_{\text{mix}}^{\mathrm{s}}+D_{\mathrm{M}}^{\mathrm{s}}\right)\cdot \mathrm{\Phi }\cdot {\displaystyle \frac{C_{\mathrm{d}}^{\mathrm{w}}-C_{\mathrm{d}}^{\text{pw}}}{\mathrm{\Delta }z}}\end{array}$$

where $D_{\text{mix}}^{\mathrm{s}}$ (m² s⁻¹) is the effective surface mixing coefficient accounting for sediment mixing processes such as bioturbation and resuspension, $D_{\mathrm{M}}^{\mathrm{s}}$ (m² s⁻¹) is the surface molecular coefficient, $C_{\mathrm{d}}^{\text{pw}}$ (mmol m⁻³) is the concentration of dissolved matter in the pore water of fluff layer, Δz (m) is the thickness of the diffusive boundary layer at sediment-water interface.

At the bottom of the model domain, the slow-decayed PON is allowed to move out as the burial process with a flux (B_s) as (Radtke et al., 2019):

$$\begin{array}{}\text{(15)}& {\displaystyle }{\displaystyle }{B}_{\mathrm{s}}=\mathit{\omega }\cdot (\mathrm{1}-\mathrm{\Phi })\cdot C_{\mathrm{s}}^{\mathrm{b}}\end{array}$$

where superscript of b represents the deepest layer of the model domain.

The bottom fluxes (B_d) of ${\mathrm{NH}}_{\mathrm{4}}^{+}$, ${\mathrm{NO}}_{\mathrm{2}}^{-}$, and ODU are expressed as:

$$\begin{array}{}\text{(16)}& {\displaystyle }{\displaystyle }{B}_{\mathrm{d}}=-\left(D_{\mathrm{B}}^{\mathrm{b}}+D_{\mathrm{M}}^{\mathrm{b}}\right)\cdot \mathrm{\Phi }\cdot {\displaystyle \frac{C_{\mathrm{d}}^{\mathrm{b}}-C_{\mathrm{d}}^{\text{bc}}}{\mathrm{\Delta }{z}^{\mathrm{b}}}}\end{array}$$

where $C_{\mathrm{d}}^{\text{bc}}$ are the concentrations in the pore water under the deepest layer of the model domain (mmol m⁻³), whose values are given from observations.

2.4 Model configurations

The model configuration, including all initial and boundary conditions used in the simulations, is described in this section.

Since this study focuses on seasonal changes at the sediment–water interface, the sediment depth in the model was set to 10 cm and uniformly divided into 100 layers.

The initial sediment concentration profile of PON was obtained from core observations in Harima Nada in July 2020, while those of ${\mathrm{NH}}_{\mathrm{4}}^{+}$, ${\mathrm{NO}}_{\mathrm{2}}^{-}$, and ${\mathrm{NO}}_{\mathrm{3}}^{-}$ were derived from observations conducted in April (Nakakuni et al., 2024). The DO profile was adopted from previous observations in the Seto Inland Sea (Sayama et al., 2002). The ODU concentration was initially set to 100 mmol m⁻³ and allowed to evolve dynamically within the model. In addition, the porosity profile was prescribed as a constant input condition throughout the simulation (Fig. S4 and Sect. S2 in the Supplement).

Monthly observations of bottom water temperature and concentrations of ${\mathrm{NH}}_{\mathrm{4}}^{+}$, ${\mathrm{NO}}_x^{-}$ (${\mathrm{NO}}_{\mathrm{2}}^{-}+{\mathrm{NO}}_{\mathrm{3}}^{-}$) and DO from April 2020 to April 2021 were used to prescribe the upper boundary conditions, while ${\mathrm{NH}}_{\mathrm{4}}^{+}$ and ${\mathrm{NO}}_{\mathrm{2}}^{-}$ concentration beneath the model domain from April 2020 to May 2021 served as the bottom boundary conditions (Nakakuni et al., 2024). The concentration of PON and ${\mathrm{NO}}_{\mathrm{2}}^{-}$ in the bottom water were collected from the seasonal dataset of Japanese Ministry of the Environment from January 2020 to May 2021. The bottom stress was derived from a hydrodynamic model for the Seto Inland Sea (Zhu et al., 2019) (Fig. 3). The ODU concentration in the bottom water was set to 0 mmol m⁻³, while its concentration in the deeper sediment was assumed to be the initial value of 100 mmol m⁻³ (Soetaert et al., 1996). Forced by these boundary conditions, the sediment model simulated the nitrogen cycle changes in the sediments of Harima Nada from April 2020 to March 2021. Additionally, the model was calculated for a long time (more than 170 years) to reach a steady-state condition, wherein the input and output fluxes of materials become balanced.

Figure 3 Model boundary conditions for the simulation period spanning April 2020–March 2021: (a) PON concentration, (b) ${\mathrm{NH}}_{\mathrm{4}}^{+}$ concentration, (c) ${\mathrm{NO}}_{\mathrm{2}}^{-}$ concentration, (d) ${\mathrm{NO}}_{\mathrm{3}}^{-}$ concentration, (e) DO concentration in the bottom water, (f) bottom water temperature, (g) bottom stress, (h) ${\mathrm{NH}}_{\mathrm{4}}^{+}$ concentration and (c) ${\mathrm{NO}}_{\mathrm{2}}^{-}$ concentration in the layer below model domain (10–12 cm). The extrapolation shown in panels (a) and (c) utilizes the observed values from January 2020 and May 2021 for interpolation, which are not displayed in this figure. For panels (h) and (i), monthly observations are smoothed using a three-point moving average to reduce intra-seasonal variability.

3 Results

3.1 Seasonal variations of the particle and dissolved nitrogen concentrations in the sediment

The monthly concentrations of PON, ${\mathrm{NH}}_{\mathrm{4}}^{+}$, DO, ${\mathrm{NO}}_{\mathrm{2}}^{-}$, and ${\mathrm{NO}}_{\mathrm{3}}^{-}$ at different layers from 0 to 10 cm in the model domain are shown in Fig. 4, alongside the corresponding observed data. The comparison between the model results and the observations, especially the seasonal variations, suggests a reliable performance of the model.

Figure 4 The seasonal variations of concentrations of (a) PON, (b) ${\mathrm{NH}}_{\mathrm{4}}^{+}$, (c) DO, (d) ${\mathrm{NO}}_x^{-}$ (${\mathrm{NO}}_{\mathrm{2}}^{-}+{\mathrm{NO}}_{\mathrm{3}}^{-}$), (e) ${\mathrm{NO}}_{\mathrm{3}}^{-}$ and (f) ${\mathrm{NO}}_{\mathrm{2}}^{-}$ in different layers in the sediment. Colours represent different depth. In all the panels, lines are model results and circles are the observation data.

The concentration of PON in the sediment was consistently greater in the surface layer compared to the underlying layers (Fig. 4a). Its seasonal variation was apparent in the upper 2 cm of sediment and almost disappears in the layers below 2 cm. In the 0–1 cm layer, the PON concentration decreased from a value of 2.43 mg g⁻¹ in April to a minimum value of 2.37 mg g⁻¹ in December, followed by an increase back to 2.43 mg g⁻¹ in March. In the deep layer of the sediment (8–10 cm), the PON concentration remained around 2.12 mg g⁻¹ throughout the year.

The ${\mathrm{NH}}_{\mathrm{4}}^{+}$ concentration had a minimum value at the surface and a maximum value at the deepest layer (Fig. 4b), which was contrary to the vertical distribution of PON concentrations. The ${\mathrm{NH}}_{\mathrm{4}}^{+}$ concentrations presented significant seasonal variations across entire sediment layers, with the range of variation increasing with depth. For example, the ${\mathrm{NH}}_{\mathrm{4}}^{+}$ concentration in the 0–1 cm layer presented the lowest value of 14.65 mmol m⁻³ in February and the highest value of 43.79 mmol m⁻³ in September; however, the ${\mathrm{NH}}_{\mathrm{4}}^{+}$ concentration in the deep layer of the sediment (8–10 cm) presented the lowest value of 106.44 mmol m⁻³ in April and the highest value of 222.35 mmol m⁻³ in October. Except for the episodic high values observed in the deep sediment layer during autumn, the model successfully reproduced the seasonal variation of ${\mathrm{NH}}_{\mathrm{4}}^{+}$ concentration across all sediment layers, which was characterized by lower levels in winter and spring and elevated levels in summer and autumn.

The DO concentration profiles exhibited a clear vertical gradient, with higher values at the surface and lower values in deeper layers (Fig. 4c). Surface DO reach a maximum of 174.51 mmol m⁻³ in January and decreased to a minimum of 8.70 mmol m⁻³ in August. The oxygen penetration depth (OPD) remained below 1 cm from June to October, and extended below 2 cm only during January and February.

${\mathrm{NO}}_{\mathrm{3}}^{-}$ was the dominant component of ${\mathrm{NO}}_x^{-}$ (Fig. 4d and e). From December to April, a subsurface maximum appeared at the layer of 1–2 cm, whereas during the rest of the year the highest concentrations occurred at the surface layer. Within the 0–4 cm layer, ${\mathrm{NO}}_{\mathrm{3}}^{-}$ exhibited a pronounced seasonal pattern, with lower concentrations in summer and higher concentrations in winter. Specifically, concentrations ranged from 1.82 to 4.76 mmol m⁻³ in the 0–1 cm layer, 0.49–8.18 mmol m⁻³ in the 1–2 cm layer, and 0.08–3.58 mmol m⁻³ in the 2–4 cm layer. Below 4 cm, the ${\mathrm{NO}}_{\mathrm{3}}^{-}$ concentrations remained close to zero throughout the year.

${\mathrm{NO}}_{\mathrm{2}}^{-}$ also showed the similar seasonal pattern, with lower concentrations in summer and higher concentrations in winter (Fig. 4f). The maximum value of 1.84 mmol m⁻³ occurred in January at the deepest layer. Vertically, the ${\mathrm{NO}}_{\mathrm{2}}^{-}$ maxima alternated between the 0–2 cm layer and the 8–10 cm layer. Notably, during January to February, the maximum concentration in the 8–10 cm layer exceeded 1 mmol m⁻³, whereas in other months the peak values remained around 1 mmol m⁻³.

3.2 The flux of PON, DIN and N-loss (N₂ and N₂O) between bottom water and sediment

In this section, we introduce the fluxes across the sediment–water interface, including the net sinking flux of PON and the exchange fluxes of solutes driven by molecular diffusion and bioturbation. Bioirrigation is not included in these fluxes because it is represented as a source–sink term distributed in the sediment column rather than as an interface flux. Its contribution is also substantially smaller than the calculated sediment–water fluxes as shown in Sect. 3.3.

At the sediment-water interface, PON sank from the bottom water to the sediment throughout the year (Fig. 5a). The PON flux exhibited a bimodal seasonal pattern, with two peaks in July (2.14 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$) and February (1.96 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$), and a minimum in November (1.50 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$).

Figure 5 The fluxes of (a) PON, (b) DIN, (c) DO, (d) N-loss (N₂ and N₂O) at the sediment-water interface in the model, and (e) the ratios of DIN and N-loss fluxes to PON flux. The red circles in (b) are the values given by observations. DIN in (b) are sum of ${\mathrm{NH}}_{\mathrm{4}}^{+}$, ${\mathrm{NO}}_{\mathrm{2}}^{-}$, and ${\mathrm{NO}}_{\mathrm{3}}^{-}$, whose fluxes are presented in different colour in the same panel. Positive values indicate upward fluxes from the sediment toward the overlying water, whereas negative values indicate downward fluxes into the sediment. In (d), ra represents the proportion of the anammox in N-loss flux.

DIN in the Harima Nada was released from the sediment into the bottom water throughout the year, exhibiting the highest flux of 1.37 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ in September and the lowest flux of 0.58 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ in February (Fig. 5b). The DIN flux and its seasonal variation given by the model exhibit a general agreement with observations (Nakakuni et al., 2024). The ${\mathrm{NH}}_{\mathrm{4}}^{+}$ flux dominated the seasonal variability of DIN flux, reaching a maximum of 1.44 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ in September and a minimum of 0.52 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ in February. In contrast, the ${\mathrm{NO}}_{\mathrm{3}}^{-}$ fluxes were substantially smaller. From June to November, ${\mathrm{NO}}_{\mathrm{3}}^{-}$ was transported from the bottom water into the sediment, with a maximum influx of 0.26 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ in August. From December to April, ${\mathrm{NO}}_{\mathrm{3}}^{-}$ was released from the sediment into the overlying water, with a peak of 0.06 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ in January. The ${\mathrm{NO}}_{\mathrm{2}}^{-}$ fluxes were negligible. A slight influx into the sediment occurred in January (0.002 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$), whereas in all other months ${\mathrm{NO}}_{\mathrm{2}}^{-}$ was released from the sediment, with a maximum of 0.02 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ in August. The DO flux continuously supplies DO from the bottom water to the sediment and is the only source of DO for the sediment. It reached the highest value in September at 10.06 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ and its lowest in January at 5.37 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ (Fig. 5c).

N-loss is produced through denitrification and anammox inside the sediment, generating N₂ and N₂O that are released to the overlying water. The N-loss flux across the sediment-water interface reached its peak of 0.14 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ in August and the lowest value of 0.03 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ in March (Fig. 5d). Denitrification-derived N₂ and N₂O constituted the dominant components of N-loss flux and controlled its seasonal variability, with values ranging from 0.02 to 0.13 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$. In contrast, the N₂ production via anammox showed weak seasonal variation, remaining approximately constant at 0.01 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$. Consequently, the relative contribution of anammox to N-loss flux (ra %) was lowest in August (8 %) and highest in March (41 %). Although N₂ was the main product of denitrification, the flux of N₂O, a potent greenhouse gas, was not negligible. It varied from 0.19 to 1.27 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$, with a mean value of 0.44 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$.

The ratio of DIN flux to PON flux (R_DIN/PON) reflects the fraction of nitrogen that sinks into the sediment and later releases back to the seawater (Fig. 5e). R_DIN/PON increased from 39 % in April to a peak of 69 % in September, with the annual average of 48 %. The ratio of N-loss flux to PON flux (R_N-loss/PON) represented the proportion of nitrogen removed from the sediment. The highest value of R_N-loss/PON was 7 % in August while the lowest value was 2 % in May. The sum of R_DIN/PON and R_N-loss/PON ranged from 32 % in February to 73 % in September, with an annual average of 51 %. It suggested that around 50 % of the PON flux sinking from seawater was retained in the sediment on average, with this retention rate potentially reaching nearly 70 % during certain months.

3.3 The budgets of PON, ${\mathrm{NH}}_{\mathrm{4}}^{+}$, and ${\mathrm{NO}}_{\mathrm{3}}^{-}$ in the sediment

The annual budgets for the four forms of nitrogen in the sediment of Harima Nada are estimated from the model results (Fig. 6). To facilitate comparison with the fluxes presented in the preceding sections, all fluxes and reaction rates are expressed as annual mean daily values ($\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$). The sinking of PON from bottom water was the source of PON in the sediment, quantified at 1.86 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$. Since the primary production in the Harima Nada in 2020 was estimated to be 8.02 $\mathrm{mmol}\mathrm{N}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ (Tada, 2021), it could be interpreted that about 23 % of the particles produced by primary production sank into the sediment. Within the deeper layers of the sediment, the burial rate of PON was 0.72 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$, representing 39 % of the PON flux into the sediment and 9 % of the primary production. The rest 1.12 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ of PON flux into the sediment was mineralized to ${\mathrm{NH}}_{\mathrm{4}}^{+}$, which became the main source of ${\mathrm{NH}}_{\mathrm{4}}^{+}$ in the sediment.

Figure 6 Nitrogen budgets in the sediment model of Harima Nada. All fluxes and reaction rates are expressed as annual mean daily values ($\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$).

In addition to mineralization, ${\mathrm{NH}}_{\mathrm{4}}^{+}$ was also produced via denitrification and DNRA in the model, but the contributions from these two processes were considerably lower. Additionally, a flux of 0.06 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ of ${\mathrm{NH}}_{\mathrm{4}}^{+}$ was transferred to the deeper sediment located beneath the model domain. Furthermore, ${\mathrm{NH}}_{\mathrm{4}}^{+}$ has four sinks: the nitrification of ${\mathrm{NH}}_{\mathrm{4}}^{+}$ to ${\mathrm{NO}}_{\mathrm{2}}^{-}$ had a rate of 0.07 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$, the anammox of ${\mathrm{NH}}_{\mathrm{4}}^{+}$ to N₂ had a rate of 0.01 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$, the upward ${\mathrm{NH}}_{\mathrm{4}}^{+}$ flux through the sediment-water interface was quantified at 0.93 and 0.07 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ of ${\mathrm{NH}}_{\mathrm{4}}^{+}$ was released from the sediment by bioirrigation.

At the sediment–water interface, ${\mathrm{NO}}_{\mathrm{2}}^{-}$ was released into the overlying water with a flux of 0.01 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$, while bioirrigation contributed less than 0.01 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$. Nitrification converted 0.06 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ of ${\mathrm{NO}}_{\mathrm{2}}^{-}$ into ${\mathrm{NO}}_{\mathrm{3}}^{-}$. ${\mathrm{NO}}_{\mathrm{3}}^{-}$ inputs from the bottom water to the sediment occurred via diffusive exchange and bioirrigation, amounting to 0.04 and less than 0.01 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$, respectively. DNRA consumed less than 0.01 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ of ${\mathrm{NO}}_{\mathrm{3}}^{-}$, whereas denitrification produced 0.09 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ of N₂ and N₂O. In total, nitrogen removal processes (denitrification and anammox) generated 0.10 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ of N₂ and N₂O, corresponding to approximately 6 % of the PON sinking flux. Overall, bioirrigation is much weaker than the diffusive flux of solutes between sediments and water, and their seasonal variations shown in Fig. S5 and Sect. S3 in the Supplement.

4 Discussion

4.1 Factors affecting the seasonal variations of nitrogen inventories in the sediment

The model results show that the concentrations of PON, ${\mathrm{NH}}_{\mathrm{4}}^{+}$, ${\mathrm{NO}}_{\mathrm{2}}^{-}$, ${\mathrm{NO}}_{\mathrm{3}}^{-}$ and DO in the sediment exhibit apparent seasonal variations. To clarify the factors controlling these seasonal variations, we carried out seven sensitivity experiments, in which we removed the seasonal variation of one of the following seven variables by substituting its monthly means with its annual mean: PON concentration (PON^bw), ${\mathrm{NH}}_{\mathrm{4}}^{+}$ concentration (${{\mathrm{NH}}_{\mathrm{4}}^{+}}^{\text{bw}}$), ${\mathrm{NO}}_{\mathrm{3}}^{-}$ concentration (${{\mathrm{NO}}_{\mathrm{3}}^{-}}^{\text{bw}}$), ${\mathrm{NO}}_{\mathrm{2}}^{-}$ concentration (${{\mathrm{NO}}_{\mathrm{2}}^{-}}^{\text{bw}}$), DO concentration (DO^bw), and water temperature in the bottom water, as well as bottom stress. Figure 7 shows the calculation results of the seven sensitivity experiments alongside the results described in Sect. 3 (referred to “Control” in Fig. 7).

Figure 7 Comparisons of the model results among eight calculations (seven sensitivity experiments and one control run). The sensitivity experiments include the cases removing the seasonal variations of PON concentration (dark green line, “Mean PON”), ${\mathrm{NH}}_{\mathrm{4}}^{+}$ concentration (light green line, “Mean ${\mathrm{NH}}_{\mathrm{4}}^{+}$”), ${\mathrm{NO}}_{\mathrm{3}}^{-}$ concentration (dark blue line, “Mean ${\mathrm{NO}}_{\mathrm{3}}^{-}$”), ${\mathrm{NO}}_{\mathrm{2}}^{-}$ concentration (light blue line, “Mean ${\mathrm{NO}}_{\mathrm{2}}^{-}$”), DO concentration (pink line, “Mean DO”) and water temperature (red line, “Mean temperature”) in the overlying water, as well as of the bottom stress (orange line, “Mean bottom stress”). (a) PON inventory; (b) ${\mathrm{NH}}_{\mathrm{4}}^{+}$ inventory; (c) ${\mathrm{NO}}_{\mathrm{3}}^{-}$ inventory; (d) ${\mathrm{NO}}_{\mathrm{2}}^{-}$ inventory; (e) DO inventory.

The concentration of PON^bw, bottom stress and the bottom water temperature significantly influence the seasonal variation of PON inventory in the sediment (Fig. 7a). Sinking of PON^bw serves as the source of PON in the sediment, and its temporal variations are controlled by those in the concentration of PON^bw and bottom stress (Fig. 8a). The variations in bottom water temperature can change the mineralization rate in the sediment (Fig. S6a and Sect. S4 in the Supplement). Although these seasonal variations can directly affect the PON inventory, the corresponding seasonal fluctuation in PON inventory represents only a very small fraction of the total sedimentary PON pool.

Among the seven factors examined, bottom water temperature is the dominant driver of the seasonal variability in the ${\mathrm{NH}}_{\mathrm{4}}^{+}$ inventory (Fig. 7b) because it directly regulates the mineralization rate that produces ${\mathrm{NH}}_{\mathrm{4}}^{+}$ in the sediment (Fig. S6 and Sect. S4 in the Supplement). Although mineralization depends on both organic matter availability and temperature, seasonal variations in PON supply induce only minor perturbations to the sedimentary PON reservoir in the model domain. Consequently, the amount of organic matter available for mineralization remains relatively stable on seasonal timescales. In contrast, seasonal changes in temperature substantially alter mineralization rates and therefore ${\mathrm{NH}}_{\mathrm{4}}^{+}$ production in the whole sediment domain. Since mineralization is the primary source of ${\mathrm{NH}}_{\mathrm{4}}^{+}$ in the sediment (Fig. 6), bottom water temperature ultimately exerts the strongest control on the seasonal variability of ${\mathrm{NH}}_{\mathrm{4}}^{+}$ inventory.

The ${\mathrm{NO}}_{\mathrm{3}}^{-}$ inventory is primarily controlled by bottom water temperature and DO^bw concentration, while the influence of ${{\mathrm{NO}}_{\mathrm{3}}^{-}}^{\text{bw}}$ concentration is comparatively weaker (Fig. 7c). During summer, high temperatures enhance denitrification and increase ${\mathrm{NO}}_{\mathrm{3}}^{-}$ consumption (Fig. S6d in the Supplement). At the same time, reduced DO supply from the overlying water suppresses nitrification (Fig. S6b and c in the Supplement), limiting the replenishment of sedimentary ${\mathrm{NO}}_{\mathrm{3}}^{-}$. The combined effects of enhanced consumption and reduced production result in relatively low ${\mathrm{NO}}_{\mathrm{3}}^{-}$ concentrations during summer. In winter, lower temperatures weaken denitrification, while higher DO availability promotes nitrification, allowing ${\mathrm{NO}}_{\mathrm{3}}^{-}$ to accumulate within the sediment. Variations in ${{\mathrm{NO}}_{\mathrm{3}}^{-}}^{\text{bw}}$ concentration can further modify sedimentary ${\mathrm{NO}}_{\mathrm{3}}^{-}$ levels through molecular diffusion across the sediment–water interface, although their influence is smaller than that of temperature and oxygen availability. Therefore, unlike ${\mathrm{NH}}_{\mathrm{4}}^{+}$, which is primarily regulated by mineralization, the ${\mathrm{NO}}_{\mathrm{3}}^{-}$ inventory is jointly controlled by nitrification and denitrification and is consequently more sensitive to changes in sediment redox conditions. As an intermediate product of nitrification, ${\mathrm{NO}}_{\mathrm{2}}^{-}$ is primarily controlled by DO^bw concentration (Fig. 7d).

Since sedimentary DO is supplied exclusively from the overlying water, DO concentrations within the sediment depend largely on DO^bw(Fig. 7e). In winter, lower temperatures reduce the consumption of DO through mineralization, allowing more DO to be preserved within the sediment. This also explains the deeper OPD developed in winter.

The different controlling factors identified for PON, ${\mathrm{NH}}_{\mathrm{4}}^{+}$, ${\mathrm{NO}}_{\mathrm{3}}^{-}$ and DO inventories reflect the hierarchical structure of nitrogen cycling within the sediment. PON serves as the primary nitrogen reservoir and provides the substrate for mineralization. However, because the sedimentary PON inventory is substantially larger than the seasonal variation in PON deposition, fluctuations in external PON supply only weakly affect the amount of organic matter available for decomposition. Consequently, the seasonal variability of ${\mathrm{NH}}_{\mathrm{4}}^{+}$ inventory is governed mainly by temperature-dependent mineralization rather than by changes in PON input. ${\mathrm{NH}}_{\mathrm{4}}^{+}$ produced through mineralization subsequently fuels nitrification. As an intermediate product of nitrification, ${\mathrm{NO}}_{\mathrm{2}}^{-}$ remains at relatively low concentrations but responds sensitively to oxygen availability because its production depends directly on nitrification activity. ${\mathrm{NO}}_{\mathrm{3}}^{-}$ is further generated through nitrification and consumed through denitrification. Therefore, the ${\mathrm{NO}}_{\mathrm{3}}^{-}$ inventory is controlled by the balance between nitrification and denitrification and is particularly sensitive to temperature and oxygen availability. These results indicate that the seasonal variability of sedimentary nitrogen cycling is regulated through a cascade of coupled biogeochemical processes, in which different environmental factors affect different stages of the nitrogen cycle with varying strengths. Among them, temperature exerts the strongest influence because it simultaneously affects mineralization, nitrification, denitrification and molecular diffusion, thereby propagating through the entire nitrogen cycle. In contrast, variations in external PON supply primarily influence the sedimentary organic matter pool and have limited effects on downstream nitrogen transformations because the large sedimentary PON reservoir buffers short-term fluctuations in organic matter deposition.

4.2 Factors affecting the seasonal variations in the sediment-water fluxes of nitrogen

The seasonal variation in PON flux across the sediment–water interface is primarily controlled by the PON^bw concentration and bottom stress (Fig. 8a). Elevated PON^bw concentrations enhance the downward settling flux of particulate organic matter into the sediment. In contrast, an inverse seasonal relationship is observed between the PON flux and bottom stress (Fig. 3h), indicating that stronger hydrodynamic forcing during winter and spring promotes sediment resuspension and reduces the net deposition of PON. These results suggest that the seasonal variability of PON flux is governed mainly by physical processes occurring at the sediment–water interface.

Figure 8 Comparisons of the model results among eight calculations (seven sensitivity experiments and one control run). The sensitivity experiments include the cases removing the seasonal variations of PON concentration (dark green line, “Mean PON”), ${\mathrm{NH}}_{\mathrm{4}}^{+}$ concentration (light green line, “Mean ${\mathrm{NH}}_{\mathrm{4}}^{+}$”), ${\mathrm{NO}}_{\mathrm{3}}^{-}$ concentration (dark blue line, “Mean ${\mathrm{NO}}_{\mathrm{3}}^{-}$”), ${\mathrm{NO}}_{\mathrm{2}}^{-}$ concentration (light blue line, “Mean ${\mathrm{NO}}_{\mathrm{2}}^{-}$”), DO concentration (pink line, “Mean DO”) and water temperature (red line, “Mean temperature”) in the overlying water, as well as of the bottom stress (orange line, “Mean bottom stress”). (a) PON flux; (b) ${\mathrm{NH}}_{\mathrm{4}}^{+}$ flux; (c) ${\mathrm{NO}}_{\mathrm{3}}^{-}$ flux; (d) ${\mathrm{NO}}_{\mathrm{2}}^{-}$ flux; (e) DO flux; (f) N-loss flux; (g) ratio of DIN flux to PON flux; (h) ratio of N-loss to PON flux.

In contrast to the PON flux, the seasonal variability of ${\mathrm{NH}}_{\mathrm{4}}^{+}$ flux is primarily controlled by bottom water temperature (Fig. 8b). Because the ${\mathrm{NH}}_{\mathrm{4}}^{+}$ concentrations within the sediment are substantially higher than those in the bottom water, seasonal changes in sedimentary ${\mathrm{NH}}_{\mathrm{4}}^{+}$ inventory directly determine the magnitude of ${\mathrm{NH}}_{\mathrm{4}}^{+}$ release flux (Figs. 3b and 4b). As demonstrated in Sect. 4.1, temperature strongly regulates mineralization rates and therefore controls ${\mathrm{NH}}_{\mathrm{4}}^{+}$ production within the sediment. Consequently, temperature-driven changes in sedimentary ${\mathrm{NH}}_{\mathrm{4}}^{+}$ inventory dominate the seasonal variability of ${\mathrm{NH}}_{\mathrm{4}}^{+}$ flux.

The seasonal variability of ${\mathrm{NO}}_{\mathrm{3}}^{-}$ flux is governed by a combination of bottom water temperature, DO^bw, and ${{\mathrm{NO}}_{\mathrm{3}}^{-}}^{\text{bw}}$ concentration (Fig. 8c). Higher temperatures during summer enhance mineralization and increase ${\mathrm{NH}}_{\mathrm{4}}^{+}$ production, thereby supplying substrate for nitrification. Elevated temperatures also accelerate molecular diffusion and strengthen solute exchange across the sediment–water interface. However, oxygen depletion within the sediment during summer suppresses nitrification and maintains low sedimentary ${\mathrm{NO}}_{\mathrm{3}}^{-}$ concentrations. This favors the diffusion of ${\mathrm{NO}}_{\mathrm{3}}^{-}$ from the overlying water into the sediment. The relatively high ${{\mathrm{NO}}_{\mathrm{3}}^{-}}^{\text{bw}}$ concentration during summer further enhances this downward flux. Compared with ${\mathrm{NH}}_{\mathrm{4}}^{+}$ and ${\mathrm{NO}}_{\mathrm{3}}^{-}$, the ${\mathrm{NO}}_{\mathrm{2}}^{-}$ fluxes are much smaller but exhibit high sensitivity to bottom water temperature as well as variations in ${{\mathrm{NO}}_{\mathrm{3}}^{-}}^{\text{bw}}$, ${{\mathrm{NO}}_{\mathrm{2}}^{-}}^{\text{bw}}$, and DO^bw concentrations (Fig. 8d), reflecting the role of ${\mathrm{NO}}_{\mathrm{2}}^{-}$ as a transient intermediate in coupled nitrification–denitrification processes.

The contrasting controls on PON and DIN fluxes reveal an important characteristic of sedimentary nitrogen cycling. Although PON deposition supplies organic matter to the sediment, seasonal changes in PON flux do not lead to proportional changes in DIN release. The depositional amount of PON on monthly timescales is several orders of magnitude smaller than the existing sedimentary PON inventory and therefore represents only a minor perturbation to the sediment organic matter reservoir. This large internal reservoir buffers the influence of short-term fluctuations in organic matter deposition, causing sedimentary mineralization rates to respond much more strongly to temperature rather than to seasonal variations in PON input. As a result, sediment DIN fluxes are controlled primarily by internal biogeochemical processing rather than by immediate organic matter supply. Since the ${\mathrm{NH}}_{\mathrm{4}}^{+}$ flux constitutes the dominant component of DIN flux, temperature-driven mineralization ultimately becomes the principal driver of seasonal DIN release. These findings suggest that approaches estimating sediment DIN fluxes solely from PON deposition or near-bottom PON concentrations may substantially misrepresent seasonal variability if sediment biogeochemical transformations and their temperature dependence are not explicitly considered.

In addition to temperature, oxygen availability also plays a critical role in regulating nitrogen transformations within the sediment. The DO flux is also largely controlled by bottom water temperature (Fig. 8e). Temperature influences the DO flux through two pathways. First, it regulates molecular diffusion and therefore affects DO transport across the sediment–water interface. Second, it controls oxygen consumption within the sediment through oxic mineralization and nitrification, which together account for most of the sediment oxygen demand. Seasonal changes in temperature therefore simultaneously affect both DO transport and consumption. In addition, the decline in DO^bw concentrations associated with summer stratification in Harima Nada reduces the concentration gradient across the sediment–water interface and further limits oxygen uptake by the sediment.

N-loss from the sediment is dominated by denitrification (Fig. 6). However, low sedimentary ${\mathrm{NO}}_{\mathrm{3}}^{-}$ concentrations during summer limit the availability of electron acceptors for denitrification, making external ${\mathrm{NO}}_{\mathrm{3}}^{-}$ supply increasingly important. Consequently, the enhanced downward flux of ${\mathrm{NO}}_{\mathrm{3}}^{-}$ from the bottom water contributes directly to the summer maximum in denitrification and N-loss flux. The seasonal variability of N-loss therefore reflects controls similar to those governing the ${\mathrm{NO}}_{\mathrm{3}}^{-}$ flux, namely bottom water temperature, DO^bw, and ${{\mathrm{NO}}_{\mathrm{3}}^{-}}^{\text{bw}}$ concentration (Fig. 8f). The ratios R_DIN/PON and R_N-loss/PON provide useful indicators of nitrogen recycling and removal efficiency in semi-enclosed inland seas. The seasonal variation of R_DIN/PON is mainly regulated by bottom water temperature and bottom stress (Fig. 8g), whereas R_N-loss/PON is controlled primarily by bottom water temperature, DO^bw, and ${{\mathrm{NO}}_{\mathrm{3}}^{-}}^{\text{bw}}$ concentration (Fig. 8h).

Overall, the seasonal dynamics of nitrogen cycling at the sediment–water interface exhibit a transition from physically controlled particulate nitrogen deposition to biogeochemically controlled dissolved nitrogen recycling and nitrogen removal. These results demonstrate that the seasonal variability of sedimentary nitrogen recycling is governed primarily by internal biogeochemical processes, particularly temperature-dependent mineralization and coupled nitrification–denitrification dynamics, rather than by short-term fluctuations in organic matter deposition.

4.3 Comparisons with other coastal seas

Numerous studies have quantified nitrogen fluxes across the sediment–water interface in coastal and shelf environments, revealing substantial variability in both the magnitude and direction of individual nitrogen fluxes (Zhou et al., 2017; Sun et al., 2021; Canion et al., 2014; Lin et al., 2017; De Vittor et al., 2012; Mu et al., 2017; Rich et al., 2020; Dale et al., 2022; Bohlen et al., 2011; Ratmaya et al., 2022; Zhou et al., 2022). Despite this variability, the reported fluxes can be interpreted within a common mechanistic framework. Specifically, the ${\mathrm{NH}}_{\mathrm{4}}^{+}$ flux is primarily controlled by the balance between sedimentary ${\mathrm{NH}}_{\mathrm{4}}^{+}$ production and ${\mathrm{NH}}_{\mathrm{4}}^{+}$ availability in bottom waters, the ${\mathrm{NO}}_{\mathrm{3}}^{-}$ flux by the balance among nitrification, denitrification, and ${\mathrm{NO}}_{\mathrm{3}}^{-}$ supply from the overlying water, and the N-loss flux by the availability of reactive nitrogen, particularly ${\mathrm{NO}}_{\mathrm{3}}^{-}$, to support denitrification. The fluxes simulated for Harima Nada are discussed below within this framework and compared with observations from other coastal environments (Table 3).

Table 3 Benthic fluxes ($\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$) reported in some regions; positive values represent the release from sediment to water.

${\mathrm{NH}}_{\mathrm{4}}^{+}$ is consistently released from the sediment to the overlying water in Harima Nada. This upward flux is driven by continuous mineralization of organic matter, which maintains sedimentary ${\mathrm{NH}}_{\mathrm{4}}^{+}$ concentrations substantially higher than those in bottom waters and thereby sustains an upward concentration gradient across the sediment–water interface. Similar upward ${\mathrm{NH}}_{\mathrm{4}}^{+}$ fluxes have been reported in the northern Adriatic Sea, Vilaine Bay, the Mauritanian upwelling region, and the Peruvian upwelling region (De Vittor et al., 2012; Ratmaya et al., 2022; Dale et al., 2022; Bohlen et al., 2011). In contrast, downward ${\mathrm{NH}}_{\mathrm{4}}^{+}$ fluxes were observed in Bohai Bay, the Yellow Sea and East China Sea, and areas off eastern Taiwan and have been associated with elevated ${\mathrm{NH}}_{\mathrm{4}}^{+}$ concentrations in bottom waters or relatively weak sedimentary ${\mathrm{NH}}_{\mathrm{4}}^{+}$ accumulation, which reduce or reverse the concentration gradient across the sediment–water interface (Mu et al., 2017; Zhou et al., 2017; Zhou et al., 2022). These observations suggest that the direction of ${\mathrm{NH}}_{\mathrm{4}}^{+}$ flux is generally determined by the balance between sedimentary ${\mathrm{NH}}_{\mathrm{4}}^{+}$ production and ${\mathrm{NH}}_{\mathrm{4}}^{+}$ concentrations in the overlying water.

The ${\mathrm{NO}}_{\mathrm{3}}^{-}$ flux exhibits pronounced seasonal variability in both magnitude and direction. In Harima Nada, this variability reflects shifts in the balance among nitrification, denitrification, and ${\mathrm{NO}}_{\mathrm{3}}^{-}$ supply from the overlying water. During summer, strong stratification promotes hypoxia and suppresses nitrification, thereby reducing in situ ${\mathrm{NO}}_{\mathrm{3}}^{-}$ production. Meanwhile, denitrification continues to consume ${\mathrm{NO}}_{\mathrm{3}}^{-}$ and bottom-water ${\mathrm{NO}}_{\mathrm{3}}^{-}$ concentrations reach seasonal maxima. As a result, sedimentary ${\mathrm{NO}}_{\mathrm{3}}^{-}$ concentrations become lower than those in the overlying water, producing a downward ${\mathrm{NO}}_{\mathrm{3}}^{-}$ flux. In winter, improved oxygen conditions enhance nitrification and reduce ${\mathrm{NO}}_{\mathrm{3}}^{-}$ depletion within the sediment, leading to an upward ${\mathrm{NO}}_{\mathrm{3}}^{-}$ flux. Similar seasonal reversals have been reported in Vilaine Bay, where temporal changes in oxygen conditions and nitrogen transformation processes lead to alternating periods of ${\mathrm{NO}}_{\mathrm{3}}^{-}$ release and uptake (Ratmaya et al., 2022). In contrast, systems such as the Mauritanian upwelling region, the Peruvian upwelling region, and areas off eastern Taiwan exhibit persistent sedimentary uptake of ${\mathrm{NO}}_{\mathrm{3}}^{-}$ because continuous ${\mathrm{NO}}_{\mathrm{3}}^{-}$ consumption within sediments, together with relatively high bottom-water ${\mathrm{NO}}_{\mathrm{3}}^{-}$ concentrations, maintains a persistent concentration gradient that drives ${\mathrm{NO}}_{\mathrm{3}}^{-}$ diffusion from the water column into the sediment (Bohlen et al., 2011; Zhou et al., 2017; Dale et al., 2022). Conversely, predominantly upward ${\mathrm{NO}}_{\mathrm{3}}^{-}$ fluxes reported in Bohai Bay and Taranto Gulf suggest that sedimentary ${\mathrm{NO}}_{\mathrm{3}}^{-}$ production exceeds consumption, allowing porewater ${\mathrm{NO}}_{\mathrm{3}}^{-}$ concentrations to be higher than those in bottom waters (Mu et al., 2017; De Vittor et al., 2015). These comparisons indicate that differences among coastal systems can largely be interpreted as shifts in the balance between sedimentary ${\mathrm{NO}}_{\mathrm{3}}^{-}$ production, sedimentary ${\mathrm{NO}}_{\mathrm{3}}^{-}$ consumption, and external ${\mathrm{NO}}_{\mathrm{3}}^{-}$ supply.

The N-loss flux in Harima Nada accounts for approximately 6 % of the sinking PON flux and is substantially lower than values reported in many other coastal and shelf seas, including the East China Sea, Chinese marginal seas, and the Peru margin (Sun et al., 2021; Lin et al., 2017; Rich et al., 2020). Our simulations indicate that denitrification is the dominant pathway of nitrogen removal, consistent with numerous previous studies. Consequently, the relatively low N-loss flux in Harima Nada can be largely attributed to limited ${\mathrm{NO}}_{\mathrm{3}}^{-}$ availability. Because both sedimentary and bottom-water ${\mathrm{NO}}_{\mathrm{3}}^{-}$ concentrations remain comparatively low, denitrification is constrained by the supply of its principal reactant, resulting in reduced N₂ production. In contrast, previous studies have reported substantially higher N-loss fluxes in the sediment environments characterized by abundant organic matter inputs and greater ${\mathrm{NO}}_{\mathrm{3}}^{-}$ availability, which support more intensive denitrification activity (Lin et al., 2017; Rich et al., 2020; Sun et al., 2021). These comparisons suggest that the magnitude of N-loss reflects the combined effects of organic matter supply, oxygen availability, and ${\mathrm{NO}}_{\mathrm{3}}^{-}$ availability. In Harima Nada, however, the relatively low ${\mathrm{NO}}_{\mathrm{3}}^{-}$ concentrations in both sediments and bottom waters appear to be the primary factor constraining denitrification and limiting N₂ production. In addition, the ${\mathrm{NO}}_{\mathrm{2}}^{-}$ concentrations remain low in the sediment and are insufficient to sustain high rates of anammox. Consequently, ${\mathrm{NO}}_{\mathrm{2}}^{-}$ availability primarily controls the seasonal variability of the relative contribution of anammox to total N-loss. During summer hypoxia, reduced ${\mathrm{NO}}_{\mathrm{2}}^{-}$ production suppresses anammox but favours denitrification, resulting in the lowest relative contribution of anammox during this period.

Taken together, the above comparisons suggest that, despite large differences in environmental settings and flux magnitudes among coastal systems, a common set of mechanisms governs benthic nitrogen cycling. The ${\mathrm{NH}}_{\mathrm{4}}^{+}$ flux is primarily controlled by the balance between organic matter mineralization and ${\mathrm{NH}}_{\mathrm{4}}^{+}$ accumulation in bottom waters, whereas the ${\mathrm{NO}}_{\mathrm{3}}^{-}$ flux reflects the balance among nitrification, denitrification, and external ${\mathrm{NO}}_{\mathrm{3}}^{-}$ supply. In turn, the magnitude of N-loss is strongly constrained by ${\mathrm{NO}}_{\mathrm{3}}^{-}$ availability because denitrification constitutes the dominant nitrogen removal pathway. These mechanisms consistently explain both the seasonal variability observed in Harima Nada and the differences reported among other coastal environments, indicating that they may represent general controls on sedimentary nitrogen cycling in coastal ecosystems.

While these results provide insights into the mechanisms controlling sedimentary nitrogen cycling in Harima Nada, several limitations should be acknowledged. First, the model is one-dimensional and therefore does not resolve spatial heterogeneity in sediment properties, organic matter deposition, and benthic biogeochemical processes across Harima Nada. Second, direct observations of individual benthic nitrogen fluxes, including ${\mathrm{NH}}_{\mathrm{4}}^{+}$, ${\mathrm{NO}}_{\mathrm{2}}^{-}$, ${\mathrm{NO}}_{\mathrm{3}}^{-}$, and N-loss fluxes, are not available for independent validation. Consequently, uncertainties remain in the quantitative estimates of individual transformation pathways despite the generally good agreement between simulations and available observations. Nevertheless, the model provides a process-based framework linking organic matter supply, oxygen availability, and bottom-water nutrient conditions to sediment nitrogen fluxes. Future studies combining direct benthic flux measurements with spatially resolved sediment biogeochemical models would help evaluate the generality of the controlling mechanisms identified here and further constrain regional nitrogen budgets.

5 Conclusions

This study quantified the seasonal variations in nitrogen fluxes across the sediment–water interface in Harima Nada using a one-dimensional sediment nitrogen cycling model and identified the key environmental factors controlling these variations. The results showed that PON deposition consistently represented the dominant pathway of nitrogen input to the sediment, greatly exceeding the contribution from downward ${\mathrm{NO}}_{\mathrm{3}}^{-}$ fluxes. The seasonal variability of PON deposition was primarily controlled by bottom-water PON concentration and bottom stress.

Although PON deposition supplied the organic matter supporting nitrogen mineralization, seasonal variations in DIN fluxes were largely decoupled from seasonal variations in PON deposition. The large sedimentary PON reservoir buffered short-term fluctuations in organic matter input, which results in only minor changes in sedimentary nitrogen inventories. Instead, temperature-dependent biogeochemical processes exerted stronger control on sediment–water DIN exchange. Continuous mineralization maintained an upward ${\mathrm{NH}}_{\mathrm{4}}^{+}$ flux throughout the year, and its seasonal variability was primarily regulated by bottom-water temperature through its influence on mineralization rates and sedimentary ${\mathrm{NH}}_{\mathrm{4}}^{+}$ inventory. In contrast, the magnitude and direction of ${\mathrm{NO}}_{\mathrm{3}}^{-}$ flux varied seasonally under the combined influence of temperature-dependent nitrogen transformation processes and bottom-water ${\mathrm{NO}}_{\mathrm{3}}^{-}$ concentrations. These findings suggest that the approaches estimating sediment-water DIN exchange flux solely from PON deposition, while neglecting temperature-dependent sediment processes and bottom-water nitrogen concentrations, may substantially misrepresent the seasonal dynamics of DIN fluxes.

Denitrification was identified as the dominant pathway of nitrogen removal and was primarily responsible for the seasonal variability of N-loss. Although ${\mathrm{NO}}_{\mathrm{3}}^{-}$ was continuously regenerated within the sediment through nitrification, sedimentary ${\mathrm{NO}}_{\mathrm{3}}^{-}$ concentrations remained relatively low under the oligotrophic conditions of Harima Nada and were frequently supplemented by ${\mathrm{NO}}_{\mathrm{3}}^{-}$ supplied from the overlying water. Consequently, denitrification was constrained by limited ${\mathrm{NO}}_{\mathrm{3}}^{-}$ availability, resulting in relatively low N-loss rates. This finding highlights the importance of ${\mathrm{NO}}_{\mathrm{3}}^{-}$ availability as a key regulator of nitrogen removal in nitrogen-limited coastal environments.

Collectively, the results demonstrate that seasonal benthic nitrogen fluxes cannot be inferred solely from organic matter deposition. Instead, sedimentary nitrogen cycling is governed by the interaction among organic matter supply, oxygen availability, and bottom-water nutrient concentrations. Within this framework, organic matter deposition controls nitrogen inputs to the sediment, oxygen conditions regulate internal nitrogen transformation pathways, and bottom-water nutrient concentrations influence the magnitude and direction of sediment–water nutrient fluxes. The balance among these factors ultimately determines the seasonal variability of benthic nitrogen fluxes and nitrogen removal.

Several limitations should be acknowledged. The present model is one-dimensional and therefore does not resolve spatial heterogeneity in sediment properties and benthic biogeochemical processes. In addition, direct observations of individual benthic nitrogen fluxes and nitrogen removal rates are still lacking in Harima Nada, introducing uncertainty into the quantitative estimates of specific nitrogen transformation pathways. Furthermore, the model was developed primarily for muddy sediments and does not explicitly account for porewater advection, which may become an important transport process in highly permeable sandy sediments. Application of the model to such environments would therefore require appropriate correction of porewater flow velocity and associated solute transport processes. Nevertheless, the model successfully reproduces the observed seasonal variations in sediment nitrogen concentrations and DIN fluxes, providing a process-based framework for understanding sedimentary nitrogen cycling in semi-enclosed inland seas dominated by muddy sediments.

Future integration of this modelling framework with long-term environmental observations will provide a valuable tool for assessing the impacts of eutrophication, de-eutrophication, and climate-driven warming on sediment nitrogen cycling, sediment–water nutrient exchange, and the capacity of coastal sediments to remove reactive nitrogen, thereby improving our understanding of long-term biogeochemical responses to environmental change in semi-enclosed coastal seas.