A 1-Dimensional Sympagic–Pelagic–Benthic Transport Model (SPBM): Coupled Simulation of Ice, Water Column, and Sediment Biogeochemistry, Suitable for Arctic Applications

Abstract

Marine biogeochemical processes can strongly interact with processes occurring in adjacent ice and sediments. This is especially likely in areas with shallow water and frequent ice cover, both of which are common in the Arctic. Modeling tools are therefore required to simulate coupled biogeochemical systems in ice, water, and sediment domains. We developed a 1D sympagic–pelagic–benthic transport model (SPBM) which uses input from physical model simulations to describe hydrodynamics and ice growth and modules from the Framework for Aquatic Biogeochemical Models (FABM) to construct a user-defined biogeochemical model. SPBM coupled with a biogeochemical model simulates the processes of vertical diffusion, sinking/burial, and biogeochemical transformations within and between the three domains. The potential utility of SPBM is demonstrated herein with two test runs using modules from the European regional seas ecosystem model (ERSEM) and the bottom-redox model biogeochemistry (BROM-biogeochemistry). The first run simulates multiple phytoplankton functional groups inhabiting the ice and water domains, while the second simulates detailed redox biogeochemistry in the ice, water, and sediments. SPBM is a flexible tool for integrated simulation of ice, water, and sediment biogeochemistry, and as such may help in producing well-parameterized biogeochemical models for regions with strong sympagic–pelagic–benthic interactions.

1. Introduction

Arctic marine ecosystems have undergone drastic changes and the most important changes are climatically driven [1,2,3,4,5]. The Coupled Model Intercomparison Project and the community climate system model studies have projected atmospheric warming in the Arctic of 1.5–4.5 times the mean global warming, and the Arctic marine environment is expected to be strongly impacted by a loss of ice cover, increasing light exposure, ocean warming, freshening, acidification, and deoxygenation [6]. Modeling simulations are needed for the analysis of present conditions and the projection of long-term impacts on Arctic marine biogeochemistry.

A biogeochemical model suitable for the Arctic should take into account the specific conditions of this region, such as the seasonal to permanent ice cover and the presence of shelf areas. Thus, the model should preferably combine processes occurring in three domains: ice, water column, and sediments. Each of these domains has some specific features and modeling challenges:

Ice. The Arctic ice-algal primary production is a significant part of the total primary production of the Arctic region [7]. Photosynthetic microorganisms extend the production season, provide a winter and early spring food source, and contribute to organic carbon export to depth [8]. A modeling study [9] estimated an average Arctic ice-algal primary production of 21.7 Tg C year ⁻¹, which equates to roughly 5% of total pelagic primary production [10] for this area. Other authors [7] estimated sea ice-algal production accounting for 5–10% of total Arctic and Southern Ocean primary productivity. Another modeling study [11] suggested that under a mild climate change scenario the sea ice community around Greenland may become generally more productive while pelagic phytoplankton productivity may decrease. It is therefore desirable to include the ice domain in biogeochemical modeling studies of the Arctic region. There are three main approaches to implement ice algae behavior according to the place where algae live in the ice column [12,13]: in the bottom layer of an ice column with fixed thickness, in the bottom layer of an ice column with variable thickness, or in any layer of an ice column. Recent research suggests that ice-algal models should resolve the ice vertically to avoid biases that may result from either assuming that ice algae are solely present at the bottom layer or that they have a homogeneous vertical distribution [10].

Water column. In the Arctic, global climate change is causing seawater acidification, accompanied by local changes in productivity and oxygen depletion [14,15]. It follows that the carbon cycle can be an important component of multidecadal-scale biogeochemical models. Oxygen dynamics and redox process parameterization can also be useful in areas affected by oxygen depletion (often in estuaries and fjords). To improve the representation of near-bottom processes the benthic boundary layer (BBL) should be resolved within the water column domain. The BBL is “the part of the marine environment that is directly influenced by the presence of the interface between the bed and its overlying water” [16]. For the Arctic, this layer is especially important since ice melting and permafrost thawing can drive strong fluxes of ungrazed organic material to the BBL [17].

Sediments. Sinking fluxes from the water column can provide sources of new energy for the benthic community. Also, it has been shown [18] that benthic, as well as pelagic, activity can be an important factor for annual pH variability in coastal areas. Sediment layers in models should therefore respond accurately to sinking fluxes and provide accurate remineralization rates. Redox processes occurring in sediments can be highly structured in the vertical direction [19], suggesting a need for explicit vertical resolution in sediment models.

In view of these features and challenges, we aimed to develop a flexible 1D vertical transport model that, when coupled with a biogeochemical model, can provide integrated simulation of biogeochemical processes in ice, water column, and sediment domains, with a vertically-resolved grid for each. The resulting sympagic–pelagic–benthic model (SPBM) uses NetCDF file inputs from hydrodynamic/ice models to describe an “offline” physical environment, and the Framework for Aquatic Biogeochemical Models (FABM) [20] to provide biogeochemical source-minus-sink terms and vertical sinking velocities. FABM is “a Fortran 2003 programming framework for biogeochemical models of marine and freshwater systems. FABM enables complex biogeochemical models to be developed as sets of stand-alone, process-specific modules.” (FABM wiki). The FABM coupling allows the user to construct their own biogeochemical model using existing modules in the FABM library plus any new modules written by the user (SPBM does not itself provide any new biogeochemical modules). The FABM library is rapidly expanding and presently includes modules from some of the most detailed published biogeochemical models, e.g., The European regional seas ecosystem model (ERSEM) [21], the bottom-redox model biogeochemistry (BROM-biogeochemistry) [22], the PCLake aquatic ecosystem model [23], and the model for adaptive ecosystems in coastal seas (MAECS) [24,25]. As with FABM, SPBM transport code is written in FORTRAN.

The paper is structured as follows: Section 1—Introduction (this part); Section 2—Description of the SPBM routines; Section 3—Results from two test simulations to demonstrate SPBM’s capabilities and its relevance to Arctic biogeochemical modeling; Section 4—A discussion of SPBM capabilities and limitations; Section 5—Conclusions.

2. Methods—A 1D Transport Model

SPBM is a 1D advection–diffusion–reaction solver that uses FABM to define an arbitrary biogeochemical model structure and to calculate reaction terms, sinking speeds within the water domain, and various optional biogeochemical diagnostics. FABM distinguishes three types of model variables: state variables, diagnostic variables, and dependencies. State variables are the basic elements for which the rates of changes must be provided (e.g., nitrate, chlorophyll concentrations). Diagnostic variables are calculated within FABM according to the values of the state variables and dependencies at each time step (e.g., pH, nitrification rate). Dependencies are the physical environment variables and interconnections within FABM (e.g., temperature, salinity). SPBM sends dependencies to FABM and updates the state variables over each time step using various advection/diffusion algorithms and the FABM-calculated reaction terms. SPBM outputs all necessary state and diagnostic variables in NetCDF files. Within SPBM, state variables are considered as solute or particulate concentrations.

2.1. Formulation and Numerical Integration

SPBM solves a system of 1-D transport equations in Cartesian coordinates for all three domains (ice, water column, and sediments). The dynamics are

$$\frac{\partial {\mathrm{C}}_{\mathrm{i}}}{\partial \mathrm{t}}=\frac{\partial }{\partial \mathrm{z}}{\mathrm{A}}_{\mathrm{f}}\mathrm{D}\frac{\partial {\mathrm{C}}_{\mathrm{i}}{\mathrm{P}}_{\mathrm{f}}}{\partial \mathrm{z}}-\frac{\partial }{\partial \mathrm{z}}{\mathrm{uC}}_{\mathrm{i}}{+\mathrm{R}}_{\mathrm{i}}$$

(1)

