2.4.1 Turbulent mixing

Vertical exchange due to turbulent mixing in the water column is prescribed externally⁴ by a turbulent diffusivity. In our case, it is taken from a three-dimensional MOM5 model run (Neumann et al., 2017). In this model set-up, turbulent vertical mixing is estimated by the KPP turbulence scheme (Large et al., 1994), which considers both local mixing and, in the case of unstable stratification, (non-local) convection. We only take into account the local part of the mixing and apply it to all tracers in the water column.

2.4.2 Particle sinking

In our model, suspended particulate matter sinks at a constant rate through the water column. We choose 4.5 m day⁻¹ for detritus, 1 m day⁻¹ for manganese and iron oxides, including the phosphate adsorbed by them, and 0.5 m day⁻¹ for large-cell phytoplankton and particulate organic carbon. In contrast, cyanobacteria are not sinking but, due to their positive buoyancy, they show an upward movement of 0.1 m day⁻¹. In reality, the sinking rate differs among individual particles; the currently chosen average values are a result of fitting the previous ERGOM model with the simplified sediment representation to observations.

2.4.3 Sedimentation and resuspension

Shear stress at the bottom determines whether erosion or sedimentation takes place. We apply the combined shear stress of currents and waves calculated by the same MOM5 model as the turbulent mixing. If this shear stress τ is below a critical value of τ_c=0.016 N m⁻² (Christiansen et al., 2002), the sinking suspended matter accumulates in the fluff layer compartment. If it is exceeded, the fluff layer material is resuspended into the lowest water cell at a constant relative rate r_ero=6 day⁻¹.

In our model, no material will ever be resuspended from the sediment itself, which starts below the fluff layer. This means that our model is incapable of realistically capturing extreme events like storms or bottom trawling which winnow the upper layers of the sediment, removing organic material, which has a lower sinking velocity, by separating it from the heavier mineral components (Bale and Morris, 1998). It also neglects a washout, which is the removal of organic matter from the sediment pores by advective transport of pore water by strong bottom currents (Rusch et al., 2001). In our model, sediment reworking by currents and waves is not explicitly represented, but rather parameterised together with the bioturbation process. This process allows for a bi-directional exchange of particulate material between the sediment and the fluff layer; see Sect. 2.4.5. The upward component of the transport represents winnowing of sediments (Bale and Morris, 1998).

2.4.4 Bioerosion

In environments with oxic bottom waters, we assume that in addition to waves and currents, macrofaunal animals or demersal fish can resuspend organic material from the fluff layer by active movements (Graf and Rosenberg, 1997). Therefore, under oxic conditions, we assume that r_biores=3 % day⁻¹ of the fluff material is resuspended independently from the shear stress conditions. This number was estimated from the calibration of a three-dimensional Baltic Sea ecosystem model (Neumann and Schernewski, 2008) in which the process proved to be critical for transporting organic matter to the deep basins below a depth of approx. 60 m. In these depths, a resuspension due to wave-induced shear stress is no longer possible.

2.4.5 Bioturbation

Bioturbation describes the movement and mixing of particles inside the sediment caused by the zoobenthos.⁵ In fact, it is difficult to discriminate what causes the vertical mixing of particles; physical effects like bottom shear may also have the same effect. We therefore include them in our “bioturbation” process.

We consider bioturbation to act as a vertical diffusivity D_B,solids(z) on the concentrations of the different solid species in the sediment. This implies that we exclude non-local mixing processes, even if they may be important in nature (Soetaert et al., 1996b), and try to represent them by local mixing. We only take intraphase mixing into account, which means we assume that the porosity Φ(z) remains constant over time.

The diffusivity D_B,solids(z) is also applied to describe the transport between the uppermost sediment layer and the fluff, which is caused by benthic organisms. In reality, the fluff layer may strongly differ in its compaction (porosity) depending on the turbulence conditions. However, we assume it to be perfectly compacted (ϕ=0) to be able to apply the above equation to describe the exchange process and therefore assume a thickness of 3 mm. This is not a physical assumption but rather a numerical trick which we use to transport the fluff material into the sediments. In reality, the fluff layer may be up to a few centimetres thick, and the incorporation of organic matter is done by macrofaunal activities (van de Bund et al., 2001).

The value 3 mm describes a volume estimate of SPM (suspended particulate matter) taken from this region: typical SPM concentrations in the lowermost 40 cm of the water column are about 8 mg L⁻¹ higher compared to the value 5 m above the sea floor (Christiansen et al., 2002). As the density of these particles is just slightly higher than that of the surrounding water, we can estimate their volume at approximately 3 L m⁻², which gives 3 mm of height if perfectly compacted. We see this explicit treatment of the fluff layer as a major advantage compared to the deposition of sinking particles directly into the surface sediments. We regard it as essential for the application of the model in a three-dimensional setting.

The vertical structure of bioturbation intensity, D_B,solids(z), is parameterised vertically as follows.

In the uppermost part of the sediment, we assume a constant bioturbation rate. Below that, it decays exponentially with depth until it reaches a maximum depth, which may be below the bottom of our model. So, we externally prescribe (a) the maximum mixing intensity⁶ and (b) three length scales describing the vertical structure of bioturbation⁷, which are the depth down to which the maximum mixing rate is applied (z_full), the length scale of the exponential decay of the mixing rate below this depth (z_decay) and the maximum depth of mixing (z_max).

The present formulation of the model has no explicit dependence of bioturbation depth on the availability of oxidants, i.e. bioturbation will take place in oxic as well as in sulfidic environments; adding this dependence should be essential if the model is applied to sulfidic areas.

2.4.6 Bioirrigation

Bioirrigation describes the mixing of solutes within the pore water and the exchange with the bottom water. We describe it as a mixing intensity D_B,liquids(z). The vertical profile of bioirrigation intensity is assumed identical to that of bioturbation. The maximum bioirrigation rate is assumed constant in time and prescribed externally⁸.

2.4.7 Molecular diffusion

Molecular diffusion in the sediment can be described by the equation

(Boudreau, 1997). Here, D₀ describes the molecular diffusivity in a particle-free solution, which is effectively reduced by the effect of hydrodynamic tortuosity θ. This describes the effect that the solutes need to travel a longer path as the direct way may be obstructed by solid particles. It is estimated from porosity by ${\mathit{\theta }}^{\mathrm{2}}=\mathrm{1}-\mathrm{2.02}\ln \left(\mathit{\varphi }\right)$ (Boudreau, 1997).

A diffusive exchange between the pore water and the overlying bottom water is controlled by the thickness of a diffusive boundary layer. While in reality this relates to the viscous sublayer thickness and is therefore inversely related to the velocity of the bottom water (Boudreau, 1997), for simplicity we assume a constant diffusive boundary layer thickness of 3 mm.

