2.9. Biogeochemistry

2.9.1. Aerosols

2.9.1.1. Basic Aerosols

Aerosols may be deposited on the ice and gradually work their way through it until the ice melts and they are passed into the ocean. They are defined as ice and snow volume tracers (Eq. 15 and 16 in CICE.v5 documentation), with the snow and ice each having two tracers for each aerosol species, one in the surface scattering layer (delta-Eddington SSL) and one in the snow or ice interior below the SSL.

Rather than updating aerosols for each change to ice/snow thickness due to evaporation, melting, snow-ice formation, etc., during the thermodynamics calculation, these changes are deduced from the diagnostic variables (melts, meltb, snoice, etc) in icepack_aerosol.F90. Three processes change the volume of ice or snow but do not change the total amount of aerosol, thus causing the aerosol concentration (the value of the tracer itself) to increase: evaporation, snow deposition and basal ice growth. Basal and lateral melting remove all aerosols in the melted portion. Surface ice and snow melt leave a significant fraction of the aerosols behind, but they do “scavenge” a fraction of them given by the parameter kscav = [0.03, 0.2, 0.02, 0.02, 0.01, 0.01] (only the first 3 are used in CESM, for their 3 aerosol species). Scavenging also applies to snow-ice formation. When sea ice ridges, a fraction of the snow on the ridging ice is thrown into the ocean, and any aerosols in that fraction are also lost to the ocean.

As upper SSL or interior layers disappear from the snow or ice, aerosols are transferred to the next lower layer, or into the ocean when no ice remains. The atmospheric flux faero_atm contains the rates of aerosol deposition for each species, while faero_ocn has the rate at which the aerosols are transferred to the ocean.

The aerosol tracer flag tr_aero must be set to true in icepack_in, and the number of aerosol species is set in icepack.settings; CESM uses 3.

2.9.1.2. Z-Aerosols

An alternate scheme for aerosols in sea ice is available using the brine motion based transport scheme of the biogeochemical tracers. All vertically resolved biogeochemical tracers (z-tracers), including aerosols, have the potential to be atmospherically deposited onto the snow or ice, scavenged during snow melt, and passed into the brine. The mobile fraction (discussed in Mobile and stationary phases) is then transported via brine drainage processes (Eq. (17)) while a stationary fraction (discussed in Mobile and stationary phases) adheres to the ice crystals. Snow deposition and the process of scavenging aerosols during snow melt is consistent with the basic aerosol scheme, though parameters have been generalized to accomodate potential atmospheric deposition for all z-tracers. For an example, see the scavenging parameter kscavz for z-tracers defined in icepack_zbgc_shared.F90.

Within the snow, z-tracers are defined as concentrations in the snow surface layer (\(h_{ssl}\)) and the snow interior (\(h_s-h_{ssl}\)). The total snow content of z-tracers per ice area per grid cell area, \(C_{snow}\) is

\[C_{snow} = C_{ssl}h_{ssl} + C_{sint}(h_{s}-h_{ssl})\]

One major difference in how the two schemes model snow aerosol transport is that the fraction scavenged from snow melt in the z-tracer scheme is not immediately fluxed into the ocean, but rather, enters the ice as a source of low salinity but potentially tracer-rich brine. The snow melt source is included as a surface flux condition in icepack_algae.F90.

All the z-aerosols are nonreactive with the exception of the dust aerosols. We assume that a small fraction of the dust flux into the ice has soluble iron (dustFe_sol in icepack_in) and so is passed to the dissolved iron tracer. The remaining dust passes through the ice without reactions.

To use z-aerosols, tr_zaero must be set to true in icepack_in, and the number of z-aerosol species is set in icepack.settings, TRZAERO. Note, the basic tracers tr_aero must be false and NTRAERO in icepack.settings should be 0. In addition, z-tracers and the brine height tracer must also be active. These are set in icepack_in with tr_brine and z_tracer set to true. In addition, to turn on the radiative coupling between the aerosols and the Delta-Eddington radiative scheme, shortwave must equal ’dEdd’ and dEdd_algae must be true in icepack_in.

2.9.2. Water Isotope

Water isotopes may be deposited on the ice from above or below, and gradually work their way through it until the ice melts and they are passed into the ocean. They are defined as ice and snow volume tracers (Eq. 15 and 16 in CICE.v5 documentation), with the snow and ice each having one tracer for each water isotope species.

Rather than updating water isotopes for each change to ice/snow thickness due to evaporation, melting, snow-ice formation, etc., during the thermodynamics calculation, these changes are deduced from the diagnostic variables (melts, meltb, snoice, etc) in icepack_isotope.F90. The water isotopes follow the “real” water in the sense that all of the mass budget changes that impact fresh water equally affect the water isotopes. The sources of water isotopes are precipitation (snow accumulation), condensation, and sea ice growth, including both frazil and congelation formation that take up water isotopes from the ocean. The sinks are evaporation and melting snow and sea ice. Isotopic fractionation occurs for vapor condensation and new sea ice growth, in which H2_16O (regular water, ${H_2}O$) is treated differently than H2_18O and HDO.

More information can be found in [7].

2.9.3. Brine height

The brine height, \(h_b\), is the distance from the ice-ocean interface to the brine surface. When tr_brine is set true in icepack_in and TRBRI is set equal to 1 in icepack.settings, the brine surface can move relative to the ice surface. Physically, this occurs when the ice is permeable and there is a nonzero pressure head: the difference between the brine height and the equilibrium sea surface. Brine height motion is computed in icepack_brine.F90 from thermodynamic variables and the ice microstructural state deduced from internal bulk salinities and temperature. This tracer is required for the transport of vertically resolved biogeochemical tracers.

Vertical transport processes are, generally, a result of the brine motion. Therefore the vertical transport equations for biogeochemical tracers will be defined only where brine is present. This region, from the ice-ocean interface to the brine height, defines the domain of the vertical bio-grid. The resolution of the bio-grid is specified in icepack.settings by setting the variable NBGCLYR. A detailed description of the bio-grid is given in section Grid and boundary conditions. The ice microstructural state, determined in icepack_brine.F90, is computed from sea ice salinities and temperatures linearly interpolated to the bio-grid. When \(h_b > h_i\), the upper surface brine is assumed to have the same temperature as the ice surface.

Brine height is transported horizontally as the fraction \(f_{bri} = h_b/h_i\), a volume conserved tracer. Note that unlike the sea ice porosity, brine height fraction may be greater than 1 when \(h_b > h_i\).

Changes to \(h_b\) occur from ice and snow melt, ice bottom boundary changes, and from pressure adjustments. The computation of \(h_b\) at \(t+\Delta t\) is a two step process. First, \(h_b\) is updated from changes in ice and snow thickness, ie.