where ${\mathrm{C}}_{\mathrm{i}}$ is the i-th state variable in units provided by the biogeochemical model through FABM, ($\mathrm{mmol} {\mathrm{m}}^{-3} \mathrm{total} \mathrm{volume}$) or ($\mathrm{mg} {\mathrm{m}}^{-3} \mathrm{total} \mathrm{volume}$) (here total volume refers to a representative control volume including both liquid and solid); $\mathrm{t}$ is the time step, ($\mathrm{s}$); $\mathrm{z}$ is the depth, ($\mathrm{m}$); ${\mathrm{A}}_{\mathrm{f}}$ is the porosity-related area restriction factor for fluxes, dimensionless; $\mathrm{D}$ is the total diffusivity, ${(\mathrm{m}}^2 {\mathrm{s}}^{-1})$; ${\mathrm{P}}_{\mathrm{f}}$ is the porosity factor, dimensionless; $\mathrm{u}$ is the sinking velocity (advection/burial in the sediments), ($\mathrm{m} {\mathrm{s}}^{-1}$). ${\mathrm{R}}_{\mathrm{i}}$ is the combined sources minus sinks of the i-th state variable provided by the biogeochemical model through FABM, ($\mathrm{mmol} {\mathrm{m}}^{-3} \mathrm{total} \mathrm{volume} {\mathrm{s}}^{-1})$ or ($\mathrm{mg} {\mathrm{m}}^{-3} \mathrm{total} \mathrm{volume} {\mathrm{s}}^{-1}$). The porosity factor ${\mathrm{P}}_{\mathrm{f}}$ is used to calculate the volume concentration in brine (in the ice column) or in pore water/solid matrix in the sediments. Exchange within the ice and sediment layers occurs through brine channels and through pores or solid matrix, so the area restriction factor ${\mathrm{A}}_{\mathrm{f}}$ is included to limit fluxes within the respective phases (intraphase mixing). The values of ${\mathrm{A}}_{\mathrm{f}}$, ${\mathrm{P}}_{\mathrm{f}}$, $\mathrm{D},$ and $\mathrm{u}$ depend on whether these parameters are calculated in ice, water column, or sediment domains and whether the state variable is solute or particulate.

In the ice domain:

For particulates, it is assumed that the concentration is the same in both the brine channels and ice matrix, hence ${\mathrm{P}}_{\mathrm{f}}=1$. However, vertical fluxes are assumed to be restricted to the brine channels where the particulates are mobilised in suspension, hence ${\mathrm{A}}_{\mathrm{f} }= \mathsf{\phi }\left(\mathrm{z}\right)$. Here, the dimensionless porosity $\mathsf{\phi }\left(\mathrm{z}\right)$ is equal to the relative volume of the brine channels in the ice [26], which can be obtained from an ice thermodynamic model or using empirical relationships (see Appendix A). Solutes are assumed to be excluded from the ice matrix, hence ${\mathrm{P}}_{\mathrm{f}}=\frac{1}{\mathsf{\phi }\left(\mathrm{z}\right)}$, and fluxes are again restricted to the brine channels, hence ${\mathrm{A}}_{\mathrm{f}}=\mathsf{\phi }\left(\mathrm{z}\right)$. The total diffusivity $\mathrm{D}$ in the ice brine channels is a sum of the molecular diffusivity ${\mathrm{D}}_{\mathrm{m}}\left(\mathrm{s}\right)$ (${\mathrm{m}}^2 {\mathrm{s}}^{-1}$) on the ice–water interface (applied only to solutes), the gravity drainage diffusivity ${\mathrm{D}}_{\mathrm{gd}}\left(\mathrm{z}\right)$ (${\mathrm{m}}^2 {\mathrm{s}}^{-1})$ at depths $\mathrm{z}$ within the ice, and the diffusivity caused by convection that occurs in the bottom layer of the growing ice ${\mathrm{D}}_{\mathrm{gi}}\left(\mathrm{s}\right)$ (${\mathrm{m}}^2 {\mathrm{s}}^{-1})$ [26]:

$$\begin{array}{l}{\mathrm{D}=\mathrm{D}}_{\mathrm{m}}{\left(\mathrm{s}\right)+\mathrm{D}}_{\mathrm{gd}}{\left(\mathrm{z}\right)+\mathrm{D}}_{\mathrm{gi}}\left(\mathrm{s}\right)\\ {\mathrm{D}}_{\mathrm{gd}}{\left(\mathrm{z}\right)=\mathrm{F}}_{\mathrm{vb}}{\mathrm{z}}_{\mathrm{b}}\\ {\mathrm{D}}_{\mathrm{gi}}{\left(\mathrm{s}\right)=10}^{-2}{\mathrm{z}}_{\mathrm{s}}(9.667\cdot {10}^{-9}+4.49\cdot {10}^{-6}\mathrm{IceGrowth}-1.39\cdot {10}^{-7}{\mathrm{IceGrowth}}^2)\end{array}$$

where $\mathrm{s}$ means that the value of the parameter is determined only on the interface between the bottom (skeletal) layer of ice and surface water layer; ${\mathrm{F}}_{\mathrm{vb}}$ is a constant mean flux volume rate from the brine channels, ($\mathrm{m} {\mathrm{s}}^{-1}$); ${\mathrm{z}}_{\mathrm{b}}$ is the vertical distance over which the ice column is influenced by brine tube convection (depths where ${\mathsf{\phi }\left(\mathrm{z}\right) > \mathsf{\phi }}_{\min }$), ($\mathrm{m}$); $\mathrm{IceGrowth}$ is the total ice growth rate (${\mathrm{cm} \mathrm{s}}^{-1}$); ${\mathrm{z}}_{\mathrm{s}}$ is the thickness of the ice layer, ($\mathrm{m}$). ${\mathrm{D}}_{\mathrm{gi}}\left(\mathrm{s}\right)$ is not equal to zero only during the period of ice build-up and only on the interface between water and ice. Alternatively, the total diffusivity $\mathrm{D}$ can be read from an input file generated by e.g., an ice thermodynamic model.

The sinking velocity $\mathrm{u}$ is non-zero only for particulate variables in the layers where ${\mathsf{\phi }\left(\mathrm{z}\right) > \mathsf{\phi }}_{\min }$ (if $\mathsf{\phi }\left(\mathrm{z}\right) \le { \mathsf{\phi }}_{\min }$ sea ice brine pockets are not interconnected) and is generally determined at each time step by the biogeochemical model through FABM. For all diatoms living in the ice column, to represent their ability to maintain their vertical position relative to the skeletal layer [26], $\mathrm{u}$ is set to a constant but possibly layer-dependent value within the ice column and zero on the ice–water interface between ice and water domains, while the total diffusivity $\mathrm{D}$ is set to zero.

In the water column domain:

Here ${\mathrm{P}}_{\mathrm{f}}=1$ and ${\mathrm{A}}_{\mathrm{f}}=1$ at all depths for both solutes and particulates, since there is only one phase to consider.

The total diffusivity $\mathrm{D}$ is composed of the molecular diffusivity ${\mathrm{D}}_0$ (${\mathrm{m}}^2 {\mathrm{s}}^{-1})$ (applied only to solutes) and the turbulence diffusivity ${\mathrm{D}}_{\mathrm{t}}\left(\mathrm{z}\right)$ (${\mathrm{m}}^2 {\mathrm{s}}^{-1})$:

$${\mathrm{D}=\mathrm{D}}_0{+\mathrm{D}}_{\mathrm{t}}\left(\mathrm{z}\right)$$