In reality, the diffusive boundary layer thickness is on the order of 1 mm at low-bottom-shear situations and becomes even shallower if the bottom shear increases (Gundersen and Jorgensen, 1990). We choose a larger value because we need to account for the transport through the fluff layer as well. A future model version might include a dependence of this parameter on the bottom shear stress.

Molecular diffusivities for the different solute species are calculated from water viscosity following Boudreau (1997). The water viscosity is determined from salinity and temperature (assumed to be identical to that in the bottom water cell). A problem occurs with the combined tracers DIC and total alkalinity, as they do not represent a specific ion but rather a set of different species with different molecular diffusivities. For simplicity, we approximate DIC diffusivity to be that of the ${\mathrm{HCO}}_{\mathrm{3}}^{-}$ ion, the most common one at the pH values we expect. For total alkalinity, we take a two-step approach: in the first step, we also take the diffusivity of the ${\mathrm{HCO}}_{\mathrm{3}}^{-}$ ion. But this is an underestimate, especially for the OH⁻ ions, which increase in concentration as the solution becomes alkaline. To take their higher diffusivity into account, we introduce two additional tracers, t_ohm_slowdiff and t_ohm_quickdiff. Before the molecular diffusion is applied during a model time step, they are both set equal to the OH⁻ concentrations. During the diffusion time step, the former diffuses with the reduced ${\mathrm{HCO}}_{\mathrm{3}}^{-}$ diffusion rate, the latter with the OH⁻ diffusivity. So afterwards, total alkalinity is corrected by adding the difference of the two, t_ohm_quickdiff-t_ohm_slowdiff. This results in a smoothed alkalinity profile.

2.4.8 Sediment accumulation

In nature, sediments grow upwards as new particulate matter is deposited onto them. In our model, this process is taken into account, but represented as the downward advection of material in the sediment. So, our coordinate system moves upward with the sediment surface. We assume that the sediment growth is supplied by terrigenous, bioinert material and prescribe⁹ a growth rate from the literature for the mud stations only (Table 7). We do not assume sediment growth for the sand and silt stations.

We use a simple Euler-forward advection to move the material from each grid cell into the cell below. Material leaving the model through the lower boundary is lost. Except for organic carbon, we assume that a part of it is mineralised, as will be explained in Sect. 2.7.1. In the top cell, new organic material from the fluff layer enters by sediment growth.

2.4.9 Parameterisation of lateral transport

The Baltic Sea sediments can be classified as accumulation, transport and erosion bottoms (Jonsson et al., 1990). The lateral transport of matter is characterised by the advection of fluff layer material from the transport and erosion bottoms in the shallower areas to the accumulation bottoms in the deep basins (Christiansen et al., 2002). As this process is not represented in our 1-D model set-ups, we need to parameterise it.

For the sandy and silty sediments, we assume transport away from the site. This is described by a constant removal rate for all material deposited in the fluff layer. For the mud stations, we assume transport of organic material towards the site. This is described by a constant input of detritus. Our model contains six detritus classes which degrade at different rates, as will be explained later in Sect. 2.6.4. We assume that the quickest-degradable part of the detritus is already mineralised in the shallow coastal areas before its lateral migration to the mud stations and therefore exclude the first two classes from this artificial input.

In the 3-D version of the model, these processes are no longer required, as the material is dynamically removed from the shallow sites and transported to deeper ones by advection.

2.5.1 Equations of motion

In this subsection, we will describe the equations of motion solved by the model. The equations in the water column can be derived from the assumption that the vertical (upward) flux of a tracer can be described by an advective and a diffusive flux, which follows Fick's law:

where c^wat(z,t) denotes the tracer concentration and D^wat is the turbulent diffusivity given as external forcing¹⁰. For particulate matter, the constant w describes its vertical velocity relative to the water, which is negative if the particles are sinking. For dissolved tracers, w is set to zero. We further assume that the water itself does not move vertically. In this case, conservation of mass yields an advection–diffusion equation:

where $q_{\mathrm{c}}^{\mathrm{wat}}(z,t)$ describes the biogeochemical sources minus sinks of the considered state variable.

The equations in the sediment are different because we need to take porosity into account and treat dissolved tracers (in the pore water) and solid tracers differently. For the pore water tracers, the upward flux is given by

where ϕ(z) is the porosity of the sediment (the ratio between pore water volume and total volume), which we assume as constant in time. The concentration c^pw(z,t) relates to the pore water volume only. The effective diffusivity D^pw is the sum of two contributions, the effective molecular diffusivity $\frac{D_{\mathrm{0}}}{{\mathit{\theta }}^{\mathrm{2}}}$ and the effective (bio)irrigation diffusivity D_B,liquids(z). The advection–diffusion equation is then given by

which is a well-known early diagenetic equation (Boudreau, 1997). For the solid-state tracers, their concentration c^sed(z,t) relates to the volume of the solids only, and the flux is given by

where w(z) is the velocity for virtual vertical downward transport. It results from sediment growth due to the deposition of particulate material, but as we keep the sediment–water interface at a constant position in our model, we need to describe the increasing depth in which we find individual sediment particles as downward advection. Volume conservation of the particulate material requires that we write w(z) as

such that the vertical velocity gets smaller in depths at which the sediment is more compacted, and w₀ describes a theoretical velocity which would occur at perfect compaction (ϕ=0)¹¹. The advection–diffusion equation then reads

Practically, we do not store the concentration c^sed(z,t) (mol m⁻³) as a state variable but rather the quantity of the tracer per area in a specific layer, C^sed(k,t) (mol m⁻²), where k is a vertical index. The transformation is straightforward:

For particulate tracers, we also consider storage in the fluff layer, C^fluff(t), which is measured in mol m⁻². The equation for C^fluff(t) is derived in the following subsection.

2.5.2 Boundary conditions

Boundary conditions are required for the partial differential equations given above. We give two boundary conditions for the water column concentrations: one at the sea surface, z_surf, and one at the sediment–water interface, z₀. We also give two boundary conditions for the sediment concentrations: one at the sediment–water interface, z₀, and one at the lower model boundary, z_bot. We start describing the boundary conditions from bottom to top for the dissolved tracers and then continue describing them from top to bottom for the particulate and solid-phase state variables.

The pore water tracers have a zero-flux boundary condition at the bottom of the model:

An exception to the zero-flux boundary is the parameterisation of sulfide production in the deep, which will be discussed later.

At the sediment–water interface, we assume that the dissolved tracers are exchanged between pore water and the water column via a diffusive boundary layer of a depth Δz_bbl. So, our upper boundary condition for the pore water tracers is given by