(1)\[h_b' = h_b(t) + \Delta h_b|_{h_i,h_s} .\]

Second, pressure driven adjustments arising from meltwater flushing and snow loading are applied to \(h'_b\). Brine flow due to pressure forces are governed by Darcy’s equation

(2)\[w = -\frac{\Pi^* \bar{\rho} g}{\mu}\frac{h_p}{h_i}.\]

The vertical component of the net permeability tensor \(\Pi^*\) is computed as

(3)\[\Pi^* = \left(\frac{1}{h}\sum_{i=1}^N{\frac{\Delta z_i}{\Pi_i}}\right)^{-1}\]

where the sea ice is composed of \(N\) vertical layers with \(i\)th layer thickness \(\Delta z_i\) and permeability \(\Pi_i\). The average sea ice density is \(\bar{\rho}\) specified in icepack_zbgc_shared.F90. The hydraulic head is \(h_p = h_b - h_{sl}\) where \(h_{sl}\) is the sea level given by

(4)\[h_{sl} = \frac{\bar{\rho}}{\rho_w}h_i + \frac{\rho_s}{\rho_w}h_s .\]

Assuming constant \(h_i\) and \(h_s\) during Darcy flow, the rate of change of \(h_b\) is

(5)\[\frac{\partial h_b}{\partial t} = -w h_p\]

where \(w_o = \Pi^* \bar{\rho} g/(h_i\mu\phi_{top})\) and \(\phi_{top}\) is the upper surface porosity. When the Darcy flow is downward into the ice (\(w_o < 0\)), then \(\phi_{top}\) equals the sea ice porosity in the uppermost layer. However, when the flow is upwards into the snow, then \(\phi_{top}\) equals the snow porosity phi_snow specified in icepack_in. If a negative number is specified for phi_snow, then the default value is used: phi_snow \(=1 - \rho_s/\rho_w\).

Since \(h_{sl}\) remains relatively unchanged during Darcy flow, (5) has the approximate solution

(6)\[\begin{aligned} h_b(t+\Delta t) \approx h_{sl}(t+\Delta t) + [h'_b - h_{sl}(t+\Delta t)]\exp\left\{-w \Delta t\right\}.\end{aligned}\]

The contribution \(\Delta h_b|_{h_i,h_s}\) arises from snow and ice melt and bottom ice changes. Since the ice and brine bottom boundaries coincide, changes in the ice bottom from growth or melt, \((\Delta h_i)_{bot}\), equal the bottom brine boundary changes. The surface contribution from ice and snow melt, however, is opposite in sign. The ice contribution is as follows. If \(h_i > h_b\) and the ice surface is melting, ie. \((\Delta h_i)_{top} < 0\)), then meltwater increases the brine height:

(7)\[\begin{split}\begin{aligned} \left(\Delta h_b\right)_{top} = \frac{\rho_i}{\rho_o} \cdot \left\{ \begin{array}{ll} -(\Delta h_i)_{top} & \mbox{if } |(\Delta h_i)_{top}| < h_i-h_b \\ h_i-h_b & \mbox{otherwise.} \end{array} \right. \end{aligned}\end{split}\]

For snow melt (\(\Delta h_s < 0\)), it is assumed that all snow meltwater contributes a source of surface brine. The total change from snow melt and ice thickness changes is

(8)\[\Delta h_b|_{h_i,h_s} = \left( \Delta h_b\right)_{top} -\left(\Delta h_i\right)_{bot} -\frac{\rho_s}{\rho_o}\Delta h_s.\]

The above brine height calculation is used only when \(h_i\) and \(h_b\) exceed a minimum thickness, thinS, specified in icepack_zbgc_shared.F90. Otherwise

(9)\[h_b(t+\Delta t) = h_b(t) + \Delta h_i\]

provided that \(|h_{sl}-h_b| \leq 0.001\). This formulation ensures small Darcy velocities when \(h_b\) first exceeds thinS.

Both the volume fraction \(f_{bri}\) and the area-weighted brine height \(h_b\) are available for output.

(10)\[{{\sum f_{bri} v_i} \over {\sum v_i}},\]

while hbri is comparable to hi (\(h_i\))

(11)\[{{\sum f_{bri} h_i a_i} \over {\sum a_i}},\]

where the sums are taken over thickness categories.

2.9.4. Sea ice ecosystem

There are two options for modeling biogeochemistry in sea ice: 1) a skeletal layer or bottom layer model that assumes biology and biological molecules are restricted to a single layer at the base of the sea ice; and 2) a vertically resolved model (zbgc) that allows for biogeochemical processes throughout the ice column. The two models may be run with the same suite of biogeochemical tracers and use the same module algal_dyn in icepack_algae.F90 to determine the biochemical reaction terms for the tracers at each vertical grid level. In the case of the skeletal-layer model this is a single layer, while for zbgc there are NBGCLYR\(+1\) vertical layers. The primary difference between the two schemes is in the vertical transport assumptions for each biogeochemical tracer. This includes the parameterizations of fluxes between ocean and ice.

In order to run with the skeletal-layer model, the code must be built with the following options in icepack.settings:

setenv TRBGCS 1   # set to 1 for skeletal layer tracers
setenv TRBGCZ 0   # set to 1 for zbgc tracers

For zbgc with 8 vertical layers:

setenv TRBRI  1   # set to 1 for brine height tracer
setenv TRBGCS 0   # set to 1 for skeletal layer tracers
setenv TRBGCZ 1   # set to 1 for zbgc tracers
setenv NBGCLYR 7  # number of zbgc layers

There are also environmental variables in icepack.settings that, in part, specify the complexity of the ecosystem and are used for both zbgc and the skeletal-layer model. These are 1) TRALG, the number of algal species; 2) TRDOC, the number of dissolved organic carbon groups, 3) TRDIC, the number of dissolved inorganic carbon groups (this is currently not yet implemented and should be set to 0); 4) TRDON, the number of dissolved organic nitrogen groups, 5) TRFEP, the number of particulate iron groups; and 6) TRFED, the number of dissolved iron groups. The current version of algal_dyn biochemistry has parameters for up to 3 algal species (diatoms, small phytoplankton and Phaeocystis sp, respectively), 2 DOC tracers (polysaccharids and lipids, respectively), 0 DIC tracers, 1 DON tracer (proteins/amino acids), 1 particulate iron tracer and 1 dissolved iron tracer. Note, for tracers with multiple species/groups, the order is important. For example, specifying TRALG = 1 will compute reaction terms using parameters specific to ice diatoms. However, many of these parameters can be modified in icepack_in.

The complexity of the algal ecosystem must be specified in both icepack.settings during the build and in the namelist, icepack_in. The procedure is equivalent for both the skeletal-layer model and zbgc. The namelist specification is described in detail in section Vertical BGC (‘’zbgc’’)

Biogeochemical upper ocean concentrations are initialized in the subroutine icepack_init_ocean_conc in icepack_zbgc.F90 unless coupled to the ocean biogeochemistry. Silicate and nitrate may be read from a file. This option is specified in the namelist by setting the variables bgc_data_type to ISPOL or NICE. The location of forcing files is specified in data_dir and the filename is also in namelist, bgc_data_file.

2.9.4.1. Skeletal Layer BGC

In the skeletal layer model, biogeochemical processing is modelled as a single layer of reactive tracers attached to the sea ice bottom. Optional settings are available via the zbgc_nml namelist in icepack_in. In particular, skl_bgc must be true and z_tracers and solve_zbgc must both be false.

Skeletal tracers \(T_b\) are ice area conserved and follow the horizontal transport Equation (1). For each horizontal grid point, local biogeochemical tracer equations are solved in icepack_algae.F90. There are two types of ice-ocean tracer flux formulations: 1) ‘Jin2006’ modeled after the growth rate dependent piston velocity and 2) ‘constant’ modeled after a constant piston velocity. The formulation is specified in icepack_in by setting bgc_flux_type equal to ‘Jin2006’ or ‘constant’.

In addition to horizontal advection and transport among thickness categories, biogeochemical tracers (\(T_b\) where \(b = 1,\ldots, N_b\)) satisfy a set of local coupled equations of the form

(12)\[\frac{d T_b}{dt} = w_b \frac{\Delta T_b}{\Delta z} + R_b({T_j : j = 1,\ldots,N_b})\]

where \(R_b\) represents the nonlinear biochemical reaction terms (described in section Reaction Equations) and \(\Delta z\) is a length scale representing the molecular sublayer of the ice-ocean interface. Its value is absorbed in the piston velocity parameters. The piston velocity \(w_b\) depends on the particular tracer and the flux formulation.

For ‘Jin2006’, the piston velocity is a function of ice growth and melt rates. All tracers (algae included) flux with the same piston velocity during ice growth, \(dh/dt > 0\):

(13)\[\begin{aligned} w_b & = & - p_g\left|m_1 + m_2 \frac{dh}{dt} - m_3 \left(\frac{dh}{dt} \right)^2\right|\end{aligned}\]

with parameters \(m_1\), \(m_2\), \(m_3\) and \(p_g\) defined in skl_biogeochemistry in icepack_algae.F90. For ice melt, \(dh/dt < 0\), all tracers with the exception of ice algae flux with

(14)\[\begin{aligned} w_b & = & p_m\left|m_2 \frac{dh}{dt} - m_3 \left(\frac{dh}{dt} \right)^2\right| \end{aligned}\]