where ${\mathrm{D}}_{\mathrm{t}}\left(\mathrm{z}\right)$ is taken from the hydrophysical model as input data. The water column domain contains the structure that could be called the BBL. It is located in the lower part of the water column. Turbulent diffusivity for each layer ${\mathrm{z}}_{\mathrm{i}}$ within the BBL is linearly decreasing from the deepest non-zero value of the diffusivity ${\mathrm{D}}_{\mathrm{t}}{(\mathrm{z}}_{\mathrm{d}})$ as follows:

$${\mathrm{D}}_{\mathrm{t}}{(\mathrm{z}}_{\mathrm{i}})=\frac{{\mathrm{D}}_{\mathrm{t}}{(\mathrm{z}}_{\mathrm{d}})}{{\mathrm{z}}_{\mathrm{d}}{-\mathrm{z}}_0}{(\mathrm{z}}_{\mathrm{i}}{-\mathrm{z}}_0)$$

where ${\mathrm{z}}_{\mathrm{d}}$ ($\mathrm{m}$) is the deepest depth with non-zero value of ${\mathrm{D}}_{\mathrm{t}}\left(\mathrm{z}\right)$ and ${\mathrm{z}}_0$ ($\mathrm{m}$) is the bottom depth.

The sinking velocity $\mathrm{u}$ is taken from the biogeochemical model through FABM for all particulates and is zero for all solutes.

In the sediment domain:

Within the sediments, particulate variables are confined to the sediment matrix and solutes are confined to the pore water. So, for solid particulates the porosity factor ${\mathrm{P}}_{\mathrm{f}}=\frac{1}{1-\mathsf{\phi }\left(\mathrm{z}\right)}$ and the area restriction factor ${\mathrm{A}}_{\mathrm{f}}=1-\mathsf{\phi }\left(\mathrm{z}\right)$ at depths $\mathrm{z}$. For solutes ${\mathrm{P}}_{\mathrm{f}}=\frac{1}{\mathsf{\phi }\left(\mathrm{z}\right)}$ and ${\mathrm{A}}_{\mathrm{f}}=\mathsf{\phi }\left(\mathrm{z}\right)$. There is no adsorption in the present version.

A time-independent porosity $\mathsf{\phi }\left(\mathrm{z}\right)$ at depths $\mathrm{z}$ through the entire sediment domain is described following [27]:

$${\mathsf{\phi }\left(\mathrm{z}\right)=\mathsf{\phi }(\mathrm{z}}_{\infty })+(\mathsf{\phi }({\mathrm{z}}_0)-\mathsf{\phi }({\mathrm{z}}_{\infty }{\left)\right)\mathrm{e}}^{-\frac{{(\mathrm{z}-\mathrm{z}}_{\mathrm{swi}})}{{\mathrm{k}}_{\mathsf{\phi }}}}$$

where ${\mathsf{\phi }(\mathrm{z}}_{\infty })$ is the deep porosity, dimensionless; ${\mathsf{\phi }(\mathrm{z}}_0)$ is the porosity at the sediment–water interface (SWI), dimensionless; ${\mathrm{z}}_{\mathrm{swi}}$ is the depth of the SWI, ($\mathrm{m}$); ${\mathrm{k}}_{\mathsf{\phi }}$ is the coefficient for exponential porosity change, ($\mathrm{m}$).

The total diffusivity $\mathrm{D}$ is a sum of the molecular diffusivity ${\mathrm{D}}_{\mathrm{m}}\left(\mathrm{z}\right)$ (${\mathrm{m}}^2 {\mathrm{s}}^{-1})$ (applied only to solutes) and the bioturbation diffusivity ${\mathrm{D}}_{\mathrm{b}}\left(\mathrm{z}\right)$ (${\mathrm{m}}^2 {\mathrm{s}}^{-1})$ [28]:

$$\begin{array}{l}{\mathrm{D}=\mathrm{D}}_{\mathrm{m}}{\left(\mathrm{z}\right)+\mathrm{D}}_{\mathrm{b}}\left(\mathrm{z}\right)\\ {\mathrm{D}}_{\mathrm{m}}\left(\mathrm{z}\right){=\mathrm{D}}_0\frac{1}{1-2\ln \mathsf{\phi }\left(\mathrm{z}\right)}{\mathsf{\mu }}_{\mathrm{d}}\\ {\mathrm{D}}_{\mathrm{b}}{\left(\mathrm{z}\right)=\mathrm{D}}_{\mathrm{bo}}\left(\mathrm{z}\right)\frac{{\mathrm{O}}_2}{{\mathrm{O}}_{2 }{+ \mathrm{K}}_{{\mathrm{O}}_2}}\end{array}$$

where ${\mathrm{D}}_0$ is the infinite-dilution molecular diffusivity, (${\mathrm{m}}^2 {\mathrm{s}}^{-1});$ ${\mathsf{\mu }}_{\mathrm{d}}$ is the relative dynamic viscosity, dimensionless; ${\mathrm{O}}_2$ is the oxygen concentration in the bottom layer of the water column, ($\mathrm{mmol} {\mathrm{m}}^{-3}$); ${\mathrm{K}}_{{\mathrm{O}}_2}$ is the half-saturation constant, ($\mathrm{mmol} {\mathrm{m}}^{-3}$). The oxygen-saturated bioturbation diffusivity [22] ${\mathrm{D}}_{\mathrm{bo}}\left(\mathrm{z}\right)$ (${\mathrm{m}}^2 {\mathrm{s}}^{-1})$ depends on the distance ${\mathrm{z}}_{\mathrm{db}}\left(\mathrm{z}\right)$ ($\mathrm{m}$) between the interface depth $\mathrm{z}$ and the depth with a constant bioturbation activity as follows:

$$\begin{array}{cc}\hfill & {\mathrm{z}}_{\mathrm{db}}{\left(\mathrm{z}\right)=\mathrm{z}-(\mathrm{z}}_{\mathrm{swi}}{+\mathrm{z}}_{\mathrm{cb}})\hfill \\ \hfill {\mathrm{if} \mathrm{z} < \mathrm{z}}_{\mathrm{swi}}{+\mathrm{z}}_{\mathrm{cb}}: & {\mathrm{D}}_{\mathrm{bo}}{\left(\mathrm{z}\right)=\mathrm{D}}_{\mathrm{bm}}\hfill \\ \hfill {\mathrm{if} \mathrm{z} > \mathrm{z}}_{\mathrm{swi}}{+\mathrm{z}}_{\mathrm{cb}}: & {\mathrm{D}}_{\mathrm{bo}}{\left(\mathrm{z}\right)=\mathrm{D}}_{\mathrm{bm}}{\mathrm{e}}^{-\frac{{\mathrm{z}}_{\mathrm{db}}\left(\mathrm{z}\right)}{{\mathrm{F}}_{\mathrm{d}}}}\hfill \end{array}$$

where ${\mathrm{z}}_{\mathrm{swi}}$ is the depth at the SWI, ($\mathrm{m}$); ${\mathrm{z}}_{\mathrm{cb}}$ is the constant bioturbation activity layer thickness, ($\mathrm{m}$); ${\mathrm{D}}_{\mathrm{bm}}$ is the maximum bioturbation diffusivity, (${\mathrm{m}}^2 {\mathrm{s}}^{-1})$; ${\mathrm{F}}_{\mathrm{d}}$ is the bioturbation decay scale, ($\mathrm{m}$).

On the SWI it is assumed that the bioturbation diffusivity mixes concentrations in units ($\mathrm{mmol} {\mathrm{m}}^{-3} \mathrm{total} \mathrm{volume}$) instead of ($\mathrm{mmol} {\mathrm{m}}^{-3} \mathrm{solids}/\mathrm{solutes}$) (interphase mixing). Therefore, special values of ${\mathrm{P}}_{\mathrm{f}}$ are needed for the layers immediately above and below the SWI (see Appendix B):