This flux can be directed into or out of the sediment, depending on where the concentration is larger.

To satisfy mass conservation, the vertical flux applied as the lower boundary condition for the dissolved species concentrations in the water column depends on the upward flux from the sediment:

The additional term ${\stackrel{\mathrm{̃}}{Q}}^{\mathrm{fluff}}\left(t\right)$ represents the sources minus sinks of the dissolved state variable, which are caused by biogeochemical transformations of the fluff layer material. At the sea surface, we apply a zero-flux condition, both for dissolved and particulate state variables:

An exception is only made for tracers which are modified by gas exchange with the atmosphere, e.g. oxygen.

Now the boundary conditions for the particulate state variables are different. The reason is that the water column and the sediment do not directly interact, but we consider the fluff layer as an intermediate layer between the two. Particulate material which sinks to the bottom is deposited in the fluff layer, from which it is incorporated into the sediments.

At the bottom of the water column, there can be two possible situations.

• If the bottom shear stress is lower than the critical shear stress, we assume a deposition of particulate material. This sinking material (w<0) vanishes from the water column because of sedimentation. It appears in the fluff layer.

• If the bottom shear stress exceeds the critical shear stress, particulate material from the fluff layer is eroded and enters the water column.

In both cases, we additionally consider the bioresuspension process which was described above in Sect. 2.4.4. We can therefore formulate the boundary condition for particulate material as

The fluff interacts with the surface sediment layer in two ways. Firstly, sediment growth means an incorporation of fluff layer material into the surface sediments. Secondly, bioturbation, which is considered diffusion–analogue mixing, leads to an exchange of particulate material between the fluff layer and surface sediment. So, the boundary condition for solids at the sediment surface is given by

Here, Δz_fluff represents a virtual thickness of the fluff layer assuming it was perfectly compacted; see the discussion in Sect. 2.4.5. In this way, the benthofaunal processes of incorporating fluff layer material into the surface sediments can be simply described as a diffusion–analogue flux of particulates. The opposite processes which cause removal of fine-grained material from the sediments, winnowing or washout, can be described in the same way as a diffusion process, in this case upward. This occurs in the model, especially when the fluff layer material is resuspended during periods of high bottom shear and the concentration C^fluff(t) is correspondingly low.

The concentration change in the fluff layer is then defined by mass conservation and is simply given by

for all particulate state variables. Here, $Q_{\mathrm{c}}^{\mathrm{pw}}\left(t\right)$ describes the sources minus sinks term from the biogeochemical transformations of the considered state variable.

Table 2Reaction network table for the processes in the water column. See Table A1 for the composition of state variables. Processes marked with * also take place in the pore water.

Download Print Version | Download XLSX

Finally, the burial of particulate material at the lower model boundary can be described by the following boundary condition:

So, we assume the particulate material to be buried forever when it leaves the model domain. An exception, as mentioned before, is the parameterisation of deep sulfide formation, which is described in Sect. 2.7.

2.6.1 Primary production and phytoplankton growth

There are three classes of phytoplankton in the model, representing large-cell and small-cell microalgae as well as diazotroph cyanobacteria. Their growth is determined by a class-specific maximum growth rate, but contains two limiting factors for nutrients and light. The light limitation is a saturation function with optimal growth at a class-specific optimum level or at 50 % of the surface radiation. The shortwave light flux at the surface is taken from a dynamically downscaled ERA40 atmospheric forcing (Uppala et al., 2005) using the regional Rossby Centre Atmosphere model (RCA). Nutrient limitation is a quadratic Michaelis–Menten term for DIN (nitrate + ammonium) or phosphate, depending on which one is limiting, based on Redfield stoichiometry. Diazotroph cyanobacteria are only limited by phosphate and not by DIN, but they are only allowed to grow in a specific salinity range. Cyanobacteria and small-cell algae also require a minimum temperature to grow (Andersson et al., 1994; Wasmund, 1997).

However, according to Engel (2002), although nutrients are limiting an enhanced polysaccharide exudation could be the result of a cellular carbon overflow whenever nutrient acquisition limits biomass production but not photosynthesis. These transparent exopolymers are included in our model, and they are assumed to have a constant sinking velocity.

2.6.2 Phytoplankton respiration and mortality

We assume a constant respiration of phytoplankton which is proportional to its biomass. As the model maintains the Redfield ratio, the degradation of biomass (catabolism) goes along with the excretion of ammonium and phosphate. This simplified description of phytoplankton growth does not describe day–night metabolism or temperature dependence. A small fraction of the nitrogen is released as dissolved organic nitrogen (DON). In the model, this represents the DON fraction which is less utilisable by phytoplankton, while the fraction with high bioavailability is considered to be part of the ammonium state variable.

Due to simplification, in our model phytoplankton experiences a constant background mortality, although we know this is far away from reality in which it is species specific and depends on abiotic (e.g. nutrient, light, etc.) and biotic conditions. An additional mortality is generated by the grazing of zooplankton as described next.

2.6.3 Zooplankton processes

Zooplankton is only represented as one bulk state variable. It grows by assimilating any type of phytoplankton; however, it has a smaller food preference for the cyanobacteria class compared to the other classes. The uptake becomes limited by a Michaelis–Menten function if the zooplankton's food approaches a saturation concentration. Feeding can only take place in oxic waters and is temperature dependent. It shows a maximum at an optimum temperature and a double-exponential decrease when this temperature is exceeded.

Both zooplankton respiration and mortality represent a closure term for the model. They are meant to include the respiration and mortality of the higher trophic levels (fish) which feed on zooplankton, and therefore we use a quadratic closure. Mortality is additionally enhanced under anoxic conditions, which do not occur in our study area.

2.6.4 Mineralisation processes

The description of detritus¹² differs from the previous ERGOM versions. We have split the detritus into six classes, depending on its degradability. This degradability is described as a decay rate constant, which ranges from 0.065 day⁻¹ for the first class to $\mathrm{1.6}\times {\mathrm{10}}^{-\mathrm{5}}$ day⁻¹ for the fifth class, while the last one is assumed to be completely bioinert. This type of model is known as a “multi-G model” (Westrich and Berner, 1984).

Details on the specific choice of the classes are given in Appendix B.

The mineralisation is, however, temperature dependent by a Q₁₀ rule (Sawicka et al., 2012; Thamdrup et al., 1998), as it is realised by microbial processes; the values given above are valid at 0 ^∘C. The 0 ^∘C choice is somewhat arbitrary. Actually, the model is not very sensitive to this choice, as an enhanced baseline temperature, meaning a lower decomposition rate of each class, would be compensated for by a shift in the class composition, leaving higher concentrations of quickly degradable detritus classes, which overall means a very similar total decomposition rate; see Appendix B.

