A 1-Dimensional Ice-Pelagic-Benthic transport model (IPBM) v0.1: Coupled simulation of ice, water column, and sediment biogeochemistry

Aquatic biogeochemical processes can strongly interact, especially in polar regions, with processes occurring in adjacent ice and sediments layers, yet there are few modelling tools to simulate these systems in a fully coupled manner. We developed a 1D Ice-Pelagic-Benthic transport model (IPBM) for coupled simulation of ice, water column, and upper sediments biogeochemistry. IPBM describes the processes of diffusion and particle sinking in both ice and water, as well as sedimentation and bioturbation processes in the sediments. To describe ice, pelagic, and benthic biogeochemical dynamics (reaction terms), 5 IPBM was partly coupled to the European Regional Seas Ecosystem Model (ERSEM) and partly to the Bottom RedOx Model biogeochemistry module (BROM-biogeochemistry) using the Framework for Aquatic Biogeochemical Models (FABM). To test the coupled system, hydrophysical forcing for a site in the Kara Sea area from a Regional Oceanic Modeling System (ROMS) simulation was used. The test run showed reasonable results for all main variables. IPBM reproduces the ice algae bloom in July followed by a pelagic phytoplankton bloom in August September, as well as seasonal variability of nutrients in 10 the water column.


Introduction
Arctic ecosystems have undergone drastic changes and the most important changes are climatically induced (Schofield et al., 2010;Bellerby et al., 2005Bellerby et al., , 2012;;Denman et al., 2011;Silyakova et al., 2013).Community Climate System Model studies have projected atmospheric warming in the Arctic of 1.5 -4.5 times the global mean warming (Holland and Bitz, 2003) and the Arctic marine environment is expected to be strongly impacted by loss of ice cover, increased light exposure, ocean warming, freshening, acidification, and deoxygenation.Modeling estimates are needed for the analysis of the present conditions and projection of future changes.
A biogeochemical model suitable for the Arctic should take into account the specific conditions of this region.Since most of the Arctic is seasonally covered by ice and located in a shelf zone the model preferably should combine processes occurring in 3 domains: ice, water column, and sediments.Also each of the domains needs to implement some specific features.
Ice.The Arctic ice-algal primary production is a significant part of the total primary production of the Arctic region.Photosynthetic microorganisms extend the production season, provide a winter and early spring food source, and contribute to organic carbon export to depth (Vancoppenolle et al., 2013).A recent modeling study (Jin et al., 2012) estimated an average Arctic ice-algal primary production of 21.7 Tg C year −1 , which equates to roughly 5% of total pelagic primary production (Duarte et al., 2015) for this area.(Arrigo et al., 2010) estimates sea ice production from 5 % to 10 % of total Arctic and Southern Ocean productivity.Another modelling study suggests that under a mild climate change scenario the sea ice community is projected to be generally more productive while pelagic phytoplankton productivity may decrease (Tedesco et al., 2012).It is therefore desirable to include the ice domain in biogeochemical modeling studies of the Arctic region.Recent research suggests that ice-algal models should resolve vertically the ice to avoid biases that may result from either assuming that ice algae are solely at the bottom layer or that they are homogeneously distributed vertically (Duarte et al., 2015).
Water column.In the Arctic, global change is causing pervasive seawater acidification, accompanied by more localized changes in productivity and oxygen depletion (Bopp et al., 2013;Henson et al., 2017).Therefore the most important feature to be included in the biogeochemical model is the carbon cycle.On top of this, oxygen dynamics and redox process parameterization can be useful in areas of the Arctic affected by oxygen depletion (often in estuaries and fjords).
To improve the representation of near-bottom processes the Benthic Boundary Layer (BBL) should be incorporated into 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" (Dade et al., 2001).For the Arctic this layer is especially important since ice melting and permafrost thawing drive strong fluxes of ungrazed organic material to the bottom (Lønne, 1999).
Sediments.Sinking fluxes from the water column can provide sources of new energy for the benthic community.Also it has been shown that benthic, as well as pelagic, activity can be an important factor for annual pH variability in coastal areas (Blackford and Gilbert, 2007).Therefore sediments must respond accurately to sinking fluxes and provide the correct remineralization rates over the correct timescales.Redox processes occurring in sediments can be highly structured in the vertical, suggesting a need for explicit vertical resolution in sediment models.
Understanding and projecting the consequences of CO 2 emissions and warming on the Arctic Ocean requires analysis of the processes occurring in ice, water column, sediments, and their interactions.The goal of this work was to develop a model capable of simulating transport and biogeochemical processes in all three domains simultaneously.
As a modeling tool we have employed the Framework for Aquatic Biogeochemical Models (FABM) (Bruggeman and Bolding, 2014) with its variety of available ecosystem and biogeochemical models (such as The European Regional Seas Ecosystem Model (ERSEM) (Butenschön et al., 2016), Bottom RedOx Model -biogeochemistry (BROM-biogeochemistry) (Yakushev et al., 2017), PCLake aquatic ecosystem model (Hu et al., 2016) and etc.) FABM provides a flexible coupling of any ecosystem model within FABM with a selection of hydrodynamic models i.e. 3D General Estuarine Transport Model (GETM) (Stips et al., 2004), Regional Oceanic Modeling System (ROMS) (Shchepetkin and McWilliams, 2005), Modular Ocean Model (MOM), Nucleus for European Modeling of the Ocean (NEMO) (Madec, 2008), Finite Volume Community Ocean Model (FVCOM) (Chen et al., 2006), Semi-implicit Cross-scale Hydroscience Intergrated System Model (SCHISM) (Zhang et al., 2016) as well as a 1D General Ocean Turbulence Model (GOTM) (Umlauf et al., 2005;Burchard et al., 2006).However, none of the existing FABM hydrodynamic models provide vertically resolved structure for both ice and sediment layers.We have developed a 1D Ice-Pelagic-Benthic transport model (IPBM) to describe simultaneously the processes in ice, water column, and sediments.The model is currently coupled offline with ROMS (see Subsect. 2.6), but online with both ERSEM and BROM-biogeochemistry using FABM.The program is written using Fortran language.