$$\begin{array}{cc}\hfill \mathrm{for} \mathrm{solutes}: & {\mathrm{P}}_{\mathrm{f}}{(\mathrm{z}}_{\mathrm{a},\mathrm{b}})=\frac{\frac{{\mathsf{\phi }}_{\mathrm{swi}}}{{\mathsf{\phi }}_{\mathrm{a},\mathrm{b}}}{\mathrm{D}}_{\mathrm{m}}{(\mathrm{z}}_{\mathrm{swi}}{)+\mathrm{D}}_{\mathrm{b}}{(\mathrm{z}}_{\mathrm{swi}})}{{\mathsf{\phi }}_{\mathrm{swi}}{(\mathrm{D}}_{\mathrm{m}}{(\mathrm{z}}_{\mathrm{swi}}{)+\mathrm{D}}_{\mathrm{b}}{(\mathrm{z}}_{\mathrm{swi}}\left)\right)}\hfill \\ \hfill \mathrm{for} \mathrm{solids}: & {\mathrm{P}}_{\mathrm{f}}{(\mathrm{z}}_{\mathrm{a},\mathrm{b}})=\frac{1}{{1-\mathsf{\phi }}_{\mathrm{swi}}}\hfill \end{array}$$

where the subscripts $\mathrm{a}$, $\mathrm{b}$ and $\mathrm{swi}$ determine a location of the corresponding variables: $\mathrm{a}$ means the layer above, $\mathrm{b}$ the layer below, $\mathrm{swi}$ on the SWI.

The advection/burial velocities $\mathrm{u}\left(\mathrm{z}\right)$ are described following [22]:

$$\mathrm{for} \mathrm{solutes}: \mathrm{u}\left(\mathrm{z}\right)=\frac{{\mathsf{\phi }(\mathrm{z}}_{\infty })}{\mathsf{\phi }\left(\mathrm{z}\right)}{\mathrm{u}}_{\mathrm{b}}+\frac{1}{\mathsf{\phi }\left(\mathrm{z}\right)}{\mathrm{D}}_{\mathrm{b}}^{\mathrm{inter}}\frac{\partial \mathsf{\phi }\left(\mathrm{z}\right)}{\partial \mathrm{z}}$$

(2)

$$\mathrm{for} \mathrm{particulates}: \mathrm{u}\left(\mathrm{z}\right)=\frac{{1-\mathsf{\phi }(\mathrm{z}}_{\infty })}{1-\mathsf{\phi }\left(\mathrm{z}\right)}{\mathrm{u}}_{\mathrm{b}}-\frac{1}{1-\mathsf{\phi }\left(\mathrm{z}\right)}{\mathrm{D}}_{\mathrm{b}}^{\mathrm{inter}}\frac{\partial \mathsf{\phi }\left(\mathrm{z}\right)}{\partial \mathrm{z}}$$

(3)

where ${\mathrm{u}}_{\mathrm{b}}$ is the deep burial velocity, ($\mathrm{m} {\mathrm{s}}^{-1}$); ${\mathrm{D}}_{\mathrm{b}}^{\mathrm{inter}}$ is the interphase component of the total bioturbation diffusivity ${\mathrm{D}}_{\mathrm{b}}{=\mathrm{D}}_{\mathrm{b}}^{\mathrm{intra}}{+\mathrm{D}}_{\mathrm{b}}^{\mathrm{inter}}$, and ${\mathrm{D}}_{\mathrm{b}}^{\mathrm{inter}}$ is nonzero only on the SWI where ${\mathrm{D}}_{\mathrm{b}}^{\mathrm{inter}}{=\mathrm{D}}_{\mathrm{b}}$ (${\mathrm{m}}^2 {\mathrm{s}}^{-1}$). Note that although a non-zero ${\mathrm{D}}_{\mathrm{b}}^{\mathrm{inter}}$ beyond the SWI would alter the computed advection/burial velocities $\mathrm{u}\left(\mathrm{z}\right)$ via Equations (2) and (3), the net transport of biogeochemical tracers would not be affected because corresponding interphase components $\frac{1}{\mathsf{\phi }}{\mathrm{C}}_{\mathrm{i}}{\mathrm{D}}_{\mathrm{b}}^{\mathrm{inter}}\frac{\partial \mathsf{\phi }}{\partial \mathrm{z}}$ and $\frac{1}{1-\mathsf{\phi }}{\mathrm{C}}_{\mathrm{i}}{\mathrm{D}}_{\mathrm{b}}^{\mathrm{inter}}\frac{\partial (1-\mathsf{\phi })}{\partial \mathrm{z}}$ would need to be added to the diffusive fluxes in Equation (1), and these will exactly cancel the contributions of ${\mathrm{D}}_{\mathrm{b}}^{\mathrm{inter}}$ to advection/burial. In other words, when the porosity profile is specified and used to compute advection/burial velocities under steady state compaction, the tracer advection/diffusion depends only on the total bioturbation diffusivity, and the intraphase form assumed by SPBM for diffusion inside the sediments is correct irrespective of the relative contribution of inter- vs. intraphase mixing, see [29]. However, there can be no intraphase component of a Fickian particulate diffusion across the SWI because by definition there is no solid matrix above the SWI [22].

Equation 1 is integrated numerically over a single combined (ice, water column, sediments) grid, using a constant model time step. The coupling method follows an operator splitting approach [30]: concentrations are successively updated by contributions over one time step of diffusion, reaction, and sinking/advection/burial, in that order. Diffusive updates are calculated by a semi-implicit central-space algorithm adapted from a routine in BROM-transport [22] which in turn was adapted from the general ocean turbulence model (GOTM) [31]. Sinking/advection/burial updates are calculated using a first-order upwind differencing scheme. Reaction updates are calculated from forward Euler time steps.

2.2. The Grid

SPBM uses a fixed grid structure for the water column and sediments, and a time-dependent grid for the ice column. The number of grid points inside the ice column can vary with time but the spacing is fixed (see Figure 1). Water column layer depths ($\mathrm{m}$) are taken as input from a hydrophysical model (distances between layers can be unequal) and extra layers are incorporated in the lower part in order to fully resolve the BBL. Total ice thickness ($\mathrm{m}$) for every day of simulation is also taken as input from a hydrophysical model, and the ice column is constructed using a fixed layer thickness ($\mathrm{m}$) as an input parameter. Therefore, the ice column is discretized into layers of strictly constant thickness ${\mathrm{z}}_{\mathrm{s}}$, and when the ice column grows or melts its total thickness can change only by multiples of ${\mathrm{z}}_{\mathrm{s}}$. This simplification facilitates recalculation of the variable concentrations during melting and freezing.

2.3. Irradiance Formulation

FABM biogeochemical models generally need to know the photosynthetically active radiation (PAR), e.g., ($\mathrm{mol} \mathrm{photons} {\mathrm{m}}^{-2} {\mathrm{day}}^{-1}$), in each layer of the model grid. Some FABM models compute water column PAR given only surface PAR, but they do not assume the existence of the ice column and consider all grid points to be located within the water column. SPBM therefore provides the following simple approach to calculate PAR in both ice and water column domains.

PAR on the surface of water or ice ${\mathrm{P}}_{\mathrm{s}}$ can be calculated from the surface shortwave radiative flux ${\mathrm{F}}_{\mathrm{surf}}$ ($\mathrm{W} {\mathrm{m}}^{-2}$), depending on the solar declination ${\mathrm{k}}_{\mathrm{decl}}$ ($\mathrm{degrees}$):

$$\begin{array}{l}{\mathrm{k}}_{\mathrm{decl}}=23.5\cdot \sin \frac{2\mathsf{\pi }\mathrm{JulianDay}-81}{365}\\ {\mathrm{F}}_{\mathrm{surf}}{=\mathrm{I}}_{\mathrm{m}}\cos \frac{{\mathsf{\pi }(\mathrm{latitude}-\mathrm{k}}_{\mathrm{decl}})}{180}\\ {\mathrm{P}}_{\mathrm{s}}{=\mathrm{k}}_{\mathrm{f}}{\mathrm{F}}_{\mathrm{surf}}\end{array}$$