When organic detritus is created by plankton mortality, it is partitioned into the different classes in a constant ratio. This ratio was determined from a fit of the multi-G model to an empirical relation between detritus age and its relative decay rate, which was proposed by Middelburg (1989). The fraction of non-decaying detritus was estimated from empirically determined carbon burial rates in the Baltic Sea (Leipe et al., 2011).

The chemical composition of detritus is, in contrast to phytoplankton and zooplankton, not determined by the Redfield ratio. It is enriched in carbon and phosphorus by 50 % such that it has a C : N : P ratio of 159 : 16 : 1.5. This resembles detritus compositions as they were determined in sediment traps and by investigating fluffy layer material in the Baltic Sea (Emeis et al., 2000, 2002; Heiskanen and Leppänen, 1995; Struck et al., 2004).

In the water column, detritus can be mineralised by three different oxidants: oxygen, nitrate and sulfate. They are utilised in this order; if the preferential oxidant's concentration declines, the specific pathway is reduced by a Michaelis–Menten limiter and the next pathway takes over such that the total mineralisation is held constant. In all pathways, DIC, ammonium and phosphate are released. Nitrate reduction also produces molecular nitrogen (heterotrophic denitrification), while sulfate reduction generates hydrogen sulfide.

Mineralisation of particulate organic carbon in transparent exopolymers takes place via the same pathways, but only releases DIC. DON is also mineralised after some time and decays to ammonium (which may represent the transformation to bioavailable DON compounds).

2.6.5 Reoxidation of reduced substances

In the presence of oxygen, ammonium is nitrified to nitrate (Guisasola et al., 2005). The intermediate step, the formation of nitrite, is omitted in the model. Hydrogen sulfide can be reoxidised by oxygen or by nitrate (chemolithoautotrophic denitrification) (Bruckner et al., 2013). This takes place as a two-step process via the formation of elemental sulfur (Jørgensen, 2006). All reoxidation processes exponentially increase their rates with temperature.

In the sediments, we additionally assume that Fe²⁺ can be produced as a reduced substance. If it is released from the sediments and enters the water column, it can be reoxidised by oxygen, creating suspended iron oxyhydroxides.

2.6.6 Adsorption and desorption reactions

Dissolved phosphate can be adsorbed to iron oxyhydroxide particles suspended in the water column. In the same way, phosphate adsorbed to iron oxyhydroxide particles can be released if the ambient concentration of phosphate is low. The process is identical to the one in the sediments and is discussed in Sect. 2.7.5 in detail.

2.7.1 Mineralisation in general

The mineralisation of detritus is the dominant biogeochemical process in the sediments, as the oxidation of the carbon therein is the major supply of chemical energy for microbes.

As in the water column, oxidants are utilised in a specific order, and a smooth transition to the next mineralisation pathway occurs when the preferred one gets exhausted. However, the number of possible oxidants is increased in the sediment, as here solid components may also act as electron acceptors. The order in which they are utilised is (Boudreau, 1997)

1. oxygen,

2. nitrate,

3. manganese oxide,

4. iron oxyhydroxide,

5. Fe (III) contained in clay minerals and

6. sulfate.

After sulfate is exhausted, typically the formation of methane would start. This process is omitted in the current model, as we designed our model for the top 22 cm of the south-western part of the Baltic Sea, where we do not expect sulfate to be limiting. This depth restriction is based on the limited length of the sediment cores taken in the empirical part of our research project. We do, however, describe the process implicitly, since we assume that a part of the organic carbon which leaves the model domain through the lower boundary will be transformed to methane, which as it diffuses upward will be oxidised by sulfate and generate H₂S. Therefore, we parameterise this process by a conversion from sulfate to hydrogen sulfide at the lower boundary.

As in the water column, we distinguish six different classes of detritus with different basic mineralisation rates.

Details on the specific choice of the classes are given in Appendix B.

These rates are only controlled by temperature, not by the specific oxidant which is available. There is an ongoing controversy as to what determines the rate of sedimentary carbon decay and whether it is the oxidant (and therefore the accessible energy per mole of carbon) or the degradability of the detrital carbon itself (Arndt et al., 2013; Kristensen et al., 1995). In leaving out the explicit dependence of the oxidant, we do not favour the latter theory; we chose to adopt the decay rates proposed by Middelburg (1989), which may implicitly take the effect of the oxidant into account¹³.

Sedimentary organic phosphorus (OP) may degrade faster than the corresponding nitrate and carbon, an effect known as preferential P mineralisation (Ingall and Jahnke, 1997). We include this by introducing additional state variables t_detp_n for each class n of detritus, describing the OP concentration, as well as a constant factor pref_remin_p, which describes a redox-dependent ratio between the mineralisation speeds of OP and organic carbon and nitrogen. This factor is set equal to 1 under oxic conditions and greater than 1 under anoxic conditions (Jilbert et al., 2011). This approach follows Reed et al. (2011).

2.7.2 Specific mineralisation processes

Here, we describe the implementation of the primary redox reactions, indicated by the red numbers in Fig. 3.

Oxic mineralisation and heterotrophic denitrification are formulated in the same way as in the water column; see Sect. 2.6.4.

The next pathway is the reduction of Mn (IV) to Mn (II), which produces dissolved manganese.

The reduction of iron oxyhydroxides should produce dissolved Fe (II). This, however, may precipitate very quickly, especially where hydrogen sulfide is present. So for numerical reasons, we combine these reactions, and the reduced Fe (III) is directly converted into iron monosulfide or considered as adsorbed by clay minerals, as we describe below in Sect. 2.7.3.

Some clay minerals, especially sheet silicates which are abundant in the German part of the Baltic Sea (Belmans et al., 1993), contain structural iron which is available for redox reactions (Jaisi et al., 2007). We prescribe a station-specific content of these minerals given in Table 7 and assume that they contain a small amount (0.1 mass-%) of reducible iron; this is because a particle analysis of sheet silicates from the area of interest (Leipe, unpublished data) showed slightly lower iron contents in the sulfidic zone compared to the surface area.

The primary redox reaction follows process 32 in Table 3; we describe it in detail since it is a new process added to our model. Mineralisation of organic carbon under the reduction of structural iron in sheet silicates takes place at a rate of

where t_sed_k is the amount of detritus per area in a specific sediment layer, r_k is a basic reactivity for this class (see Appendix B), τ=0.15 K⁻¹ is a temperature sensitivity constant for the mineralisation and T is the temperature in ^∘C. The limitation functions l_? are of Ivlev type, e.g.

they have a value close to 1 at high concentrations of the corresponding oxidant and become zero if this oxidant is exhausted. Here, O_2,min denotes a threshold concentration at which the limitation occurs, but the Ivlev function generates a soft transition between the different oxidation pathways, so they can occur simultaneously. The product of the last five factors in Eq. (19) means this process will run at a substantial rate only if