IPBM: A 1D transport model
In general FABM provides three types of model variables: state variables, diagnostic variables, and dependencies.State variables are the basic elements for which the rates of changes are provided.Diagnostic variables are calculated by the biogeochemical models according to the values of the state variables at each time step.Dependencies are the physical environment parameters and interconnections within FABM.IPBM operates with state variables, sends dependencies to FABM, and outputs all necessary state and diagnostic variables in NetCDF files.Within IPBM, state variables are considered as solute or solid concentrations.

Formulation and numerical integration
IPBM 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 concentration of the ith state variable in units provided by the biogeochemical model through FABM, [mmol m −3 total volum or [mg m −3 total volume]; t is the time step, [s]; z is the depth, [m]; R i is the combined sources minus sinks of the ith state variable provided by the biogeochemical model through FABM, [mmol m −3 total volume s −1 ] or [mg m −3 total volume s −1 ]; A f is the porosity-related area restriction factor for fluxes, dimensionless; P f is the porosity factor, dimensionless; D is the total diffusivity, [m 2 s −1 ]; u is the sinking velocity for particles, [m 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 between in 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 to within the respective phases (intraphase mixing).
The values of A f , P f , D and u depend on whether these parameters are calculated in the ice, water column, or sediment domains and whether state variable can be considered as solute or solid.
In the ice domain: For solids, 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 solids are mobilised in suspension, hence where the dimensionless porosity ϕ(z) is equal to the relative volume of the brine channels in ice (Arrigo et al., 1993) (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) 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 the brine tube convection (depths where ϕ(z) > 72 ppt), [m]; IceGrowth is the total ice increase [cm sec −1 ]; z s is the thickness of the ice layer, [m].D gi (s) is equal to zero everywhere and at all times except during the period of ice build-up and only on the interface between water and ice.
The sinking velocity for particles u is non-zero only for particulate variables and generally determined at each time step by the biogeochemical model through FABM.To represent the ability of sea ice diatoms to maintain their vertical position relative to the skeletal layer (Arrigo et al., 1993) the diatom state variable sedimentation velocity is set to 3 cm day −1 within the ice column and zero on the interface between ice and water domains.
In the water column domain: Here P f = 1 and A f = 1 at all depths for both solutes and solids, since there is only one phase to consider.
The total diffusivity D is composed of the molecular diffusivity D 0 [m 2 s −1 ] and the turbulence diffusivity D t (z) [m 2 s −1 ]: 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 (see Subsect. 2.2).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 for particles u is taken from the biogeochemical model through FABM for all particulate state variables.
For all solute state variables u = 0.
In the sediments domain: Here the particulate variables become part of the sediments.So for solid state variables the porosity factor P f = 1 1−ϕ(z) and the area restriction factor A f = 1 − ϕ(z) at depths z.For solute state variables P f = 1 ϕ(z) and A f = ϕ(z).A time-independent porosity ϕ(z) at depths z through the entire sediments domain is described following (Soetaert et al., 1996): where ϕ(z ∞ ) is the porosity at the infinite sediment depth, dimensionless; ϕ(z 0 ) is the porosity at the sediment-water  , 1997): for solutes: where z swi is the depth at the SWI, [m]; z cb is the constant bioturbation activity layer width, [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 solids: where the subscripts a, b and swi determinate a location of the corresponding variables: a means the layer above, b -the layer below, swi -on the SWI.The sinking velocity u(z) is described following (Yakushev et al., 2017): for solids: where ϕ(z ∞ ) is the deep porosity, dimensionless; u b is the deep burial velocity, [m s −1 ].
Equation ( 1) is integrated numerically over a single combined (ice, water column, sediments) grid (see Subsect. 2.2) and using a constant model time step.Time stepping follows an operator splitting approach (Butenschön et al., 2012): concentrations are successively updated by contributions over one time step of diffusion, reaction, and sedimentation, in that order.Diffusive updates are calculated by a semi-implicit central-space algorithm adapted from a routine in BROM-transport (Yakushev et al., 2017) which in turn was adapted from the General Ocean Turbulence Model (GOTM) (Umlauf et al., 2005).Sedimentation updates are calculated using a first-order upwind differencing scheme.Reaction updates are calculated from forward Euler time steps.Figure 1 shows the IPBM grid structure with the constant parameter values used for the test case (these can all be redefined).For the test case, a time-dependent total ice thickness and time-independent water column structure were derived from a ROMS simulation, while the BBL was inserted with the following parameters: width = 50 cm, resolution = 10 cm.The grid in the sediments domain was continued for another 10 cm with resolution 2 cm.

Ice column structure
The total ice thickness is taken from hydrophysical model output and a fine multilayer structure is imposed for better representation of ice algal behavior.IPBM recalculates all state variables concentrations during freezing and melting following algorithms 1 and 2, where IceGrowth is the number of freezing (positive value) or melting (negative value) layers; C is the state variable concentration; i is the layer index; s is the surface water layer index; b is the bottom ice layer index.The layers are enumerated from the bottom upwards.It can be available from some hydrophysical model output.Also FABM has an ability to calculate water column PAR(z) at any depths z which is based on PAR at the water surface, but at present there is no FABM routine to propagate PAR through ice.

Algorithm 1 State variable concentration recalculation during ice freezing
IPBM 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 is calculated according to the surface radiative flux F surf [W m −2 ].This is is calculated according to the solar declination k decl [degrees]: If there is ice the PAR after considering albedo influence P a becomes (Light et al., 2008): if snow depth ≤ 5 mm: if snow depth > 5 mm: 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 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 given by: if there is an ice: P (z water ) = P (z IceBottom ) exp (−k water z water ) (26 if there is no ice: P (z water ) = P s exp (−k water z water )

Initial and boundary conditions
Initial conditions for all state variable concentrations are provided through FABM using its YAML type configuration file (Bruggeman and Bolding, 2014).
Boundary conditions for the state variables can be presented for upper and lower boundaries as: no diffusive flux at the bottom or surface boundary; diffusive fluxes of O 2 and CO 2 at the surface boundary are provided by biogeochemical models through FABM (only for ice free periods).The test case (see Sect. 4) uses forcing data from a ROMS simulation undertaken by Norwegian Meteorological Institute under the project BaSIC.The BaSIC project provides a number of the parameters available, but only depth, turbulence, temperature, salinity, ice thickness, snow thickness, and ice surface temperature are used.The period of the available data extends from 1980 till 2010.

FABM models
All of the biogeochemical models available through FABM can be used by IPBM.Most of the original FABM models were developed for use only within the water column, but IPBM can extend their applicability to the ice and sediments domains where it is possible.The FABM allows different state variables from different models to be combined in a modular fashion (using its YAML type configuration file); however care is still needed to ensure compatibility between state variables.In the test case we combine components from two biogeochemical models using the FABM: 1) the European Regional Seas Ecosystem Model (ERSEM, (Butenschön et al., 2016)), and 2) the biogeochemical modules of the Bottom RedOx Model (BROM-biogeochemistry (Yakushev et al., 2017)).The general coupling scheme is illustrated in Fig. 2.
ERSEM modules are used to model the non-bacterial lower trophic levels (plankton).These are based on the conception of the functional groups (Baretta et al., 1995;Vichi et al., 2007) which describe macroscopic role of the main ecosystem compounds.In the present implementation we use, for simplicity, only two phytoplankton groups (diatoms and microphytoplankton) and one zooplankton grooup (microzooplankton) but the configuration can be easily extended resolve more functional groups.
BROM-biogeochemistry modules are used to represent bacteria and biogeochemistry within the ice, water and sediments domains.The code of the model was significantly modified to make it work with ERSEM but the processes are the same.The former model (Yakushev et al., 2017) consists of three biogeochemical modules: ecological, redox processes, and carbonate system.To avoid duplicating ERSEM state variables, original BROM-biogeochemistry modules were split into separate modules for: equilibrium constants, carbon, calcium, pH, phytoplankton and zooplankton, bacterial, nitrogen cycle, methane, sulfur cycle, iron, manganese, and particulate silicon.

Results
The test case modelling point is located at 78.3 • N 78.8 • E i.e. the position of the ROMS grid point.Here we present IPBM output for the year 1984.The initial condition is derived from two spin-ups, repeating the first day 100 times and then repeating the first year 10 times.The main parameter values and forcing properties are provided in Appendix C (common parameters in Table C1, ice parameters in Table C2, sediments parameters in Table C3, irradiance parameters in Table C4, and forcing properties in Table C5).For demonstration purposes the snow depth was set to zero since the ROMS values were too high to allow ice algae growth during the melting season.Sediment Water Interface (c) 0.00 × 10 0 2.00 × 10 1 4.00 × 10 1 6.00 × 10 1 8.00 × 10 1 5.00 × 10 1 1.00 × 10 0 1.50 × 10 0 2.00 × 10 0 2.50 × 10 0 5.00 × 10 7 1.00 × 10 6 1.50 × 10 6 2.00 × 10 6 pH in this model is a diagnostic parameter and its value in the sea ice corresponds to the pH of the brine channels with very high salinity (greater than 200 psu).The pH shown in (Fig. 5) was calculated using the carbonate system constant given for the oceanic water.In the water column pH is characterized by a maximum (up to 7.98) connected with the algae bloom.Near the bottom pH decreases to 7.90.Ranges and pattern of the seasonal variability of pH are typical for the Arctic Ocean.

Conclusions
We have developed a 1D transport model that allows simultaneous simulation of the biogeochemistry of 3 different media: ice, water, and sediments.Description of transportation processes in ice, water, and sediments for both solutes and solids was   IPBM can be used to evaluate biogeochemical processes and their interconnections occurring in ice, water column, and sediments using biogeochemical models available through FABM.IPBM algorithms could also be used to add transport processes for ice and sediments into existing hydrophysical models coupled with FABM, and thereby extend their applicability.
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 ] and the bioturbation diffusivity D b (z) ) 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 O2 is the half-saturation constant, [mmol m −3 ].The oxygen limited bioturbation diffusivity (Yakushev et al., 2017) 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: Geosci.Model Dev.Discuss., https://doi.org/10.5194/gmd-2017-299Manuscript under review for journal Geosci.Model Dev. Discussion started: 5 December 2017 c Author(s) 2017.CC BY 4.0 License.