where ${\mathrm{I}}_{\mathrm{m}}$ is the theoretical maximum of 24-h average surface downwelling shortwave irradiance in air, ($\mathrm{W} {\mathrm{m}}^{-2}$); ${\mathrm{k}}_{\mathrm{f}}$ is the factor to convert downwelling shortwave irradiance in air to scalar PAR in water, ($\mathrm{mol} \mathrm{photons} {\mathrm{day}}^{-1} {\mathrm{W}}^{-1}$) [32]. Alternatively, ${\mathrm{P}}_{\mathrm{s}}$ (or ${\mathrm{F}}_{\mathrm{surf}}$) can be read from an input file.

In the presence of ice, PAR after considering albedo influence ${\mathrm{P}}_{\mathrm{a}}$ becomes [33]:

$$\begin{array}{cc}\hfill \mathrm{if} \mathrm{snow} \mathrm{depth} \le 5 \mathrm{mm}: & {\mathrm{P}}_{\mathrm{a}}{=\mathrm{P}}_{\mathrm{s}}{\mathrm{k}}_{\mathrm{scatter}}{(1-\mathrm{A}}_{\mathrm{ice}})\hfill \\ \hfill \mathrm{if} \mathrm{snow} \mathrm{depth} > 5 \mathrm{mm}: & {\mathrm{P}}_{\mathrm{a}}{=\mathrm{P}}_{\mathrm{s}}{\mathrm{k}}_{\mathrm{scatter}}{(1-\mathrm{A}}_{\mathrm{snow}}{)\mathrm{e}}^{{-\mathrm{k}}_{\mathrm{snow}}{\mathrm{z}}_{\mathrm{snow}}}\hfill \end{array}$$

where ${\mathrm{k}}_{\mathrm{scatter}}$ is the fraction of radiation transmitted through the highly scattering surface of the ice, dimensionless; ${\mathrm{A}}_{\mathrm{ice}}$ is the ice albedo for visible light, dimensionless; ${\mathrm{A}}_{\mathrm{snow}}$ is the snow albedo for visible light, dimensionless; ${\mathrm{k}}_{\mathrm{snow}}$ is the snow light extinction coefficient, (${\mathrm{m}}^{-1}$); ${\mathrm{z}}_{\mathrm{snow}}$ is the snow depth, ($\mathrm{m}$).

PAR at any depth in the ice ${\mathrm{P}(\mathrm{z}}_{\mathrm{ice}})$ is given by:

$${\mathrm{P}(\mathrm{z}}_{\mathrm{ice}}{)=\mathrm{P}}_{\mathrm{a}}{\mathrm{e}}^{{-\mathrm{k}}_{\mathrm{ice}}{\mathrm{z}}_{\mathrm{ice}}}$$

where ${\mathrm{k}}_{\mathrm{ice}}$ is the ice light extinction coefficient, (${\mathrm{m}}^{-1}$); ${\mathrm{z}}_{\mathrm{ice}}$ is the ice depth, ($\mathrm{m}$).

PAR in the water column ${\mathrm{P}(\mathrm{z}}_{\mathrm{water}})$ is calculated according to the Beer–Lambert formulation:

$$\begin{array}{cc}\hfill \mathrm{if} \mathrm{there} \mathrm{is} \mathrm{ice}: & {\mathrm{P}(\mathrm{z}}_{\mathrm{water}}{)=\mathrm{P}(\mathrm{z}}_{\mathrm{IceBottom}}{)\mathrm{e}}^{{{\displaystyle \int }}_0^{\mathrm{z}}{\mathrm{K}}_{\mathrm{water}}{(\mathrm{z}}_{\mathrm{water}}{)\mathrm{dz}}_{\mathrm{water}}}\hfill \\ \hfill \mathrm{if} \mathrm{there} \mathrm{is} \mathrm{no} \mathrm{ice}: & {\mathrm{P}(\mathrm{z}}_{\mathrm{water}}{)=\mathrm{P}}_{\mathrm{s}}{\mathrm{e}}^{{{\displaystyle \int }}_0^{\mathrm{z}}{\mathrm{K}}_{\mathrm{water}}{(\mathrm{z}}_{\mathrm{water}}{)\mathrm{dz}}_{\mathrm{water}}}\hfill \end{array}$$

where ${\mathrm{P}(\mathrm{z}}_{\mathrm{IceBottom}})$ is the PAR at the ice bottom layer; ${\mathrm{K}}_{\mathrm{water}}$ is the vertically varying water light extinction coefficient provided by the FABM models, describing attenuation due to living and non-living optically-active substances, (${\mathrm{m}}^{-1}$); ${\mathrm{z}}_{\mathrm{water}}$ is the water layer depth, ($\mathrm{m}$).

In the sediment domain, PAR equals zero in all layers.

2.4. Initial and Boundary Conditions; the Forcing Data

Initial conditions for all state variable concentrations are provided through FABM using its YAML type configuration file [20]. By default, zero gradient boundary conditions are used at upper and lower boundaries for all state variables except O₂ and CO₂. Diffusive fluxes of O₂ and CO₂ are provided by the biogeochemical model through FABM at the surface boundary (only for ice-free periods) and are set to zero at the lower boundary. It is possible to change both boundary conditions according to the user’s needs.

SPBM requires time-dependent input forcing for the entire period of simulation for the water column (turbulent diffusivity (${\mathrm{m}}^2 {\mathrm{s}}^{-1}$) on layer interfaces; temperature ($\mathrm{C}$) and salinity ($\mathrm{psu}$) on layer centers) and for the ice column (total thickness ($\mathrm{m}$), snow thickness ($\mathrm{m}$), and surface temperature ($\mathrm{C}$)). Additional forcings may be required depending on the FABM biogeochemical models. Downwelling shortwave radiation and PAR can be read from an input file instead of using the formulae provided in Section 2.3. Other optional input forcing includes brine volumes and diffusion coefficients in the ice, input fluxes at the water surface, and horizontal mixing fluxes at any depth. Input fluxes are based on concentrations $\mathrm{C}$ which can be provided in three ways: read from text or NetCDF file; set as fixed sinusoidal variation in time defined by a maximum value $\mathrm{M}$ and $\mathrm{Phase}$ parameters (${\mathrm{C}=2}^{-1}{\mathrm{M}+4}^{-1}{\mathrm{M}(1+\sin (365}^{-1}\cdot 2\mathsf{\pi }(\mathrm{JulianDay}-\mathrm{Phase})\left)\right)$); set as fixed constant value. $\mathrm{M}$ and the boundary concentrations $\mathrm{C}$ should be in units corresponding to the state variables of the appropriate FABM model, $\mathrm{Phase}$ is in ($\mathrm{days})$. SPBM uses input data files in NetCDF and text formats.

3. Results—Test Runs

The purpose of the test runs is only to demonstrate the flexibility of SPBM and its relevance to Arctic marine modeling. Rigorous, site-specific adaptation and skill assessment of particular SPBM- ‘biogeochemical model’ configurations are not within our present remit. SPBM itself does not require validation since it is based on a standard advection–diffusion solver (with a possibility to solve within the ice, water, and sediment domains simultaneously). However, mass conservation of state variables has been checked.

The test runs use forcing data from a regional oceanic modeling system (ROMS) simulation to provide a hydrodynamic scenario for the Laptev Sea. The whole period of this simulation spans a period from 1980 till 2010. A time-dependent total ice thickness and a time-independent water column structure were derived from this simulation, while the BBL was inserted with the following parameters: width = 50 $\mathrm{cm}$, resolution = 10 $\mathrm{cm}$. The grid in the sediment domain was continued for another 10 $\mathrm{cm}$ with resolution 2 $\mathrm{cm}$.