• oxygen concentration is low,

• nitrate concentration is low,

• manganese oxides in the sediment are depleted and

• iron oxyhydroxides in the sediment are depleted, but

• reducible structural Fe (III) in the clay minerals is still abundant.

Sulfate reduction produces hydrogen sulfide. As discussed above, it represents the terminal mineralisation process in our model. This process, described by processes 33 and 34 in Table 3, follows the formulation by Reed et al. (2011) in our model. Organic material leaving the lower boundary of our model because of sediment growth will also be mineralised by this process. We assume that a fraction of the buried material will be mineralised by either sulfate reduction or methanogenesis, the rest being buried. For the methane produced, we assume that it will not enter the model domain but rather be oxidised by sulfate, producing H₂S below the model domain. We assume for simplicity that all these reactions happen instantaneously, which results in the same net reaction as the sulfate reduction.

2.7.3 Precipitation and dissolution reactions

Solids can precipitate from a solution when it becomes supersaturated. This happens in an aqueous solution when the actual ion activity product exceeds the respective solubility product and a critical degree of supersaturation is reached (Böttcher and Dietzel, 2010; Sunagawa, 1994).

Diagenetic models often simplify the calculation by multiplying the concentrations rather than the activities (van Cappellen and Wang, 1996). The resulting product is then proportional to the actual solubility product as long as the ionic strength of the solution does not change. As the ionic strength of seawater is almost completely defined by its salinity (Millero and Leung, 1976), this assumption is well justified for most marine environments. The Baltic Sea, however, is a brackish sea with strong spatial and temporal changes in bottom salinity, especially in the western part (Leppäranta and Myrberg, 2009). For this reason, we take the activity coefficients, which transform concentrations to activities, into account. This is done by using the Davies equation (Davies, 1938), which determines the individual activity coefficient a_i as (Stumm and Morgan, 2012)

where I is the ionic strength expressed in mol L⁻¹, z_i is the ion charge and A is the Davies parameter calculated after Kalka (2018):

Here, ϵ is the dielectric constant of water calculated after Gadani et al. (2012)¹⁴, and T_K is the temperature in K.

Ca–rhodochrosite precipitates at elevated concentrations of manganese and carbonate. Its solubility product is composition dependent, as the Ca : Mn ratio varies (Böttcher, 1997; Böttcher and Dietzel, 2010). For Baltic Sea muds in which ratios around 0.6 occur, an effective solubility product (including the effect of oversaturation) of 10^−9.5 to 10⁻⁹ M² can be deduced from Jakobsen and Postma (1989). In our model, the reaction follows process 53 in Table 5. 10 % of the dissolved manganese will precipitate per day if the solution is undersaturated. Saturation is calculated by the formula (Jakobsen and Postma, 1989)

where $\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]$ is the concentration of carbonate ions, which depends on DIC concentration and pH; see the description of the carbon cycle in Sect. 2.8. The term a₂ is the activity coefficient for ions with a charge of 2; see Eq. (21). The term in the denominator is the solubility product for rhodochrosite (Jakobsen and Postma, 1989).

If the solution becomes undersaturated, rhodochrosite will be dissolved again. Then, process 53 is reversed, and a fixed amount of 10⁻⁶ mol kg⁻¹ day⁻¹ of manganese (II) is released until saturation is reached.

Iron monosulfide precipitates on contact with dissolved Fe (II) and sulfide, depending on pH, with a solubility product taken from Morse et al. (1987). But, as stated in Sect. 2.7.2, we assume for numerical reasons that this process takes place directly after Fe (III) reduction. The solubility product is then used in an inverse way to determine the equilibrium concentration of dissolved Fe (II) at the current pH, sulfide concentration and salinity:

where a_i represents the Davies activity coefficients (see Eq. 21), and 10^−2.95 mol L⁻¹ is the solubility product for iron monosulfide (Morse et al., 1987; Theberge and Iii, 1997).

We then assume a precipitation or dissolution of iron monosulfide, which relaxes the present concentration of Fe (II) against this equilibrium. This is in agreement with a pore water chemistry model for the central Baltic Sea (Kulik et al., 2000), which states that dissolved iron concentrations in the pore water are buffered by iron sulfides (mackinawite and greigite). The dissolution of iron monosulfide also takes place if clay minerals in the same grid cell are capable of adsorbing additional Fe (II). This process is described in Sect. 2.7.5.

As a simplification, we neglect the change in porosity which would be caused by precipitation (or dissolution) of any solids.

2.7.4 Pyrite formation

Pyrite (FeS₂) is a crystalline compound formed in early diagenesis (Rickard and Luther, 2007). Its formation from iron monosulfide is included in most early diagenetic models. This process is not a simple precipitation process, but rather a redox process. While both sulfide and iron monosulfide contain sulfur of oxidation state −2, the redox state of S in pyrite is −1. This implies that an electron acceptor is required to create pyrite. A generally accepted mechanism for pyrite creation is the use of zero-valent sulfur from polysulfides; this may be created by the oxidation of sulfate with Fe (III). However, this process alone cannot explain the high degrees of pyritisation in Baltic deep sediments (Boesen and Postma, 1988).

An additional pathway which does not rely on elemental sulfur, but instead reduces hydrogen sulfide to hydrogen gas, has been proposed by Drobner et al. (1990), Rickard and Luther (1997), and Rickard (1997). Similar to how it was done in early diagenetic models (Wijsman et al., 2002), we include this pathway and therefore assume that whenever iron monosulfide and H₂S are present, pyrite is formed from them. The generated H₂ will be consumed by sulfate-reducing bacteria (Stephenson and Stickland, 1931), so in the net reaction, sulfate acts as the electron acceptor.

In our model, the reaction therefore follows process 39 in Table 4. The speed of the transformation process is determined by

where 𝚝_𝚒𝚖𝚜 is the concentration of iron monosulfide, 𝚝_𝚑𝟸𝚜 is the concentration of hydrogen sulfide and k_pyr is a kinetic constant for conversion of iron monosulfide to pyrite. Its value of 8.9 kg mol⁻¹ day⁻¹ was adopted from Wijsman et al. (2002) and Rickard and Luther (1997), which use 8.9 L mol⁻¹ day⁻¹.

2.7.5 Adsorption balances

Adsorption in our model takes place on the surfaces of two particle types: iron oxyhydroxides and clay minerals. Adsorption on silicate particles is not explicitly represented in the model, but parameterised by a reduction of the effective diffusivity of phosphate and ammonium, following Boudreau (1997).