with \(p_m\) defined in skl_biogeochemistry. The ‘Jin2006’ formulation also requires that for both expressions, \(|w_b| \leq 0.9 h_{sk}/\Delta t\). The concentration difference at the ice-ocean boundary for each tracer, \(\Delta T_b\), depends on the sign of \(w_b\). For growing ice, \(w_b <0\), \(\Delta T_b = T_b/h_{sk} - T_{io}\), where \(T_{io}\) is the ocean concentration of tracer \(i\). For melting ice, \(w_b > 0\), \(\Delta T_b = T_b/h_{sk}\).

In ‘Jin2006’, the algal tracer (\(N_a\)) responds to ice melt in the same manner as the other tracers (14). However, this is not the case for ice growth. Unlike dissolved nutrients, algae are able to cling to the ice matrix and resist expulsion during desalination. For this reason, algal tracers do not flux between ice and ocean during ice growth unless the ice algal brine concentration is less than the ocean algal concentration (\(N_o\)). Then the ocean seeds the sea ice concentration according to

(15)\[\begin{aligned} w_b \frac{\Delta N_a}{\Delta z} = \frac{N_oh_{sk}/\phi_{sk} - N_a}{\Delta t}\end{aligned}\]

The ‘constant’ formulation uses a fixed piston velocity (PVc) for positive ice growth rates for all tracers except \(N_a\). As in ‘Jin2006’, congelation ice growth seeds the sea ice algal population according to (15) when \(N_a < N_o h_{sk}/\phi_{sk}\). For bottom ice melt, all tracers follow the prescription

(16)\[\begin{split}\begin{aligned} w_b \frac{\Delta T_b}{\Delta z} & = & \left\{ \begin{array}{ll} T_b |dh_i/dt|/h_{sk} \ \ \ \ \ & \mbox{if } |dh_i/dt|\Delta t/h_{sk} < 1 \\ T_b/\Delta t & \mbox{otherwise.} \end{array} \right. \end{aligned}\end{split}\]

A detailed description of the biogeochemistry reaction terms is given in section Reaction Equations.

2.9.4.2. Vertical BGC (‘’zbgc’’)

In order to solve for the vertically resolved biogeochemistry, several flags in icepack_in must be true: a) tr_brine, b) z_tracers, and c) solve_zbgc.

• tr_brine = true turns on the dynamic brine height tracer, \(h_b\), which defines the vertical domain of the biogeochemical tracers. z-Tracer horizontal transport is conserved on ice volume\(\times\)brine height fraction.

• z_tracers = true indicates use of vertically resolved biogeochemical and z-aerosol tracers. This flag alone turns on the vertical transport scheme but not the biochemistry.

• solve_zbgc = true turns on the biochemistry for the vertically resolved tracers and automatically turns on the algal nitrogen tracer flag tr_bgc_N. If false, tr_bgc_N is set false and any other biogeochemical tracers in use are transported as passive tracers. This is appropriate for the black carbon and dust aerosols specified by tr_zaero true.

With the above flags, the default biochemistry is a simple algal-nitrate system: tr_bgc_N and tr_bgc_Nit are true. Options exist in icepack_in to use a more complicated ecosystem which includes up to three algal classes, two DOC groups, one DON pool, limitation by nitrate, silicate and dissolved iron, sulfur chemistry plus refractory humic material.

The icepack_in namelist options are described in the Tables of Namelist Options.

Vertically resolved z-tracers are brine- volume conserved and thus depend on both the ice volume and the brine height fraction tracer (\(v_{in}f_b\)). These tracers follow the conservation equations for multiply dependent tracers (see, for example Equation (17) where \(a_{pnd}\) is a tracer on \(a_{lvl}a_{i}\))

The following sections describe the vertical transport equation for mobile tracers, the partitioning of tracers into mobile and stationary fractions and the biochemical reaction equations. The vertical bio-grid is described in the Grid and boundary conditions section.

2.9.4.2.1. Mobile and stationary phases

Purely mobile tracers are tracers which move with the brine and thus, in the absence of biochemical reactions, evolve like salinity. For vertical tracer transport of purely mobile tracers, the flux conserved quantity is the bulk tracer concentration multiplied by the ice thickness, i.e. \(C = h\phi [c]\), where \(h\) is the ice thickness, \(\phi\) is the porosity, and \([c]\) is the tracer concentration in the brine. \(\phi\), \([c]\) and \(C\) are defined on the interface bio grid (igrid):

\[\mbox{igrid}(k) = \Delta (k-1) \ \ \ \mbox{for }k = 1:n_b+1 \ \ \mbox{and }\Delta = 1/n_b.\]

The biogeochemical module solves the following equation:

(17)\[\begin{aligned} \frac{\partial C}{\partial t} & =& \frac{\partial }{\partial x}\left\{ \left( \frac{v}{h} + \frac{w_f}{h\phi} - \frac{\tilde{D}}{h^2\phi^2}\frac{\partial \phi}{\partial x} \right) C + \frac{\tilde{D}}{h^2\phi}\frac{\partial C}{\partial x} \right\} + h\phi R([c])\end{aligned}\]

where \(D_{in} = \tilde{D}/h^2 = (D + \phi D_m)/h^2\) and \(R([c])\) is the nonlinear biogeochemical interaction term (see [29]).

The solution to (17) is flux-corrected and positive definite. This is accomplished using a finite element Galerkin discretization. Details are in Flux-corrected, positive definite transport scheme.

In addition to purely mobile tracers, some tracers may also adsorb or otherwise adhere to the ice crystals. These tracers exist in both the mobile and stationary phases. In this case, their total brine concentration is a sum \(c_m + c_s\) where \(c_m\) is the mobile fraction transported by equation (17) and \(c_s\) is fixed vertically in the ice matrix. The algae are an exception, however. We assume that algae in the stationary phase resist brine motion, but rather than being fixed vertically, these tracers maintain their relative position in the ice. Algae that adhere to the ice interior (bottom, surface), remain in the ice interior (bottom, surface) until release to the mobile phase.

In order to model the transfer between these fractions, we assume that tracers adhere (are retained) to the crystals with a time-constant of \(\tau_{ret}\), and release with a time constant \(\tau_{rel}\), i.e.

\[\begin{split}\begin{aligned} \frac{\partial c_m}{\partial t} & = & -\frac{c_m}{\tau_{ret}} + \frac{c_s}{\tau_{rel}} \\ \frac{\partial c_s}{\partial t} & = &\frac{c_m}{\tau_{ret}} - \frac{c_s}{\tau_{rel}}\end{aligned}\end{split}\]

We use the exponential form of these equations:

\[\begin{aligned} c_m^{t+dt} & = & c_m^t\exp\left(-\frac{dt}{\tau_{ret}}\right) + c^t_s\left(1-\exp\left[-\frac{dt}{\tau_{rel}}\right]\right) \end{aligned}\]

\[\begin{aligned} c_s^{t+dt} & = & c_s^t\exp\left(-\frac{dt}{\tau_{rel}}\right) + c_m^t\left(1-\exp\left[-\frac{dt}{\tau_{ret}}\right]\right) \end{aligned}\]

The time constants are functions of the ice growth and melt rates (\(dh/dt\)). All tracers except algal nitrogen diatoms follow the simple case: when \(dh/dt \geq 0\), then \(\tau_{rel} \rightarrow \infty\) and \(\tau_{ret}\) is finite. For \(dh/dt < 0\), then \(\tau_{ret} \rightarrow \infty\) and \(\tau_{rel}\) is finite. In other words, ice growth promotes transitions to the stationary phase and ice melt enables transitions to the mobile phase.

The exception is the diatom pool. We assume that diatoms, the first algal nitrogen group, can actively maintain their relative position within the ice, i.e. bottom (interior, upper) algae remain in the bottom (interior, upper) ice, unless melt rates exceed a threshold. The namelist parameter algal_vel sets this threshold.