For the test simulations, we use FABM to combine components from two published biogeochemical models. Here we will explain only the most basic aspects of these models; the reader can find detailed descriptions in the provided references. We remind that SPBM calculates only the transport terms in Equation (1), while the FABM biogeochemical modules provide the combined sources-minus-sinks terms ${\mathrm{R}}_{\mathrm{i}}$ and the sinking velocities $\mathrm{u}$ in the water column. FABM model formulations and parameter values were derived from existing parameterizations with some limited adaptation to the Arctic scenario.

The first biogeochemical model is the European regional seas ecosystem model, ERSEM [21]. Originally a coastal ecosystem model for the North Sea, ERSEM has evolved into a generic tool for ecosystem simulations from shelf seas to the global ocean. Model dynamics within each functional group describe processes occurring inside a ’standard organism’ [34,35]. ERSEM accounts for flexible elemental stoichiometry in planktonic processes by allowing decoupled fluxes of carbon, nitrogen, phosphorus, silicate, and chlorophyll a. This requires multiple state variables to describe each functional group biomass (e.g., diatom carbon, diatom nitrogen, etc.) and results in a relatively complex model.

The second biogeochemical model is the bottom-redox model (BROM) biogeochemistry module [22]. This model represents key biogeochemical processes in the water and upper sediments, with a focus on oxygen dynamics and redox biogeochemistry. Compared to ERSEM, it simulates the coupled cycles of more elements (N, P, Si, C, O, S, Mn, and Fe), resolves more structure in the bacterial community (four functional groups vs. one in ERSEM), and calculates carbonate chemistry in more detail; however it assumes fixed stoichiometry for all forms of organic matter (nitrogen currency), resolves only one functional group each for phytoplankton and zooplankton, and does not resolve dissolved organic matter into different lability classes (in the present version).

The general coupling scheme is illustrated in Figure 2. A quasi-stationary solution is derived from two spin-ups, repeating the first day 100 times and then repeating the first year 10 times. Along with SPBM requirements there is additional forcing required by ERSEM and BROM: wind speed ($\mathrm{m} {\mathrm{s}}^{-1}$) and concentration of ${\mathrm{CO}}_2$ in air ($\mathrm{ppm}$). Surface shortwave radiation at sea-surface level is provided by ROMS and is read from an input file.

We present two test cases with the same hydrophysical forcing but different FABM model configurations. The first demonstrates the simulation of multiple primary producer functional groups and shows their variability in contrasting conditions. The second test case demonstrates changing redox conditions in the three domains in response to a constant input of organic matter (OM). The main joint parameter values and forcing properties are provided in Appendix C (common parameters in

Table A1. Common parameters.

Parameter	Description	Value	Unit	Reference
$\mathrm{t}$	Time step	300	$\mathrm{s}$
${\mathrm{D}}_0$	Infinite-dilution molecular diffusivity	10−9	${\mathrm{m}}^2{\mathrm{s}}^{-1}$	[28]
${\mathrm{V}}_{\mathrm{ws}}$	Wind speed	5	${\mathrm{ms}}^{-1}$
${\mathrm{CO}}_2$(g)	Concentration of ${\mathrm{CO}}_2$ in air	380	ppm

, ice parameters in

Table A2. Ice parameters.

Parameter	Description	Value	Unit	Reference
${\mathrm{D}}_{\mathrm{m}}\left(\mathrm{s}\right)$	Diffusivity on sea-water interface	10−5	${\mathrm{m}}^2{\mathrm{s}}^{-1}$	[51]
${\mathrm{u}}_{\mathrm{d}}$	Diatoms vertical movement velocity	3	${\mathrm{cm} \mathrm{d}}^{-1}$
${\mathrm{F}}_{\mathrm{vb}}$	Flux rate from the brine channels	10−8	${\mathrm{ms}}^{-1}$	[26]
${\mathrm{z}}_{\mathrm{s}}$	Thickness of the ice layer	0.06	m
${\mathsf{\phi }}_{\min }$	Minimum porosity to enable brine convection	0.072 (0.12)	-	[26]

, sediments parameters in

Table A3. Sediments parameters.

Parameter	Description	Value	Unit	Reference
${\mathsf{\phi }(\mathrm{z}}_{\infty })$	Porosity at the infinite sediments depth	0.8	-	[27]
${\mathsf{\phi }(\mathrm{z}}_0)$	Porosity at the sediments–water interface	0.95	-	[27]
${\mathrm{k}}_{\mathsf{\phi }}$	Coefficient for exponential porosity change	0.04	m	[27]
${\mathsf{\mu }}_{\mathrm{d}}$	Relative dynamic viscosity	0.94	-	[28]
${\mathrm{K}}_{\mathrm{O}}{}_2$	Oxygen half-saturation constant	5	${\mathrm{mmol} \mathrm{m}}^{-3}$	[22]
${\mathrm{z}}_{\mathrm{cb}}$	Constant bioturbation activity layer width	0.02	m	[28]
${\mathrm{D}}_{\mathrm{bm}}$	Maximum bioturbation diffusivity	10−11	${\mathrm{m}}^2{\mathrm{s}}^{-1}$	[28]
${\mathrm{F}}_{\mathrm{d}}$	Bioturbation decay scale	0.01	m	[28]
${\mathrm{u}}_{\mathrm{b}}$	Deep burial velocity	10−10	${\mathrm{ms}}^{-1}$	[28]

, irradiance parameters are provided in

Table A4. Irradiance parameters.

Parameter	Description	Value	Unit	Reference
${\mathrm{k}}_{\mathrm{f}}$	Factor converting irradiance to PAR	0.5	${\mathrm{mol} \mathrm{photons} \mathrm{d}}^{-1}{\mathrm{W}}^{-1}$	[32]
${\mathrm{k}}_{\mathrm{scatter}}$	Fraction of transmitted radiation	0.97	-	[33]
${\mathrm{A}}_{\mathrm{ice}}$	Ice albedo	0.744	-	[33]
${\mathrm{A}}_{\mathrm{snow}}$	Snow albedo	0.9	-	[52]
${\mathrm{k}}_{\mathrm{snow}}$	Snow light extinction coefficient	4.3	${\mathrm{m}}^{-1}$	[52]
${\mathrm{k}}_{\mathrm{ice}}$	Ice light extinction coefficient	0.93	${\mathrm{m}}^{-1}$	[33]

, and forcing properties in

Table A5. Forcing properties were chosen according to the values provided by [53].

Parameter	Position	Type	Value
DIC	Water surface layer	Constant	${1930 \mathrm{mmol} \mathrm{m}}^{-3} \mathrm{total} \mathrm{volume}$
TA	Water surface layer	Constant	${2000 \mathrm{mmol} \mathrm{m}}^{-3} \mathrm{total} \mathrm{volume}$
PO4	Water surface layer	Sinusoidal	M = 0.4, Phase = 130
NO3	Water surface layer	Sinusoidal	M = 3, Phase = 130
Si	Water surface layer	Sinusoidal	M = 10, Phase = 130
O2	Water surface layer	Sinusoidal	M = 330, Phase = 130
DON	Water surface layer	Sinusoidal	M = 12.5, Phase = 130
DIC	Water bottom layer	Constant	${2280 \mathrm{mmol} \mathrm{m}}^{-3} \mathrm{total} \mathrm{volume}$
TA	Water bottom layer	Constant	${2350 \mathrm{mmol} \mathrm{m}}^{-3} \mathrm{total} \mathrm{volume}$

).

3.1. Test Case 1