Iron oxyhydroxides adsorb dissolved phosphate. This is a well-known process responsible for the sedimentary retention of phosphate derived from mineralisation processes (Sundby et al., 1992). As both phosphate and hydroxide ions can occupy the adsorption sites at the surface, adsorption is less efficient in alkaline environments. In our model, we use a formula from Lijklema (1980) which describes the adsorbed P : Fe ratio at a given phosphate concentration and pH.

Here, DIP gives the dissolved phosphate concentration. But we use it inversely. We calculate an equilibrium concentration for dissolved phosphate at the current P : Fe ratio and pH.

If the current concentration of dissolved phosphate is above this equilibrium concentration, adsorption takes place and PO₄ in the pore water is decreased. If it is below the equilibrium concentration, desorption takes place. The maximum function is added to treat situations when both pH and the amount of adsorbed phosphate get so low that the formula by Lijklema (1980) gives no real solution for DIP. In this case, we assume that all currently dissolved phosphate will become adsorbed. The model processes p_po4_ads_ips and p_ips_diss_po4 will change the phosphate concentration in the pore water by

where k_ips_dissolution = 0.1 day⁻¹ is a reaction rate constant we assume. We chose this probably unrealistically low value for reasons of numerical stability.

Following Reed et al. (2011), we define two classes of iron oxyhydroxides. The first one is fresh, amorphous and adsorbs phosphate. The second one is a more crystalline phase, for which we assume no adsorption. The first phase is transformed to the second one with a constant rate in time, implying continuous phosphate release.

Clay minerals, due to their large surface area, can also adsorb pore water species. We allow for an adsorption of phosphate, ammonium and Fe (II). For simplicity, we assume that the ratio of adsorbed species to clay mass is proportional to the pore water concentration until a saturation threshold is exceeded. For Fe (II), this proportionality constant is derived from Jaisi et al. (2007), for ammonium from Raaphorst and Malschaert (1996), and for phosphate from Edzwald et al. (1976). In all three cases, we calculate a pore water concentration which is in equilibrium with the current ratio of adsorbed species to clay mineral mass. Then adsorption or release processes take place to relax the present pore water concentration towards the equilibrium value.

To calculate ${{\mathrm{Fe}}^{\mathrm{2}+}}_{\mathrm{eq},\mathrm{clay}}$, we first determine the mass of clay minerals present in a specific model layer per square metre. This is done by the formula

Here, ${\mathit{\rho }}_{\mathrm{clay}}=\mathrm{2.7}\times {\mathrm{10}}^{\mathrm{3}}$ kg m⁻³ gives the density of montmorillonite (Osipov, 2012), Δz (m) gives the thickness of the layer, (1−ϕ) gives the ratio between the volume of the solids and the total volume of the sediments, and r_clay is the volume fraction of clay minerals among the total volume of solids. So, m_clay has a unit of kg m⁻². In the next step, we find out how much iron gets adsorbed to clay depending on the dissolved concentration. For Fe (II) concentrations much smaller than 1 mmol L⁻¹ as we observe in our sediments, we can linearise the adsorption isotherms given by Jaisi et al. (2007) and obtain

where $q_{max}^{\mathrm{mass}}$ is a mass-specific sorption capacity (mol kg⁻¹), α is a binding energy constant (L mol⁻¹) and [Fe²⁺] is the concentration of iron in the pore water (mol L⁻¹). We can rearrange the equation to obtain the equilibrium concentration of dissolved Fe (II):

For the product $\mathit{\alpha }\cdot q_{max}^{\mathrm{mass}}$, Jaisi et al. (2007) find values between 500 and 3000 L kg⁻¹ for different types of clay, which means 1 kg of clay added to 0.5 to 3 m³ of water would adsorb the same amount of Fe²⁺ as would remain in the solution. We adopt a value of 1000 L kg⁻¹ for our model.

For numerical reasons, we allow for an immediate precipitation of the desorbed Fe (II) as iron monosulfide in the case of oversaturation, leaving out the intermediate transformation to dissolved Fe (II). The inverse is also true: if iron monosulfide is dissolved, the released Fe (II) may directly be adsorbed by the clay minerals instead of being released to the pore water first.

This is described by process 52 in Table 5.

For the adsorption isotherm of phosphate on clay minerals, we follow the study by Edzwald et al. (1976). They give maximum adsorption capacities $m_{\mathrm{P},\mathrm{ads},\max }/m_{\mathrm{clay}}$ between 0.09 mg g⁻¹ for P on kaolinite and 2.58 mg g⁻¹ for illite. These values were obtained at a pH close to 7.5, and the pH dependence of adsorption differs among the different clay minerals. Since the composition of the clay minerals is unknown to us, we choose a conservative value of 0.2 mg g⁻¹; this could be adapted when such data are available. In contrast to Fe (II) adsorption, a half-saturation of P adsorption is already reached at concentrations around 1 mg L⁻¹, which corresponds to approx. 0.03 mmol L⁻¹. We model this saturation in a very simple way by a linear dependency of dissolved and adsorbed phosphate below a threshold concentration of dissolved phosphate and a constant amount of adsorbed phosphate if the threshold is exceeded:

where M_P=31 g mol⁻¹ is the molar mass of phosphate, m_clay is the mass of clay per square metre in the given grid cell (see Eq. 29), $m_{\mathrm{P},\mathrm{ads},\max }/m_{\mathrm{clay}}=\mathrm{0.2}$ mg g⁻¹ is the assumed maximum adsorption capacity and ${\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-,\mathrm{sat}}=\mathrm{0.03}$ mmol L⁻¹ is the concentration of dissolved phosphate at which we assume this saturation is reached. We then define an adsorption and a desorption reaction following process 49 in Table 5. The adsorption process is assumed to happen instantaneously, but for numerical reasons we limit the process rate by demanding that at maximum (a) 10 % of the dissolved phosphate is removed per day or (b) 10 % of the lack of adsorbed phosphate with reference to the equilibrium concentration is precipitated. This artificial deceleration of the precipitation process had to be included to avoid numerical difficulties. The desorption process works in a similar way. At maximum, 10 % of the adsorbed phosphate which exceeds the equilibrium concentration is released per day or 10 % of the saturation concentration ${\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-,\mathrm{sat}}$, whichever is less.

For ammonium adsorption to clay minerals, the processes are in principle identical to those of phosphate. Since the adsorption is weak compared to that of phosphorus (Raaphorst and Malschaert, 1996), we neglect the effect by setting the maximum amount of adsorbable ammonium to zero in our present set-up. So while the model is able to include the dynamics of ammonium adsorption to clay minerals, we make no use of it in the present application.

2.7.6 Reoxidation of reduced substances

Reduced substances can be reoxidised if the appropriate oxidant is present in a sufficient concentration. Table 6 gives a summary of the redox reactions implemented in our model, which will be described one by one in this section.

