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

: 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-deﬁned biogeochemical model. SPBM coupled with a biogeochemical model simulates the processes of vertical di ﬀ usion, 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 ﬁrst 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 ﬂexible 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.


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.

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 where C i is the i-th state variable in units provided by the biogeochemical model through FABM, (mmol m −3 total volume) or (mg m −3 total volume) (here total volume refers to a representative control volume including both liquid and solid); t is the time step, (s); z is the depth, (m); A f is the porosity-related area restriction factor for fluxes, dimensionless; D is the total diffusivity, (m 2 s −1 ; P f is the porosity factor, dimensionless; u is the sinking velocity (advection/burial in the sediments), (m s −1 ). R i is the combined sources minus sinks of the i-th state variable provided by the biogeochemical model through FABM, ( mmol m −3 total volume s −1 or (mg m −3 total volume s −1 ). The porosity factor P 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 A f is included to limit fluxes within the respective phases (intraphase mixing). The values of A f , P f , D, and 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 P f = 1. However, vertical fluxes are assumed to be restricted to the brine channels where the particulates are mobilised in suspension, hence A f = ϕ(z). Here, the dimensionless porosity ϕ(z) 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 P f = 1 ϕ(z) , and fluxes are again restricted to the brine channels, hence A f = ϕ(z). The total diffusivity D in the ice brine channels is a sum of the molecular diffusivity D m (s) (m 2 s −1 ) on the ice-water interface (applied only to solutes), the gravity drainage diffusivity D gd (z) ( m 2 s −1 at depths z within the ice, and the diffusivity caused by convection that occurs in the bottom layer of the growing ice D gi (s) ( m 2 s −1 [26]: D = D m (s) + D gd (z) + D gi (s) D gd (z) = F vb z b D gi (s) = 10 −2 z s (9.667 · 10 −9 +4.49 · 10 −6 IceGrowth − 1.39 · 10 −7 IceGrowth 2 ) where 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; F vb is a constant mean flux volume rate from the brine channels, (m s −1 ); z b is the vertical distance over which the ice column is influenced by brine tube convection (depths where ϕ(z) > ϕ min ), (m); IceGrowth is the total ice growth rate (cm s −1 ); z s is the thickness of the ice layer, (m). D gi (s) 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 D can be read from an input file generated by e.g., an ice thermodynamic model.
The sinking velocity u is non-zero only for particulate variables in the layers where ϕ(z) > ϕ min (if ϕ(z) ≤ ϕ 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], 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 D is set to zero.
In the water column domain: Here P f = 1 and A f = 1 at all depths for both solutes and particulates, since there is only one phase to consider.
The total diffusivity D is composed of the molecular diffusivity D 0 ( m 2 s −1 (applied only to solutes) and the turbulence diffusivity D t (z) ( m 2 s −1 : where D t (z) 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 z i within the BBL is linearly decreasing from the deepest non-zero value of the diffusivity D t (z d ) as follows: where z d (m) is the deepest depth with non-zero value of D t (z) and z 0 (m) is the bottom depth. The sinking velocity 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 P f = 1 1−ϕ(z) and the area restriction factor A f = 1 − ϕ(z) at depths z. For solutes P f = 1 ϕ(z) and A f = ϕ(z). There is no adsorption in the present version.
A time-independent porosity ϕ(z) at depths z through the entire sediment domain is described following [27]: where ϕ(z ∞ ) is the deep porosity, dimensionless; ϕ(z 0 ) is the porosity at the sediment-water interface (SWI), dimensionless; z swi is the depth of the SWI, (m); k ϕ is the coefficient for exponential porosity change, (m).
The total diffusivity D is a sum of the molecular diffusivity D m (z) ( m 2 s −1 (applied only to solutes) and the bioturbation diffusivity D b (z) ( m 2 s −1 [28]: where D 0 is the infinite-dilution molecular diffusivity, (m 2 s −1 ); µ d is the relative dynamic viscosity, dimensionless; O 2 is the oxygen concentration in the bottom layer of the water column, (mmol m −3 ); K O 2 is the half-saturation constant, (mmol m −3 ). The oxygen-saturated bioturbation diffusivity [22] D bo (z) ( m 2 s −1 depends on the distance z db (z) (m) between the interface depth z and the depth with a constant bioturbation activity as follows: where z swi is the depth at the SWI, (m); z cb is the constant bioturbation activity layer thickness, (m); D bm is the maximum bioturbation diffusivity, ( m 2 s −1 ; F d is the bioturbation decay scale, (m). On the SWI it is assumed that the bioturbation diffusivity mixes concentrations in units (mmol m −3 total volume) instead of (mmol m −3 solids/solutes) (interphase mixing). Therefore, special values of P f are needed for the layers immediately above and below the SWI (see Appendix B): for solutes : for solids : P f (z a,b ) = 1 1 − ϕ swi where the subscripts a, b and swi determine a location of the corresponding variables: a means the layer above, b the layer below, swi on the SWI.
The advection/burial velocities u(z) are described following [22]: for particulates :  (1), and these will exactly cancel the contributions of D inter b 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.

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 (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 (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 (m) as an input parameter. Therefore, the ice column is discretized into layers of strictly constant thickness z s , and when the ice column grows or melts its total thickness can change only by multiples of z s . This simplification facilitates recalculation of the variable concentrations during melting and freezing.
Water 2019, 11, x FOR PEER REVIEW 6 of 23 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 ( 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 (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 (m) as an input parameter. Therefore, the ice column is discretized into layers of strictly constant thickness z s , and when the ice column grows or melts its total thickness can change only by multiples of z s . This simplification facilitates recalculation of the variable concentrations during melting and freezing.

Irradiance Formulation
FABM biogeochemical models generally need to know the photosynthetically active radiation (PAR), e.g., (mol photons m -2 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 P s can be calculated from the surface shortwave radiative flux F surf (W m -2 ), depending on the solar declination k decl (degrees): 2πJulianDay -81 365

Irradiance Formulation
FABM biogeochemical models generally need to know the photosynthetically active radiation (PAR), e.g., (mol photons m −2 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.
where I m is the theoretical maximum of 24-h average surface downwelling shortwave irradiance in air, (W m −2 ); k f is the factor to convert downwelling shortwave irradiance in air to scalar PAR in water, (mol photons day −1 W −1 ) [32]. Alternatively, P s (or F surf ) can be read from an input file.
In the presence of ice, PAR after considering albedo influence P a becomes [33]: where k scatter is the fraction of radiation transmitted through the highly scattering surface of the ice, dimensionless; A ice is the ice albedo for visible light, dimensionless; A snow is the snow albedo for visible light, dimensionless; k snow is the snow light extinction coefficient, (m −1 ); z snow is the snow depth, (m). PAR at any depth in the ice P(z ice ) is given by: where k ice is the ice light extinction coefficient, (m −1 ); z ice is the ice depth, (m). PAR in the water column P(z water ) is calculated according to the Beer-Lambert formulation: if there is ice : P(z water ) = P(z IceBottom )e z 0 K water (z water )dz water if there is no ice : P(z water ) = P s e z 0 K water (z water )dz water where P(z IceBottom ) is the PAR at the ice bottom layer; K 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, (m −1 ); z water is the water layer depth, (m).
In the sediment domain, PAR equals zero in all layers.

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 2 and CO 2 . Diffusive fluxes of O 2 and CO 2 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 (m 2 s −1 ) on layer interfaces; temperature (C) and salinity (psu) on layer centers) and for the ice column (total thickness (m), snow thickness (m), and surface temperature (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 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 M and Phase parameters (C = 2 −1 M + 4 −1 M(1 + sin(365 −1 · 2π(JulianDay − Phase)))); set as fixed constant value. M and the boundary concentrations C should be in units corresponding to the state variables of the appropriate FABM model, Phase is in ( days). SPBM uses input data files in NetCDF and text formats.

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 cm, resolution = 10 cm. The grid in the sediment domain was continued for another 10 cm with resolution 2 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 R i and the sinking velocities 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 (m s −1 ) and concentration of CO 2 in air (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, ice parameters in Table A2, sediments parameters in Table A3, irradiance parameters are  provided in Table A4, and forcing properties in Table A5). 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 (m s -1 ) and concentration of CO 2 in air (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, ice parameters in Table A2, sediments parameters in Table A3, irradiance parameters are provided in Table A4, and forcing properties in Table A5).

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

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 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). 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 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 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). 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;   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. 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.  . Variability of ERSEM primary producer functional groups during the maximum production year (1983) and its preceding year: diatoms chl a, ice diatoms chl a, nanophytoplankton chl a, picophytoplankton chl a, and microphytoplankton chl a (SPBM-ERSEM-BROM coupled system output).

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. Figures 5-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   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 Figures 6-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 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 (NO − 2 , 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., 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 Figures 6 and 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 Mn 2+ and dissolved 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.

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 (Figures 3-7) and in the sediments (Figures 5-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].

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.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
Porosityϕ(z) at depth z in the ice column is considered as relative volume of brine channels in ice [26]: Brine salinity, S b (z) (ppt) [26] and corresponding sea ice temperature (degrees Celsius), T i (z): where α 0 , α 1 , α 2 and α 3 are different for 3 ranges of temperatures: where Z p is the ratio between the distance from the ice surface and ice thickness. : where ρ 0 = 912 × 10 3 g m −3 is the density of pure ice.

Appendix B
The molecular diffusivity D m (m 2 s −1 ) mixes concentrations of the solutes in units (mmol m −3 solutes). While the bioturbation diffusivity D b (m 2 s −1 ) mixes concentration of the both solutes and solids in units (mmol m −3 total volume). So there is a flux for solutes on the SWI: In the 1D model the flux is calculated in the form where the porosity factor P f (z a,b ) should be determined: Comparing Equation (A1) and Equation (A2): For solids since D m = 0 and 1−ϕ swi instead of ϕ swi : where C is the concentration of the variable, (mmol m −3 total volume); ϕ is porosity, dimensionless. Subscripts a, b, and swi determine a location of the corresponding variables: a means the layer above, b-the layer below, swi-on the SWI.  In ERSEM there are 4 adjustable parameters that influence gross production without affecting nutrient or temperature limitation: maximum specific productivity at reference temperature (g), initial slope of PI-curve (α), photoinhibition parameter (β), and maximum effective chlorophyll to carbon photosynthesis ratio (φ). Tuning these parameters one can get different photosynthesis-irradiance curves which would represent different irradiance requirements [21] (Figure 13, p. 1333). We wanted our new primary producer functional group to be more tolerant to lower PAR conditions, but we did not want to change significantly its behavior. So, we adjusted only g, α, and β in a way that increased the initial slope of the photosynthesis-irradiance curve but preserved the area under this curve from 0 to 20 Wm −2 . The parameters for the original ERSEM diatom and the derived parameters for the new ice diatom are presented in Table A6.