Figure 1 .
Figure 1.The test case grid structure Layers are enumerated from bottom to top.The ice column is discretized into layers of strictly constant thickness z s , and when 6 Geosci.Model Dev.Discuss., https://doi.org/10.5194/gmd-2017-299Manuscript under review for journal Geosci.Model Dev. Discussion started: 5 December 2017 c Author(s) 2017.CC BY 4.0 License. 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 (see Subsect. 2.3).

Require: |IceGrowth| > 0
for all i in the ice domain do Ci = C i−|IceGrowth| end for for all i of the new ice layers do Ci = Cs Cs = Cs − Cs IceLayerW idth Surf aceLayerW idth end for 2.4 Irradiance formulation FABM biogeochemical models generally need to know the surface photosynthetically active radiation (PAR) [mol photons m −2 day −1 ].
Dev. Discuss., https://doi.org/10.5194/gmd-2017-299Manuscript under review for journal Geosci.Model Dev. Discussion started: 5 December 2017 c Author(s) 2017.CC BY 4.0 License.Algorithm 2 State variable concentration recalculation during ice melting Require: |IceGrowth| > 0 if there is a diatom state variable and not the last layer has melted then for all i of the melted ice layers do C b = C b + Ci end for for all i in the ice domain except the bottom ice layer do Ci = C i+|IceGrowth| end for else for all i of the melted ice layers do Cs = Cs + Ci IceLayerW idth Surf aceLayerW idth end for for all i in the ice domain do Ci = C i+|IceGrowth| end for end if where I m is the theoretical maximum of 24-hour 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 ] (Mobley and Boss, 2012).
IPBM requires time-dependent input forcings for the water column (turbulent diffusivity [m 2 s −1 ] on layer interfaces and temperature [C • ], salinity [psu] on layer centres) 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 (for ERSEM/BROM we require wind speed [m s −1 ] and concentration of CO 2 in air [ppm]).For better estimates of surface shortwave or photosynthetically active radiation it can be read from an input file instead of using the formulas provided in Subsect.2.4.IPBM uses NetCDF format data files for input.Other important forcing functions for the model are the input fluxes of some state variables in the layers of choice.Their estimates are based on concentrations C which can be provided in 3 ways: read from text file; fixed sinusoidal variation in time defined by amplitude A, mean value M ean, and P hase parameters (C = A(M ean+sin 2π JulianDay−P hase 365 )); fixed constant value.A, M ean, and boundary concentrations C should be in units corresponding to the state variables of the appropriate FABM model, P hase is in [days].

Figure 2 .
Figure 2. The general scheme of models coupling

Figures 3 -
Figures3-5show IPBM output for the state variables chlorophyll-a and dissolved oxygen and for the diagnostic variable pH (total scale).All state variable concentrations are provided per total volume, but diagnostic variable concentrations are provided per volume of solutes/solids.

Figure 3 .
Figure 3. Seasonal variability of diatom chlorophyll-a in the ice (a), in the water column and BBL (b), and in the BBL and sediments (c).

Figure 4 .
Figure 4. Seasonal variability of dissolved oxygen in the ice (a), in the water column and BBL (b), and in the BBL and sediments (c).

Figure 5 .
Figure 5. Seasonal variability of pH in the ice (a), in the water column and BBL (b), and in the BBL and sediments (c).