Ammonium is oxidised to nitrate in the presence of oxygen. The rate of this process is proportional to both the ammonium and the oxygen concentration and, as in the water column, increases exponentially with temperature.

Dissolved manganese (II) will be oxidised in the presence of oxygen and precipitates as manganese oxide. This is also assumed to be a second-order process proportional to both precursor concentrations.

Dissolved Fe (II) is oxidised by oxygen in a pH-dependent way. The rate of this process is proportional to the Fe (II) and oxygen concentration, as well as to the square of the hydronium ion concentration. It is also influenced by temperature and ionic strength, as described by Millero et al. (1987). For numerical reasons, we also allow for the direct oxidation of Fe (II) adsorbed to clay minerals. Alternatively, dissolved Fe (II) can be oxidised by reducing manganese. This process follows Reed et al. (2011). The generated Fe (III) immediately precipitates as iron oxyhydroxide.

Structural iron in clay minerals can be reoxidised as well. We only allow this process in the presence of oxygen, when it transforms back to Fe (III), which is kept bound in the clay minerals.

This reaction follows process 43 in Table 4. The process runs at the speed of

The oxidation occurs only until the maximum amount of reducible Fe (III) in the clay material, i3i_max, is reached. It occurs at a rate of r_i2i_ox = 0.1 day⁻¹, a somewhat arbitrary value indicating that the process is typically fast compared to the vertical transport of clay minerals. It is Ivlev-limited by a factor $l_{\mathrm{o}\mathrm{2}}=\mathrm{1}-\exp \left(\frac{-\left[{\mathrm{O}}_{\mathrm{2}}\right]}{{\mathrm{O}}_{\mathrm{2},\min }}\right)$ with ${\mathrm{O}}_{\mathrm{2},\min }=\mathrm{2.0}\times {\mathrm{10}}^{-\mathrm{5}}$ mol kg⁻¹, consistent with the concentration at which carbon oxidation becomes limited.

Table 6Reaction network of secondary redox reactions in the sediment, giving the possible reoxidation processes in the presence of the oxidants listed in the first row.

Download Print Version | Download XLSX

Hydrogen sulfide can reduce any of the previously mentioned oxidants being converted to sulfate. The reaction with oxygen or nitrate is carried out as a two-step reaction. The intermediate species formed in these reactions is elemental sulfur, which can be further oxidised to sulfate. These processes follow the same kinetics as in the water column; see Sect. 2.6.5. Hydrogen sulfide can alternatively react with manganese oxides or iron oxyhydroxides, producing dissolved Mn (II) or Fe (II). For the generated Fe (II), however, we assume either an immediate precipitation to iron monosulfide or an immediate absorption to clay minerals, whichever is more favourable. We assume these reactions to be proportional to both the concentration of sulfide and the metal oxides. Hydrogen sulfide can also reduce structural Fe (III) in the clay minerals, and the Fe (II) will in this case remain in the clay.

This reaction follows process 47 in Table 4. The process runs at a speed of

The model parameter r_i2i_ox describes a relative speed of 0.1 day⁻¹ at which H₂S is reoxidised by this process, a somewhat arbitrary value just expressing our assumption that the process is fast compared to vertical transport of the clay minerals. The model shows low sensitivity to this rate parameter, as shown in the Supplement. The process is Ivlev-limited by the factor l_i3i.

Iron monosulfide is typically not directly oxidised but dissolves at low sulfide concentrations. However, if it is exposed to oxygen, we assume complete oxidation to Fe (III) and sulfate.

Finally, pyrite can be oxidised in the presence of oxygen or manganese (IV), but in marine environments not by Fe (III) (Schippers and Jørgensen, 2002). We assume complete oxidation to sulfate and iron oxyhydroxides.

2.9.1 Numerics of physical processes

Vertical diffusion is done explicitly by multiplying each vertical tracer vector by a diffusion matrix. This includes turbulent mixing in the water column as well as pore water diffusion, bioturbation (faunal solid transport) and bioirrigation (faunal solute transport). This diffusion matrix is tridiagonal, and for a small time step, which is in our case limited by the thin layers at the top of the sediment, a Euler-forward method can be applied. Larger time steps could be split into smaller Euler-forward steps, which means a repeated multiplication by the tridiagonal matrix. We instead use an efficient algorithm to calculate powers of the tridiagonal matrix (Al-Hassan, 2012) and perform the multiplication only once.

2.9.2 Numerics of biogeochemical processes

The sources and sinks for the different tracers are calculated from the process rates. These include not only biogeochemistry, but also parameterisations for lateral transport processes as well as sedimentation and resuspension.

To calculate the changes in a tracer concentration with time, we form the sum of the processes consuming or producing it (Radtke and Burchard, 2015):

Here, T_i represents the concentration of tracer i, p_k is the rate at which process k runs and q_ik (and s_ik) is the stoichiometric ratio in which process k produces (or consumes) tracer i. In order to ensure both non-negativity of the tracer concentrations and mass conservation, we apply the positive Euler-forward method from Radtke and Burchard (2015). It is a clipping method which, in the case that a tracer concentration becomes negative during one Euler-forward time step, first executes a partial time step until this tracer is zero. Then the rest of the time step is continued without the processes consuming this tracer, i.e. they are switched off. More than two partial time steps may be needed if more than one tracer is exhausted.

4.3.1 Penalty function

The first step in such an optimisation is to define a metric or a penalty function quantifying the misfit between model and observations. Our aim is then to minimise this function.

We chose to penalise the relative deviation between model and measurement and define the penalty function by

Here, i_max is the number of state variables we compare and j_max(i) is the number of observations of this state variable at all depths and at all stations. The expression r_i,j is a measure for the relative deviation between model value m_i,j and observation o_i,j:

The term Δ_i¹⁵ is included to avoid huge relative errors between values which are close to zero. It denotes the random deviation in this parameter, quantified from duplicate or triplicate measurements of the same parameter in different sediment cores of the same sampling. Obviously r_i,j becomes zero if the model and observations match. Finally, the weight w_i assigned to each comparison in Eq. (48) is defined as the ratio

between the average observed value of the variable and the random deviation. The weight is applied to make sure that fitting the most certain variables has the highest priority.

The pore water species fitted are ammonium, phosphate, silicate, sulfide, iron, manganese, the total alkalinity and the relative sulfate deficit¹⁶.

4.3.2 Optimisation strategies

After we defined a penalty function, the second step is to choose an algorithm to minimise it. Several such algorithms exist; however, our choice of methods was restricted by the relatively long runtime of a single model iteration. Since it took about 8 min to run a single station for 100 years, we had to choose methods which