The variable bgc_tracer_type determines the mobile to stationary transition timescales for each z-tracer. It is multi-dimensional with a value for each z-tracer. For bgc_tracer_type``(k) equal to -1, the kth tracer remains solely in the mobile phase. For ``bgc_tracer_type equal to 1, the tracer has maximal rates in the retention phase and minimal in the release. For bgc_tracer_type equal to 0, the tracer has maximal rates in the release phase and minimal in the retention. Finally for bgc_tracer_type equal to 0.5, minimum timescales are used for both transitions. Table Types of Mobile and Stationary Transitions summarizes the transition types. The tracer types are: algaltype_diatoms, algaltype_sp (small plankton), algaltype_phaeo (phaeocystis), nitratetype, ammoniumtype, silicatetype, dmspptype, dmspdtype, humtype, doctype_s (saccharids), doctype_l (lipids), dontype_protein, fedtype_1, feptype_1, zaerotype_bc1 (black carbon class 1), zaerotype_bc2 (black carbon class 2), and four dust classes, zaerotype_dustj, where j takes values 1 to 4. These may be modified to increase or decrease retention. Another option is to alter the minimum tau_min and maximum tau_max timescales which would impact all the z-tracers.

Types of Mobile and Stationary Transitions

bgc_tracer_type	\(\tau_{ret}\)	\(\tau_{rel}\)	Description
-1.0	\(\infty\)	0	entirely in the mobile phase
0.0	min	max	retention dominated
1.0	max	min	release dominated
0.5	min	min	equal but rapid exchange
2.0	max	max	equal but slow exchange

The fraction of a given tracer in the mobile phase is independent of ice depth and stored in the tracer variable zbgc_frac. The horizontal transport of this tracer is conserved on brine volume and so is dependent on two tracers: brine height fraction (\(f_b\)) and ice volume (\(v_{in}\)). The conservation equations are given by

\[{\partial\over\partial t} (f_{b}v_{in}) + \nabla \cdot (f_{b}v_{in} {\bf u}) = 0.\]

The tracer, zbgc_frac, is initialized to 1 during new ice formation, because all z-tracers are initially in the purely mobile phase. Similarly, as the ice melts, z-tracers return to the mobile phase. Very large release timescales will prevent this transition and could result in an unphysically large accumulation during the melt season.

2.9.4.2.2. Flux-corrected, positive definite transport scheme

Numerical solution of the vertical tracer transport equation is accomplished using the finite element Galerkin discretization. Multiply (17) by “w” and integrate by parts

\[\begin{split}\begin{aligned} & & \int_{h}\left[ w\frac{\partial C}{\partial t} - \frac{\partial w}{\partial x} \left(-\left[\frac{v}{h} + \frac{w_f}{h\phi}\right]C + \frac{D_{in}}{\phi^2}\frac{\partial \phi}{\partial x}C - \frac{D_{in}}{\phi}\frac{\partial C}{\partial x} \right) \right]dx \nonumber \\ & + & w\left.\left( -\left[\frac{1}{h}\frac{dh_b}{dt}+ \frac{w_f}{h\phi}\right]C + \frac{D_{in}}{\phi^2}\frac{\partial \phi}{\partial x}C -\frac{D_{in}}{\phi}\frac{\partial C}{\partial x}\right)\right|_{bottom} + w\left[\frac{1}{h}\frac{dh_t}{dt} + \frac{w_f}{h\phi}\right]C|_{top} = 0\end{aligned}\end{split}\]

The bottom boundary condition indicated by \(|_{bottom}\) satisfies

\[\begin{split}\begin{aligned} -w\left.\left( -\left[\frac{1}{h}\frac{dh_b}{dt}+ \frac{w_f}{h\phi}\right]C + \frac{D_{in}}{\phi^2}\frac{\partial \phi}{\partial x}C -\frac{D_{in}}{\phi}\frac{\partial C}{\partial x}\right)\right|_{bottom} & = & \nonumber \\ w\left[\frac{1}{h}\frac{dh_b}{dt} + \frac{w_f}{h \phi_{N+1}}\right](C_{N+2} \mbox{ or }C_{N+1}) - w\frac{D_{in}}{\phi_{N+1}(\Delta h+g_o)}\left(C_{N+1} - C_{N+2}\right) & & \end{aligned}\end{split}\]

where \(C_{N+2} = h\phi_{N+1}[c]_{ocean}\) and \(w = 1\) at the bottom boundary and the top. The component \(C_{N+2}\) or \(C_{N+1}\) depends on the sign of the advection boundary term. If \(dh_b + w_f/\phi > 0\) then use \(C_{N+2}\) otherwise \(C_{N+1}\).

Define basis functions as linear piecewise, with two nodes (boundary nodes) in each element. Then for \(i > 1\) and \(i < N+1\)