According to the [36], there are few tools yet developed to simulate different groups of ice algae. Therefore, we constructed a test case to simulate different primary producers in different ice conditions. We used the ERSEM primary producer functional group formulation with parameterizations from [21] to simulate diatoms, nanophytoplankton, picophytoplankton, and microphytoplankton. One more algal group called ice diatoms was added, based on an original diatom primary producer parametrization which we adapted to improve growth in low irradiance conditions (see Appendix D). Thus, in total we used five groups. All groups were given the same initial conditions (prior to spin-up) in both ice, water, and sediment domains; hence differences in the steady state abundances were determined by the environment and the growth parameters/sinking velocities of the different functional groups. To calculate pH a corresponding module from BROM was connected. All configuration files are available in the supplementary material.

Figure 3 shows a SPBM–ERSEM–BROM coupled system output (for chlorophyll a, pH, oxygen, and nutrients) in the ice and upper water column layers during the period with maximum ice algae chlorophyll a concentration. This maximum is a result of thin ice (<50 $\mathrm{cm}$) during at least one month and favorable irradiance conditions. Chlorophyll a, oxygen, and nutrients are all state variables and therefore output as concentrations per unit total volume; pH is a diagnostic variable and is output as the negative logarithm of hydrogen ion concentrations of brine or seawater. The modeled values were compared with observed concentrations of biogeochemical tracers in sea ice during spring in the Arctic provided by [8] (p. 210). Most model ranges fall inside the observational ranges (see

Table 1. Comparing SPBM–ERSEM–BROM output during ice melting period with concentrations of biogeochemical tracers in sea ice and surface seawater observed for typical spring conditions in the Arctic according to [8] (p. 210). For the observed values in water no value range is provided by the authors.

Parameter	Si µM	PO4 µM	NO3 µM	$\mathbf{Chl} \mathbf{a} \mathbf{mg} {\mathbf{m}}^{-3}$	O2 µM
Observed in ice	-	0–0.7	0–1	1–100	50–250
Modeled in ice	0.4–1.6	0.01–0.07	0–0.5	1–120	50–80
Observed in water	-	≈1.2	≈7	≈1	≈380
Modeled in water	8.5–9.5	0.34–0.38	2.5–2.9	0.1–0.7	300–320

). Also, Figure 3 shows that the modeled vertical distributions in sea ice reproduce some commonly observed features [8]: during ice melt chlorophyll a concentrations are highest in the bottom layer; during freezing the pH increases in the upper ice layers, reaching values higher than 10 during winter (see supplementary material) in accordance with observations [37]; nutrients have maximum values on the lowest ice layer, phosphates and nitrates are almost depleted in all ice layers except the bottom one; oxygen profiles in ice have a complicated structure with two maxima, in the middle and the top ice layers. Modeled oxygen are somewhat low compared to the observational range in [8,38] (see Table 1). Observed values include oxygen from gas bubbles that are incorporated into the ice and which are not available for biogeochemical reactions in brine channels, while the modeled ice oxygen is only the oxygen dissolved in ice brines.

Figure 4 shows ice/water biomass concentrations of the five primary producer functional groups during the year of maximum primary production and the preceding year. The year of maximum primary production (1983) shows the springtime migration of the modeled ice diatoms and diatoms from the top to the bottom ice layer. The nanophytoplankton and picophytoplankton remain in the uppermost ice layers, while the microphytoplankton migrate downwards during the final stage of the melting season. All modeled primary producers start growing from the surface ice layers, where they were frozen during previous year ice buildup. In years with high irradiance and low summertime ice cover the concentrations of nanophytoplankton and picophytoplankton in the water column are much higher (see supplementary material) and can exceed the diatom concentration.

3.2. Test Case 2

Test case 2 demonstrates the potential utility of SPBM for studying anaerobic processes in the sea ice, water, and sediment columns. Here we use BROM to simulate biogeochemical processes occurring in low oxygen environments. The BROM configuration file is provided following the link in the code availability section. To facilitate suboxic conditions, we forced a supplementary flux of particulate OM to the upper level of the water column (see along with other forcing properties in Table A5). For presentation, we chose the same years as for test case 1. Figure 5, Figure 6 and Figure 7 show some of the available BROM variables in the ice, water, and sediment domains. Concentrations in ice are strongly driven by surface water concentrations, which in turn are strongly influenced by hydrophysical conditions. For variables not involved in redox reactions, vertical distributions in the ice column mainly reflect temporal distributions in the upper water layer during freezing.

The simulation reproduces some general features of sea ice–water–sediment seasonal biogeochemistry connected with the seasonal production and decomposition of OM. Phytoplankton (one state variable in BROM) starts to bloom in the ice and below the ice during ice melting (see supplementary material). Phytoplankton blooms in the upper water column lead to a seasonal increase in dissolved organic matter (DON, Figure 5A, in nitrogen units), dead particulate organic matter (PON, Figure 5B, in nitrogen units), and dissolved oxygen (Figure 5C), followed by oxygen consumption in the lower water column/BBL and the generation of reduced forms of nitrogen (${\mathrm{NO}}_2^{-}$, ${\mathrm{NH}}_4^{+}$) at depth (Figure 6C). DON is incorporated into the ice during freezing, but since BROM only models the labile fraction at present, the model DON is rapidly oxidized and contributes little to OM in the upper ice layers (Figure 5A). Here, the main reduction agent is PON, which occurs in ice due to decomposition of zooplankton (Figure 5B). The modeled oxygen in ice is almost depleted (Figure 5C). In the upper sediments, oxygen is depleted almost year-round, indicating active redox processes in this domain.

In contrast to nitrogen species, silicate does not participate in redox reactions and therefore its distribution in the ice mainly reflects the surface water concentrations during ice buildup, which are mainly driven by phytoplankton uptake and mixing with meltwater and deep water (Figure 6A). Similarly, silicate in the sediments mainly reflects bottom water concentrations, which are mainly driven by remineralization, mixing with surface water, and transport into the sediments. In the case of no additional supply of OM, and in absence of low oxygen conditions, mineral forms of nitrogen (e.g., ${\mathrm{NO}}_3^{-}$) would have very similar profiles to silicate. As it is, the low wintertime concentrations of oxygen in the surface water (Figure 5C), and in the latter case bacteria can use nitrate as an oxidizing agent. Therefore, during ice freezing in some ice layers (where there is more oxygen) nitrification occurs and within other layers (where there is less oxygen) denitrification and anammox processes prevail (see supplementary material). By the start of the melting season there are almost no active nitrogen redox transformations in the ice (see supplementary material), and the distributions of nitrates, nitrites, and ammonium in the ice are mainly determined by melting processes.

In Figure 6 and Figure 7 it is demonstrated that a long ice melting season leads to a strong stratification in the water column that prevents oxygen supply from the surface layer downwards. OM produced during the phytoplankton bloom leads to oxygen consumption in the subsurface layers, thereby triggering the process of denitrification and the consumption of nitrate and nitrite for OM decomposition. Variability of nitrogen species illustrates this as a temporal decrease of concentrations of nitrate and an increase of ammonia in the water column (Figure 6). Before and after this nitrate minimum, nitrate maxima are formed. For this oxygen depleted period the model also predicts an increase in content of dissolved ${\mathrm{Mn}}^{2+}$ and dissolved ${\mathrm{Fe}}^{2+}$ (Figure 7B,C), connected with the reduction of oxidized forms of these metals. Finally, trace concentrations of hydrogen sulfide appear in the bottom water (Figure 7A) due to the process of sulphate reduction that starts when oxidized forms of nitrogen, manganese, and iron are depleted. This sequence of changes of electron acceptors corresponds to the theory [39] and to the classical features of the structure of the water column redox interfaces, e.g., in the Black Sea or the Baltic Sea [40,41,42]. Similar temporal changes are observed in places with variable redox conditions, e.g., Norwegian fjords and coastal lagoons [43,44,45]. The duration and intensity of oxygen depletion in the water column (and therefore in the sediments) are determined by the oxygen and OM supply in combination with the peculiarity of the ice regime (time and duration of the ice-melting period), that affect the distributions of chemical and biological parameters in the water column.