1. needed a relatively small number of iteration steps and therefore

2. allowed for a high degree of parallelism in the individual optimisation step in order to effectively search the 115-dimensional parameter space.

Our first choice was the Adaptive Hierarchical Recombination–Evolutionary Strategies (AHR-ES) algorithm implemented in the R package calibraR (Oliveros-Ramos and Shin, 2016). We were, however, not satisfied with the optimisation result. Possibly we just failed to find out the optimal settings for the algorithm, such as the survival rate.

Therefore, our second choice was a simple alternative algorithm: our own extension of the generalised pattern search (GPS) algorithm (Hooke and Jeeves, 1961). Every optimisation step consists of two sub-steps. The first sub-step is the most simple “grid search” step in which all 115 parameters are varied by a predefined step width. We run 230 sets of seven models in parallel such that each parameter can be both increased and decreased. In the second sub-step, 230 combinations of the most successful changes are formed. Even if no single parameter change could improve the existing solution, sometimes combinations of them can, which is the basic idea of GPS. The parameter vector with the best score which was obtained in any of the two sub-steps is then chosen as the starting point for the next step. If neither of the two sub-steps leads to an improvement of the overall fit, the step size is reduced by a factor of 2.

The optimisation converged after 30 iteration steps and reduced the error function from 6363 (the value obtained by previous manual tuning) to 4797.

The algorithm obviously does not guarantee that we reach a global optimum, which can be seen as a drawback. The automatic method was started after manual calibration of the model. Since the optimisation method is deterministic, the local optimum is defined by this initial condition. However, in a vector space with a dimension as high as ours, it is anyway difficult to find a non-local point with a better score, no matter if it is by manual optimisation or a different search algorithm.

5.1.1 Pore water profiles at mud stations

Figure 4 shows a comparison of simulated and measured pore water profiles at the three mud stations.

Figure 4Pore water concentrations of several dissolved species at the three mud stations Lübeck Bight (a), Mecklenburg Bight (b) and Arkona Basin (c). Points and horizontal lines indicate the range of measurements. For station AB, open triangles indicate the range of observations from the January 2015 cruise, which were taken at a different location in the same area 23 km apart. Curves and shading present the model results and indicate year-average concentrations and the seasonal range. Please note the different horizontal scales.

Download

In the left panels, we see that the rise of alkalinity with depth is captured well by the model, except for the AB site where observations show a higher alkalinity below 10 cm of depth. The decline in sulfate follows the lower range of the observations.

The panels in the second column show that the vertical profiles of ammonium and silicate are also represented relatively well by the model. However, especially at the Arkona Basin station, the observed range of both ammonium and sulfide shows strong variation (by an order of magnitude). The model does not capture that but rather sticks to the lower range of the observations. Most probably, the variability in the observations is not due to seasonality, but a consequence of spatial variability between sampling sites, since the samples were taken from two sites 23 km apart.

Surprisingly, the model is able to reflect the differently steep sulfide profiles between the stations LB and MB. While Lipka (2018) see the low sulfide concentrations which occur especially in March 2014 at MB as an indication for a preceding mixing event, our model cannot adopt this interpretation due to its limitation by a temporally constant vertical mixing. In contrast, our model suggests a higher deposition of iron oxyhydroxides at this site.

The right panels show that the modelled manganese concentrations match the observations quite well. The dissolved iron profiles show their maxima at the correct depths and a relatively large seasonal spread. The measurements show an even larger spread than the model. For the phosphate profiles, the model results mostly resemble the lowest of the measured values, except for station AB where we see a clear underestimation. This can be seen as an artefact of our fitting method, more precisely of the choice of our penalty function. Giving a penalty for the relative error means that the same absolute error is punished more heavily if the observation is smaller, making the model try harder to fit low values compared to high ones.

The model results for the mud stations fit quite well, considering the fact that the real pore water profiles may be shaped by very different temporal variations. These include, for example, mixing events, changing loads of organic matter, or temperature and salinity variations. Our model, not knowing the sediments' past, can only try to estimate the average conditions that might produce similar pore water concentrations.

5.1.2 Pore water profiles at sand stations

Figure 5 illustrates the model fit at the sandy stations.

Figure 5Pore water concentrations of several dissolved species at the three sand stations Stoltera (a), Darss Sill (b) and Oder Bank (c). Points and horizontal lines indicate the range of measurements. Curves and shading present the model results and indicate year-average concentrations and the seasonal range. Please note the different horizontal scales.

Download

All of the sandy stations have one major error in common: sulfide concentrations are strongly overestimated at depths below 5 cm. We suppose that the precipitation or reoxidation of sulfide is underestimated. For all other pore water species, the agreement between measured and modelled ranges is reasonable. The rise of alkalinity with depth is especially well captured by the model. The sulfate deficit in the empirical data has a large uncertainty, as it is calculated as a small difference of similarly large quantities.

In our model, the sandy sites show a more pronounced seasonal cycle in the pore water profiles compared to the muddy stations. Especially iron and manganese concentrations vary considerably due to the seasonally different supply of quickly degradable organic matter and corresponding differences in mixing intensity. While the variability in the supply of fresh organic matter is captured by the model, the variation in mixing is not. Still, the simulated ranges are supported by the variability in the observed pore water concentrations.

Figure 6Pore water concentrations of several dissolved species at the silt station Tromper Wiek. Points and horizontal lines indicate the range of measurements. Curves and shading present the model results and indicate year-average concentrations and the seasonal range.

Download

5.1.3 Pore water profiles at the silt station Tromper Wiek

For the station Tromper Wiek, we used data from two different cruises in April and June 2014. Even if the idea in the SECOS project was to repeatedly sample the same station, the locations were approximately 6 km apart for this station, and the substrate type at the station sampled in April was sand rather than silt. The amount of sulfide in the pore waters showed a large difference between the April and the June cruise, with the latter concentrations exceeding the former by a factor of 20. This reflects spatial rather than temporal variations. Some of the depth intervals were only sampled during the June cruise, which explains the different observed ranges at the different depths.

The good agreement in the profiles of ammonium and phosphate (middle and right panel in Fig. 6) suggests that total mineralisation is captured well. The left panel shows that the model estimates the sulfate deficit at the lower range of observations, but the rise of alkalinity with depth at the higher range. The model overestimates the vertical extent of Fe (II) in the pore waters. However, the model reasonably reproduced the range of iron concentrations and also the fact that dissolved manganese concentrations were always low compared to those of iron.

Figure 7Mass fractions of nitrogen (blue), sulfur (yellow) and organic carbon (black) in the dry sediment: model results (curves and shading for seasonal range) versus measurements (vertical segments). Please note that the scales on the horizontal axes differ by a factor of 40.

Download