\[\begin{split}\begin{aligned} w_i(x) & = & \left\{ \begin{array}{llll} 0 & x < x_{i-1} \\ (x - x_{i-1})/\Delta & x_{i-1}< x \leq x_{i} \\ 1 - (x-x_i)/\Delta & x_i \leq x < x_{i+1} \\ 0, & x \geq x_{i+1} \end{array} \right.\end{aligned}\end{split}\]

For \(i=1\)

\[\begin{split}\begin{aligned} w_1(x) & = & \left\{ \begin{array}{ll} 1 - x/\Delta & x < x_{2} \\ 0, & x \geq x_{2} \end{array} \right.\end{aligned}\end{split}\]

and \(i = N+1\)

\[\begin{split}\begin{aligned} w_{N+1}(x) & = & \left\{ \begin{array}{ll} 0, & x < x_{N} \\ (x-x_{N})/\Delta & x \geq x_{N} \end{array} \right.\end{aligned}\end{split}\]

Now assume a form

\[C_h = \sum_j^{N+1} c_j w_j\]

Then

\[\begin{split}\begin{aligned} \int_h C_h dx & = & c_1\int_0^{x_{2}}\left(1-\frac{x}{\Delta}\right)dx + c_{N+1}\int_{x_N}^{x_{N+1}}\frac{x-x_{N}}{\Delta}dx \nonumber \\ & + & \sum_{j=2}^{N}c_j\left\{\int_{j-1}^{j}\frac{x-x_{j-1}}{\Delta}dx + \int_{j}^{j+1}\left[1 - \frac{(x-x_j)}{\Delta}\right]dx\right\} \nonumber \\ & = & \Delta \left[\frac{c_1}{2} + \frac{c_{N+1}}{2} + \sum_{j = 2}^{N}c_{j}\right]\end{aligned}\end{split}\]

Now this approximate solution form is substituted into the variational equation with \(w = w_h \in \{w_j\}\)

\[\begin{split}\begin{aligned} 0 &= & \int_{h}\left[ w_h\frac{\partial C_h}{\partial t} - \frac{\partial w_h}{\partial x} \left(\left[-\frac{v}{h} - \frac{w_f}{h\phi}+ \frac{D_{in}}{\phi^2}\frac{\partial \phi}{\partial x}\right]C_h - \frac{D_{in}}{\phi}\frac{\partial C_h}{\partial x} \right) \right]dx \nonumber \\ & + & w_h\left.\left( -\left[\frac{1}{h}\frac{dh_b}{dt}+ \frac{w_f}{h\phi}\right]C_h + \frac{D_{in}}{\phi^2}\frac{\partial \phi}{\partial x}C -\frac{D_{in}}{\phi}\frac{\partial C_h}{\partial x}\right)\right|_{bottom} + w_h\left[\frac{1}{h}\frac{dh_t}{dt} + \frac{w_f}{h\phi}\right]C_h|_{top} \end{aligned}\end{split}\]

The result is a linear matrix equation

\[M_{jk}\frac{\partial C_k(t)}{\partial t} = [K_{jk}+S_{jk}] C_k(t) + q_{in}\]

where

\[\begin{split}\begin{aligned} M_{jk} & = & \int_h w_j(x)w_k(x)dx \nonumber \\ K_{jk} & = & \left[-\frac{v}{h} - \frac{w_f}{h\phi}+ \frac{D_{in}}{\phi^2}\frac{\partial \phi}{\partial x}\right]\int_h \frac{\partial w_j}{\partial x} w_k dx \nonumber \\ &-& \left. w_j\left(-\left[\frac{v}{h} + \frac{w_f}{h\phi_k}\right]w_k + \frac{D_{in}}{\phi^2}\frac{\partial \phi_k}{\partial x} w_k - \frac{D_{in}}{\phi_k}\frac{\partial w_k}{\partial x}\right)\right|_{bot} \nonumber \\ & = & -V_k\int_h \frac{\partial w_j}{\partial x} w_k dx - \left. w_j\left(-V_kw_k - \frac{D_{in}}{\phi_k}\frac{\partial w_k}{\partial x}\right)\right|_{bot} \nonumber \\ S_{jk} & = & -\frac{D_{in}}{\phi_k}\int_h \frac{\partial w_j}{\partial x} \cdot \frac{\partial w_k}{\partial x}dx \nonumber \\ q_{in} & = & -V C_{t} w_j(x)|_{t}\end{aligned}\end{split}\]

and \(C_{N+2} = h\phi_{N+1}[c]_{ocean}\).

For the top condition \(q_{in}\) is applied to the upper value \(C_2\) when \(VC_t < 0\), i.e. \(q_{in}\) is a source.

Compute the \(M_{jk}\) integrals:

\[\begin{split}\begin{aligned} M_{jj} & = & \int_{x_{j-1}}^{x_j}\frac{(x- x_{j-1})^2}{\Delta^2}dx + \int_{x_{j}}^{x_{j+1}}\left[ 1-\frac{(x- x_{j})}{\Delta}\right]^2dx = \frac{2\Delta}{3} \ \ \ \mbox{for }1 < j < N+1 \nonumber \\ M_{11} & = & \int_{x_{1}}^{x_{2}}\left[ 1-\frac{(x- x_{2})}{\Delta}\right]^2dx = \frac{\Delta}{3} \nonumber \\ M_{N+1,N+1} & = &\int_{x_{N}}^{x_{N+1}}\frac{(x- x_{N})^2}{\Delta^2}dx = \frac{\Delta}{3}\nonumber \end{aligned}\end{split}\]

Off diagonal components:

\[\begin{split}\begin{aligned} M_{j,j+1} & = & \int_{x_{j}}^{x_{j+1}}\left[1 - \frac{(x- x_{j})}{\Delta}\right]\left[ \frac{x-x_{j}}{\Delta}\right]dx = \frac{\Delta}{6} \ \ \ \mbox{for }j < N+1 \nonumber \\ M_{j,j-1} & = & \int_{x_{j-1}}^{x_{j}}\left[ \frac{x-x_{j-1}}{\Delta}\right]\left[1 - \frac{(x- x_{j-1})}{\Delta}\right]dx = \frac{\Delta}{6} \ \ \ \mbox{for }j > 1 \nonumber \\\end{aligned}\end{split}\]

Compute the \(K_{jk}\) integrals:

\[\begin{split}\begin{aligned} K_{jj} &=& k'_{jj}\left[ \int_{x_{j-1}}^{x_j} \frac{\partial w_j}{\partial x} w_j dx + \int_{x_{j}}^{x_{j+1}} \frac{\partial w_j}{\partial x} w_j dx \right] \nonumber \\ & = & \frac{1}{2} + -\frac{1}{2} = 0 \ \ \ \mbox{for } 1 < j < N+1 \nonumber \\ K_{11} & = & -\frac{k'_{11}}{2} = \frac{1}{2}\left[\frac{v}{h} + \frac{w_f}{h\phi}\right] \nonumber \\ K_{N+1,N+1} & = & \frac{k'_{N+1,N+1}}{2} +\min\left[0, \left(\frac{1}{h}\frac{dh_b}{dt} + \frac{w_f}{h \phi_{N+1}}\right)\right] - \frac{D_{in}}{\phi_{N+1}(g_o/h)} \nonumber \\ & = & \left[-\frac{v}{h} - \frac{w_f}{h\phi}+ \frac{D_{in}}{\phi^2}\frac{\partial \phi}{\partial x}\right]\frac{1}{2} +\min\left[0, \left(\frac{1}{h}\frac{dh_b}{dt} + \frac{w_f}{h \phi_{N+1}}\right)\right] - \frac{D_{in}}{\phi_{N+1}(g_o/h)} \end{aligned}\end{split}\]

Off diagonal components:

\[\begin{split}\begin{aligned} K_{j(j+1)} &=& k'_{j(j+1)}\int_{x_{j}}^{x_{j+1}} \frac{\partial w_j}{\partial x} w_{j+1} dx = -k'_{j(j+1)}\int_{x_{j}}^{x_{j+1}}\frac{(x-x_j)}{\Delta^2}dx \nonumber \\ & = & -\frac{k'_{j(j+1)}}{\Delta^2}\frac{\Delta^2}{2} = -\frac{k'_{j(j+1)}}{2} = p5*\left[\frac{v}{h} + \frac{w_f}{h\phi}- \frac{D_{in}}{\phi^2}\frac{\partial \phi}{\partial x}\right]_{(j+1)} \ \ \ \mbox{for } j < N+1 \nonumber \\ K_{j(j-1)} &=& k'_{j(j-1)}\int_{x_{j-1}}^{x_{j}} \frac{\partial w_j}{\partial x} w_{j-1} dx = k'_{j(j-1)}\int_{x_{j-1}}^{x_{j}}\left[1 - \frac{(x-x_{j-1})}{\Delta^2}\right] dx \nonumber \\ & = & \frac{k'_{j(j-1)}}{\Delta^2}\frac{\Delta^2}{2} = \frac{k'_{j(j-1)}}{2} = -p5*\left[\frac{v}{h} + \frac{w_f}{h\phi}- \frac{D_{in}}{\phi^2}\frac{\partial \phi}{\partial x}\right]_{(j-1)} \ \ \ \ \ \mbox{for } j > 1 \end{aligned}\end{split}\]

For \(K_{N+1,N}\), there is a boundary contribution:

\[K_{N+1,N} = \frac{k'_{N+1(N)}}{2} - \frac{D_N}{\Delta \phi_{N}}\]

The bottom condition works if \(C_{bot} = h \phi_{N+2}[c]_{ocean}\), \(\phi^2\) is \(\phi_{N+1}\phi_{N+2}\) and

\[\begin{aligned} \left. \frac{\partial \phi}{\partial x}\right|_{bot} & = & \frac{\phi_{N+2} - \phi_{N}}{2\Delta};\end{aligned}\]

then the \(D_{N+1}/\phi_{N+1}/\Delta\) cancels properly with the porosity gradient. In general

\[\begin{aligned} \left. \frac{\partial \phi}{\partial x}\right|_{k} & = & \frac{\phi_{k+2} - \phi_{k}}{2\Delta}.\end{aligned}\]

When evaluating the integrals for the diffusion term, we will assume that \(D/\phi\) is constant in an element. For \(D_{in}/i\phi\) defined on interface points, \(D_1 = 0\) and for \(j = 2,...,N\) \(D_j/b\phi_j = (D_{in}(j) + D_{in}(j+1))/(i\phi_j + i\phi_{j+1})\). Then the above integrals will be modified as follows:

Compute the \(S_{jk}\) integrals:

\[\begin{split}\begin{aligned} S_{jj} & = & - \left[\frac{D_{j-1}}{b\phi_{j-1}} \int_{x_{j-1}}^{x_j}\left( \frac{\partial w_j}{\partial x}\right)^2 dx + \frac{D_{j}}{b\phi_{j}} \int_{x_{j}}^{x_{j+1}} \left(\frac{\partial w_j}{\partial x}\right)^2 dx \right] \nonumber \\ & = & -\frac{1}{\Delta}\left[\frac{D_{j-1}}{b\phi_{j-1}} + \frac{D_{j}}{b\phi{j}}\right]\ \ \mbox{for }1 < j < N+1 \nonumber \\ S_{11} & = & \frac{s'_{11}}{\Delta} = 0 \nonumber \\ S_{N+1,N+1} & = & \frac{s'_{N+1,N+1}}{\Delta} = -\frac{(D_{N})}{b\phi_{N}\Delta}\end{aligned}\end{split}\]

Compute the off-diagonal components of \(S_{jk}\):

\[\begin{split}\begin{aligned} S_{j(j+1)} & = & s'_{j(j+1)}\int_{x_j}^{x_{j+1}}\frac{\partial w_j}{\partial x} \frac{\partial w_{j+1}}{\partial x} dx = -\frac{s'_{j(j+1)}}{\Delta} = \frac{D_{j}}{b\phi_{j}\Delta} \ \ \ \mbox{for } j < N+1 \nonumber \\ S_{j(j-1)} & = & s'_{j(j-1)}\int_{x_{j-1}}^{x_{j}}\frac{\partial w_j}{\partial x} \frac{\partial w_{j-1}}{\partial x} dx = -\frac{s'_{j(j-1)}}{\Delta} = \frac{D_{j-1}}{b\phi_{j-1}} \ \ \ \mbox{for } j > 1\end{aligned}\end{split}\]

We assume that \(D/\phi^2 \partial \phi/\partial x\) is constant in the element \(i\). If \(D/\phi_j\) is constant, and \(\partial \phi/\partial x\) is constant then both the Darcy and \(D\) terms go as \(\phi^{-1}\). Then \(\phi = (\phi_j - \phi_{j-1})(x-x_j)/\Delta + \phi_j\) and \(m = (\phi_j - \phi_{j-1})/\Delta\) and \(b = \phi_j - mx_j\).

The first integral contribution to the Darcy term is:

\[\begin{split}\begin{aligned} K^1_{jj} & = & \frac{-1}{\Delta ^2}\left(\frac{w_f}{h}-\frac{D}{\phi}\frac{\partial \phi}{\partial x}\right)\int_{j-1}^{j}(x-x_{j-1})\frac{1}{mx + b}dx \nonumber \\ & = &-\left(\frac{w_f}{h}-\frac{D}{\phi}\frac{\partial \phi}{\partial x}\right) \frac{1}{\Delta^2}\left[ \int_{j-1}^{j}\frac{x}{mx + b}dx - x_{j-1}\int_{j-1}^{j}\frac{1}{mx + b}dx \right] \nonumber \\ & = &- \left(\frac{w_f}{h}-\frac{D}{\phi}\frac{\partial \phi}{\partial x}\right) \frac{1}{\Delta^2}\left[ \frac{mx - b\log(b + mx)}{m^2} - x_{j-1}\frac{\log(b+mx)}{m}\right]^{x_j}_{x_{j-1}} \nonumber \\ & = & -\left(\frac{w_f}{h}-\frac{D}{\phi}\frac{\partial \phi}{\partial x}\right)\frac{1}{\Delta_{\phi}}\left[ 1 + \log \left( \frac{\phi_j}{\phi_{j-1}} \right) - \frac{\phi_j}{\Delta_{\phi_j}}\log \left( \frac{\phi_j}{\phi_{j-1}} \right)\right] \nonumber \\ & = & -\left(\frac{w_f}{h}-\frac{D}{\phi}\frac{\partial \phi}{\partial x}\right)\frac{1}{\Delta_{\phi}}\left[ 1 + \frac{\phi_{j-1}}{\Delta_{\phi}}\log \left( \frac{\phi_j}{\phi_{j-1}} \right) \right]\end{aligned}\end{split}\]

\[\begin{split}\begin{aligned} K^2_{jj} & = & \left(\frac{w_f}{h}-\frac{D}{\phi}\frac{\partial \phi}{\partial x}\right)\frac{1}{\Delta}\int_{x_{j}}^{x_{j+1}}\left[1 - \frac{(x-x_{j})}{\Delta}\right]\frac{1}{mx+b} dx \nonumber \\ & = & \left(\frac{w_f}{h}-\frac{D}{\phi}\frac{\partial \phi}{\partial x}\right)\frac{1}{\Delta}\left[\frac{ (b+m(x_j+\Delta))\log(b+mx)-mx}{\Delta m^2}\right]_{x_{j}}^{x_{j+1}} \nonumber \\ & = & \left(\frac{w_f}{h}-\frac{D}{\phi}\frac{\partial \phi}{\partial x}\right)\frac{1}{\Delta_{\phi}}\left[ 1 - \frac{\phi_{j+1}}{\Delta_{\phi}}\log \left( \frac{\phi_{j+1}}{\phi_{j}} \right) \right]\end{aligned}\end{split}\]

Now \(m = (\phi_{j+1}-\phi_{j})/\Delta\) and \(b = \phi_{j+1} - mx_{j+1}\).

Source terms \(q_{bot} = q_{N+1}\) and \(q_{top} = q_{1}\) (both positive)

\[\begin{split}\begin{aligned} q_{bot} & = &\max\left[0, \left(\frac{1}{h}\frac{dh_b}{dt} + \frac{w_f}{h \phi_{N+1}}\right)\right]C|_{bot} + \frac{D_{in}}{\phi_{N+1}(g_o/h)}C|_{bot} \nonumber \\ C|_{bot}&= & \phi_{N+1}[c]_{ocean}\end{aligned}\end{split}\]

where \(g_o\) is not zero.

\[\begin{split}\begin{aligned} q_{in}& = & -\min\left[ 0, \left(\frac{1}{h}\frac{dh_t}{dt} + \frac{w_f}{h\phi}\right)C|_{top} \right] \nonumber \\ C|_{top}& = & h [c]_o\phi_{min}\end{aligned}\end{split}\]

Calculating the low order solution:

1) Find the lumped mass matrix \(M_l = diag\{m_i\}\)

\[\begin{split}\begin{aligned} m_j & = & \sum_i m_{ji} = m_{j(j+1)} + m_{j(j-1)} + m_{jj} \\ & = & \frac{\Delta}{6} + \frac{\Delta}{6} + \frac{2\Delta}{3} = \Delta \ \ \ \mbox{for }1 < j < N+1 \\ m_1 & = & m_{11} + m_{12} = \frac{\Delta}{3} + \frac{\Delta}{6} = \frac{\Delta}{2} \\ m_{N+1} & = & m_{N+1,N} + m_{N+1,N+1} = \frac{\Delta}{6} + \frac{\Delta}{3} = \frac{\Delta}{2}\end{aligned}\end{split}\]

2. Define artificial diffusion \(D_a\)

\[\begin{split}\begin{aligned} d_{j,(j+1)} & = & \max\{ -k_{j(j+1)},0,-k_{(j+1)j}\} = d_{(j+1)j} \\ d_{jj} & = & -\sum_{i\neq j} d_{ji}\end{aligned}\end{split}\]

3. Add artificial diffusion to \(K\): \(L = K + D_a\).

4. Solve for the low order predictor solution:

\[(M_l -\Delta t [L+S])C^{n+1} = M_l C^{n} + \Delta t q\]

Conservations terms for the low order solution are:

\[\begin{split}\begin{aligned} \int \left[C^{n+1} -C^{n}\right] w(x)dx & = & \Delta \left[\frac{c^{n+1}_1-c^{n}_1}{2} + \frac{c^{n+1}_{N+1}-c^{n}_{N+1}}{2} + \sum_{j = 2}^{N}(c^{n+1}_{j}-c^{n}_{j})\right] \nonumber \\ & = & \Delta t \left[ q_{bot} + q_{in} + (K_{N+1,N+1}+ K_{N,N+1} )C^{n+1}_{N+1} + (K_{1,1} + K_{2,1})C^{n+1}_{1}\right] \nonumber \end{aligned}\end{split}\]

Now add the antidiffusive flux: compute the F matrix using the low order solution \(c^{n+1}\). Diagonal components are zero. For \(i \neq j\)

\[\begin{aligned} f_{ij} & = & m_{ij}\left[\frac{\Delta c_i}{\Delta t} - \frac{\Delta c_j}{\Delta t} + d_{ij}(c^{n+1}_i - c^{n+1}_j\right].\end{aligned}\]

2.9.4.3. Reaction Equations

The biogeochemical reaction terms for each biogeochemical tracer (see Table Biogeochemical Tracers for tracer definitions) are defined in icepack_algae.F90 in the subroutine algal_dyn. The same set of equations is used for the bottom layer model (when skl_bgc is true) and the multi-layer biogeochemical model (when z_tracers and solve_zbgc are true).

Biogeochemical Tracers

Text Variable	Variable in code	flag	Description	units
N (1)	Nin(1)	tr_bgc_N	diatom	\(mmol\) \(N/m^3\)
N (2)	Nin(2)	tr_bgc_N	small phytoplankton	\(mmol\) \(N/m^3\)
N (3)	Nin(3)	tr_bgc_N	Phaeocystis sp	\(mmol\) \(N/m^3\)
DOC (1)	DOCin(1)	tr_bgc_DOC	polysaccharids	\(mmol\) \(C/m^3\)
DOC (2)	DOCin(2)	tr_bgc_DOC	lipids	\(mmol\) \(C/m^3\)
DON	DONin(1)	tr_bgc_DON	proteins	\(mmol\) \(C/m^3\)
fed	Fedin(1)	tr_bgc_Fe	dissolved iron	\(\mu\) \(Fe/m^3\)
fep	Fepin(1)	tr_bgc_Fe	particulate iron	\(\mu\) \(Fe/m^3\)
\({\mbox{NO$_3$}}\)	Nitin	tr_bgc_Nit	\({\mbox{NO$_3$}}\)	\(mmol\) \(N/m^3\)
\({\mbox{NH$_4$}}\)	Amin	tr_bgc_Am	\({\mbox{NH$_4$}}\)	\(mmol\) \(N/m^3\)
\({\mbox{SiO$_3$}}\)	Silin	tr_bgc_Sil	\({\mbox{SiO$_2$}}\)	\(mmol\) \(Si/m^3\)
DMSPp	DMSPpin	tr_bgc_DMS	particulate DMSP	\(mmol\) \(S/m^3\)
DMSPd	DMSPdin	tr_bgc_DMS	dissolved DMSP	\(mmol\) \(S/m^3\)
DMS	DMSin	tr_bgc_DMS	DMS	\(mmol\) \(S/m^3\)
PON	PON \(^a\)	tr_bgc_PON	passive mobile tracer	\(mmol\) \(N/m^3\)
hum	hum \(^{ab}\)	tr_bgc_hum	passive sticky tracer	\(mmol\) \(/m^3\)
BC (1)	zaero(1) \(^a\)	tr_zaero	black carbon species 1	\(kg\) \(/m^3\)
BC (2)	zaero(2) \(^a\)	tr_zaero	black carbon species 2	\(kg\) \(/m^3\)
dust (1)	zaero(3) \(^a\)	tr_zaero	dust species 1	\(kg\) \(/m^3\)
dust (2)	zaero(4) \(^a\)	tr_zaero	dust species 2	\(kg\) \(/m^3\)
dust (3)	zaero(5) \(^a\)	tr_zaero	dust species 3	\(kg\) \(/m^3\)
dust (4)	zaero(6) \(^a\)	tr_zaero	dust species 4	\(kg\) \(/m^3\)

\(^a\) not modified in algal_dyn

\(^b\) may be in C or N units depending on the ocean concentration

The biochemical reaction term for each algal species has the form:

\[\Delta {\mbox{N}}/dt = R_{{\mbox{N}}} = \mu (1- f_{graze} - f_{res}) - M_{ort}\]

where \(\mu\) is the algal growth rate, \(M_{ort}\) is a mortality loss, \(f_{graze}\) is the fraction of algal growth that is lost to predatory grazing, and \(f_{res}\) is the fraction of algal growth lost to respiration. Algal mortality is temperature dependent and limited by a maximum loss rate fraction (\(l_{max}\)):

\[M_{ort} = \min( l_{max}[{\mbox{N}}], m_{pre} \exp\{m_{T}(T-T_{max})\}[{\mbox{N}}])\]

Note, \([\cdot]\) denotes brine concentration.

Nitrate and ammonium reaction terms are given by

\[\begin{split}\begin{aligned} \Delta{\mbox{NO$_3$}}/dt & = & R_{{\mbox{NO$_3$}}} = [{\mbox{NH$_4$}}] k_{nitr}- U^{tot}_{{\mbox{NO$_3$}}} \nonumber \\ \Delta{\mbox{NH$_4$}}/dt & = & R_{{\mbox{NH$_4$}}} = -[{\mbox{NH$_4$}}] k_{nitr} -U^{tot}_{{\mbox{NH$_4$}}} + (f_{ng}f_{graze}(1-f_{gs})+f_{res})\mu^{tot} \nonumber \\ & + & f_{nm} M_{ort} \nonumber \\ & = & -[{\mbox{NH$_4$}}]k_{nitr} -U^{tot}_{{\mbox{NH$_4$}}} + N_{remin}\end{aligned}\end{split}\]

where the uptake \(U^{tot}\) and algal growth \(\mu^{tot}\) are accumulated totals for all algal species. \(k_{nitr}\) is the nitrification rate and \(f_{ng}\) and \(f_{nm}\) are the fractions of grazing and algal mortality that are remineralized to ammonium and \(f_{gs}\) is the fraction of grazing spilled or lost. Algal growth and nutrient uptake terms are described in more detail in Algal growth and nutrient uptake.

Dissolved organic nitrogen satisfies the equation

\[\begin{aligned} \Delta {\mbox{DON}}/dt & = & R_{{\mbox{DON}}} = f_{dg}f_{gs}f_{graze}\mu^{tot} - [{\mbox{DON}}]k_{nb}\end{aligned}\]

With a loss from bacterial degration (rate \(k_{nb}\)) and a gain from spilled grazing that does not enter the \({\mbox{NH$_4$}}\) pool.

A term Z\(_{oo}\) closes the nitrogen cycle by summing all the excess nitrogen removed as zooplankton/bacteria in a timestep. This term is not a true tracer, i.e. not advected horizontally with the ice motion, but provides a diagnostic comparison of the amount of \(N\) removed biogeochemically from the ice \({\mbox{N}}\)-\({\mbox{NO$_3$}}\)-\({\mbox{NH$_4$}}\)-\({\mbox{DON}}\) cycle at each timestep.

\[\begin{aligned} \mbox{Z}_{oo} & = & [(1-f_{ng}(1-f_{gs}) - f_{dg}f_{gs}]f_{graze}\mu^{tot}dt + (1-f_{nm})M_{ort}dt + [{\mbox{DON}}]k_{nb}dt \nonumber\end{aligned}\]

Dissolved organic carbon may be divided into polysaccharids and lipids. Parameters are two dimensional (indicated by superscript \(i\)) with index 1 corresponding to polysaccharids and index 2 appropriate for lipids. The \({\mbox{DOC}}^i\) equation is:

\[\begin{aligned} \Delta {\mbox{DOC}}^i/dt & = & R_{{\mbox{DOC}}} = f^i_{cg}f_{ng}\mu^{tot} + R^i_{c:n}M_{ort}-[{\mbox{DOC}}]k^i_{cb}\end{aligned}\]

Silicate has no biochemical source terms within the ice and is lost only through algal uptake:

\[\begin{aligned} \Delta {\mbox{SiO$_3$}}/dt & = & R_{{\mbox{SiO$_3$}}} = -U_{{\mbox{SiO$_3$}}}^{tot}\end{aligned}\]

Dissolved iron has algal uptake and remineralization pathways. In addition, \({\mbox{fed}}\) may be converted to or released from the particulate iron pool depending on the dissolve iron (\({\mbox{fed}}\)) to polysaccharid (\({\mbox{DOC}}(1)\)) concentration ratio. If this ratio exceeds a maximum value \(r^{max}_{fed:doc}\) then the change in concentration for dissolved and particulate iron is

\[\begin{split}\begin{aligned} \Delta_{fe}{\mbox{fed}}/dt & = & -[{\mbox{fed}}]/\tau_{fe} \nonumber \\ \Delta_{fe}{\mbox{fep}}/dt & = & [{\mbox{fed}}]/\tau_{fe}\end{aligned}\end{split}\]

For values less than \(r^{max}_{fed:doc}\)

\[\begin{split}\begin{aligned} \Delta_{fe}{\mbox{fed}}/dt & = & [{\mbox{fep}}]/\tau_{fe} \nonumber \\ \Delta_{fe}{\mbox{fep}}/dt & = & -[{\mbox{fep}}]/\tau_{fe}\end{aligned}\end{split}\]

Very long timescales \(\tau_{fe}\) will remove this source/sink term. The default value is currently set at 3065 days to turn off this dependency (any large number will do to turn it off). 61-65 days is a more realistic option (Parekh et al., 2004).

The full equation for \({\mbox{fed}}\) including uptake and remineralization is

\[\begin{aligned} \Delta {\mbox{fed}}/dt & = & R_{{\mbox{fed}}} = -U^{tot}_{{\mbox{fed}}} + f_{fa}R_{fe:n}N_{remin} + \Delta_{fe}{\mbox{fed}}/dt\end{aligned}\]

Particulate iron also includes a source term from algal mortality and grazing that is not immediately bioavailable. The full equation for \({\mbox{fep}}\) is

\[\begin{aligned} \Delta {\mbox{fep}}/dt & = & R_{{\mbox{fep}}} = R_{fe:n}[\mbox{Z}_{oo}/dt + (1-f_{fa})]N_{remin} + \Delta_{fe}{\mbox{fep}}/dt\end{aligned}\]

The sulfur cycle includes \({\mbox{DMS}}\) and dissolved DMSP (\({\mbox{DMSPd}}\)). Particulate DMSP is assumed to be proportional to the algal concentration, i.e. \({\mbox{DMSPp}}= R^i_{s:n}{\mbox{N}}^i\) for algal species \(i\). For \({\mbox{DMSP}}\) and \({\mbox{DMS}}\),

\[\begin{split}\begin{aligned} \Delta {\mbox{DMSPd}}/dt & = & R_{{\mbox{DMSPd}}} = R_{s:n}[ f_{sr}f_{res}\mu^{tot} +f_{nm}M_{ort} ] - [{\mbox{DMSPd}}]/\tau_{dmsp} \nonumber \\ \Delta {\mbox{DMS}}/dt & = & R_{{\mbox{DMS}}} = y_{dms}[{\mbox{DMSPd}}]/\tau_{dmsp} - [{\mbox{DMS}}]/\tau_{dms}\end{aligned}\end{split}\]

See BGC Tuning Parameters for a more complete list and description of biogeochemical parameters.

2.9.4.4. Algal growth and nutrient uptake

Nutrient limitation terms are defined in the simplest ecosystem for \({\mbox{NO$_3$}}\). If the appropriate tracer flags are true, then limitation terms may also be found for \({\mbox{NH$_4$}}\), \({\mbox{SiO$_3$}}\), and \({\mbox{fed}}\)

\[\begin{split}\begin{aligned} {\mbox{NO$_3$}}_{lim} & = & \frac{[{\mbox{NO$_3$}}]}{[{\mbox{NO$_3$}}] + K_{{\mbox{NO$_3$}}}} \nonumber \\ {\mbox{NH$_4$}}_{lim} & = & \frac{[{\mbox{NH$_4$}}]}{[{\mbox{NH$_4$}}] + K_{{\mbox{NH$_4$}}}}\nonumber \\ N_{lim} & = &\mbox{min}(1,{\mbox{NO$_3$}}_{lim} + {\mbox{NH$_4$}}_{lim}) \nonumber \\ {\mbox{SiO$_3$}}_{lim} & = & \frac{[{\mbox{SiO$_3$}}]}{[{\mbox{SiO$_3$}}] + K_{{\mbox{SiO$_3$}}}} \nonumber \\ {\mbox{fed}}_{lim} & = & \frac{[{\mbox{fed}}]}{[{\mbox{fed}}] + K_{{\mbox{fed}}}} \end{aligned}\end{split}\]

Light limitation \(L_{lim}\) is defined in the following way: \(I _{sw}(z)\) (in \(W/m^2\)) is the shortwave radiation at the ice level and the optical depth is proportional to the chlorophyll concentration, \(op_{dep} =\) chlabs [Chl*a*]. If ( \(op_{dep} > op_{min}\)) then

\[I_{avg} = I_{sw}(1- \exp(-op_{dep}))/op_{dep}\]

otherwise \(I_{avg} = I_{sw}\).

\[L_{lim} = (1 - \exp(-\alpha I_{avg}))\exp(-\beta I_{avg})\]

The maximal algal growth rate before limitation is

\[\begin{split}\begin{aligned} \mu_o & = & \mu_{max}\exp(\mu_T\Delta T)f_{sal}[{\mbox{N}}] \\ \mu' & = & min(L_{lim},N_{lim},{\mbox{SiO$_3$}}_{lim},{\mbox{fed}}_{lim}) \mu_o\end{aligned}\end{split}\]

where \(\mu'\) is the initial estimate of algal growth rate for a given algal species and \(\Delta T\) is the difference between the local tempurature and the maximum (in this case T\(_{max} = 0^o\)C).

The initial estimate of the uptake rate for silicate and iron is

\[\begin{split}\begin{aligned} \tilde{U}_{{\mbox{SiO$_3$}}} & = & R_{si:n}\mu' \\ \tilde{U}_{{\mbox{fed}}} & = & R_{fe:n}\mu' \end{aligned}\end{split}\]

For nitrogen uptake, we assume that ammonium is preferentially acquired by algae. To determine the nitrogen uptake needed for each algal growth rate of \(\mu\), first determine the “potential” uptake rate of ammonium:

\[U'_{{\mbox{NH$_4$}}} = {\mbox{NH$_4$}}_{lim}\mu_o\]

Then

\[\begin{split}\begin{aligned} \tilde{U}_{{\mbox{NH$_4$}}} & = & \min(\mu', U'_{{\mbox{NH$_4$}}}) \\ \tilde{U}_{{\mbox{NO$_3$}}} & = & \mu' - \tilde{U}_{{\mbox{NH$_4$}}}\end{aligned}\end{split}\]

We require that each rate not exceed a maximum loss rate \(l_{max}/dt\). This is particularly important when multiple species are present. In this case, the accumulated uptake rate for each nutrient is found and the fraction (\(fU^i\)) of uptake due to algal species \(i\) is saved. Then the total uptake rate is compared with the maximum loss condition. For example, the net uptake of nitrate when there are three algal species is

\[\tilde{U}^{tot}_{{\mbox{NO$_3$}}} = \sum^{3}_{i=1}\tilde{U}^i_{{\mbox{NO$_3$}}}\ \ \ .\]

Then the uptake fraction for species \(i\) and the adjusted total uptake is

\[\begin{split}\begin{aligned} fU^i_{{\mbox{NO$_3$}}} & = & \frac{\tilde{U}^i_{{\mbox{NO$_3$}}}}{\tilde{U}^{tot}_{{\mbox{NO$_3$}}}} \nonumber \\ U^{tot}_{{\mbox{NO$_3$}}} & = & \min(\tilde{U}^{tot}_{{\mbox{NO$_3$}}}, l_{max}[{\mbox{NO$_3$}}]/dt)\end{aligned}\end{split}\]

Now, for each algal species the nitrate uptake is

\[U^i_{{\mbox{NO$_3$}}} = fU^i_{{\mbox{NO$_3$}}} U^{tot}_{{\mbox{NO$_3$}}}\]

Similar expressions are found for all potentially limiting nutrients. Then the true growth rate for each algal species \(i\) is

\[\begin{aligned} \mu^i & = & \min(U^i_{{\mbox{SiO$_3$}}}/R_{si:n}, U^i_{{\mbox{NO$_3$}}} + U^i_{{\mbox{NH$_4$}}}, U^i_{{\mbox{fed}}}/R_{fe:n})\end{aligned}\]

Preferential ammonium uptake is assumed once again and the remaining nitrogen is taken from the nitrate pool.