4. Discussion

Overall, the test simulations demonstrate potentially important interactions between ice, water, and sediment biogeochemistry, as well as how distinct vertical structure can emerge in the ice (Figure 3, Figure 4, Figure 5, Figure 6 and Figure 7) and in the sediments (Figure 5, Figure 6 and Figure 7). This suggests that SPBM is a potentially useful tool for marine biogeochemical modeling in the Arctic. SPBM is not restricted to a particular biogeochemical model. Instead, it can calculate the transport of variables provided by any model (or models) already available in the FABM library or written according to the FABM application programming interface. SPBM has been developed initially to study vertical interactions between ice, water column, and sediments. However, by setting the number of ice or sediment layers to zero a user can choose domains of interest. Also, SPBM can be applied to test newly developed biogeochemical models since the user can vary the SPBM grid from a box to a multilayer model. Compared to 3D models, 1D tools are less suitable for forecasts but can be used as “complex calculators” for processes investigations. In this regard, SPBM provides an important ability to quantify fluxes of biogeochemical elements between the ice, water, and sediment domains. SPBM thus contains all necessary domains and is an ideally suited instrument for studying polar biogeochemistry especially in shallow waters. However, there are some important limitations of the present version that the user should consider.

First, as a 1D model, SPBM does not explicitly account for horizontal transports (advection/diffusion) which may be significant in the water column. Horizontal mixing fluxes can however be implemented, with depth/time-varying mixing concentrations and mixing rate coefficients. Of course, this latter is more likely to give reasonable results if the real horizontal transports are of a mixing/exchange character, rather than e.g., a persistent advective flux divergence. An alternative that could be implemented in future versions is to allow arbitrary depth/time-varying horizontal transport contributions or advective timescales to be read from file, perhaps based on a 3D biogeochemical model simulation e.g., [46].

Second, the present ice parameterization in SPBM is most suitable for one-year-old sea ice since it is largely based on formulations from [26]. If this is not adequate, the ice brine volume and diffusion coefficients can alternatively be taken from an ice thermodynamic model (if available). Also, the present SPBM implementation of gases in ice takes into account only the dissolved part of them and does not include bubbles. Since most of the biogeochemistry processes in ice occur in brine channels it is not crucial in the context of oxygen availability for redox reaction representation. But the fact that this process is not included can result in overestimation of initial values for dissolved gases incorporated in an ice core (e.g., [38] estimate the bubbles contribute roughly a third of the total oxygen content in ice). This will also be addressed in further work.

A third potential weakness of SPBM is in the parameterization of transport in the sediments (porosity, diffusion, and burial velocities). Equations (2) and (3) assume a fixed (time-invariant) porosity profile, a fixed deep burial velocity for solutes and particulates, and no net contribution of biogeochemical transformations to the total particulate volume fraction (see [22] Appendix B). This in turn implies a fixed total particulate volume flux or “sedimentation rate” at the SWI. Future versions might allow some temporal variability in this total volume flux, perhaps using input files and/or including an explicit contribution from the seasonal sinking flux of SPBM-modeled particulates (e.g., PON). A subtlety with the latter approach is that if, as in SPBM, the bottom layer of the water column is considered a “fluff layer” with particles entering at sinking velocity and leaving at burial velocity, then additional assumptions are required to determine the flux of modeled particulates that enters the sediments and becomes part of the sediment matrix, rather than remaining as “fluff” on the SWI [22]. A related issue here is the lack of explicit erosion/resuspension processes in the present SPBM model. Within the sediments, the neglect of solute adsorption in the present version may also be an issue in some applications.

Fourth: the ability of FABM to combine state variables from different models in a modular fashion (as specified in the configuration file) should accelerate the development of well-parameterized sympagic–pelagic–benthic biogeochemical modules within SPBM, and will also ensure that the same module code is used in any subsequent 3D simulations as long as the 3D model is also FABM-coupled (e.g., ROMS, FVCOM, GETM, MOM, NEMO). However, care is needed to ensure compatibility between state variables, and differences in model structure and currencies may present obstacles. For example, in test case 1 we combined state variables from ERSEM and BROM, seeking to combine the pelagic process resolution in ERSEM with the BROM resolution of redox biogeochemistry and sedimentary nutrient recycling. However, it was not possible to fully couple these models because they use different currencies (BROM uses nitrogen units, while ERSEM uses carbon, nitrogen, phosphorus, and silicon). Complete coupling of ERSEM and BROM may ultimately require some recoding of BROM state variable modules to allow for flexible elemental stoichiometry. Furthermore, while FABM allows modules to be repurposed to describe domain-specific variables (e.g., ice diatoms from the ERSEM primary producer module in test case 1) the user must exercise caution to ensure that parameter values are suitably adapted and that the module has sufficient flexibility to describe the domain-specific variable (if not, a new FABM module may need to be written). Also, a structural approximation that may be adequate in one domain (e.g., the use of one lability class for DON in BROM) may not be adequate in another (e.g., for DON in ice, as in test case 2).

Finally, the SPBM approach of simulating all variables in all domains implies some computationally inefficiency, e.g., calculating the transport of minute quantities of phytoplankton in the sediments. Future versions could implement domain-specific screening switches to avoid this, although in a 1D context this may not be necessary. Also, a base assumption that “everything is everywhere, but the environment selects” [47] can be insightful. For example, in test case 1, SPBM simulated some limited growth of phytoplankton groups in the ice domain, and this may be important in seeding the phytoplankton blooms in the water column following ice melt [48].

5. Conclusions

We aimed to develop a flexible 1D vertical transport model to allow simultaneous simulation of the marine biogeochemistry of 3 different media: ice, water, and sediments. The resulting sympagic–pelagic–benthic model (SPBM) includes vertically-resolved ice and sediment domains, and allows fine resolution of the benthic boundary layer. SPBM reads input file data on ice growth and water column physics (and optionally also brine volumes and ice diffusivity) and uses the Framework for Aquatic Biogeochemical Models (FABM) to provide a user-defined model for biogeochemical transformations and water column sinking velocities, based on published models in the FABM library and possible new modules written by the user. Two test simulations demonstrated the potential utility of SPBM for modeling systems with strong interactions between ice, water column, and sediment biogeochemistry, as are often found in the Arctic Ocean and shelf seas. In the first test case, the FABM coupling was used to combine modules from two complex biogeochemical models (ERSEM and BROM) and to adapt an existing ERSEM diatom parameterization to simulate a new sea ice diatom group in combination with the other four ERSEM phytoplankton groups. The simulation demonstrated a strong interaction between water column and ice domains with respect to algal blooms, with the ice providing seed populations of phytoplankton and the water column providing an income of nutrients. It also demonstrated that different groups of primary producers have different spatial and temporal variabilities both in the ice and water domains due to different requirements and limitations. A second test case demonstrated strong interactions between ice, water, and sediment domains, with spatial variability of nutrients in sea water during sea ice congelation season determining the processes occurring in the ice core in the following winter, and the melting season features determining the redox reactions occurring in the sediments. Although there are some notable limitations of the present SPBM version, the results herein suggest that SPBM can already provide a useful tool for tuning existing biogeochemical models, accelerating the development of new biogeochemical models for regions with strong interactions between ice, water column, and sediment, and for investigating the potential importance of such interactions in determining the response of Arctic ecosystems to local and global anthropogenic drivers.