CSIRO Environmental Modelling Suite (EMS): scientific description of the optical and biogeochemical models (vB3p0)

Since the mid-1990s, Australia's Commonwealth Science Industry and Research Organisation (CSIRO) has been developing a biogeochemical (BGC) model for coupling with a hydrodynamic and sediment model for application in estuaries, coastal waters and shelf seas. The suite of coupled models is referred to as the CSIRO Environmental Modelling Suite (EMS) and has been applied at tens of locations around the Australian continent. At a mature point in the BGC model's development, this paper presents a full mathematical description, as well as links to the freely available code and user guide. The mathematical description is structured into processes so that the details of new parameterisations can be easily identified, along with their derivation. In EMS, the underwater light field is simulated by a spectrally resolved optical model that calculates vertical light attenuation from the scattering and absorption of 20+ optically active constituents. The BGC model itself cycles carbon, nitrogen, phosphorous and oxygen through multiple phytoplankton, zooplankton, detritus and dissolved organic and inorganic forms in multiple water column and sediment layers. The water column is dynamically coupled to the sediment to resolve deposition, resuspension and benthic–pelagic biogeochemical fluxes. With a focus on shallow waters, the model also includes detailed representations of benthic plants such as seagrass, macroalgae and coral polyps. A second focus has been on, where possible, the use of geometric derivations of physical limits to constrain ecological rates. This geometric approach generally requires population-based rates to be derived from initially considering the size and shape of individuals. For example, zooplankton grazing considers encounter rates of one predator on a prey field based on summing relative motion of the predator with the prey individuals and the search area; chlorophyll synthesis includes a geometrically derived self-shading term; and the bottom coverage of benthic plants is calculated from their biomass using an exponential form derived from geometric arguments. This geometric approach has led to a more algebraically complicated set of equations when compared to empirical biogeochemical model formulations based on populations. But while being algebraically complicated, the model has fewer unconstrained parameters and is therefore simpler to move between applications than it would otherwise be. The version of EMS described here is implemented in the eReefs project that delivers a near-real-time coupled hydrodynamic, sediment and biogeochemical simulation of the Great Barrier Reef, northeast Australia, and its formulation provides an example of the application of geometric reasoning in the formulation of aquatic ecological processes.

1 Introduction

The first model of marine biogeochemistry was developed more than 70 years ago to explain phytoplankton blooms (Riley, 1947). Today, the modelling of estuarine, coastal and global biogeochemical systems has been used for a wide variety of applications including coastal eutrophication (Madden and Kemp, 1996; Baird et al., 2003), shelf carbon and nutrient dynamics (Yool and Fasham, 2001; Dietze et al., 2009), plankton ecosystem diversity (Follows et al., 2007), ocean acidification (Orr et al., 2005), impact of local developments such as fish farms and sewerage treatment plants (Wild-Allen et al., 2010), fishery production (Stock et al., 2008) and operational forecasting (Fennel et al., 2019), to name a few. As a result of these varied applications, a diverse range of biogeochemical models has emerged, with some models developed over decades and being capable of investigating a suite of biogeochemical phenomena (Butenschön et al., 2016). With model capabilities typically dependent on the history of applications for which a particular model has been developed, and perhaps even the backgrounds and interests of the developers themselves, significant differences exist between models. Thus, it is vital that biogeochemical models are accurately described in full (e.g. Butenschön et al., 2016; Aumont et al., 2015 and Dutkiewicz et al., 2015), so that model differences can be understood and, where useful, innovations shared between modelling teams.

Estuarine, coastal and shelf modelling projects undertaken over the past 20+ years by Australia's national science agency, the Commonwealth Science Industry and Research Organisation (CSIRO), have led to the development of the CSIRO Environmental Modelling Suite (EMS). EMS contains a suite of hydrodynamic, transport, sediment, optical and biogeochemical models that can be run coupled or offline. The EMS biogeochemical model, the subject of this paper, has been applied around the Australian coastline (Fig. 1) leading to characteristics of the model which have been tailored to the Australian environment and its challenges.

Figure 1 Model domains of the CSIRO EMS hydrodynamic and biogeochemical applications from 1996 onwards. Additionally, EMS was used for the nationwide Simple Estuarine Response Model (SERM) that was applied generically around Australia's 1000+ estuaries (Baird et al., 2003). Brackets refer to specific funding bodies. EMS has also been applied in the Los Lagos region of Chile. A full list of past and current applications and funding bodies is available at https://research.csiro.au/cem/projects/ (last access: 22 September 2020).

Australian shelf waters range from tropical to temperate, micro- to macro-tidal, with shallow waters containing coral-, seagrass- or algae-dominated benthic communities. With generally narrow continental shelves, and being surrounded by two poleward-flowing boundary currents (Thompson et al., 2009), primary production in Australian coastal environments is generally limited by dissolved nitrogen in marine environments, phosphorus in freshwaters, and unlimited by silica and iron. The episodic nature of rainfall on the Australian continent, especially in the tropics, and a lack of snow cover, delivers intermittent but occasionally extreme river flows to coastal waters. With a low population density, continent-wide levels of human impacts are small relative to other continents but can be significant locally, often due to large isolated developments such as dams, irrigation schemes, mines and ports. Global changes such as ocean warming and acidification affect all regions. The EMS biogeochemical (BGC) model has many structural features similar to other models (e.g. multiple plankton functional types, nutrient and detrital pools, an increasing emphasis on optical and carbon chemistry components). Nonetheless, the geographical characteristics of, and anthropogenic influences on, the Australian continent have shaped the development of EMS and led to a BGC model with many unique features.

As the national science body, CSIRO needed to develop a numerical modelling system that could be deployed across the broad range of Australian coastal environments and capable of resolving multiple anthropogenic impacts. With a long coastline (60 000+ km by one measure), containing over 1000 estuaries, an Australian-wide configuration has insufficient resolution to be used for many applied environmental challenges. Thus, in 1999, the EMS biogeochemical model development was targeted to increase its applicability across a range of ecosystems. In particular, given limited resources to model a large number of environments/ecosystems, developments aimed to minimise the need for re-parameterisation of biogeochemical processes for each application. Two innovations arose from this imperative: (1) the software development of a process-based modelling architecture, such that model processes could be included, or excluded, while using the same executable file; and (2) the use, where possible, of geometric descriptions of physical limits to ecological processes as a means of reducing parameter uncertainty (Baird et al., 2003). It is the use of these geometric descriptions that has led to the greatest differences between EMS and other aquatic biogeochemical models.

In the aquatic sciences, there has been a long history of experimental and process studies that use geometric arguments to quantify ecological processes, but these derivations have rarely been applied in biogeochemical models, with notable exceptions (microalgal light absorption and plankton sinking rates generally, surface-area-to-volume considerations, Reynolds, 1984, size-focused trait-based modelling, Litchman and Klausmeier, 2008). By prioritising geometric arguments, EMS has included a number of previously published geometric forms including diffusion limitation of microalgae nutrient uptake (Hill and Whittingham, 1955), absorption cross-sections of microalgae (Fig. 2c; Duysens, 1956; Kirk, 1975; Morel and Bricaud, 1981), diffusion limits to macroalgae and coral nutrient uptake (Munk and Riley, 1952; Atkinson and Bilger, 1992; Zhang et al., 2011) and encounter-rate limitation of grazing rates (Fig. 2b; Jackson, 1995).

Figure 2 Examples of geometric descriptions of ecological processes. (a) The relative difference in the 2-D experience to nutrient and light fields of leaves compared to the 3-D experience of cells, as typified by the ratio of surface area (coloured) to projected area (hashed area). (b) The encounter rate of prey per individual predator as a function of the radius of encounter (the sum of the predator and prey radii) and the relative motion and prey concentration following Jackson (1995). (c) The use of ray tracing and the mass-specific absorption coefficient to calculate an absorption cross-section for a randomly oriented spheroid following Kirk (1975). (d) Fraction of the bottom covered as seen from above as a result of increasing the number of randomly placed leaves (Baird et al., 2016a), based on the assumption that leaves are randomly placed; the cover reaches $\mathrm{1}-\exp (-\mathrm{1})=\mathrm{0.63}$ when the sum of the shaded areas induced by all individual leaves equals the ground area (i.e. a leaf area index of 1).

Perhaps the most important consequence of using geometric constraints in the BGC model is the representation of benthic flora as two-dimensional surfaces, while plankton are represented as three-dimensional suspended objects (Baird and Middleton, 2004). Thus, leafy benthic plants such as macroalgae take up nutrients and absorb light on a 2-D surface. In contrast, nutrient uptake to microalgae occurs through a 3-D field, while light uptake of the 3-D cell is limited by the 2-D projected area (Fig. 2a). These contrasting geometric properties, from which the model equations are derived, generate greater potential light absorption relative to nutrient uptake of benthic communities relative to the same potential light absorption relative to nutrient uptake in unicellular algae (Baird et al., 2004). In the most simple terms, this can be related to the surface-area-to-projected-area ratio of a leaf being 0.25 times that of a microalgae cell (Fig. 2a). Thus, the competition for nutrients, ultimately being driven by light absorption and its rate compared to nutrient uptake, is explicitly determined by the contrasting geometries of cells and leaves.

In addition to geometric constraints derived by others, a number of novel geometric descriptions have been introduced into the EMS BGC model, including

1. geometric derivation of the relationship between biomass, B, and fraction of the bottom covered, $A_{\mathrm{eff}}=\mathrm{1}-\exp (-\mathrm{\Omega }B)$, where Ω is the nitrogen-specific leaf area (Sect. 6);

2. impact of self-shading on chlorophyll synthesis quantified by the incremental increase in absorption with the increase in pigment content (Sect. 5.1.3);

3. mass-specific absorption coefficients of photosynthetic pigments better utilised to determine phytoplankton absorption cross-sections (Duysens, 1956; Kirk, 1975; Morel and Bricaud, 1981) through the availability of a library of mass-specific absorption coefficients (Clementson and Wojtasiewicz, 2019) and their wavelength correction using the refractive index of the solvent used in the laboratory determinations (Fig. 5);

4. the space limitation of zooxanthellae within coral polyps using zooxanthellae projected areas in a two-layer gastrodermal cell anatomy (Sect. 6.3.1); and

5. preferential ammonium uptake, which is often calculated using different half-saturation coefficients of nitrate and ammonium uptake (Lee et al., 2002), determined by allowing ammonium uptake to proceed up to the diffusion limit. Should this diffusion limit not meet the required demand, nitrate uptake supplements the ammonium uptake. This representation has the benefit that no additional parameters are required to assign preference, with the same approach applied for both microalgae and benthic plants (Sect. 9.1).

To be clear, these geometric definitions have their own set of assumptions (e.g. a single cell size for a population) and simplifications (e.g. spherical shape). Nonetheless, the effort to apply geometric descriptions of physical limits across the BGC model appears to have been beneficial, as measured by the minimal amount of re-parameterisation that has been required to apply the model to contrasting environments. Of the above-mentioned new formulations, the most useful and easily applied is the bottom cover calculation (Fig. 2d). In fact, it is so simple, and such a clear improvement on empirical forms as demonstrated in Baird et al. (2016a), that it is likely to have been applied in other ecological/biogeochemical models, although we are unaware of any other implementation.

The geometrically constrained relationship between bottom cover and seagrass biomass, B, is $\mathrm{cover}=\mathrm{1}-\exp (-\mathrm{\Omega }B)$ and can be used to illustrate how geometric arguments can produce model equations with tightly constrained parameters. This geometric relationship contains only one parameter, Ω, that is the initial slope between cover and biomass. At low biomass, there is no overlapping of leaves, so the Ω is the area of leaves per unit of biomass (or nitrogen-specific leaf area) and has been determined by many authors in hundreds of seagrass studies. Comparison with data is shown in Appendix A of Baird et al. (2016a) and Fig. 2d. Thus, by using geometric arguments in developing the equation, the form contains only one parameter which has a physical meaning that is tightly constrained.

In addition to using geometric descriptions, there are a few other features unique to the EMS BGC model including

1. calculation of scalar irradiance from downwelling irradiance, vertical attenuation and a photon balance within a layer (Sect. 4.1.2);

2. an oxygen balance achieved through use of biological and chemical oxygen demand tracers (Sect. 10.3.2); and

3. the stoichiometric link of excess photons to reactive oxygen production in zooxanthellae.

1.1 Paper outline

This document provides a summary of the biogeochemical processes included in the model (Sect. 2), a summary of the transport model that integrates the advection–diffusion and sinking terms (Sect. 3), and full descriptions of the optical (Sect. 4) and ecological (Sect. 5–9) model equations. The description of both the optical and biogeochemical models is divided into the primary environmental zones: pelagic, epibenthic and sediment. Sect. 9 details parameterisations that are common across numerous ecological processes, such as temperature dependence, and Sect. 10 provides details of the numerical integration techniques. Further sections detail the model evaluation (Sect. 11) and test case generation (Sect. 12). The discussion (Sect. 13) details how past and present applications have influenced the development of the EMS BGC model and anticipates some future developments. Code availability is detailed at the end of the paper. Finally, the Supplement provides tables of processes (S1), state variables (S3) and parameter values (S4), with both mathematical and numerical code details (S2), and additional model evaluation (S5), from the Great Barrier Reef configuration.

2 Overview of the EMS optical and biogeochemical models

The optical and biogeochemical models are linked by the dependence of scattering and absorption on the state of optically active biogeochemical quantities. The optical model undertakes calculations at distinct wavelengths of light (say 395, 405, 415, … 705 nm) representative of individual wavebands (say 400–410, 410–420 nm, etc.) of the vertically resolved downwelling and scalar irradiance that are used by the biogeochemical model to drive photosynthesis. The optical model includes the effect of Earth–Sun distance, Sun angle, atmospheric transmission, surface albedo and refraction on the downwelling surface irradiance. In the water column, the model attenuates light based on the spectrally resolved total absorption and scattering of microalgae, detritus, dissolved organic matter, inorganic particles and the water itself (Fig. 3). The light reaching the bottom is further attenuated by macroalgae, seagrass, corals and benthic microalgae.

Figure 3 Schematic of the CSIRO Environmental Modelling Suite illustrating the biogeochemical processes in the water column, epipelagic and sediment zones, as well as the carbon chemistry and gas exchange used in vB3p0 for the Great Barrier Reef application. Orange labels represent components that scatter or absorb light.

The biogeochemical model is organised into three zones: pelagic, epibenthic and sediment. Depending on the grid formulation, the pelagic zone may have one or many layers of similar or varying thickness. The epibenthic zone overlaps with the lowest pelagic layer and the top sediment layer and shares the same dissolved and suspended particulate material fields. The sediment is modelled in multiple layers with a thin layer of easily resuspendable material overlying thicker layers of more consolidated sediment. Each sediment layer contains both particles and porewater (Margvelashvili, 2009).

Dissolved and particulate biogeochemical tracers are advected and diffused throughout the model domain in an identical fashion to temperature and salinity. Additionally, biogeochemical particulate substances sink and are resuspended in the same way as sediment particles. Biogeochemical processes are organised into pelagic processes of phytoplankton and zooplankton growth and mortality, detritus remineralisation and fluxes of dissolved oxygen, nitrogen and phosphorus; epibenthic processes of growth and mortality of macroalgae, seagrass and corals, and sediment-based processes of plankton mortality, microphytobenthos growth, detrital remineralisation and fluxes of dissolved substances (Fig. 3).

The biogeochemical model considers four groups of microalgae (small and large phytoplankton representing the functionality of photosynthetic cyanobacteria and diatoms, respectively, microphytobenthos and Trichodesmium), four macrophytes types (seagrass types corresponding to Zostera, Halophila, deep Halophila and macroalgae) and coral polyps. For temperate system applications of the EMS, dinoflagellates, Nodularia and multiple macroalgal species have also been characterised (Wild-Allen et al., 2013; Hadley et al., 2015a).

Photosynthetic growth is determined by concentrations of dissolved nutrients (nitrogen and phosphate) and photosynthetically active radiation. Autotrophs take up dissolved ammonium, nitrate, phosphate and inorganic carbon. Microalgae incorporate carbon (C), nitrogen (N) and phosphorus (P) at the Redfield ratio ($\mathrm{106}\mathrm{C}:\mathrm{16}\mathrm{N}:\mathrm{1}\mathrm{P}$), while macrophytes do so at the Atkinson ratio ($\mathrm{550}\mathrm{C}:\mathrm{30}\mathrm{N}:\mathrm{1}\mathrm{P}$). Microalgae contain chlorophyll a and a suite of accessories pigments, and have variable carbon : pigment ratios determined using a photo-adaptation model.

Micro- and mesozooplankton graze on small and large phytoplankton, respectively, at rates determined by particle encounter rates and maximum ingestion rates. Additionally, large zooplankton consume small zooplankton. Of the grazed material that is not incorporated into zooplankton biomass, a fraction is released as dissolved and particulate carbon, nitrogen and phosphate, with the remainder forming detritus. Additional detritus accumulates by mortality. Detritus and dissolved organic substances are remineralised into inorganic carbon, nitrogen and phosphate with labile detritus transformed most rapidly (days), refractory detritus slower (months) and dissolved organic material transformed over the longest timescales (years). The production (by photosynthesis) and consumption (by respiration and remineralisation) of dissolved oxygen is also included in the model and, depending on prevailing concentrations, facilitates or inhibits the oxidation of ammonium to nitrate and its subsequent denitrification to dinitrogen gas which is then lost from the system.

Additional water column chemistry calculations are undertaken to solve for the equilibrium carbon chemistry ion concentrations necessary to undertake ocean acidification (OA) studies and to consider sea–air fluxes of oxygen and carbon dioxide. The adsorption and desorption of phosphorus onto inorganic particles as a function of the oxic state of the water are also considered.

In the sediment porewaters, similar remineralisation processes occur as in the water column (Fig. 4). Additionally, nitrogen is denitrified and lost as N₂ gas, while phosphorus can become adsorbed onto inorganic particles and become permanently immobilised in sediments.

Schematics of sediment nitrogen (top) and phosphorus (bottom) pools and fluxes

Figure 4 Schematics of sediment nitrogen (top) and phosphorus (bottom) pools and fluxes. Processes represented include phytoplankton mortality, detrital decomposition, denitrification (nitrogen only), phosphorus adsorption (phosphorus only) and microphytobenthic growth. Grey boxes are particles.

2.1 Structure of the model description

The biogeochemical model presented in this paper is process based. That is, the rate of change of each ecological state variable is determined by a mathematical representation of each process that moves mass between one variable and another, conserving total mass. For dissolved inorganic phosphorus, the equation in the bottom water column layer (excluding advection, diffusion and particle sinking) could be written as

$$\begin{array}{ll}{\displaystyle \frac{\mathrm{d}P}{\mathrm{d}t}}=& {\displaystyle }-\sum ^{\mathrm{4}}\mathrm{microaglae}\mathrm{uptake}-\sum ^{\mathrm{3}}\mathrm{seagrass}\mathrm{uptake}\\ {\displaystyle }& {\displaystyle }-\mathrm{macroalgae}\mathrm{uptake}-\mathrm{zooxanthellae}\mathrm{uptake}\\ {\displaystyle }& {\displaystyle }-\mathrm{water}\mathrm{column}/\mathrm{sediment}\mathrm{porewater}\mathrm{exchange}\\ {\displaystyle }& {\displaystyle }-\mathrm{phosphorus}\mathrm{adsorption}/\mathrm{desorption}\\ {\displaystyle }& {\displaystyle }+\sum ^{\mathrm{4}}\mathrm{microalgae}\mathrm{respiration}+\sum ^{\mathrm{3}}\mathrm{seagrass}\mathrm{respiration}\\ {\displaystyle }& {\displaystyle }+\mathrm{macroalgae}\mathrm{respiration}+\mathrm{zooxanthellae}\mathrm{respiration}\\ {\displaystyle }& {\displaystyle }+\sum ^{\mathrm{2}}\mathrm{zooplankton}\mathrm{sloppy}\mathrm{feeding}\\ {\displaystyle }& {\displaystyle }+\sum ^{\mathrm{2}}\mathrm{zooplankton}\mathrm{respiration}\\ {\displaystyle }& {\displaystyle }+\mathrm{remineralisation}\mathrm{of}\mathrm{labile}\mathrm{detritus}\\ {\displaystyle }& {\displaystyle }+\mathrm{remineralisation}\mathrm{of}\mathrm{refractory}\mathrm{detritus}\\ {\displaystyle }& {\displaystyle }+\mathrm{remineralisation}\mathrm{of}\mathrm{dissolved}\mathrm{organic}\mathrm{matter}.\end{array}$$

As the number of processes in the model has grown, the representation of all the terms affecting one variable has become unworkable. Thus, instead of presenting the full equation for each state variable, we present the full set of equations for each process.

2.1.1 Presentation of process equations

In Sects. 5–9, descriptions are sorted by processes, such as microalgae growth, coral growth and food web interactions. This organisation allows the model to be explained, with individual notation, in self-contained chunks. For each process, the complete set of model equations, parameter values and state variables is given in tables. Within each process, the equations are required to conserve mass of oxygen, carbon, nitrogen and phosphorus. Furthermore, each process description is independent of any other processes in the model. As the code itself allows the inclusion/exclusion of processes at runtime, the process-based structuring of the scientific description aligns with the architecture of the numerical code.

2.1.2 Model stoichiometry

The model contains state variables that quantify the mass of carbon, nitrogen, phosphorus and oxygen, as well as state variables that contain stoichiometrically constant combinations of carbon, nitrogen, phosphorus ($\mathrm{O}:\mathrm{C}:\mathrm{N}:\mathrm{P}$ of $\mathrm{110}:\mathrm{106}:\mathrm{16}:\mathrm{1}$ for plankton and animals; $\mathrm{554}:\mathrm{550}:\mathrm{30}:\mathrm{1}$ for benthic plants). While the majority of the state variables and parameters are specified in units of nitrogen, the model could equally be specified by carbon or phosphorus. Furthermore, while the structural material of microalgae (including benthic microalgae and zooxanthellae) is at the Redfield ratio, changing reserves in microalgae of fixed carbon, nitrogen and phosphorus mean that the microalgae have a variable stoichiometry. The model has separate state variables for refractory detrital carbon, nitrogen and phosphorus, so the sum of all refractory detritus components has a variable stoichiometry. As explained later, we represent stoichiometric coefficients in the model equations as integers, a simple approximation (Table 41) to improve the readability of the mathematical equations.

3 Transport model

The local rate of change of concentration, C, of each dissolved and particulate constituent contains sink/source terms, S_C, which are described in length in the process descriptions of this document, and the advection, diffusion and sinking terms:

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\partial C}{\partial t}}+\mathbf{v}\cdot {\mathrm{\nabla }}^{\mathrm{2}}C=\mathrm{\nabla }\cdot \left(K\mathrm{\nabla }C\right)+w_{\mathrm{sink}}{\displaystyle \frac{\partial C}{\partial z}}+S_C,\end{array}$$

where the symbol $\mathrm{\nabla }=\left(\frac{\partial }{\partial x},\frac{\partial }{\partial y},\frac{\partial }{\partial z}\right)$, v is the velocity field, K is the eddy diffusion coefficient which varies in space and time, and w_sink is the local sinking rate (positive downwards) with the z coordinate positive upwards. The calculation of v and K is described in the hydrodynamic model (Herzfeld, 2006; Gillibrand and Herzfeld, 2016). The advection–diffusion terms of Eq. (1), based on the continuum hypothesis for a fluid (Vichi et al., 2007), are solved by either an in-line advection scheme with the baroclinic time step of the hydrodynamic model or an offline transport scheme using a potentially much longer time step (Gillibrand and Herzfeld, 2016). Options for advection and transport schemes in EMS include mass conservative Lagrangian and flux-form approaches described in Herzfeld (2006) and Gillibrand and Herzfeld (2016).

The microalgae are particulates that contain internal concentrations of dissolved nutrients (C, N, P) and pigments that are specified on a per cell basis. To conserve mass, the local rate of change of the concentration of microalgae, B, multiplied by the content (or reserve) of the cell, R, is given by

$$\begin{array}{}\text{(2)}& \begin{array}{rl}{\displaystyle \frac{\partial \left(BR\right)}{\partial t}}+\mathbf{v}\cdot {\mathrm{\nabla }}^{\mathrm{2}}\left(BR\right)=& \mathrm{\nabla }\cdot \left(K\mathrm{\nabla }\left(BR\right)\right)\\ & +w_C{\displaystyle \frac{\partial \left(BR\right)}{\partial z}}+S_{BR}.\end{array}\end{array}$$

For more information, see Sect. 5.1.6 and Sect. 3.1 of Baird et al. (2004) which describe the coupling of the plankton component of the biogeochemical model to the Princeton Ocean Model.

4 Optical model

The optical model calculates the spectrally resolved light field in each vertical column and uses it to drive the photosynthesis of phytoplankton and benthic plants in the biogeochemical model. Following the terminology of aquatic optics (Mobley, 1994), we divide the description of the model into calculations of inherent optical properties (IOPs) followed by apparent optical properties (AOPs). IOPs are properties of the medium (e.g. scattering and absorption) and do not depend on the ambient light field. The optical model uses the value of the optically active state variables, and their mass-specific absorption and scattering properties, to calculate the total absorption and scattering. AOPs are those properties that depend both on the medium (the IOPs) and on the surface light field. AOPs include downwelling and scalar irradiance. Thus, the optical model uses the vertical distribution of IOPs, and the surface light field, to determine the vertical distribution of the AOPs.

4.1 Water column optical model

4.1.1 Inherent optical properties (IOPs)

Phytoplankton absorption. The model contains four phytoplankton types (small and large phytoplankton, benthic microalgae and Trichodesmium), each with a unique ratio of internal concentration of accessory photosynthetic pigments to chlorophyll a. To calculate the absorption due to each pigment, we use a database of spectrally resolved mass-specific absorption coefficients (Clementson and Wojtasiewicz, 2019). As it can be assumed that accessory pigments stay in a constant ratio to chlorophyll a, the model needs only a state variable for chlorophyll a for each phytoplankton type. The model then calculates the chlorophyll a specific absorption coefficient due to all pigments by using a state variable quantifying the chlorophyll concentration of the population, the number of cells in the population and the ratio of concentration of the accessory pigment to chlorophyll a, and the mass-specific absorption coefficient of each of the accessory pigments. Thus, the chlorophyll a specific absorption coefficient due to all photosynthetic pigments for small phytoplankton at wavelength λ, γ_small,λ, is given by

$$\begin{array}{}\text{(3)}& \begin{array}{rl}{\mathit{\gamma }}_{\mathrm{small},\mathit{\lambda }}=& \mathrm{1.0}{\mathit{\gamma }}_{\mathrm{Chl}a,\mathit{\lambda }}+\mathrm{0.35}{\mathit{\gamma }}_{\mathrm{Zea},\mathit{\lambda }}+\mathrm{0.05}{\mathit{\gamma }}_{\mathrm{Echi},\mathit{\lambda }}\\ & +\mathrm{0.1}{\mathit{\gamma }}_{\mathit{\beta }-\mathrm{car},\mathit{\lambda }}+\mathrm{2}{\mathit{\gamma }}_{\mathrm{PE},\mathit{\lambda }}+\mathrm{1.72}{\mathit{\gamma }}_{\mathrm{PC},\mathit{\lambda }},\end{array}\end{array}$$

where Chl a is the pigment chlorophyll a, Zea is zeaxanthin, Echi is echinenone, β-car is β carotene, PE is phycoerythrin, and PC is phycocyanin, and the ratios of chlorophyll a to accessory pigment concentration are determined from Wojtasiewicz and Stoń-Egiert (2016). Note that the coefficient in Eq. (3) for Chl a is 1.0 because the ratio of chlorophyll a to chlorophyll a is 1. The resulting chlorophyll a specific absorption coefficient for picoplankton is shown in Fig. 5.

Figure 5 Pigment-specific absorption coefficients for the dominant pigments found in small phytoplankton determined using laboratory standards in solvent in a 1 cm vial. Green and red lines are photosynthetic pigments constructed from 563 measured wavelengths. Circles represent the wavelengths at which the optical properties are calculated in the simulations. The black line represents the weighted sum of the photosynthetic pigments (Eq. 3), with the weighting calculated from the ratio of each pigment concentration to chlorophyll a. The spectra are wavelength shifted from their raw measurement by the ratio of the refractive index of the solvent to the refractive index of water (1.352 for acetone used with chlorophyll a and β carotene; 1.361 for ethanol used with zeaxanthin and echinenone; 1.330 for water used with phycoerythrin and phycocyanin).

Similarly, for large phytoplankton and microphytobenthos (Wright et al., 1996),

$$\begin{array}{}\text{(4)}& {\mathit{\gamma }}_{\mathrm{large},\mathit{\lambda }}=\mathrm{1.0}{\mathit{\gamma }}_{\mathrm{Chl}a,\mathit{\lambda }}+\mathrm{0.6}{\mathit{\gamma }}_{\mathrm{Fuco},\mathit{\lambda }},\end{array}$$

where Fuco is fucoxanthin (Fig. 6) and for Trichodesmium (Carpenter et al., 1993),

$$\begin{array}{}\text{(5)}& \begin{array}{rl}{\mathit{\gamma }}_{\mathrm{Tricho},\mathit{\lambda }}=& \mathrm{1.0}{\mathit{\gamma }}_{\mathrm{Chl}a,\mathit{\lambda }}+\mathrm{0.1}{\mathit{\gamma }}_{\mathrm{Zea},\mathit{\lambda }}+\mathrm{0.02}{\mathit{\gamma }}_{\mathrm{Myxo},\mathit{\lambda }}\\ & +\mathrm{0.09}{\mathit{\gamma }}_{\mathit{\beta }-\mathrm{car},\mathit{\lambda }}+\mathrm{2.5}{\mathit{\gamma }}_{\mathrm{PE},\mathit{\lambda }},\end{array}\end{array}$$

where Myxo is myxoxanthophyll (Fig. 7).

Figure 6 Pigment-specific absorption coefficients for the dominant pigments found in large phytoplankton and microphytobenthos determined using laboratory standards in solvent in a 1 cm vial. The aqua line represents the weighted sum of the photosynthetic pigments (Eq. 4), with the weighting calculated from the ratio of each pigment concentration to chlorophyll a. See Fig. 5 for more details. Fucoxanthin was dissolved in ethanol.

Figure 7 Pigment-specific absorption coefficients for the dominant pigments found in Trichodesmium determined using laboratory standards in solvent in a 1 cm vial. The aqua line represents the weighted sum of the photosynthetic pigments (Eq. 5), with the weighting calculated from the ratio of each pigment concentration to chlorophyll a. See Fig. 5 for more details. Myxoxanthophyll was dissolved in acetone.

The absorption cross-section at wavelength λ (α_λ) of a spherical cell of radius (r), chlorophyll a specific absorption coefficient (γ_λ) and homogeneous intracellular chlorophyll a concentration (c_i) can be calculated using geometric optics (i.e. ray tracing without considering internal scattering) and is given by (Duysens, 1956; Kirk, 1975)

$$\begin{array}{}\text{(6)}& {\mathit{\alpha }}_{\mathit{\lambda }}=\mathit{\pi }{r}^{\mathrm{2}}\left(\mathrm{1}-{\displaystyle \frac{\mathrm{2}(\mathrm{1}-(\mathrm{1}+\mathrm{2}{\mathit{\gamma }}_{\mathit{\lambda }}{c}_{\mathrm{i}}r\left)e^{-\mathrm{2}{\mathit{\gamma }}_{\mathit{\lambda }}{c}_{\mathrm{i}}r}\right)}{(\mathrm{2}{\mathit{\gamma }}_{\mathit{\lambda }}{c}_{\mathrm{i}}r{)}^{\mathrm{2}}}}\right),\end{array}$$

where πr² is the projected area of a sphere, and the bracketed term is 0 for no absorption (γc_ir=0) and approaches 1 as the cell becomes fully opaque (γc_ir→∞). Note that the bracketed term in Eq. (6) is mathematically equivalent to the dimensionless efficiency factor for absorption, Q_a (used in Morel and Bricaud, 1981, Finkel, 2001 and Bohren and Huffman, 1983), of homogeneous spherical cells with an index of refraction close to that of the surrounding water. Note that the intracellular chlorophyll concentration, c_i, changes as a result of chlorophyll synthesis (described later in Eq. 34).

The use of an absorption cross-section of an individual cell has two significant advantages. Firstly, the same model parameters used here to calculate absorption in the water column are used to determine photosynthesis by individual cells, including the effect of packaging of pigments within cells. Secondly, the dynamic chlorophyll concentration determined later can be explicitly included in the calculation of phytoplankton absorption. Thus, the absorption of a population of n cell m⁻³ is given by nα m⁻¹, while an individual cell absorbs αE_o light, where E_o is the scalar irradiance.

Coloured dissolved organic matter (CDOM) absorption. Two equations for CDOM absorption are presently being trialled. The two schemes are as follows: Scheme 1. The absorption of CDOM, a_CDOM,λ, is determined from a relationship with salinity in the region (Schroeder et al., 2012):

$$\begin{array}{}\text{(7)}& a_{\mathrm{CDOM},\mathrm{443}}=-\mathrm{0.0332}S+\mathrm{1.2336},\end{array}$$

where S is the salinity. In order to avoid unrealistic extrapolation, the salinity used in this relationship is the minimum of the model salinity and 36. In some cases, coastal salinities exceed 36 due to evaporation. The absorption due to CDOM at other wavelengths is calculated using a CDOM spectral slope for the region (Blondeau-Patissier et al., 2009):

$$\begin{array}{}\text{(8)}& a_{\mathrm{CDOM},\mathit{\lambda }}=a_{\mathrm{CDOM},\mathrm{443}}\exp \left(-S_{\mathrm{CDOM}}\left(\mathit{\lambda }-\mathrm{443.0}\right)\right),\end{array}$$

where S_CDOM is an approximate spectral slope for CDOM, with observations ranging from 0.01 to 0.02 nm⁻¹ for significant concentrations of CDOM. Lower magnitudes of the spectral slope generally occur at lower concentrations of CDOM (Blondeau-Patissier et al., 2009).

Scheme 2. The absorption of CDOM, a_CDOM,λ, is directly related to the concentration of dissolved organic carbon, D_C:

$$\begin{array}{}\text{(9)}& a_{\mathrm{CDOM},\mathit{\lambda }}=k_{\mathrm{CDOM},\mathrm{443}}^{*}{D}_{\mathrm{C}}\exp \left(-S_{\mathrm{CDOM}}\left(\mathit{\lambda }-\mathrm{443.0}\right)\right),\end{array}$$

where $k_{\mathrm{CDOM},\mathrm{443}}^{*}$ is the dissolved organic carbon-specific CDOM absorption coefficient at 443 nm.

Both schemes have drawbacks. Scheme 2, using the concentration of dissolved organic carbon, is closer to reality but is likely to be sensitive to poorly known parameters such as remineralisation rates and initial detritial concentrations. Scheme 1, a function of salinity, will be more stable but perhaps less accurate, especially in estuaries where hypersaline waters may have large estuarine loads of coloured dissolved organic matter.

Absorption due to non-algal particulate material. The waters of the Great Barrier Reef contain suspended sediments originating from various marine sources, such as the white calcium carbonate fragments generated by coral erosion and sediments derived from terrestrial sources such as granite (Soja-Woźniak et al., 2019). The model uses spectrally resolved mass-specific absorption coefficients (and also total scattering measurements) from a database of laboratory measurements conducted on either pure mineral suspensions, or mineral mixtures, at two ranges of size distributions (Fig. 8; Stramski et al., 2007). In this model version, we use the calcium carbonate sample CAL1 for CaCO₃-based particles.

Figure 8 The remote-sensing reflectance of the 21 mineral mixtures suspended in water as measured by Stramski et al. (2007). Laboratory measurements of absorption and scattering properties are used to calculated remote-sensing reflectance (Baird et al., 2016b). Line colouring corresponds to that produced by the mineral suspended in clear water as calculated using the MODIS true colour algorithm (Gumley et al., 2010). CAL1, with a median particle diameter of 2 µm, is used for Mud${}_{{\mathrm{CaCO}}_{\mathrm{3}}}$.

For the terrestrially sourced particles, we used observations from Gladstone Harbour in the central Great Barrier Reef (GBR) (Fig. 9). These IOPs gave a realistic surface colour for the Queensland river sediment plumes (Baird et al., 2016b). In the model, optically active non-algal particulates (NAPs) include the inorganic particulates (such as sand and mud; see Sect. 7.1) and detritus. We assumed the optical properties of the detritus were the same as the optical properties in Gladstone Harbour, although open-ocean studies have used a detritial absorption that is more like CDOM (Dutkiewicz et al., 2015).

IOPs from Gladstone Harbour.

Figure 9 Inherent optical properties (total absorption and total scattering) at sample sites in Gladstone Harbour on 13–19 September 2013 (Babcock et al., 2015). The line colour is rendered like that in Fig. 8. The site labelling is ordered in time, from the first sample collected during neap tides at the top to the last sample collected at spring tides on the bottom. The IOPs used for the Mud${}_{\mathrm{non}-{\mathrm{CaCO}}_{\mathrm{3}}}$ end-member are from the WIT site (see legend) at the centre of the harbour, was dominated by inorganic particles. The measured concentration of NAP at the site was 33.042 mg L⁻¹ and is used to calculate mass-specific IOPs.

The absorption due to calcite-based NAP is given by

$$\begin{array}{}\text{(10)}& a_{{\mathrm{NAP}}_{{\mathrm{CaCO}}_{\mathrm{3}},\mathit{\lambda }}}=c_{\mathrm{1}}{\mathrm{NAP}}_{{\mathrm{CaCO}}_{\mathrm{3}}},\end{array}$$

where c₁ is the spectrally resolved mass-specific absorption coefficient determined from laboratory experiments (Fig. 8). The absorption due to non-calcite NAPs, NAP${}_{\mathrm{non}-{\mathrm{CaCO}}_{\mathrm{3}}}$, combined with detritus, is given by

$$\begin{array}{}\text{(11)}& \begin{array}{rl}& a_{{\mathrm{NAP}}_{\mathrm{non}-{\mathrm{CaCO}}_{\mathrm{3}},\mathit{\lambda }}}=c_{\mathrm{2}}{\mathrm{NAP}}_{\mathrm{non}-{\mathrm{CaCO}}_{\mathrm{3}}}\\ & +\left({\displaystyle \frac{\mathrm{550}}{\mathrm{30}}}{\displaystyle \frac{\mathrm{12}}{\mathrm{14}}}{D}_{\mathrm{Atk}}+{\displaystyle \frac{\mathrm{106}}{\mathrm{16}}}{\displaystyle \frac{\mathrm{12}}{\mathrm{14}}}{D}_{\mathrm{Red}}+D_{\mathrm{C}}\right)/{\mathrm{10}}^{\mathrm{6}},\end{array}\end{array}$$

where c₂ is the spectrally resolved mass-specific absorption coefficient determined from field measurements (Fig. 9), NAP${}_{\mathrm{non}-{\mathrm{CaCO}}_{\mathrm{3}}}$ is quantified in kg m⁻³, D_Atk and D_Red are quantified in mg N m⁻³, and D_C is quantified in mg C m⁻³.

Total absorption. The total absorption, a_T,λ, is given by

$$\begin{array}{}\text{(12)}& \begin{array}{rl}{a}_{\mathrm{T},\mathit{\lambda }}=& a_{\mathrm{w},\mathit{\lambda }}+a_{{\mathrm{NAP}}_{\mathrm{non}-{\mathrm{CaCO}}_{\mathrm{3}}},\mathit{\lambda }}+a_{{\mathrm{NAP}}_{{\mathrm{CaCO}}_{\mathrm{3}}},\mathit{\lambda }}\\ & +a_{\mathrm{CDOM},\mathit{\lambda }}+\sum _{x=\mathrm{1}}^N{n}_x{\mathit{\alpha }}_{x,\mathit{\lambda }},\end{array}\end{array}$$

where a_w,λ is clear-water absorption (Fig. 10) and N is the number of phytoplankton classes (see Table 3).

Table 1 Constants and parameter values used in the optical model.

^a Kirk (1994). ^b Kirk (1991) using an average cosine of scattering of 0.924 (Mobley, 1994). ^c Blondeau-Patissier et al. (2009); see also Cherukuru et al. (2019).

Table 2 State and derived variables in the water column optical model.

Table 3 Traits of suspended microalgae.

^* At T_ref=20 ^∘C.

Figure 10 Spectrally resolved energy distribution of sunlight, clear-water absorption and clear-water scattering (Smith and Baker, 1981). The fraction of solar radiation between 400 and 700 nm for clear-sky irradiance at the particular spectral resolution is given in panel (a). The centre of each waveband used in the model simulations is identified by a cross on each curve. Panel (d) shows the pigment-specific absorption of Chl a and generic photosynthetic carotenoids (Ficek et al., 2004) that were used in earlier versions of this model (Baird et al., 2016b) before the mass-specific absorption coefficients of multiple accessory pigments were implemented (Figs. 5, 6 and 7).

Scattering. The total scattering coefficient is given by

$$\begin{array}{}\text{(13)}& \begin{array}{rl}{b}_{\mathrm{T},\mathit{\lambda }}=& b_{\mathrm{w},\mathit{\lambda }}+c_{\mathrm{1}}{\mathrm{NAP}}_{\mathrm{non}-{\mathrm{CaCO}}_{\mathrm{3}}}+c_{\mathrm{2}}{\mathrm{NAP}}_{{\mathrm{CaCO}}_{\mathrm{3}}}\\ & +b_{\mathrm{phy},\mathit{\lambda }}\sum _{x=\mathrm{1}}^N{n}_x{c}_{i,x}{V}_x,\end{array}\end{array}$$

where NAP is the concentration of non-algal particulates, b_w,λ is the scattering coefficient due to clear water (Fig. 10), c₁ and c₂ are the spectrally resolved mass-specific coefficients (Figs. 8 and 9) and phytoplankton scattering is the product of the chlorophyll-specific phytoplankton scattering coefficient, b_phy,λ and the water column chlorophyll concentration of all classes, $\sum n_x{c}_{i,x}{V}_x$ (where c_i is the chlorophyll concentration in the cell, and V is the cell volume). The value for b_phy,λ is set to 0.2 (mg Chl a m⁻²)⁻¹ for all wavelengths, a typical value for marine phytoplankton (Kirk, 1994). For more details, see Baird et al. (2007b).

4.1.2 Apparent optical properties (AOPs)

The optical model is forced with the downwelling short-wave radiation just above the sea surface, based on remotely sensed cloud fraction observations and calculations of top-of-the-atmosphere clear-sky irradiance and solar angle. The calculation of downwelling radiation and surface albedo (a function of solar elevation and cloud cover) is detailed in Sect. 9.1.1 of the hydrodynamic scientific description (https://research.csiro.au/cem/software/ems/ems-documentation/, last access: 22 September 2020).

The downwelling irradiance just above the water interface is split into wavebands using the weighting for clear-sky irradiance (Fig. 10). Snell's law is used to calculate the azimuth angle of the mean light path through the water, θ_sw, as calculated from the atmospheric azimuth angle, θ_air, and the refraction of light at the air–water interface (Kirk, 1994):

$$\begin{array}{}\text{(14)}& {\displaystyle \frac{\sin {\mathit{\theta }}_{\mathrm{air}}}{\sin {\mathit{\theta }}_{\mathrm{sw}}}}=\mathrm{1.33}.\end{array}$$

Calculation of in-water light field. Given the IOPs determined above, the exact solution for AOPs would require a radiative transfer model (Mobley, 1994) which is too computationally expensive for a complex ecosystem model such as developed here. Instead, the in-water light field is solved for using empirical approximations of the relationship between IOPs and AOPs (Kirk, 1991; Mobley, 1994).

The vertical attenuation coefficient at wavelength λ when considering absorption and scattering, K_λ, is given by

$$\begin{array}{}\text{(15)}& K_{\mathit{\lambda }}={\displaystyle \frac{a_{\mathrm{T},\mathit{\lambda }}}{\cos {\mathit{\theta }}_{\mathrm{sw}}}}\sqrt{\mathrm{1}+\left(g_{\mathrm{i}}\cos {\mathit{\theta }}_{\mathrm{sw}}-g_{\mathrm{ii}}\right){\displaystyle \frac{b_{\mathrm{T},\mathit{\lambda }}}{a_{\mathrm{T},\mathit{\lambda }}}}}.\end{array}$$

The term outside the square root quantifies the effect of absorption, where a_T,λ is the total absorption. The term within the square root of Eq. (15) represents scattering as an extended pathlength through the water column, where g_i and g_ii are empirical constants and take values of 0.402 and 0.180, respectively. The values of g_i and g_ii depend on the average cosine of scattering. For filtered water with scattering only due to water molecules, the values of g_i and g_ii are quite different to natural waters. But for waters ranging from coastal to open ocean, the average cosine of scattering varies by only a small amount (0.86–0.95; Kirk, 1991), and thus uncertainties in g_i and g_ii do not strongly affect K_λ.

The downwelling irradiance at wavelength λ at the bottom of a layer h thick, $E_{\mathrm{d},\mathit{\lambda },\mathrm{bot}}$, is given by

$$\begin{array}{}\text{(16)}& E_{\mathrm{d},\mathrm{bot},\mathit{\lambda }}=E_{\mathrm{d},\mathrm{top},\mathit{\lambda }}{e}^{-K_{\mathit{\lambda }}h},\end{array}$$

where $E_{\mathrm{d},\mathrm{top},\mathit{\lambda }}$ is the downwelling irradiance at wavelength λ at the top of the layer and K_λ is the vertical attenuation coefficient at wavelength λ, a result of both absorption and scattering processes.

Assuming a constant attenuation rate within the layer, the average downwelling irradiance at wavelength λ, E_d,λ, is given by

$$\begin{array}{}\text{(17)}& E_{\mathrm{d},\mathit{\lambda }}={\displaystyle \frac{\mathrm{1}}{h}}\underset{\mathrm{bot}}{\overset{\mathrm{top}}{\int }}{E}_{\mathrm{d},z,\mathit{\lambda }}{e}^{-K_{\mathit{\lambda }}z}\mathrm{d}z={\displaystyle \frac{E_{\mathrm{d},\mathrm{top},\mathit{\lambda }}-E_{\mathrm{d},\mathrm{bot},\mathit{\lambda }}}{K_{\mathit{\lambda }}h}}.\end{array}$$

We can now calculate the scalar irradiance, E_o, for the calculation of phytoplankton absorption, from downwelling irradiance, E_d. The light absorbed within a layer must balance the difference in downwelling irradiance from the top and bottom of the layer (since scattering in this model only increases the pathlength of light); thus,

$$\begin{array}{}\text{(18)}& E_{o,\mathit{\lambda }}{a}_{\mathrm{T},\mathit{\lambda }}h=E_{\mathrm{d},\mathrm{top},\mathit{\lambda }}-E_{\mathrm{d},\mathrm{bot},\mathit{\lambda }}=E_{\mathrm{d},\mathit{\lambda }}{K}_{\mathit{\lambda }}h.\end{array}$$

Cancelling h, the scalar irradiance as a function of downwelling irradiance is given by

$$\begin{array}{}\text{(19)}& E_{o,\mathit{\lambda }}={\displaystyle \frac{E_{\mathrm{d},\mathit{\lambda }}{K}_{\mathit{\lambda }}}{a_{\mathrm{T},\mathit{\lambda }}}}.\end{array}$$

This correction conserves photons within the layer, although it is only as good as the original approximation of the impact of scattering and azimuth angle on vertical attenuation (Eq. 15).

Vertical attenuation of heat. The vertical attenuation of heat is given by

$$\begin{array}{}\text{(20)}& K_{\mathrm{heat}}=-\int {\displaystyle \frac{\mathrm{1}}{E_{\mathrm{d},z,\mathit{\lambda }}}}{\displaystyle \frac{\partial E_{\mathrm{d},z,\mathit{\lambda }}}{\partial z}}\mathrm{d}\mathit{\lambda },\end{array}$$

and the local heating by

$$\begin{array}{}\text{(21)}& {\displaystyle \frac{\partial T}{\partial t}}=-{\displaystyle \frac{\mathrm{1}}{\mathit{\rho }{c}_p}}\int {\displaystyle \frac{\partial E_{\mathrm{d},\mathit{\lambda }}}{\partial z}}\mathrm{d}\mathit{\lambda },\end{array}$$

where T is temperature, ρ is the density of water, and c_p is the specific heat of water.

4.2 Epibenthic optical model

The spectrally resolved light field at the base of the water column is attenuated, from top to bottom, by macroalgae and seagrass (Zostera, then shallow and then deep forms of Halophila), followed by the zooxanthellae in corals. The downwelling irradiance at wavelength λ after passing through each macroalgae and seagrass species is given by E_below,λ:

$$\begin{array}{}\text{(22)}& E_{\mathrm{below},\mathit{\lambda }}=E_{\mathrm{d},\mathrm{above},\mathit{\lambda }}{e}^{-A_{\mathit{\lambda }}{\mathrm{\Omega }}_{X}X},\end{array}$$

where E_above,λ for macroalgae is $E_{\mathrm{d},\mathrm{bot},\mathit{\lambda }}$, the downwelling irradiance of the bottom water column layer, A_λ is the leaf-specific absorptance, Ω is the nitrogen-specific leaf area, and X is the leaf nitrogen biomass.

The light absorbed by corals is assumed to be entirely due to zooxanthellae and is given by

$$\begin{array}{}\text{(23)}& E_{\mathrm{below},\mathit{\lambda }}=E_{\mathrm{above},\mathit{\lambda }}{e}^{-n{\mathit{\alpha }}_{\mathit{\lambda }}},\end{array}$$

where $n=\mathrm{CS}/m_{N,\mathrm{CS}}$ is the areal density of zooxanthellae cells and α_λ is the absorption cross-section of a cell a result of the absorption of multiple pigment types.

4.3 Sediment optical model

The optical model in the sediment only concerns the benthic microalgae growing in the porewaters of the top sediment layer. The calculation of light absorption by benthic microalgae assumes they are the only attenuating component in a layer that lies on top layer of sediment, with a perfectly absorbing layer below and no scattering by any other components in the layer. Thus, no light penetrates through to the second sediment layer where benthic microalgae also reside. Thus, the downwelling irradiance at wavelength λ at the bottom of a layer, $E_{\mathrm{d},\mathit{\lambda },\mathrm{bot}}$, is given by

$$\begin{array}{}\text{(24)}& E_{\mathrm{d},\mathrm{bot},\mathit{\lambda }}=E_{\mathrm{d},\mathrm{top},\mathit{\lambda }}{e}^{-n{\mathit{\alpha }}_{\mathit{\lambda }}h},\end{array}$$

where $E_{\mathrm{d},\mathrm{top},\mathit{\lambda }}$ is the downwelling irradiance at wavelength λ at the top of the layer and α_λ is the absorption cross-section of the cell at wavelength λ, and n is the concentration of cells in the layer. The layer thickness used here, h, is the thickness of the top sediment layer, so as to convert the concentration of cells in that layer, n, into the areal concentration of cells in the biofilm, n h.

Given no scattering in the cell, and that the vertical attenuation coefficient is independent of azimuth angle, the scalar irradiance that the benthic microalgae is exposed to in the surface biofilm is given by

$$\begin{array}{}\text{(25)}& E_{o,\mathit{\lambda }}=\left(E_{\mathrm{d},\mathrm{top},\mathit{\lambda }}-E_{\mathrm{d},\mathrm{bot},\mathit{\lambda }}\right)/\left(n{\mathit{\alpha }}_{\mathit{\lambda }}h\right).\end{array}$$

The photons captured by each cell, and the microalgae process, follow the same equations as for the water column (Sect. 5.1.3).

5 Pelagic processes

5.1 Microalgae

The model contains four functional groups of suspended microalgae: small and large phytoplankton, microphytobenthos and Trichodesmium. The growth from internal reserves for each of the functional groups is identical and explained below. The differences in the ecological interactions of the four functional groups are summarised in Table 3. Trichodesmium, a nitrogen fixer, contains additional processes (Sect. 5.2).

5.1.1 Microalgal growth

The growth of microalgae has been modelled in many ways, from simple exponential growth and logistic growth curves, to single- and multiple-nutrient-based curves, through to equations that contain a state variable for the physiological state of the cell (variously described as stores, quotas, reserves, etc.) and to consider the complex processing of photons in the microalgae photosystem. It is now common for complex biogeochemical models to contain state variables for the physiological state of each of the potentially limiting nutrients (Baretta-Bekker et al., 1997; Vichi et al., 2007) and include adaptation to photosystems (Geider et al., 1998). In the context of many different microalgae models, the model that is described here has taken another path again. As articulated above, we chose to base nutrient uptake and light absorption on using geometric constraints. This meant that any growth model needed to be formulated around the maximum rate of supply of each of the limiting nutrients (and light) (see Fig. 2 of Baird et al., 2006).

In the microalgae model (most fully described in Baird et al., 2001), the uptake of nutrients and light absorption increases the reserves of nutrients and light, as quantified by a reserve, R, which has units of mass per cell. In the equations, we often use a normalised reserve, R^*, which is a quantity between 0 and 1 (Table 4). The reserves are in turn consumed to generate structural material. Thus, the total content of nitrogen in the microalgae is equal to the sum of the structural material and the reserves.

Table 4 State and derived variables for the microalgae growth model. DIN is given by the sum of nitrate and ammonium concentrations, [NO₃]+[NH₄].

The model considers the diffusion-limited supply of dissolved inorganic nutrients (N and P) and the absorption of light, delivering N, P and fixed C to the internal reserves of the cell (Fig. 11). Nitrogen and phosphorus are taken directly into the reserves, but carbon is first fixed through photosynthesis (Kirk, 1994):

$$\begin{array}{}\text{(26)}& \begin{array}{rl}\mathrm{106}{\mathrm{CO}}_{\mathrm{2}}+& \mathrm{212}{\mathrm{H}}_{\mathrm{2}}\mathrm{O}\stackrel{\mathrm{1060}\mathrm{photons}}{⟶}\mathrm{106}{\mathrm{CH}}_{\mathrm{2}}\mathrm{O}\\ & +\mathrm{106}{\mathrm{H}}_{\mathrm{2}}\mathrm{O}+\mathrm{106}{\mathrm{O}}_{\mathrm{2}}.\end{array}\end{array}$$

The internal reserves of C, N and P are consumed to form structural material at the Redfield ratio (Redfield et al., 1963):

$$\begin{array}{}\text{(27)}& \begin{array}{rl}\mathrm{106}{\mathrm{CH}}_{\mathrm{2}}\mathrm{O}& +\mathrm{16}{\mathrm{NH}}_{\mathrm{4}}^{+}+{\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}+\mathrm{16}{\mathrm{H}}_{\mathrm{2}}\mathrm{O}\\ ⟶& \left({\mathrm{CH}}_{\mathrm{2}}\mathrm{O}{)}_{\mathrm{106}}\right({\mathrm{NH}}_{\mathrm{3}}{)}_{\mathrm{16}}{\mathrm{H}}_{\mathrm{3}}{\mathrm{PO}}_{\mathrm{4}}+\mathrm{13}{\mathrm{H}}^{+},\end{array}\end{array}$$

where we have represented nitrogen as ammonium (NH₄) in Eq. (27). When the nitrogen source to the cell is nitrate, NO₃, it is assumed to lose its oxygen at the cell wall (Sect. 9.1). The growth rate of microalgae is given by the maximum growth rate, μ^max, multiplied by the normalised reserves, R^*, of each of C, N and P:

$$\begin{array}{}\text{(28)}& \mathit{\mu }={\mathit{\mu }}^{max}{R}_{\mathrm{C}}^{*}{R}_{\mathrm{N}}^{*}{R}_{\mathrm{P}}^{*}.\end{array}$$

The mass of the reserves (and therefore the total $\mathrm{C}:\mathrm{N}:\mathrm{P}:\mathrm{Chl}$ a ratio) of the cell depends on the interaction of the supply and consumption rates (Fig. 11). When consumption exceeds supply, and the supply rates are non-Redfield, the normalised internal reserves of the non-limiting nutrients approach 1, while the limiting nutrient becomes depleted. Thus, the model behaves like a “law of the minimum” growth model, except during fast changes in nutrient supply rates.

Figure 11 Schematic of the process of microalgae growth from internal reserves. Blue circle: structural material; red pie: nitrogen reserves; purple pie: phosphorus reserves; yellow pie: carbon reserves; green pie: pigment content. Here, a circular pie has a value of 1, representing the normalised reserve (a value between 0 and 1). The box shows that generating structural material for an additional cell requires the equivalent of 100 % internal reserves of carbon, nitrogen and phosphorus of one cell. This figure shows the discrete growth of two cells to three, requiring both the generation of new structural material from reserves and the reserves being diluted as a result of the number of cells in which they are divided increasing from two to three. Thus, the internal reserves for nitrogen after the population increases from two to three is given by two from the initial two cells, minus one for building structural material of the new cell, shared across the three offspring, to give one-third. The same logic applies to carbon and phosphorus reserves, with phosphorus reserves being reduced to one-sixth and carbon reserves being exhausted. In contrast, pigment is not required for structural material so the only reduction is through dilution; the three-fourths content of two cells is shared among three cells to equal one-half in the three cells. This schematic shows one limitation of a population-style model whereby reserves are “shared” across the population (as opposed to individual-based modelling; Beckmann and Hense, 2004). A proof of the conservation of mass for this scheme, including under mixing of populations of suspended microalgae, is given in Baird et al. (2004). The model equations also include terms affecting internal reserves through nutrient uptake, light absorption, respiration and mortality that are not shown in this simple schematic.

The molar ratio of a cell, the addition of structural material and reserves, is given by

$$\begin{array}{}\text{(29)}& \mathrm{C}:\mathrm{N}:\mathrm{P}=\mathrm{106}(\mathrm{1}+R_{\mathrm{C}}^{*}):\mathrm{16}(\mathrm{1}+R_{\mathrm{N}}^{*}):\mathrm{1}+R_{\mathrm{P}}^{*}.\end{array}$$

5.1.2 Nutrient uptake

The diffusion-limited nutrient uptake to a single phytoplankton cell, J, is given by

$$\begin{array}{}\text{(30)}& J=\mathit{\psi }D\left(C_{\mathrm{b}}-C_{\mathrm{w}}\right),\end{array}$$

where ψ is the diffusion shape factor (=4πr for a sphere), D is the molecular diffusivity of the nutrient, C_b is the average extracellular nutrient concentration, and C_w is the concentration at the wall of the cell. The diffusion shape factor is determined by equating the divergence of the gradient of the concentration field in the vicinity of the cell to zero (∇²C=0). A semi-empirical correction to Eq. (30), to account for fluid motion around the cell, and the calculation of non-spherical diffusion shape factors, has been applied in earlier work (Baird and Emsley, 1999). For the purposes of biogeochemical modelling, these uncertain corrections for small-scale turbulence and non-spherical shapes are not quantitatively important and have not been pursued here.

Numerous studies have considered diffusion-limited transport to the cell surface at low nutrient concentrations saturating to a physiologically limited nutrient uptake at higher concentrations (Hill and Whittingham, 1955; Pasciak and Gavis, 1975; Mann and Lazier, 2006). The physiological limitation is typically considered using a Michaelis–Menten-type equation. Here, we simply consider the diffusion-limited uptake to be saturated by the filling up of reserves, $\left(\mathrm{1}-R^{*}\right)$. Thus, nutrient uptake is given by

$$\begin{array}{}\text{(31)}& J=\mathit{\psi }D{C}_{\mathrm{b}}\left(\mathrm{1}-R^{*}\right),\end{array}$$

where R^* is the normalised reserve of the nutrient being considered. As shown later when considering preferential ammonium uptake, under extreme limitation relative to other nutrients, R^* approaches 0, and uptake approaches the diffusion limitation.

5.1.3 Light capture and chlorophyll synthesis

Light absorption by microalgae cells has already been considered above (Eq. 6). The same absorption cross-section, α, is used to calculate the capture of photons per cell:

$$\begin{array}{}\text{(32)}& {\displaystyle \frac{\partial R_{\mathrm{C}}}{\partial t}}=\left(\mathrm{1}-R_{\mathrm{C}}^{*}\right){\displaystyle \frac{({\mathrm{10}}^{\mathrm{9}}hc{)}^{-\mathrm{1}}}{A_V}}\int {\mathit{\alpha }}_{\mathit{\lambda }}{E}_{o,\mathit{\lambda }}\mathit{\lambda }\mathrm{d}\mathit{\lambda },\end{array}$$

where $\left(\mathrm{1}-R_{\mathrm{C}}^{*}\right)$ accounts for the reduced capture of photons as the reserves becomes saturated, and $\frac{({\mathrm{10}}^{\mathrm{9}}hc{)}^{-\mathrm{1}}}{A_V}$ converts from energy to photons. The absorption cross-section is a function of intracellular pigment concentration, which is a dynamic variable determined below. While a drop-off of photosynthesis occurs as the carbon reserves become replete, this formulation does not consider photo-inhibition due to photo-oxidative stress, although it has been considered elsewhere for zooxanthellae (Baird et al., 2018).

The rate of synthesis of pigment is based on the incremental benefit of adding pigment to the rate of photosynthesis. This calculation includes both the reduced benefit when carbon reserves are replete, $(\mathrm{1}-R_{\mathrm{C}}^{*})$, and the reduced benefit due to self-shading, χ. The factor χ is calculated for the derivative of the absorption cross-section per unit projected area (see Eq. 6), α∕PA, with non-dimensional group ρ=γc_ir. For a sphere of radius r (Baird et al., 2013),

$$\begin{array}{}\text{(33)}& {\displaystyle \frac{\mathrm{1}}{\mathrm{PA}}}{\displaystyle \frac{\partial \mathit{\alpha }}{\partial \mathit{\rho }}}={\displaystyle \frac{\mathrm{1}-e^{-\mathrm{2}\mathit{\rho }}(\mathrm{2}{\mathit{\rho }}^{\mathrm{2}}+\mathrm{2}\mathit{\rho }+\mathrm{1})}{{\mathit{\rho }}^{\mathrm{3}}}}=\mathit{\chi },\end{array}$$

where χ represents the area-specific incremental rate of change of absorption with ρ. The rate of chlorophyll synthesis is given by

$$\begin{array}{}\text{(34)}& {\displaystyle \frac{\partial c_{\mathrm{i}}}{\partial t}}=k_{\mathrm{Chl}}^{\max }(\mathrm{1}-R_{\mathrm{C}}^{*})\stackrel{\mathrm{‾}}{\mathit{\chi }}\mathrm{if}\mathrm{C}:\mathrm{Chl}>{\mathit{\theta }}_{min},\end{array}$$

where $k_{\mathrm{Chl}}^{\max }$ is the maximum rate of synthesis, θ_min is the minimum C:Chl ratio, and the line in $\stackrel{\mathrm{‾}}{\mathit{\chi }}$ signifies the mean over the photosynthetically available radiation. At θ_min, pigment synthesis is zero. Both self-shading and the rate of photosynthesis itself are based on photon absorption rather than energy absorption (Table 5), as experimentally shown in Nielsen and Sakshaug (1993).

Table 5 Microalgae growth model equations. The term $B/m_{B,\mathrm{N}}$ is the concentration of cells. The equation for organic matter formation gives the stoichiometric constants: 12 g C mol C⁻¹; 32 g O mol O₂⁻¹. The equations are for scalar irradiance specified as an energy flux.

For each phytoplankton type, the model considers multiple pigments with distinct absorption spectra. The model needs to represent all photo-absorbing pigments as the C:Chl model calculates the pigment concentration based on that required to maximise photosynthesis. If only Chl a was represented, the model would predict a Chl a concentration that was accounting for the absorption of Chl a and the auxiliary pigments, thus overpredicting the Chl a concentration when compared to observations.

5.1.4 Carbon fixation/respiration

When photons are captured, there is an increase in reserves of carbon, $k_{\mathrm{I}}(\mathrm{1}-R_{\mathrm{C}}^{*})$ (Eq. 46), and an accompanying uptake of dissolved inorganic carbon, $\frac{\mathrm{106}}{\mathrm{1060}}\mathrm{12}{k}_{\mathrm{I}}(\mathrm{1}-R_{\mathrm{C}}^{*})$ (Eq. 42), and release of oxygen, $\frac{\mathrm{106}}{\mathrm{1060}}\mathrm{32}{k}_{\mathrm{I}}(\mathrm{1}-R_{\mathrm{C}}^{*})$ (Eq. 43), per cell to the water column (Table 5). Additionally, there is a basal respiration, representing a constant cost of cell maintenance. The loss of internal reserves, ${\mathit{\mu }}_B^{\max }{m}_{B,\mathrm{C}}\mathit{\varphi }{R}_{\mathrm{C}}^{*}$, results in a gain of water column dissolved inorganic carbon per cell, $\frac{\mathrm{106}}{\mathrm{1060}}\frac{\mathrm{12}}{\mathrm{14}}{\mathit{\mu }}_B^{\max }\mathit{\varphi }{R}_{\mathrm{C}}^{*}$, as well as an uptake of dissolved oxygen per cell, $\frac{\mathrm{106}}{\mathrm{1060}}\frac{\mathrm{32}}{\mathrm{14}}{\mathit{\mu }}_B^{\max }\mathit{\varphi }{R}_{\mathrm{C}}^{*}$ (Table 5). The loss in water column dissolved oxygen per cell represents an instantaneous respiration of the fixed carbon reserves. Basal respiration decreases internal reserves, and therefore growth rate, but does not directly lead to cell mortality at zero carbon reserves. Implicit in this scheme is that the basal cost is higher when the cell has more carbon reserves, $R_{\mathrm{C}}^{*}$.

A linear mortality term, resulting in the loss of structural material and carbon reserves, is considered later.

5.1.5 Application of single cell rates to a population

As mentioned above, the nutrient uptake and light absorption rates are calculated on a per cell basis. This has allowed geometric considerations to be explicitly used and contrasts with most biogeochemical models that formulate planktonic rates based on population interactions. However, the state variables for microalgae (and zooplankton) are for the population. Therefore, rates per cell need to be multiplied by the number of cells to obtain population rates. In the case of microalgae, the number of cells n is given by $B/m_{B,\mathrm{N}}$. This assumes all cells in the population are identical and that the state variable for the population, B, is quantifying only the nitrogen (or oxygen, carbon and phosphorus) associated with the structural material. It should also be noted that all cells in a population have the same state of their reserves.

5.1.6 Conservation of mass of microalgae model

The conservation of mass of cells containing structural material and reserves during transport, growth and mortality is established in Baird et al. (2004). Briefly, for microalgal growth, total concentration of nitrogen in microalgae cells is given by $B+B{R}_{\mathrm{N}}^{*}$. For conservation of mass, the time derivatives must equate to zero:

$$\begin{array}{}\text{(35)}& {\displaystyle \frac{\partial B}{\partial t}}+{\displaystyle \frac{\partial \left(R_{\mathrm{N}}B/R_{\mathrm{N}}^{max}\right)}{\partial t}}=\mathrm{0},\end{array}$$

using the product rule to differentiate the second term on the left-hand side:

$$\begin{array}{}\text{(36)}& {\displaystyle \frac{\partial B}{\partial t}}+{\displaystyle \frac{\partial B}{\partial t}}{\displaystyle \frac{R_{\mathrm{N}}}{R_{\mathrm{N}}^{max}}}+{\displaystyle \frac{B}{R_{\mathrm{N}}^{max}}}{\displaystyle \frac{\partial R_{\mathrm{N}}}{\partial t}}=\mathrm{0},\end{array}$$

where

$$\begin{array}{}\text{(37)}& {\displaystyle }& {\displaystyle \frac{\partial B}{\partial t}}=+{\mathit{\mu }}_B^{max}{R}_{\mathrm{C}}^{*}{R}_{\mathrm{N}}^{*}{R}_{\mathrm{P}}^{*}B\text{(38)}& {\displaystyle }& {\displaystyle \frac{\partial B}{\partial t}}{\displaystyle \frac{R_{\mathrm{N}}}{R_{\mathrm{N}}^{max}}}=+{\mathit{\mu }}_B^{max}{R}_{\mathrm{C}}^{*}{R}_{\mathrm{N}}^{*}{R}_{\mathrm{P}}^{*}B{\displaystyle \frac{R_{\mathrm{N}}}{R_{\mathrm{N}}^{max}}}\text{(39)}& {\displaystyle }& {\displaystyle \frac{B}{R_{\mathrm{N}}^{max}}}{\displaystyle \frac{\partial R_{\mathrm{N}}}{\partial t}}=-B(\mathrm{1}+R_{\mathrm{N}}^{*}){\mathit{\mu }}_B^{max}{R}_{\mathrm{C}}^{*}{R}_{\mathrm{N}}^{*}{R}_{\mathrm{P}}^{*}{\displaystyle \frac{R_{\mathrm{N}}}{R_{\mathrm{N}}^{max}}},\end{array}$$

thus demonstrating conservation of mass when $m_{B,\mathrm{N}}=R_{\mathrm{N}}^{max}$, as used here.

The state variables, equations and parameter values for microalgae growth are listed in Tables 4, 5 and 6, respectively. The equations in Table 5 described nitrogen uptake from the DIN pool. The partitioning of uptake between nitrate and ammonium due to preferential ammonium uptake is described in Sect. 9.1. Earlier published versions of the microalgae model are described with multiple nutrient limitation (Baird et al., 2001), with variable C:N ratios (Wild-Allen et al., 2010) and variable C:Chl ratios (Baird et al., 2013). Further, demonstration of the conservation of mass during transport is given in Baird et al. (2004).

Table 6 Constants and parameter values used in the microalgae model. V is cell volume in µm³.

^a Figs. 5, 6 and 7. ^b Straile (1997). ^c Finkel (2001), Sathyendranath et al. (2009) using high-performance liquid chromatography (HPLC) determination which isolates Chl a. ^d Li and Gregory (1974).

5.2 Nitrogen-fixing Trichodesmium

The growth of Trichodesmium follows the microalgae growth and C:Chl model above, with the following additional processes of nitrogen fixation and physiological-dependent buoyancy adjustment, as described in Robson et al. (2013). Additional parameter values for Trichodesmium are given in Table 7.

Table 7 Parameter values used in the Trichodesmium model (Robson et al., 2013).

5.2.1 Nitrogen fixation

Nitrogen fixation occurs when the DIN concentration falls below a critical concentration, DIN_crit, typically 0.3 to 1.6 µmol L⁻¹ (i.e. 4 to 20 mg N m⁻³; Robson et al., 2013), at which point Trichodesmium produces nitrogenase to allow fixation of N₂. It is assumed that nitrogenase becomes available whenever ambient DIN falls below the value of DIN_crit and carbon and phosphorus are available to support nitrogen uptake. The rate of change of internal reserves of nitrogen, R_N, due to nitrogen fixation if DIN<DIN_crit is given by

$$\begin{array}{}\text{(55)}& \begin{array}{rl}{N}_{\mathrm{fix}}=& {\displaystyle \frac{\partial R_{\mathrm{N}}}{\partial t}}{|}_{N_{\mathrm{fix}}}=\max \left(\mathrm{4}\mathit{\pi }r{D}_{{\mathrm{NO}}_{\mathrm{3}}}{\mathrm{DIN}}_{\mathrm{crit}}{R}_{\mathrm{P}}^{*}{R}_{\mathrm{C}}^{*}\right.\\ & (\mathrm{1}-R_{\mathrm{N}}^{*})-\mathrm{4}\mathit{\pi }r{D}_{{\mathrm{NO}}_{\mathrm{3}}}\left[{\mathrm{NO}}_{\mathrm{3}}+{\mathrm{NH}}_{\mathrm{4}}\right]\\ & \left.(\mathrm{1}-R_{\mathrm{N}}^{*}),\mathrm{0}\right),\end{array}\end{array}$$

where N_fix is the rate of nitrogen fixation per cell and r is the radius of the individual cell. Using this formulation, Trichodesmium is able to maintain its nitrogen uptake rate at that achieved through diffusion-limited uptake at DIN_crit even when DIN drops below DIN_crit, provided phosphorus and carbon reserves, $R_{\mathrm{P}}^{*}$ and $R_{\mathrm{C}}^{*}$, respectively, are available.

The energetic cost of nitrogen fixation is represented as a fixed proportion of carbon fixation, f_Nfix, equivalent to a reduction in quantum efficiency, and as a proportion, f_nitrogenase, of the nitrogen fixed:

$$\begin{array}{}\text{(56)}& {\displaystyle \frac{\partial R_{\mathrm{C}}}{\partial t}}=-(\mathrm{1}-f_{\mathrm{Nfix}})(\mathrm{1}-f_{\mathrm{nitrogenase}})k_{\mathrm{I}},\end{array}$$

where k_I is the rate of photon absorption per cell obtain from the microalgal growth model (Table 5).

5.2.2 Buoyancy adjustment

The rate of change of Trichodesmium biomass, B, as a result of density difference between the cell and the water, is approximated by Stokes' law:

$$\begin{array}{}\text{(57)}& {\displaystyle \frac{\partial B}{\partial t}}=-{\displaystyle \frac{\mathrm{2}}{\mathrm{9}}}{\displaystyle \frac{g{r}_{\mathrm{col}}^{\mathrm{2}}}{\mathit{\mu }}}\left(\mathit{\rho }-{\mathit{\rho }}_{\mathrm{w}}\right){\displaystyle \frac{\partial B}{\partial z}},\end{array}$$

where z is the distance in the vertical (positive up), μ is the dynamic viscosity of water, g is acceleration due to gravity, r_col is the equivalent spherical radius of the sinking mass representing a colony radius, ρ_w is the density of water, and ρ is the cell density given by

$$\begin{array}{}\text{(58)}& \mathit{\rho }={\mathit{\rho }}_{min}+R_{\mathrm{C}}^{*}\left({\mathit{\rho }}_{max}-{\mathit{\rho }}_{min}\right),\end{array}$$

where $R_{\mathrm{C}}^{*}$ is the normalised carbon reserves of the cell (see above), and ρ_min and ρ_max are the densities of the cell when there is no carbon reserves and full carbon reserves, respectively. Thus, when light reserves are depleted, the cell is more buoyant, facilitating the retention of Trichodesmium in the surface waters.

5.3 Water column inorganic chemistry

5.3.1 Carbon chemistry

The major pools of dissolved inorganic carbon species in the ocean are ${\mathrm{HCO}}_{\mathrm{3}}^{-}$, ${\mathrm{CO}}_{\mathrm{3}}^{-}$ and dissolved CO₂, which influence the speciation of H⁺ and OH⁻ ions, and therefore pH. The interaction of these ions reaches an equilibrium in seawater within a few tens of seconds (Zeebe and Wolf-Gladrow, 2001). In the BGC model here, where calculation time steps are of the order of tens of minutes, it is reasonable to assume that the carbon chemistry system is at equilibrium.

The Ocean Carbon-Cycle Model Intercomparison Project (OCMIP) has developed numerical methods to quantify air–sea carbon fluxes and carbon dioxide system equilibria (Najjar and Orr, 1999). Here, we use a modified version of the OCMIP-2 Fortran code developed for MOM4 (Geophysical Fluid Dynamics Laboratory Modular Ocean Model version 4; Griffies et al., 2004). The OCMIP procedures quantify the state of the carbon dioxide (CO₂) system using two prognostic variables, the concentration of dissolved inorganic carbon, DIC, and total alkalinity, A_T. The value of these prognostic variables, along with salinity and temperature, are used to calculate the pH and partial pressure of carbon dioxide, pCO₂, in the surface waters using a set of governing chemical equations which are solved using a Newton–Raphson method (Najjar and Orr, 1999).

We altered the OCMIP scheme by increasing the search space for the iterative scheme from ±0.5 pH units (appropriate for global models) to ±2.5. With this change, the scheme converges over a broad range of DIC and A_T values (Munhoven, 2013). For more details, see Mongin and Baird (2014) and Mongin et al. (2016b).

5.3.2 Nitrification

Nitrification is the oxidation of ammonium to form nitrite followed by the rapid oxidation of nitrite to nitrate. This is represented in a one-step process, with the rate of nitrification given by

$$\begin{array}{}\text{(65)}& {\displaystyle \frac{\partial \left[{\mathrm{NH}}_{\mathrm{4}}\right]}{\partial t}}=-{\mathit{\tau }}_{\mathrm{nit},\mathrm{wc}}\left[{\mathrm{NH}}_{\mathrm{4}}\right]{\displaystyle \frac{\left[{\mathrm{O}}_{\mathrm{2}}\right]}{K_{\mathrm{nit},\mathrm{O}}+\left[{\mathrm{O}}_{\mathrm{2}}\right]}},\end{array}$$

where the equations and parameter values are defined in Tables 9 and 10.

Table 8 State and derived variables for the water column inorganic chemistry model.

Table 9 Equations for the water column inorganic chemistry.

Table 10 Constants and parameter values used in the water column inorganic chemistry.

5.3.3 Phosphorus absorption–desorption

The rate of phosphorus desorption from particulates is given by

$$\begin{array}{}\text{(66)}& \begin{array}{rl}{\displaystyle \frac{\partial P}{\partial t}}=& {\mathit{\tau }}_{\mathrm{Pabs}}\left({\displaystyle \frac{\mathrm{PIP}}{k_{\mathrm{Pads},\mathrm{wc}}\mathrm{NAP}}}-{\displaystyle \frac{\left[{\mathrm{O}}_{\mathrm{2}}\right]P}{K_{{\mathrm{O}}_{\mathrm{2}},\mathrm{abs}}+\left[{\mathrm{O}}_{\mathrm{2}}\right]}}\right)\\ & =-{\displaystyle \frac{\partial \mathrm{PIP}}{\partial t}},\end{array}\end{array}$$

where [O₂] is the concentration of oxygen, P is the concentration of dissolved inorganic phosphorus, PIP is the concentration of particulate inorganic phosphorus, NAP is the sum of the non-algal inorganic particulate concentrations, and τ_Pabs, k_Pads,wc and $K_{{\mathrm{O}}_{\mathrm{2}},\mathrm{abs}}$ are model parameters described in Table 10.

At steady state, the PIP concentration is given by

$$\begin{array}{}\text{(67)}& \mathrm{PIP}=k_{\mathrm{Pads},\mathrm{wc}}P{\displaystyle \frac{\left[{\mathrm{O}}_{\mathrm{2}}\right]}{K_{{\mathrm{O}}_{\mathrm{2}},\mathrm{abs}}+\left[{\mathrm{O}}_{\mathrm{2}}\right]}}\mathrm{NAP}.\end{array}$$

As an example for rivers flowing into the GBR configuration, [O₂]=7411 mg m⁻³ (90 % saturation at T=25, S=0), NAP =0.231 kg m⁻³, $k_{\mathrm{Pads},\mathrm{wc}}=\mathrm{30}$ kg NAP⁻¹, $K_{{\mathrm{O}}_{\mathrm{2}},\mathrm{abs}}=\mathrm{74}$ mg O m⁻³, and P=4.2 mg m⁻³; thus, the ratio PIP ∕ DIP =6.86 (see Fig. 12). Limited available observations of absorption–desorption include from the Johnstone River (Pailles and Moody, 1992) and the GBR (Monbet et al., 2007).

Figure 12 Phosphorus adsorption–desorption equilibria, $K_{{\mathrm{O}}_{\mathrm{2}},\mathrm{abs}}=\mathrm{74}$ mg O m⁻³ at two values of K_Pads,wc.

5.4 Zooplankton herbivory

In the simple food web of the model, herbivory involves small zooplankton consuming small phytoplankton, and large zooplankton consuming large phytoplankton, microphytobenthos and Trichodesmium. For simplicity, the state variables and equations are only given for small plankton grazing (Tables 11, 13), but the parameters are given for all grazing terms (Table 12).

Table 11 State and derived variables for the zooplankton grazing. Zooplankton cell mass, $m_{\mathrm{Z}}=\mathrm{16}\mathrm{000}\times \mathrm{14.01}\times \mathrm{10.5}{V}_{\mathrm{Z}}$ mg N cell⁻¹, where V_Z is the volume of zooplankton (Hansen et al., 1997).

Table 12 Constants and parameter values used for zooplankton grazing. Dissipation rate of turbulent kinetic energy (TKE) is considered constant.

Table 13 Equations for zooplankton grazing. The terms represent a predator Z consuming a phytoplankton B. (1) If the zooplankton diet contains multiple phytoplankton classes, and grazing is prey saturated, then phytoplankton loss must be reduced to account for the saturation by other types of microalgae. (2) $\frac{Z}{m_{\mathrm{Z}}}$ is the number of individual zooplankton. (3) Phytoplankton pigment is lost to water column without being conserved. Chl a has a chemical formulae of C₅₅H₇₂O₅N₄Mg and a molecular weight of 893.49 g mol⁻¹. The uptake (and subsequent remineralisation) of molecules for chlorophyll synthesis could make up a maximum (at $\mathrm{C}:\mathrm{Chl}=\mathrm{20}$) of $(\mathrm{660}/\mathrm{893})/\mathrm{20}$ and $(\mathrm{56}/\mathrm{893})/\mathrm{20}\times (\mathrm{16}/\mathrm{106})\times (\mathrm{14}/\mathrm{12}))$, or ∼4 % and ∼2 % of the exchange of C and N between the cell and water column, and will cancel out over the lifetime of a cell. Thus, the error in ignoring chlorophyll loss to the water column is small.

The rate of zooplankton grazing is determined by the encounter rate of the predator and all its prey up until the point at which it saturates the growth of the zooplankton (Eq. 75), and then it is constant. This is a Holling type I grazing response (Gentleman, 2002). Under the condition of multiple prey types, there is no preferential grazing other than that determined by the chance of encounter. The encounter rate is the result of the relative motion of individuals brought about by diffusion (Eq. 77), swimming (Eq. 78) and shear (Eq. 79) determined relative velocities (Eq. 80) (Jackson, 1995; Baird, 2003). One particular advantage of formulating the encounter rate on individuals is that should the number of populations considered in the model change (i.e. an additional phytoplankton class is added), there is no need for empirical coefficients in the model to change. More recent uses of encounter-based grazing functions are described in Flynn and Mitra (2016).

Unlike the microalgae, zooplankton does not contain reserves of nutrients and fixed carbon, and therefore has a fixed stoichiometry of the Redfield ratio. As the zooplankton are grazing on the phytoplankton that contain internal reserves of nutrients, an additional flux of dissolved inorganic nutrients ($g{R}_{\mathrm{N}}^{*}$ for nitrogen) is returned to the water column (for more details, see Sect. 5.4.1).

5.4.1 Conservation of mass in zooplankton grazing

It is important to note that the microalgae model presented above represents internal reserves of nutrients, carbon and chlorophyll as a per cell quantity. Using this representation, there are no losses of internal quantities with either grazing or mortality. However, the implication of their presence is represented in the $\left(g{R}_{\mathrm{N}}^{*}\right)$ terms (Table 13) that return the reserves to the water column. These terms represent the fast return of a fraction of phytoplankton nitrogen due to processes like “sloppy eating”.

An alternative and equivalent formulation would be to consider total concentration of microalgal reserves in the water column, then the change in water column concentration of reserves due to mortality (either grazing or natural mortality) must be considered. This alternate representation will not be undertaken here as the above considered equations are fully consistent, but it is worth noting that the numerical solution of the model within the code represents total water column concentrations of internal reserves and therefore must include the appropriate loss terms due to mortality.

5.5 Zooplankton carnivory

Large zooplankton consume small zooplankton. This process uses similar encounter rate and consumption rate limitations calculated for zooplankton herbivory (Table 13). As zooplankton contain no internal reserves, the equations are simplified from the herbivory case to those listed in Table 14. Assuming that the efficiency of herbivory, γ, is equal to that of carnivory, and therefore assigned the same parameter, the additional process of carnivory adds no new parameters to the biogeochemical model.

Table 14 Equations for zooplankton carnivory, representing large zooplankton Z_L consuming small zooplankton Z_S. The parameters values and symbols are given in Tables 11 and 12.

5.6 Zooplankton respiration

In the model, there is no change in water column oxygen concentration if organic material is exchanged between pools with the same elemental ratio. Thus, when zooplankton consume phytoplankton, no oxygen is consumed due to the consumption of phytoplankton structural material (B_P). However, the excess carbon reserves represent a pool of fixed carbon, which when released from the phytoplankton must consume oxygen. Further, zooplankton mortality and growth inefficiency result in detritial production, which when remineralised consumes oxygen. Additionally, carbon released to the dissolved inorganic pool during inefficiency grazing on phytoplankton structural material also consumes oxygen. Thus, zooplankton respiration is implicitly captured in these associated processes.

5.7 Non-grazing plankton mortality

The rate of change of plankton biomass, B, as a result of natural mortality is given by

$$\begin{array}{}\text{(89)}& {\displaystyle \frac{\partial B}{\partial t}}=-m_{\mathrm{L}}B-m_{\mathrm{Q}}{B}^{\mathrm{2}},\end{array}$$

where m_L is the linear mortality coefficient and m_Q is the quadratic mortality coefficient.

A combination of linear and quadratic mortality rates are used in the model. When the mortality term is the sole loss term, such as zooplankton in the water column or benthic microalgae in the sediments, a quadratic term is employed to represent increasing predation/viral disease losses in dense populations. For suspended microalgae, we have used only a linear term (i.e. m_Q=0). Linear terms have been used to represent a basal respiration rate.

As described in Sect. 5.1.6, the mortality terms need to account for the internal properties of lost microalgae. For definitions of the state variables, see Tables 15 and 16.

Table 15 Constants and parameter values used for plankton mortality.

Table 16 Equations for linear phytoplankton mortality.

Table 17 Equations for the zooplankton mortality. f_Z2det is the fraction of zooplankton mortality that is remineralised and is equal to 0.5 for both small and large zooplankton.

5.8 Air–sea gas exchange

Air–sea gas exchange is calculated using wind speed (we choose a cubic relationship; Wanninkhof and McGillis, 1999), saturation state of the gas (described below) and the Schmidt number of the gas (Wanninkhof, 1992). The transfer coefficient, k, is given by

$$\begin{array}{}\text{(100)}& k={\displaystyle \frac{\mathrm{0.0283}}{\mathrm{360}\mathrm{000}}}{u}_{\mathrm{10}}^{\mathrm{3}}{\left(Sc/\mathrm{660}\right)}^{-\mathrm{1}/\mathrm{2}},\end{array}$$

where 0.0283 cm h⁻¹ is an empirically determined constant (Wanninkhof and McGillis, 1999), u₁₀ is the short-term steady wind at 10 m above the sea surface (m s⁻¹), the Schmidt number, Sc, is the ratio of the diffusivity of momentum and that of the exchanging gas, and is given by a cubic temperature relationship (Wanninkhof, 1992). Finally, a conversion factor of 360 000 m s⁻¹ (cm h⁻¹)⁻¹ is used.

In practice, the hydrodynamic model can contain thin surface layers as the surface elevation moves between z levels. Further, physical processes of advection and diffusion and gas fluxes are done sequentially, allowing concentrations to build up through a single time step. To avoid unrealistic changes in the concentration of gases in thin surface layers, the shallowest layer thicker than 20 cm receives all the surface fluxes.

5.8.1 Oxygen

The saturation state of oxygen [O₂]_sat is determined as a function of temperature and salinity following Weiss (1970). The change in concentration of oxygen in the surface layer due to a sea–air oxygen flux (positive from sea to air) is given by

$$\begin{array}{}\text{(101)}& {\displaystyle \frac{\partial \left[{\mathrm{O}}_{\mathrm{2}}\right]}{\partial t}}=k_{{\mathrm{O}}_{\mathrm{2}}}\left([{\mathrm{O}}_{\mathrm{2}}{]}_{\mathrm{sat}}-[{\mathrm{O}}_{\mathrm{2}}]\right)/h,\end{array}$$

where $k_{{\mathrm{O}}_{\mathrm{2}}}$ is the transfer coefficient for oxygen (Eq. 100), [O₂] is the dissolved oxygen concentration in the surface waters, and h is the thickness of the surface layer of the model into which sea–air flux flows.

5.8.2 Carbon dioxide

The change in surface dissolved inorganic carbon concentration, DIC, resulting from the sea–air flux (positive from sea to air) of carbon dioxide is given by

$$\begin{array}{}\text{(102)}& {\displaystyle \frac{\partial \mathrm{DIC}}{\partial t}}=k_{{\mathrm{CO}}_{\mathrm{2}}}\left([{\mathrm{CO}}_{\mathrm{2}}{]}_{\mathrm{atm}}-[{\mathrm{CO}}_{\mathrm{2}}]\right)/h,\end{array}$$

where $k_{{\mathrm{CO}}_{\mathrm{2}}}$ the transfer coefficient for carbon dioxide (Eq. 100), [CO₂] is the dissolved carbon dioxide concentration in the surface waters determined from DIC and A_T using the carbon chemistry equilibria calculations described in Sect. 5.3.1, [CO₂]_atm is the partial pressure of carbon dioxide in the atmosphere, and h is the thickness of the surface layer of the model into which sea–air flux flows.

Note the carbon dioxide flux is determined not by the gradient in DIC but the gradient in [CO₂]. At pH values around 8, [CO₂] makes up only approximately 1∕200 of DIC in seawater, significantly reducing the air–sea exchange. Counteracting this reduced gradient, note that changing DIC results in an approximately 10-fold change in [CO₂] (quantified by the Revelle factor; Zeebe and Wolf-Gladrow, 2001). Thus, the gas exchange of CO₂ is approximately $\mathrm{1}/\mathrm{200}\times \mathrm{10}=\mathrm{1}/\mathrm{20}$ of the oxygen flux for the same proportional perturbation in DIC and oxygen. At a Sc number of 524 (25 ^∘C seawater) and a wind speed of 12 m s⁻¹, 1 m of water equilibrates with CO₂ in the atmosphere with an e-folding timescale of approximately 1 d.

6 Epibenthic processes

In the model, benthic communities are quantified as a biomass per unit area or areal biomass. At low biomass, the community is composed of a few specimens spread over a small fraction of the bottom, with no interaction between the nutrient and energy acquisition of individuals. Thus, at low biomass, the areal fluxes are a linear function of the biomass.

As biomass increases, the individuals begin to cover a significant fraction of the bottom. For nutrient and light fluxes that are constant per unit area, such as downwelling irradiance and sediment releases, the flux per unit biomass decreases with increasing biomass. Some processes, such as photosynthesis in a thick seagrass meadow or nutrient uptake by a coral reef, become independent of biomass (Atkinson, 1992) as the bottom becomes completely covered. To capture the non-linear effect of biomass on benthic processes, we use an effective projected area fraction, A_eff.

To restate, at low biomass, the area on the bottom covered by the benthic community is a linear function of biomass. As the total leaf area approaches and exceeds the projected area, the projected area for the calculation of water-community exchange approaches 1 and becomes independent of biomass. This is represented using

$$\begin{array}{}\text{(103)}& A_{\mathrm{eff}}=\mathrm{1}-\exp (-{\mathrm{\Omega }}_{B}B),\end{array}$$

where A_eff is the effective projected area fraction of the benthic community (m² m⁻²), B is the biomass of the benthic community (g N m⁻²), and Ω_B is the nitrogen-specific leaf area coefficient (m² g N⁻¹). For further explanation of Ω_B, see Baird et al. (2016a).

The parameter Ω_B is critical: it provides a means of converting between biomass and fractions of the bottom covered and is used in calculating the absorption cross-section of the leaf and the nutrient uptake of corals and macroalgae. That Ω_B has a simple physical explanation and can be determined from commonly undertaken morphological measurement (see below), gives us confidence in its use throughout the model.

6.1 Macroalgae

The macroalgae model considers the diffusion-limited supply of dissolved inorganic nutrients (N and P) and the absorption of light, delivering N, P and fixed C, respectively. Unlike the microalgae model, no internal reserves are considered, implying that the macroalgae has a fixed stoichiometry that can be specified as

$$\begin{array}{}\text{(104)}& \begin{array}{rl}\mathrm{550}{\mathrm{CO}}_{\mathrm{2}}& +\mathrm{30}{\mathrm{NO}}_{\mathrm{3}}^{-}+{\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}+\mathrm{792}{\mathrm{H}}_{\mathrm{2}}\mathrm{O}\stackrel{\mathrm{5500}\mathrm{photons}}{⟶}\\ & \left({\mathrm{CH}}_{\mathrm{2}}\mathrm{O}{)}_{\mathrm{550}}\right({\mathrm{NH}}_{\mathrm{3}}{)}_{\mathrm{30}}{\mathrm{H}}_{\mathrm{3}}{\mathrm{PO}}_{\mathrm{4}}+\mathrm{716}{\mathrm{O}}_{\mathrm{2}},\end{array}\end{array}$$

where the stoichiometry is based on Atkinson and Smith (1983) (see also Baird and Middleton, 2004; Hadley et al., 2015a, b). Note that when ammonium is taken up instead of nitrate there is a slightly different O₂ balance (Sect. 9.1). In the next section, we will consider the maximum nutrient uptake and light absorption and then bring them together to determine the realised growth rate.

6.1.1 Nutrient uptake

Nutrient uptake by macroalgae is a function of nutrient concentration, water motion (Hurd, 2000) and internal physiology. The maximum flux of nutrients is specified as a mass transfer limit per projected area of macroalgae and is given by (Falter et al., 2004; Zhang et al., 2011)

$$\begin{array}{}\text{(105)}& S_x=\mathrm{2850}{\left({\displaystyle \frac{\mathrm{2}\mathit{\tau }}{\mathit{\rho }}}\right)}^{\mathrm{0.38}}S{c_x}^{-\mathrm{0.6}},S{c}_x={\displaystyle \frac{\mathit{\nu }}{D_x}},\end{array}$$

where S_x is the mass transfer rate coefficient of element x= N, P, τ is the shear stress on the bottom, ρ is the density of water, and Sc_x is the Schmidt number. The Schmidt number is the ratio of the diffusivity of momentum, ν, and mass, D_x (Table 6), and varies with temperature, salinity and nutrient species. The rate constant S can be thought of as the height of water cleared of mass per unit of time by the water–macroalgae exchange.

6.1.2 Light capture

The calculation of light capture by macroalgae involves estimating the fraction of light that is incident upon the leaves and the fraction that is absorbed. The rate of photon capture is given by

$$\begin{array}{}\text{(106)}& \begin{array}{rl}{k}_{\mathrm{I}}=& {\displaystyle \frac{{\left({\mathrm{10}}^{\mathrm{9}}hc\right)}^{-\mathrm{1}}}{A_V}}\int E_{\mathrm{d},\mathit{\lambda }}\\ & \left(\mathrm{1}-\exp \left(-A_{\mathrm{L},\mathit{\lambda }}{\mathrm{\Omega }}_{\mathrm{MA}}\mathrm{MA}\right)\right)\mathit{\lambda }\mathrm{d}\mathit{\lambda },\end{array}\end{array}$$

where h, c and A_V are fundamental constants, 10⁹ nm m⁻¹ accounts for the typical representation of wavelength, λ, in nm, and A_L,λ is the spectrally resolved leaf-specific absorptance. As shown in Eq. (103), the term $\mathrm{1}-\exp \left(-{\mathrm{\Omega }}_{\mathrm{MA}}\mathrm{MA}\right)$ gives the effective projected area fraction of the community. In the case of light absorption of macroalgae, the exponent is multiplied by the leaf-specific absorptance, A_L,λ, to account for the transparency of the leaves. At low macroalgae biomass, absorption at wavelength λ is equal to $E_{\mathrm{d},\mathit{\lambda }}{A}_{\mathrm{L},\mathit{\lambda }}{\mathrm{\Omega }}_{\mathrm{MA}}\mathrm{MA}$, increasing linearly with biomass as all leaves at low biomass are exposed to full light (i.e. there is no self-shading). At high biomass, the absorption by the community asymptotes to E_d,λ, at which point increasing biomass does not increase the absorption as all light is already absorbed.

For more details on the calculation of Ω_MA, see Baird et al. (2016a).

6.1.3 Growth

The growth rate combines nutrient, light and maximum organic matter synthesis rates following

$$\begin{array}{}\text{(107)}& \begin{array}{rl}{\mathit{\mu }}_{\mathrm{MA}}=& \min \left[{\mathit{\mu }}_{\mathrm{MA}}^{max},{\displaystyle \frac{\mathrm{30}}{\mathrm{5500}}}\mathrm{14}{\displaystyle \frac{k_{\mathrm{I}}}{\mathrm{MA}}},{\displaystyle \frac{S_{\mathrm{N}}{A}_{\mathrm{eff}}N}{\mathrm{MA}}},\right.\\ & \left.{\displaystyle \frac{\mathrm{30}}{\mathrm{1}}}{\displaystyle \frac{\mathrm{14}}{\mathrm{31}}}{\displaystyle \frac{S_{\mathrm{P}}{A}_{\mathrm{eff}}P}{\mathrm{MA}}}\right],\end{array}\end{array}$$

and the production of macroalgae is given by μ_MAMA. We have used the commonly applied multiple minimum function (von Liebig, 1840), although it is noted that others use the multiple of limitation terms (Fasham, 1993). The microalgae model described above uses dynamical reserves to determine the growth rate. The growth approximated using dynamical reserves closer approximates a multiple minimum function than a multiple of minimum terms, so it was deemed more appropriate to use a multiple minimum function for macroalgae and seagrass for which internal reserves were not resolved.

The maximum growth rate is inside the minimum operator. This allows the growth of macroalgae to be independent of temperature at low light but still have an exponential dependence at maximum growth rates (Baird et al., 2003).

6.1.4 Mortality

Mortality is defined as a simple linear function of biomass:

$$\begin{array}{}\text{(119)}& {\displaystyle \frac{\partial \mathrm{MA}}{\partial t}}=-{\mathit{\zeta }}_{\mathrm{MA}}\mathrm{MA}.\end{array}$$

A quadratic formulation is not necessary as both the nutrient and light capture rates become independent of biomass as $\mathrm{MA}\gg \mathrm{1}/{\mathrm{\Omega }}_{\mathrm{MA}}$. Thus, the steady-state biomass of macroalgae under nutrient limitation is given by

$$\begin{array}{}\text{(120)}& {\left(\mathrm{MA}\right)}_{\mathrm{SS}}={\displaystyle \frac{S_{\mathrm{N}}{A}_{\mathrm{eff}}N}{\mathit{\zeta }}},\end{array}$$

and for light-limited growth by

$$\begin{array}{}\text{(121)}& {\left(\mathrm{MA}\right)}_{\mathrm{SS}}={\displaystyle \frac{k_{\mathrm{I}}}{\mathit{\zeta }}}.\end{array}$$

The full macroalgae variables, equations and parameters are listed in Tables 18, 19 and 20.

Table 18 State and derived variables for the macroalgae model. For simplicity, in the equations, all dissolved constituents are given in grams, although elsewhere they are shown in milligrams.

Table 19 Equations for the macroalgae model. Other constants and parameters are defined in Table 20: 14 g N mol N⁻¹; 12 g C mol C⁻¹; 31 g P mol P⁻¹; 32 g O mol O₂⁻¹. Uptake shown here is for nitrate; see Sect. 9.1 for ammonium uptake.

Table 20 Constants and parameter values used to model macroalgae.

^* Spectrally resolved values.

6.2 Seagrass

Seagrasses are quantified per m² with a constant stoichiometry ($\mathrm{C}:\mathrm{N}:\mathrm{P}=\mathrm{550}:\mathrm{30}:\mathrm{1}$) for both above-ground, SG_A, and below-ground, SG_B, biomass and can translocate organic matter at this constant stoichiometry between the two stores of biomass. Growth occurs only in the above-ground biomass, but losses (grazing, decay, etc.) occur in both. Multiple seagrass varieties are represented. The varieties are modelled using the same equations for growth, respiration and mortality but with different parameter values. Here, we just list the variables, equations and parameters (Tables 21, 22 and 23) for the seagrass submodel. A description of the seagrass processes of growth, translocation between roots and leaves, and mortality has been published in Baird et al. (2016a), along with a comparison to observations from Gladstone Harbour on the northeast Australian coast.

Table 21 State and derived variables for the seagrass model. For simplicity, in the equations, all dissolved constituents are given in grams, although elsewhere they are shown in milligrams. The bottom water column thickness is spatially variable depending on bathymetry.

Table 22 Equations for the seagrass model. Other constants and parameters are defined in Table 23. The equation for organic matter formation gives the stoichiometric constants; 14 g N mol N⁻¹; 12 g C mol C⁻¹; 31 g P mol P⁻¹; 32 g O mol O₂⁻¹.

Table 23 Constants and parameter values used to model seagrass. ^a Factor of 2 for nighttime, factor of 2 for roots. ^b Zostera is calculated from leaf characteristics in Kemp et al. (1987); Hansen et al. (2000); Halophila ovalis is calculated from leaf dimensions in Vermaat et al. (1995); Ω_SG can also be determined from specific leaf area such as determined in Cambridge and Lambers (1998) for nine Australian seagrass species. ^c Spectrally resolved values in Baird et al. (2016a). ^d Duarte and Chiscano (1999). ^e Loosely based on Kaldy et al. (2013). ^f Thalassia testudinum (Gras et al., 2003). ^g Thalassia testudinum (Lee and Dunton, 1999). ^h Chartrand et al. (2012); Longstaff (2003); Chartrand et al. (2017). ⁱ Roberts (1993).

6.3 Coral polyps

The coral polyp parameterisation consists of a microalgae growth model to represent zooxanthellae growth based on Baird et al. (2013), and the parameterisation of coral–zooxanthellae interaction based on the host–symbiont model of Gustafsson et al. (2013), and the photo-adaptation, photo-inhibition and reaction centre dynamics model of Baird et al. (2018). The extra detail on the zooxanthellae photosystem is required due to its important role in thermal-stress-driven coral bleaching (Yonge, 1930; Suggett et al., 2008).

6.3.1 Coral host, symbiont and the environment

The state variables for the coral polyp model (Table 24) include the biomass of coral tissue, CH (g N m⁻²), and the structure material of the zooxanthellae cells, CS (mg N m⁻²). The structure material of the zooxanthellae, CS, in addition to nitrogen, contains carbon and phosphorus at the Redfield ratio. The zooxanthellae cells also contain reserves of nitrogen, R_N (mg N m⁻²), phosphorus, R_P (mg P m⁻²), and carbon, R_C (mg C m⁻²).

Table 24 Model state variables for the coral polyp model. Note that water column variables are three-dimensional, benthic variables are two-dimensional, and unnormalised reserves are per cell.

The zooxanthellae light absorption capability is resolved by considering the time-varying concentrations of pigments chlorophyll a, Chl, diadinoxanthin, X_p, and diatoxanthin, X_h, for which the state variable represents the areal concentration. A further three pigments, chlorophyll c₂, peridinin and β carotene, are considered in the absorption calculations, but their concentrations are in fixed ratios to chlorophyll a. Exchanges between the coral community and the overlying water can alter the water column concentrations of dissolved inorganic carbon, DIC, nitrogen, N, and phosphorus, P, as well as particulate phytoplankton, B, zooplankton, Z, and detritus, D, where multiple nitrogen, plankton and detritus types are resolved (Table 24).

The coral host is able to assimilate particulate organic nitrogen either through translocation from the zooxanthellae cells or through the capture of water column organic detritus and/or plankton. The zooxanthellae varies its intracellular pigment content depending on potential light limitation of growth, and the incremental benefit of adding pigment, allowing for the package effect (Baird et al., 2013). The coral tissue is assumed to have a Redfield $\mathrm{C}:\mathrm{N}:\mathrm{P}$ stoichiometry (Redfield et al., 1963), as shown by Muller-Parker et al. (1994). The zooxanthellae are modelled with variable $\mathrm{C}:\mathrm{N}:\mathrm{P}$ ratios (Muller-Parker et al., 1994), based on a structure material at the Redfield ratio but with variable internal reserves. The fluxes of C, N and P with the overlying water column (nutrient uptake and detritial/mucus release) can therefore vary from the Redfield ratio.

Here, we present the state variables (Table 24), derived variables (Table 25), equations (Tables 26, 27, 28 and 29) and parameters values (Tables 30 and 31) for the coral submodel. The description of the coral processes has been published in Baird et al. (2018), along with a comparison to observations from the Great Barrier Reef on the northeast Australian coast.

Table 25 Derived variables for the coral polyp model.

Table 26 Equations for the interactions of coral host, symbiont and environment excluding bleaching loss terms that appear in Table 29. The term $\mathrm{CS}/m_{B,\mathrm{N}}$ is the concentration of zooxanthellae cells. The equation for organic matter formation gives the stoichiometric constants; 12 g C mol C⁻¹; 32 g O mol O₂⁻¹.

Table 27 Equations for the coral polyp model. The term $\mathrm{CS}/m_{B,\mathrm{N}}$ is the concentration of zooxanthellae cells. The equation for organic matter formation gives the stoichiometric constants; 12 g C mol C⁻¹; 32 g O mol O₂⁻¹. Other constants and parameters are defined in Table 31.

Table 28 Equations for symbiont reaction centre dynamics. Bleaching loss terms appear in Table 29.

Table 29 Equations describing the expulsion of zooxanthellae and the resulting release of inorganic and organic molecules into the bottom water column layer.

Table 30 Constants and parameter values used to model coral polyps. V is zooxanthellae cell volume in µm³.

^a Baird et al. (2016a). ^b Straile (1997). ^c Redfield et al. (1963) and Kirk (1994). ^d Finkel (2001). ^e Ribes and Atkinson (2007); Wyatt et al. (2010). ^f, g Gustafsson et al. (2013, 2014).

Table 31 Constants and parameter values used in the coral bleaching model. “RCII” indicates the reaction centre of photosystem II.

^a In Suggett et al. (2009). ^b Ratio of constant terms in multivariate analysis in Hochberg et al. (2006). ^c Fitted parameter based on the existence of non-bleaching threshold (Suggett et al., 2009) and a comparison of observed bleaching and model output in the ∼1 km model.

6.3.2 Coral calcification

The rate of coral calcification is a function of the water column aragonite saturation, Ω_a, and the normalised reserves of fixed carbon in the symbiont, $R_{\mathrm{C}}^{*}$. The rates of change of DIC and total alkalinity, A_T, in the bottom water column layer of thickness h_wc due to calcification becomes

$$\begin{array}{}\text{(184)}& {\displaystyle }& {\displaystyle \frac{\partial \mathrm{DIC}}{\partial t}}=-\mathrm{12}g{A}_{\mathrm{eff}}/h_{\mathrm{wc}}\text{(185)}& {\displaystyle }& {\displaystyle \frac{\partial A_{\mathrm{T}}}{\partial t}}=-\mathrm{2}g{A}_{\mathrm{eff}}/h_{\mathrm{wc}}\text{(186)}& {\displaystyle }& {\displaystyle }g=k_{\mathrm{day}}\left({\mathrm{\Omega }}_{\mathrm{a}}-\mathrm{1}\right)(R_{\mathrm{C}}^{*}{)}^{\mathrm{2}}+k_{\mathrm{night}}\left({\mathrm{\Omega }}_{\mathrm{a}}-\mathrm{1}\right),\end{array}$$

where g is the rate of net calcification, k_day and k_night are defined in Table 30 with habitat-specific values (Anthony et al., 2011; Mongin and Baird, 2014). The fluxes are scaled by the effective projected area of the community, A_eff. The power of 2 for $R_{\mathrm{C}}^{*}$ ensures that generally light-replete symbionts provide the host with sufficient energy for calcification.

Table 32 Equations for coral polyp calcification and dissolution. The concentration of carbonate ions, $\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]$, is determined from equilibrium carbon chemistry as a function of A_T, DIC, temperature and salinity, and the concentration of calcium ions, [Ca²⁺], is a mean oceanic value; 12 g C mol C⁻¹. Other constants and parameters are defined in Table 30.

6.3.3 Dissolution of shelf carbonate sands

In addition to the dissolution of carbonate sands on a growing coral reef, which is captured in the net dissolution quantified above, the marine carbonates on the continental shelf dissolve (Eyre et al., 2018). Like above, the dissolution of marine carbonates is approximated as a source of DIC and alkalinity but does not affect the properties (mass, porosity, etc.) of the underlying sediments.

We assume carbonate dissolution from the sediment bed is proportional to the fraction of the total surface sediment that is composed of either sand or mud carbonates. Other components, whose fraction do not release DIC and alkalinity, include carbonate gravel and non-carbonate mineralogies. Thus, the change in DIC and A_T in the bottom water column layer is given by

$$\begin{array}{}\text{(187)}& {\displaystyle }& {\displaystyle \frac{\partial \mathrm{DIC}}{\partial t}}=-\mathrm{12}{d}_{{\mathrm{CaCO}}_{\mathrm{3}}}\left({\displaystyle \frac{{\mathrm{Mud}}_{{\mathrm{CaCO}}_{\mathrm{3}}}+{\mathrm{Sand}}_{{\mathrm{CaCO}}_{\mathrm{3}}}}{M}}\right)/h_{\mathrm{wc}}\text{(188)}& {\displaystyle }& {\displaystyle \frac{\partial A_{\mathrm{T}}}{\partial t}}=-\mathrm{2}{d}_{{\mathrm{CaCO}}_{\mathrm{3}}}\left({\displaystyle \frac{{\mathrm{Mud}}_{{\mathrm{CaCO}}_{\mathrm{3}}}+{\mathrm{Sand}}_{{\mathrm{CaCO}}_{\mathrm{3}}}}{M}}\right)/h_{\mathrm{wc}},\end{array}$$

where M is the total mass of surface layer inorganic sediments (see Sect. 7), $d_{{\mathrm{CaCO}}_{\mathrm{3}}}$ is the dissolution rate of CaCO₃ and is the reverse reaction to calcification, and h_wc is the thickness of the water column layer. The dissolution rate, $d_{{\mathrm{CaCO}}_{\mathrm{3}}}$ (mmol m⁻² d⁻¹), is assumed to be a function of Ω_a (Eyre et al., 2018):

$$\begin{array}{}\text{(189)}& d_{{\mathrm{CaCO}}_{\mathrm{3}}}=-\mathrm{11.51}{\mathrm{\Omega }}_{\mathrm{a}}+\mathrm{33.683}.\end{array}$$

7 Sediment processes

7.1 Brief summary of processes in the sediments

The EMS model contains a multi-layered sediment compartment with time- and space-varying vertical layers and the same horizontal grid as the water column and epibenthic models. All state variables that exist in the water column layers have an equivalent in the sediment layers. The dissolved tracers are given as a concentration in the porewater, while the particulate tracers are given as a concentration per unit volume (see Sect. 10.3.2).

The sediment model contains inorganic particles of different sizes (dust, mud, sand and gravel) and different mineralogies (carbonate and non-carbonate) (Table 33). The sediment model includes the processes of particulate advection and mixing in the water column, resuspension sinking and settling, as well as sediment overturning and bioturbation (Margvelashvili, 2009). These processes, along with initial conditions, determine the mass of each inorganic particulate type in the sediments.

Table 33 Characteristics of the particulate classes. Y – yes, N – no, I – initial condition, C – catchment, OM – remineralisation from organic matter, B – brown, W – white (Condie et al., 2009; Margvelashvili, 2009).

The critical shear stress for resuspension and the sinking rates are generally larger for large particles, while mineralogy only affects the optical properties. The size-class dust only has a non-carbonate form. FineSed also only has a non-carbonate form and has the same physical and optical properties as non-carbonate mud, except that it is initialised with a zero value. Dust and FineSed are particles that enter the domain from rivers during the simulation.

The organic matter classes are discussed in the Sect. 8.1. The inorganic and organic particulate classes are summarised in Table 33 and undergo resuspension, sinking, settling, sediment overturning and bioturbation in a manner similar to the inorganic particulates.

7.2 Sediment chemistry

7.2.1 Sediment nitrification–denitrification

Nitrification in the sediment is similar to the water column but with a sigmoid rather than hyperbolic relationship at low oxygen for numerical reasons (Eq. 204). Denitrification occurs only in the sediment.

7.2.2 Sediment phosphorus absorption–desorption

Sediment phosphorus absorption–desorption is similar to water column (Eq. 206). There is an additional pool of immobilised particulate inorganic phosphorus (PIPI) which accumulates in the model over time as PIP becomes immobilised and represents permanent sequestration (Eq. 207).

8 Common water, epibenthic and sediment processes

8.1 Detritus remineralisation

The non-living components of the C, N and P cycles include the particulate labile and refractory pools, and a dissolved pool (Fig. 4). The labile detritus has a pool at the Redfield ratio, D_Red, and at the Atkinson ratio, D_Atk, resulting from mortality of organic matter at these ratios (see processes above). The labile detritus from both pools then breaks down into refractory detritus and dissolved organic matter. In order to account for the mixed source of labile detritus, the refractory detritus and dissolved organic matter pools are quantified by individual elements (C, N, P). A component of the breakdown of each of these pools is returned to dissolved inorganic components. Finally, as the refractory and dissolved components are separated into C, N and P components, this introduces the possibility to have P components break down quicker than C and N. This is specified as the breakdown rate of P relative to N, ${\mathrm{\Phi }}_{{\mathrm{RD}}_{\mathrm{P}}}$ and ${\mathrm{\Phi }}_{{\mathrm{DOM}}_{\mathrm{P}}}$, respectively, for refractory and dissolved detritus. The variables, equations and parameter values for detritus remineralisation can be found in Tables 37, 38 and 39, respectively.

Table 34 State and derived variables for the sediment inorganic chemistry model.

Table 35 Constants and parameter values used in the sediment inorganic chemistry.

Table 36 Equations for the sediment inorganic chemistry.

Table 37 State and derived variables for the detritus remineralisation model in both the sediment and water column.

Table 38 Equations for detritus remineralisation in the water column and sediment.

Table 39 Constants and parameter values used in the water column detritus remineralisation model. Red is the Redfield ratio ($\mathrm{C}:\mathrm{N}:\mathrm{P}=\mathrm{106}:\mathrm{16}:\mathrm{1}$); Atk is the Atkinson ratio ($\mathrm{C}:\mathrm{N}:\mathrm{P}=\mathrm{550}:\mathrm{30}:\mathrm{1}$); see Lønborg et al. (2017). n/a – not applicable

8.1.1 Anaerobic and anoxic respiration

The processes of remineralisation, phytoplankton mortality and zooplankton grazing return carbon dioxide to the water column. In oxic conditions, these processes consume oxygen in a ratio of DIC : $\frac{\mathrm{32}}{\mathrm{12}}\left[{\mathrm{O}}_{\mathrm{2}}\right]$. At low oxygen concentrations, the oxygen consumed is

$$\begin{array}{}\text{(221)}& {\displaystyle \frac{\partial \left[{\mathrm{O}}_{\mathrm{2}}\right]}{\partial t}}=-{\displaystyle \frac{\partial \mathrm{DIC}}{\partial t}}{\displaystyle \frac{\mathrm{32}}{\mathrm{12}}}{\displaystyle \frac{[{\mathrm{O}}_{\mathrm{2}}{]}^{\mathrm{2}}}{K_{\mathrm{OA}}^{\mathrm{2}}+[{\mathrm{O}}_{\mathrm{2}}{]}^{\mathrm{2}}}},\end{array}$$

where K_OA is the half-saturation constant for anoxic respiration (Boudreau, 1996). A sigmoid saturation term is used because it is more numerically stable as the oxygen concentration approaches 0. The anoxic component of remineralisation results in an increased chemical oxygen demand (COD):

$$\begin{array}{}\text{(222)}& {\displaystyle \frac{\partial \mathrm{COD}}{\partial t}}={\displaystyle \frac{\partial \mathrm{DIC}}{\partial t}}{\displaystyle \frac{\mathrm{32}}{\mathrm{12}}}\left(\mathrm{1}-{\displaystyle \frac{[{\mathrm{O}}_{\mathrm{2}}{]}^{\mathrm{2}}}{K_{\mathrm{OA}}^{\mathrm{2}}+[{\mathrm{O}}_{\mathrm{2}}{]}^{\mathrm{2}}}}\right),\end{array}$$

where COD is a dissolved tracer, with the same units as oxygen.

When oxygen and COD coexist, they react to reduce both, following

$$\begin{array}{}\text{(223)}& {\displaystyle }& {\displaystyle \frac{\partial \left[{\mathrm{O}}_{\mathrm{2}}\right]}{\partial t}}=-{\mathit{\tau }}_{\mathrm{COD}}\min [\mathrm{COD},\mathrm{8000}]{\displaystyle \frac{\left[{\mathrm{O}}_{\mathrm{2}}\right]}{\mathrm{8000}}}\text{(224)}& {\displaystyle }& {\displaystyle \frac{\partial \mathrm{COD}}{\partial t}}=-{\mathit{\tau }}_{\mathrm{COD}}\min [\mathrm{COD},\mathrm{8000}]{\displaystyle \frac{\left[{\mathrm{O}}_{\mathrm{2}}\right]}{\mathrm{8000}}},\end{array}$$

where 8000 mg O m⁻³ is approximately the saturation concentration of oxygen in seawater, and τ_COD is the timescale of this reduction. The term min[COD,8000] is required because COD represents the end stage of anoxic reduction and can become large for long simulations. Even with this limitation, if τ_COD=1 h⁻¹, the processes in Eqs. (223) and (224) proceed faster than most of the other porewater processes.

9 Common ecological parameterisations

Most of the ecological processes contain a temperature dependence and, for those uptaking dissolved inorganic nitrogen, preferential ammonium uptake. To simplify the description of the above processes, these common parameterisations are described separately in this section. An additional processes common to all variables, and across multiples zones, is the diffusive sediment–water exchange.

9.1 Preferential uptake of ammonium

The model contains two forms of dissolved inorganic nitrogen (DIN) used by photosynthetic organisms, dissolved ammonium (NH₄) and dissolved nitrate (NO₃):

$$\begin{array}{}\text{(225)}& N=\left[{\mathrm{NH}}_{\mathrm{4}}\right]+\left[{\mathrm{NO}}_{\mathrm{3}}\right],\end{array}$$

where N is the concentration of DIN, [NH₄] is the concentration of dissolved ammonium and [NO₃] is the concentration of nitrate. In the model, the ammonium component of the DIN pool is assumed to be taken up first, followed by the nitrate, with the caveat that the uptake of ammonium cannot exceed the diffusion limit for ammonium. The underlying principle of this assumption is that photosynthetic organisms can entirely prefer ammonium but that the uptake of ammonium is still limited by diffusion to the organism's surface.

As the nitrogen uptake formulation varies for the different autotrophs, the formulation of the preference of ammonium also varies. The diffusion coefficient of ammonium and nitrate are only 3 % different, so for simplicity we have used the nitrate diffusion coefficient for both.

Thus, for microalgae (Eq. 40) and Trichodesmium (Eq. 55), which both contain internal reserves of nitrogen, the partitioning of nitrogen uptake is given by

$$\begin{array}{}\text{(226)}& {\displaystyle \frac{\partial N}{\partial t}}=-\mathit{\psi }{D}_{\mathrm{N}}N(\mathrm{1}-R_{\mathrm{N}}^{*})\left(B/m_{B,\mathrm{N}}\right)\text{(227)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\partial \left[{\mathrm{NH}}_{\mathrm{4}}\right]}{\partial t}}=& -\min \left[\mathit{\psi }{D}_{\mathrm{N}}N(\mathrm{1}-R_{\mathrm{N}}^{*}),\mathit{\psi }{D}_{\mathrm{N}}\left[{\mathrm{NH}}_{\mathrm{4}}\right]\right]\\ & \left(B/m_{B,\mathrm{N}}\right)\end{array}\text{(228)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\partial \left[{\mathrm{NO}}_{\mathrm{3}}\right]}{\partial t}}=& -\left(\mathit{\psi }{D}_{\mathrm{N}}N(\mathrm{1}-R_{\mathrm{N}}^{*})-\min \left[\mathit{\psi }{D}_{\mathrm{N}}N\right.\right.\\ & \left.\left.(\mathrm{1}-R_{\mathrm{N}}^{*}),\mathit{\psi }{D}_{\mathrm{N}}\left[{\mathrm{NH}}_{\mathrm{4}}\right]\right]\right)\left(B/m_{B,\mathrm{N}}\right).\end{array}\end{array}$$

For macroalgae (Eq. 108) and seagrass leaves (Eq. 122), which also have diffusion limits to uptake but are not represented with internal reserves of nitrogen, the terms are

$$\begin{array}{}\text{(229)}& {\displaystyle \frac{\partial N}{\partial t}}=-{\mathit{\mu }}_{\mathrm{MA}}\mathrm{MA}\text{(230)}& {\displaystyle \frac{\partial \left[{\mathrm{NH}}_{\mathrm{4}}\right]}{\partial t}}=-\min \left[S{A}_{\mathrm{eff}}\left[{\mathrm{NH}}_{\mathrm{4}}\right],{\mathit{\mu }}_{\mathrm{MA}}\mathrm{MA}\right]\text{(231)}& {\displaystyle }\begin{array}{rl}& {\displaystyle \frac{\partial \left[{\mathrm{NO}}_{\mathrm{3}}\right]}{\partial t}}=-\left({\mathit{\mu }}_{\mathrm{MA}}\mathrm{MA}\right.\\ & \left.-\min \left[S{A}_{\mathrm{eff}}\left[{\mathrm{NH}}_{\mathrm{4}}\right],{\mathit{\mu }}_{\mathrm{MA}}\mathrm{MA}\right]\right).\end{array}\end{array}$$

Zooxanthellae is a combination of the two cases above, because in the model they contain reserves like microalgae, but the uptake rate is across a 2-D surface like macroalgae.

In the case of nutrient uptake by seagrass roots (Eq. 124), which has a saturating nitrogen uptake functional form, the terms are

$$\begin{array}{}\text{(232)}& {\displaystyle \frac{\partial N_s}{\partial t}}=-{\mathit{\mu }}_{\mathrm{SG}}\mathrm{SG}\text{(233)}& {\displaystyle \frac{\partial [{\mathrm{NH}}_{\mathrm{4}}{]}_s}{\partial t}}=-\min \left[{\mathit{\mu }}_{\mathrm{SG}}\mathrm{SG},{\displaystyle \frac{{\mathit{\mu }}_{\mathrm{SG}}^{max}[{\mathrm{NH}}_{\mathrm{4}}{]}_s\mathrm{SG}}{K_N+[{\mathrm{NH}}_{\mathrm{4}}{]}_s}}\right]\text{(234)}& {\displaystyle }\begin{array}{rl}& {\displaystyle \frac{\partial [{\mathrm{NO}}_{\mathrm{3}}{]}_s}{\partial t}}=-\left({\mathit{\mu }}_{\mathrm{SG}}\mathrm{SG}-\min \right.\\ & \left.\left[{\mathit{\mu }}_{\mathrm{SG}}\mathrm{SG},{\displaystyle \frac{{\mathit{\mu }}_{\mathrm{SG}}^{max}[{\mathrm{NH}}_{\mathrm{4}}{]}_s\mathrm{SG}}{K_N+[{\mathrm{NH}}_{\mathrm{4}}{]}_s}}\right]\right),\end{array}\end{array}$$

where K_N is a function of the ratio of above-ground to below-ground biomass described in Baird et al. (2016a).

One feature worth noting is that the above formulation for preferential ammonium uptake requires no additional parameters, which is different to other classically applied formulations (Fasham et al., 1990) that require a new parameter, potentially for each autotroph. This simple use of the geometric constraint has an important role in reducing model complexity.

9.2 Oxygen release during nitrate uptake

For all autotrophs, the uptake of a nitrate ion results in the retention of the one nitrogen atom in their reserves or structural material, and the release of the three oxygen atoms into the water column or porewaters.

$$\begin{array}{}\text{(235)}& {\displaystyle \frac{\partial \left[\mathrm{O}\right]}{\partial t}}=-{\displaystyle \frac{\mathrm{48}}{\mathrm{14}}}{\displaystyle \frac{\partial \left[{\mathrm{NO}}_{\mathrm{3}}\right]}{\partial t}}.\end{array}$$

The oxygen that is part of the structural material of microalgae is assumed to have been taken up through photosynthesis.

For simplicity, in the equations for autotroph-driven changes in dissolved oxygen above, we have assumed that DIN uptake is ammonium. Thus, after partitioning on nitrogen uptake, the term in Eq. (235) needs to be added to change in oxygen in microalgae (Eq. 40), Trichodesmium (Eq. 55) and other autotrophs.

9.3 Temperature dependence of ecological rates

Physiological rate parameters (maximum growth rates, mortality rates, remineralisation rates) have a temperature dependence that is determined from

$$\begin{array}{}\text{(236)}& r_T=r_{T_{\mathrm{ref}}}{Q}_{\mathrm{10}}^{\left(T-T_{\mathrm{ref}}\right)/\mathrm{10}},\end{array}$$

where r_T is the physiological rate parameter (e.g. μ, ζ) at temperature T, T_ref is the reference temperature (nominally 20 ^∘C for GBR), $r_{T_{\mathrm{ref}}}$ the physiological rate parameter at temperature T_ref, Q₁₀ is the Q10 temperature coefficient and represents the rate of change of a biological rate as a result of increasing temperature by 10 ^∘C.

Note that while physiological rates may be temperature dependent, the ecological processes they are included in may not. For example, for extremely light-limited growth, all autotrophs capture light at a rate independent of temperature. With the reserves of nutrients replete, the steady-state realised growth rate, μ, becomes the rate of photon capture, k. This can be shown algebraically: $\mathit{\mu }={\mathit{\mu }}^{max}{R}_{\mathrm{C}}^{*}=k(\mathrm{1}-R_{\mathrm{C}}^{*})$, where $R_{\mathrm{C}}^{*}$ is the reserves of carbon. Rearranging, $R_{\mathrm{C}}^{*}=k/({\mathit{\mu }}^{max}+k)$. At $k\ll {\mathit{\mu }}^{max}$, $R_{\mathrm{C}}^{*}=k/{\mathit{\mu }}^{max}$; thus, $\mathit{\mu }={\mathit{\mu }}^{max}k/{\mathit{\mu }}^{max}=k$. This corresponds with observations of no temperature dependence of photosynthesis at low light levels (Kirk, 1994).

Similar arguments show that extremely nutrient limited autotrophs will have the same temperature dependence to that of the diffusion coefficient. Thus, the autotroph growth model has a temperature dependence that adjusts appropriately to the physiological condition of the autotroph and is a combination of constant, exponential and polynomial expressions.

Physiological rates in the model that are not temperature dependent are mass transfer rate constant for particulate grazing by corals, S_Part; net coral calcification g; maximum chlorophyll synthesis, $k_{\mathrm{Chl}}^{max}$; and rate of translocation between leaves and roots in seagrass, τ_tran.

9.4 Diffusive exchange of dissolved tracers across sediment–water interface

Due to the thin surface sediment layer, and the potentially large epibenthic drawdown of porewater dissolved tracers, the exchange of dissolved tracers between the bottom water column layer and the top sediment layer is integrated in the same numerical operation as the ecological tracers (other transport processes occurring between ecological time steps). The flux, J, is given by

$$\begin{array}{}\text{(237)}& J=k(C_{\mathrm{s}}-C),\end{array}$$

where C and C_s are the concentration in water column and sediment, respectively, and k is the transfer coefficient. For the simulations, we used $k=D/h=\mathrm{4.6}\times {\mathrm{10}}^{-\mathrm{7}}$ m s⁻¹, where $D=\mathrm{3}\times {\mathrm{10}}^{-\mathrm{9}}$ m² s⁻¹ is the diffusion coefficient and h=0.0065 mm is the thickness of the diffusive layer.

While, in reality, k would vary with water column and sediment hydrodynamics as influenced by community type, etc., these complexities have not been considered. In addition to the diffusive flux between the sediment and water column, particulate deposition entrains water column water into the sediments, and particulate resuspension releases porewaters into the water column. Sediment model details can be found at https://research.csiro.au/cem/software/ems/ems-documentation/ (last access: 22 September 2020).

10 Numerical integration

10.1 Splitting of physical and ecological integrations

The numerical solution of the time-dependent advection–diffusion–reaction equations for each of the ecological tracers is implemented through sequential solving of the partial differential equations (PDEs) for advection and diffusion, and the ordinary differential equations (ODEs) for reactions. This technique, called operator splitting, is common in geophysical science (Hundsdorfer and Verwer, 2003; Butenschön et al., 2012).

Under the sequential operator splitting technique used, first the advection–diffusion processes are solved for the period of the time step (15 min–1 h; Table 40). The values of the tracers at the end of this PDE integration, and the initial time, are then used as initial conditions for the ODE integration. After the ODE integration has run for the same time period, the values of the tracers are updated, and time is considered to have moved forward just one time step. The integration continues to operate sequentially for the whole model simulation. The errors due to operator splitting can be significant (Butenschön et al., 2012), although tests in relatively coarse (4 km) models show that reducing the time step from 60 to 30 min does not substantially change the model solution. For higher-resolution models, shorter timescales are required to resolve finer-scale motion and its interaction with ecological processes.

(Dormand and Prince, 1980)

Table 40 Integration details. Optical wavelengths are 290, 310, 330, 350, 370, 390, 410, 430, 440, 450, 470, 490, 510, 530, 550, 570, 590, 610, 630, 650, 670, 690, 710 and 800 nm.

^* For the fourth- to fifth-order integrators, the ecological derivatives are evaluated at least every approximately $\mathrm{3600}/\mathrm{4}=\mathrm{900}$ s.

The PDE solvers are described in the physical model description available at https://research.csiro.au/cem/software/ems/ems-documentation/ (last access: 22 September 2020).

The code allows fourth- to fifth-order and seventh- to eighth-order adaptive ODE solvers following Dormand and Prince (1980), as well as the Euler method and adaptive first- and second-order solvers. The preferred scheme is the adaptive fourth–fifth order (similar to ode45 in MATLAB) and implemented in numerous biogeochemical models (Yool, 1997). This requires seven function evaluations for the first step and six for each step after. A relative tolerance of 10⁻⁵ is required for each state variable for the integration step to proceed.

The solution of the ecological equations are independent for each vertical column and depend only on the layers above through which the light has propagated. For an n_wc-layer water column and n_sed-layer sediment, the integrator sequentially solves the top n_wc−1 water column layers; the nth water column layer, epibenthic and top sediment layer together; and then the n_sed−1 to bottom sediment layers.

10.2 Optical integration

The inherent and apparent optical properties are calculated between the hydrodynamic and ecological integrations. The light field used for each ecological time step is that calculated at the start time of the ecological integration. The spectral resolution of 25 wavebands has been chosen to resolve the absorption peaks associated with Chl a and to span the optical wavelengths. The wavelengths integrated have been chosen such that the lower end of one waveband and the top end of another fall on 400 and 700 nm, respectively, allowing precise calculation of photosynthetically available radiation (PAR).

10.3 Additional integration details

10.3.1 Approximation of stoichiometric coefficients

In this model description, we have chosen to explicitly include atomic mass as integer values, so that the unit conversions are more readable in the equations than if they had all been rendered as mathematical symbols. Nonetheless, these values are more precisely given in the numerical code (Table 41).

Table 41 Atomic mass of the C, N, P and O₂, both in the model description where two significant figures are used for brevity and in the numerical code, where precision is more important.

It is worth remembering that the atomic masses are approximations assuming the ratio of isotopes found in the periodic table (Atkins, 1994) based on the natural isotopic abundance of the Earth. So, for example, ¹⁴N and ¹⁵N have atomic masses of 14.00307 and 15.00011, respectively, with ¹⁴N making up 99.64 % of the abundance on Earth. Thus, the value 14.01 comes from $\mathrm{14.00307}\times \mathrm{99.64}+\mathrm{15.00011}\times \mathrm{0.36}=\mathrm{14.0067}$. If the model had state variables for ¹⁴N and ¹⁵N, then the equations would change to contain coefficients of 14 for the ¹⁴N isotope equations and 15 for the ¹⁵N isotope equations.

10.3.2 Mass conservation in water column and sediment porewaters

A check of mass conservation of total C, TC, total N, TN, total P, TP, and oxygen, [O₂], is undertaken in each grid cell at each time step. To establish mass conservation, the sum of the change in mass (of N, P, C and O) with time and the sinks/sources (such as sea–air fluxes, denitrification) must equate to zero.

The total mass and conservation equations are same for the water column and porewaters, with the caveats that (1) air–sea fluxes only affect surface layers of the water, (2) denitrification only occurs in the sediment, and (3) the porosity, ϕ, of the water column is 1. In the sediment, the concentration of particulates is given per unit volume of space, while the concentration of dissolved tracers is given per unit volume of porewater. Thus, the concentration of dissolved tracer, X, per unit space is given by ϕX.

The total carbon in a unit volume of space and its conservation are given by

$$\begin{array}{}\text{(238)}& {\displaystyle }\begin{array}{rl}\mathrm{TC}=& \mathit{\varphi }\left(\mathrm{DIC}+O_{\mathrm{C}}\right)+\left({\displaystyle \frac{\mathrm{550}}{\mathrm{30}}}{\displaystyle \frac{\mathrm{12}}{\mathrm{14}}}{D}_{\mathrm{Atk}}+D_{\mathrm{C}}\right.\\ & \left.+{\displaystyle \frac{\mathrm{106}}{\mathrm{16}}}{\displaystyle \frac{\mathrm{12}}{\mathrm{14}}}\left(D_{\mathrm{red}}+\sum B(\mathrm{1}+R_{\mathrm{C}}^{*})+\sum Z\right)\right)\end{array}\text{(239)}& {\displaystyle \frac{\partial \mathrm{TC}}{\partial t}}+\underset{\mathrm{sea}-\mathrm{air}\mathrm{flux}}{\underbrace{k_{{\mathrm{CO}}_{\mathrm{2}}}\left(\left[{\mathrm{CO}}_{\mathrm{2}}\right]-[{\mathrm{CO}}_{\mathrm{2}}{]}_{\mathrm{atm}}\right)/h}}=\mathrm{0}.\end{array}$$

The total nitrogen in a unit volume of space and its conservation are given by

$$\begin{array}{}\text{(240)}& {\displaystyle }\begin{array}{rl}\mathrm{TN}=& \mathit{\varphi }\left(\left[{\mathrm{NO}}_{\mathrm{3}}\right]+\left[{\mathrm{NH}}_{\mathrm{4}}\right]+O_{\mathrm{N}}\right)+\left(D_{\mathrm{Atk}}\right.\\ & \left.+D_{\mathrm{red}}+D_{\mathrm{N}}+\sum B(\mathrm{1}+R_{\mathrm{N}}^{*})+\sum Z\right)\end{array}\text{(241)}& {\displaystyle }\begin{array}{rl}& {\displaystyle \frac{\partial \mathrm{TN}}{\partial t}}+\left(\mathrm{denitrification}-\mathrm{nitrogen}\mathrm{fixation}\right)/\mathit{\varphi }-\\ & \mathrm{dust}\mathrm{input}/h=\mathrm{0}.\end{array}\end{array}$$

The total phosphorus in a unit volume of space and its conservation are given by

$$\begin{array}{}\text{(242)}& {\displaystyle }\begin{array}{rl}\mathrm{TP}=& \mathit{\varphi }\left(\mathrm{DIP}+O_{\mathrm{P}}\right)+\mathrm{PIP}+\mathrm{PIPI}+{\displaystyle \frac{\mathrm{1}}{\mathrm{30}}}{\displaystyle \frac{\mathrm{31}}{\mathrm{14}}}{D}_{\mathrm{Atk}}\\ & +D_{\mathrm{P}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{16}}}{\displaystyle \frac{\mathrm{31}}{\mathrm{14}}}\left(D_{\mathrm{red}}+\sum B(\mathrm{1}+R_{\mathrm{P}}^{*})+\sum Z\right)\end{array}\text{(243)}& {\displaystyle \frac{\partial \mathrm{TP}}{\partial t}}-\mathrm{dust}\mathrm{input}/h=\mathrm{0}.\end{array}$$

The concept of oxygen conservation in the model is more subtle than that of C, N and P due to the mass of oxygen in the water molecules not being considered. When photosynthesis occurs, C is transferred from the dissolved phase to reserves within the cell. With both dissolved and particulate pools considered, mass conservation of C is straightforward. In contrast to C, during photosynthesis, oxygen is drawn from the water molecules (i.e. H₂O), whose mass is not being considered, and released into the water column. Conversely, when organic matter is broken down, oxygen is consumed from the water column and released as H₂O.

In order to obtain a mass conservation for oxygen, the concept of biological oxygen demand (BOD) is used. Often BOD represents the biological demand for oxygen in say a 5 d incubation, BOD₅. Here, for the purposes of mass conservation checks, we use BOD_∞, the oxygen demand over an infinite time for breakdown. This represents the total oxygen removed from the water molecules for organic matter creation.

Anaerobic respiration reduces BOD_∞ without reducing O₂ but instead creating reduced-oxygen species. This is accounted for in the oxygen balance by the prognostic tracer COD. In other biogeochemical modelling studies, this is represented by a negative oxygen concentration. Thus, at any time point, the biogeochemical model will conserve the oxygen concentration minus BOD_∞ minus COD, plus or minus any sources and sinks such as sea–air fluxes:

$$\begin{array}{}\text{(244)}& {\displaystyle }\begin{array}{rl}& \left[{\mathrm{O}}_{\mathrm{2}}\right]+{\displaystyle \frac{\mathrm{48}}{\mathrm{14}}}\left[{\mathrm{NO}}_{\mathrm{3}}\right]-{\mathrm{BOD}}_{\mathrm{\infty }}-\mathrm{COD}=\\ & \mathit{\varphi }\left(\left[{\mathrm{O}}_{\mathrm{2}}\right]+{\displaystyle \frac{\mathrm{48}}{\mathrm{14}}}\left[{\mathrm{NO}}_{\mathrm{3}}\right]-\mathrm{COD}+{\displaystyle \frac{\mathrm{32}}{\mathrm{12}}}{O}_{\mathrm{C}}\right)\\ & -\left({\displaystyle \frac{\mathrm{550}}{\mathrm{30}}}{\displaystyle \frac{\mathrm{32}}{\mathrm{14}}}{D}_{\mathrm{Atk}}+{\displaystyle \frac{\mathrm{32}}{\mathrm{12}}}{D}_{\mathrm{C}}+{\displaystyle \frac{\mathrm{106}}{\mathrm{16}}}{\displaystyle \frac{\mathrm{32}}{\mathrm{14}}}\right.\\ & \left.\left(D_{\mathrm{red}}+\sum B_{\mathrm{N}}(\mathrm{1}+R_{\mathrm{C}}^{*})\right)\right)\end{array}\text{(245)}& {\displaystyle }\begin{array}{rl}& {\displaystyle \frac{\partial \left(\right[{\mathrm{O}}_{\mathrm{2}}]+\frac{\mathrm{48}}{\mathrm{14}}[{\mathrm{NO}}_{\mathrm{3}}]-{\mathrm{BOD}}_{\mathrm{\infty }}-\mathrm{COD})}{\partial t}}+\mathcal{R}\\ & -\stackrel{\text{sea–air flux}}{\overbrace{{\displaystyle \frac{k_{{\mathrm{O}}_{\mathrm{2}}}\left([{\mathrm{O}}_{\mathrm{2}}{]}_{\mathrm{sat}}-[{\mathrm{O}}_{\mathrm{2}}]\right)}{h}}}}\\ & \underset{\mathrm{nitrification}}{\underbrace{-\mathrm{2}{\displaystyle \frac{\mathrm{106}}{\mathrm{16}}}{\displaystyle \frac{\mathrm{32}}{\mathrm{14}}}{\mathit{\tau }}_{\mathrm{nit},\mathrm{wc}}\left[{\mathrm{NH}}_{\mathrm{4}}\right]{\displaystyle \frac{\left[{\mathrm{O}}_{\mathrm{2}}\right]}{K_{\mathrm{nit},\mathrm{O}}+\left[{\mathrm{O}}_{\mathrm{2}}\right]}}}}=\mathrm{0}\end{array},\end{array}$$

where ℛ is respiration of organic matter.

In addition to dissolved oxygen, BOD and COD, nitrate (NO₃) appears in the oxygen mass balance. This is necessary because the N associated with nitrate uptake is not taken into the autotrophs but rather released into the water column or porewater. Other entities that contain oxygen in the ocean include the water molecule (H₂O) and the phosphorus ion (PO₄). In the case of water, this oxygen reservoir is considered very large, with the small flux associated with its change balanced by BOD. In the case of PO₄, this is a small reservoir. As oxygen remains bound to P through the entire processes of uptake into reserves and incorporated into structural material and then release, it is not necessary to include it in the oxygen balance for the purposes of ensuring consistency. Nonetheless, strictly the water column and porewater oxygen reservoirs could include a term plus $\frac{\mathrm{64}}{\mathrm{31}}\left[{\mathrm{PO}}_{\mathrm{4}}\right]$, and the BOD would have similar quantities for reserves and structural material.

10.3.3 Mass conservation in the epibenthic

Mass conservation in the epibenthos requires consideration of fluxes between the water column, porewaters and the epibenthic organisms (macroalgae, seagrass and coral hosts and symbionts). The total carbon in the epibenthos, and its conservation, is given by

$$\begin{array}{}\text{(246)}& {\displaystyle }\begin{array}{rl}\mathrm{TC}=& {\displaystyle \frac{\mathrm{550}}{\mathrm{30}}}{\displaystyle \frac{\mathrm{12}}{\mathrm{14}}}\left(\mathrm{MA}+{\mathrm{SG}}_{\mathrm{A}}+{\mathrm{SG}}_{\mathrm{B}}\right)\\ & +{\displaystyle \frac{\mathrm{106}}{\mathrm{16}}}{\displaystyle \frac{\mathrm{12}}{\mathrm{14}}}\left(\mathrm{CS}\left(\mathrm{1}+R_{\mathrm{C}}^{*}\right)+\mathrm{CH}\right)\end{array}\text{(247)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\partial \mathrm{TC}}{\partial t}}& {|}_{\mathrm{epi}}+h_{\mathrm{wc}}{\displaystyle \frac{\partial \mathrm{TC}}{\partial t}}{|}_{\mathrm{wc}}+h_{\mathrm{sed}}{\displaystyle \frac{\partial \mathrm{TC}}{\partial t}}{|}_{\mathrm{sed}}\\ & +\underset{\mathrm{coral}\mathrm{calcification}-\mathrm{dissolution}}{\underbrace{\mathrm{12}\left(g{A}_{\mathrm{eff}}-d_{{\mathrm{CaCO}}_{\mathrm{3}}}\right)}}=\mathrm{0},\end{array}\end{array}$$

where h_wc and h_sed are the thickness of the bottom water column and top sediment layers, $R_{\mathrm{C}}^{*}$ is the normalised internal reserves of carbon in zooxanthellae, 12 g is the rate coral calcification per unit area of coral, A_eff is the area of the bottom covered by coral per m², and the diffusion terms between porewaters and the water column cancel, so they do not appear in the equations. Note the units of mass of CS needs to be in g N, and some configurations may have multiple seagrass and macroalgae species.

Similarly, for nitrogen, phosphorus and oxygen in the epibenthos,

$$\begin{array}{}\text{(248)}& {\displaystyle }\mathrm{TN}=\mathrm{MA}+{\mathrm{SG}}_{\mathrm{A}}+{\mathrm{SG}}_{\mathrm{B}}+\mathrm{CS}\left(\mathrm{1}+R_{\mathrm{N}}^{*}\right)+\mathrm{CH}\text{(249)}& {\displaystyle \frac{\partial \mathrm{TN}}{\partial t}}{|}_{\mathrm{epi}}+h_{\mathrm{wc}}{\displaystyle \frac{\partial \mathrm{TN}}{\partial t}}{|}_{\mathrm{wc}}+h_{\mathrm{sed}}{\displaystyle \frac{\partial \mathrm{TN}}{\partial t}}{|}_{\mathrm{sed}}=\mathrm{0}\text{(250)}& {\displaystyle }\begin{array}{rl}\mathrm{TP}=& {\displaystyle \frac{\mathrm{1}}{\mathrm{30}}}{\displaystyle \frac{\mathrm{31}}{\mathrm{14}}}\left(\mathrm{MA}+{\mathrm{SG}}_{\mathrm{A}}+{\mathrm{SG}}_{\mathrm{B}}\right)\\ & +{\displaystyle \frac{\mathrm{1}}{\mathrm{16}}}{\displaystyle \frac{\mathrm{31}}{\mathrm{14}}}\left(\mathrm{CS}\left(\mathrm{1}+R_{\mathrm{P}}^{*}\right)+\mathrm{CH}\right)\end{array}\text{(251)}& {\displaystyle \frac{\partial \mathrm{TP}}{\partial t}}{|}_{\mathrm{epi}}+h_{\mathrm{wc}}{\displaystyle \frac{\partial \mathrm{TP}}{\partial t}}{|}_{\mathrm{wc}}+h_{\mathrm{sed}}{\displaystyle \frac{\partial \mathrm{TP}}{\partial t}}{|}_{\mathrm{sed}}=\mathrm{0}\text{(252)}& {\displaystyle }\begin{array}{rl}{\mathrm{BOD}}_{\mathrm{\infty }}=& {\displaystyle \frac{\mathrm{550}}{\mathrm{30}}}{\displaystyle \frac{\mathrm{32}}{\mathrm{14}}}\left(\mathrm{MA}+{\mathrm{SG}}_{\mathrm{A}}+{\mathrm{SG}}_{\mathrm{B}}\right)\\ & +{\displaystyle \frac{\mathrm{106}}{\mathrm{16}}}{\displaystyle \frac{\mathrm{32}}{\mathrm{14}}}\left(\mathrm{CS}\left(\mathrm{1}+R_{\mathrm{C}}^{*}\right)+\mathrm{CH}\right)\end{array}\text{(253)}& {\displaystyle }\begin{array}{rl}-{\displaystyle \frac{\partial {\mathrm{BOD}}_{\mathrm{\infty }}}{\partial t}}& {|}_{\mathrm{epi}}+h_{\mathrm{wc}}{\displaystyle \frac{\partial \left(\right[{\mathrm{O}}_{\mathrm{2}}]-{\mathrm{BOD}}_{\mathrm{\infty }})}{\partial t}}{|}_{\mathrm{wc}}\\ & +h_{\mathrm{sed}}{\displaystyle \frac{\partial \left(\right[{\mathrm{O}}_{\mathrm{2}}]-{\mathrm{BOD}}_{\mathrm{\infty }})}{\partial t}}{|}_{\mathrm{sed}}=\mathrm{0},\end{array}\end{array}$$

where there is no dissolved oxygen in the epibenthos.

10.3.4 Unconditional stability

In addition to the above standard numerical techniques, a number of innovations are used to ensure model solutions are reached. Should an integration step fail in a grid cell, no increment of the state variables occurs, and the model continues with a warning flag registered (as Ecology Error). Generally, the problem does not reoccur due to the transport of tracers alleviating the stiff point in phase space of the model. Furthermore, when a water column becomes dry (the sea level drops below the seabed depth), ecological processes are turned off.

11 Model evaluation

The EMS BGC model has been deployed in a range of environments around Australia. With each deployment, a skill assessment has been undertaken (for a history of these applications, see Sect. 13). Recently, the EMS BGC model has been thoroughly assessed against remotely sensed and in situ observations on the GBR as part of the eReefs project (Schiller et al., 2014; Steven et al., 2019b). The assessment of version B1p0 of the eReefs marine model configuration of the EMS included a range of model configurations (4 km, 1 km and relocatable fine-resolution versions) (Herzfeld et al., 2016). The optical and carbon chemistry outputs were assessed in Baird et al. (2016b) and Mongin et al. (2016b), respectively.

A more recent assessment of the BGC model (vB2p0) in the GBR compared simulations against a range of in situ observations that included 24 water quality moorings, two nutrient sampling programmes (with a total of 18 stations) and time series of taxon-specific plankton abundance. In addition to providing a range of skill metrics, the assessment included analysis of seasonal plankton dynamics (Skerratt et al., 2019).

In this section, we assess version B3p0 in the 4 km GBR configuration. First, we consider the behaviour of the microalgae physiology as a means to understanding the dynamics of the microalgal growth model. Secondly, the techniques and observations used in Skerratt et al. (2019) have been applied to the model version described in this paper (vB3p0) with highlights discussed here and the full analysis appearing in the Supplement.

11.1 Analysis of microalgae growth and pigment synthesis dynamics

The microalgae growth model (Sect. 5.1.1) was derived from first quantifying the fluxes into a cell of energy (for fixing carbon), N, P and O. Each flux adds to the reserves of that element in the cell. A second process, the consumption of reserves to create microalgae structural material with a constant stoichiometry ($\mathrm{C}:\mathrm{N}:\mathrm{P}=\mathrm{106}:\mathrm{16}:\mathrm{1}$), increases the number of cells but reduces the reserves of each element both per cell and of the population. Thus, the microalgae in the growth model are generating organic matter at the Redfield ratio while being exposed to external nutrient fields at non-Redfield ratios. At the same time as the microalgae grow, the model represents the synthesis of chlorophyll based on the cell's need for more carbon fixation and the benefit of adding pigment on the rate of photon absorption (Fig. 11).

To illustrate these dynamics, we look at a vertical profile of a deep site in the Coral Sea with a 1 and 2.5 µm radius microalgae (Fig. 13). The expectation is that, for the same environmental conditions, the 1 µm cell will be less nutrient limited and more light limited than the 2.5 µm cell – a result of the 1 µm cell having a greater diffusion rate per unit volume than the larger cell. Further, that near-surface cells will be more nutrient limited, deeper cells more light limited. Finally, light-limited cells will generally have more pigment per cell.

Figure 13 Vertical profiles of physiochemical variables (a) and physiological variables from 1 µm (b, d) and 2.5 µm (c, e) radii microalgae. Panel (a) shows photosynthetically available radiation (PAR; mol m⁻² s⁻¹), dissolved inorganic phosphorus (DIP; mg P m⁻³), dissolved inorganic nitrogen (DIN; mg N m⁻³), vertical diffusivity (K_z; m³ s⁻¹) and temperature (^∘C). Panels (b) (1 µm) and (c) (2.5 µm) show biomass of structural material in nitrogen (B_N; mg N m⁻³), chlorophyll a concentration (Chl; mg Chl a m⁻³), normalised reserves of fixed carbon ($R_{\mathrm{C}}^{*}$), nitrogen ($R_{\mathrm{N}}^{*}$) and phosphorus ($R_{\mathrm{P}}^{*}$), cell opaqueness (absorption cross-section divided by projected area, α∕(πr²)) and normalised growth rate of $R_{\mathrm{C}}^{*}{R}_{\mathrm{N}}^{*}{R}_{\mathrm{P}}^{*}$. Panels (d) (1 µm) and (e) (2.5 µm) show stoichiometric ratios of carbon to nitrogen (C:N), phosphorus (C:P) and chlorophyll (C:Chl a). PAR, K_z, temperature and B_N are all scaled for plotting.

In addition to being a measure of the quantity of nutrient reserves, normalised reserves (R^*) of each nutrient is a metric of how limiting that nutrient is, with 1 being unlimited, and 0 being completely limited. At the surface at midday, the 1 µm cells have a biomass of 0.2 mg m⁻³ (Fig. 13b). The cells are strongly phosphorus limited ($R_{\mathrm{P}}^{*}=\mathrm{0.22}$), slightly nitrogen limited ($R_{\mathrm{N}}^{*}$ = 0.86) and almost light unlimited ($R_{\mathrm{C}}^{*}=\mathrm{0.99}$). The realised growth rate, as a fraction of the maximum growth rate, is $R_{\mathrm{C}}^{*}{R}_{\mathrm{N}}^{*}{R}_{\mathrm{P}}^{*}=\mathrm{0.19}$. The larger cells are more nutrient limited ($R_{\mathrm{P}}^{*}=\mathrm{0.14}$, $R_{\mathrm{N}}^{*}=\mathrm{0.54}$), and again light unlimited ($R_{\mathrm{C}}^{*}=\mathrm{0.98}$), realising a normalised growth rate of 0.07 (Fig. 13c).

The elemental ratios of the microalgae can be calculated from the reserves (in wt wt⁻¹: $\mathrm{C}:\mathrm{N}=(\mathrm{12}/\mathrm{14})(\mathrm{106}/\mathrm{16})(\mathrm{1}+R_{\mathrm{C}}^{*})/(\mathrm{1}+R_{\mathrm{N}}^{*})$; $\mathrm{C}:\mathrm{P}=(\mathrm{12}/\mathrm{31})(\mathrm{106}/\mathrm{1})(\mathrm{1}+R_{\mathrm{C}}^{*})/(\mathrm{1}+R_{\mathrm{N}}^{*})$). The C:N and C:P ratios are both higher in the nutrient-replete surface waters, with C:P varying more due to greater P limitation in the surface waters. A deep chlorophyll maximum has formed for the 1 µm microalgae at 40 m and for 2.5 µm microalgae at 60 m. Here, we will explain this distribution based on microalgae growth alone (ignoring grazing and sinking terms). The 1 µm microalgae has a growth maximum at 40 m, as nutrients have become non-limiting and fixed carbon reserves are still relatively high ($R_{\mathrm{C}}^{*}=\mathrm{0.8}$). Growth below 50 m becomes primarily light limited, so normalised growth is equal to normalised fixed carbon reserves. The 2.5 µm microalgae are nutrient limited deeper into the water column, resulting in a deep growth maximum.

The 1 and 2.5 µm cells have a slight biomass maximum at 40 and 50 m, respectively, but the chlorophyll maximum is more pronounced (Fig. 13b, c). At the surface, the 1 µm microalgae have low pigment content and are therefore relatively transparent (opaqueness, or the absorption cross-section divided by the projected area, α∕(πr²), is 0.11), with a C:Chl ratio of 100 [g g⁻¹]. At the deep chlorophyll maxima for 1 µm cells, the pigment content has increased, as shown by the C:Chl dropping to 20 [g g⁻¹], the minimum C:Chl ratio observed in the ocean (Sathyendranath et al., 2009), resulting in an opaqueness of 0.25. The 2.5 µm cells are less light limited than the smaller cells. At the surface, chlorophyll synthesis generates in 2.5 µm cells an opaqueness of 0.52. Given the larger size of the cell, this is achieved at a C:Chl ratio of 150. At the 50 m deep chlorophyll maxima, the C:Chl dropped to 20, but in the larger cells this achieves an opaqueness of 0.91 (i.e. the absorption cross-section is almost equal to the projected area).

In summary, the application of simple physical limits to uptake, a restraint of constant stoichiometric conversion to structural material, and cells synthesising chlorophyll to maximise photon absorption when light limited, generates the typical physiological properties of microalgae seen in vertical profiles in the ocean.

11.2 Model assessment of the GBR configuration

A detailed comparison of a GBR simulation against observations of Chl a concentration, dissolved inorganic carbon, nitrogen, phosphorus and ammonium, dissolved organic nitrogen and phosphorus, alkalinity, pH, aragonite saturation, mass of suspended sediments, turbidity and Secchi depth appears in the Supplement. Here, we highlight the carbon chemistry, nutrient and plankton components that are driven by the BGC model.

11.2.1 Chlorophyll dynamics

The most accurate measurements of water column chlorophyll concentrations in the GBR are obtained using high-performance liquid chromatography (HPLC) and chlorophyll extractions from water column samples. Chlorophyll extractions have been taken at 36 locations along the GBR (S5, Sect. 10; for site locations, see S5, Sect. E1). As an example, a time series at Pelorus Island in the central GBR (Fig. 14) shows large variability in both the observations and the simulations, driven by interannual climatic forcing, with 2011–2013 experiencing much greater river loads than 2014–2016, intra-annual trends driven by greater loads of nutrients during the wet season (January–June) than the remainder of the year, as well as monthly variability related to tidal movements and predator–prey oscillations. Even given this variability, comparison of the instantaneous state of the simulations against extracted chlorophyll concentrations showed the model was able to achieve across all 36 sites a bias ± root mean square error (RMSE) of $-\mathrm{0.11}\pm \mathrm{0.32}$ mg m⁻³ (Fig. 15).

Figure 14 Observed surface chlorophyll concentration from chlorophyll extraction (red dots) at the Pelorus Island Marine Monitoring Program site ($\mathrm{18}{}^{\circ }{\mathrm{33}}^{\prime }$ S, $\mathrm{146}{}^{\circ }{\mathrm{29}}^{\prime }$ E) with a comparison to configurations vB2p0 (pink line) and vB3p0 (blue line). Statistics listed include the Willmott d2 metric (Willmott et al., 1985), mean absolute percent error (MAPE) and RMSE.

Moored fluorometers are generally less accurate than chlorophyll extractions but provide a greater temporal resolution of chlorophyll dynamics. Here, we show observations from a mooring at Palm Passage (Fig. 16), away from the influence of the river discharge and dependent primary on shelf break interactions. The observed time series from 60 m depth show interannual variability in fluorescence related to primarily to an October–March maximum in intrusion events (for details of deployment details and oceanographic processes, see Benthuysen et al., 2016), which are also seen in the vB3p0 simulations. Comparison to other moored fluorometers on the GBR (S5, Sects. 19 and 21) shows a range of predictive skill.

Figure 15 Skill metrics for the comparison of chlorophyll extracts at the long-term monitoring sites against observations for model versions 2p0 and 3p0. For more information, see Fig. 14 and S5, Sect. 10.

11.2.2 Nutrient dynamics

The model represents dissolved nitrate, ammonium and phosphorus nutrients. In the surface waters of the inshore GBR, nutrients are generally at very low concentrations, with modest increases seen each wet season. At High Island in the central GBR (Fig. 17), DIP has a mean concentration of 2 mg m⁻³, with concentrations varying between 0 and 5 mg m⁻³. The concentrations of nitrate and ammonium are both very low, with occasional peaks driven by river plume exposure. The simulated nutrients generally follow the expected time-varying patterns: peaks in the wet season, larger peaks in the wettest years (2011, 2012, 2013, 2018) and extremely low concentrations in the dry seasons. Across all sites (S5, Sect. 13), vB3p0 predicted DIP with a skill of (bias ± RMSE) of $-\mathrm{0.88}\pm \mathrm{2.17}$ mg P m⁻³, nitrate of $-\mathrm{0.70}\pm \mathrm{3.86}$ mg N m⁻³ and ammonium of $-\mathrm{0.77}\pm \mathrm{1.63}$ mg N m⁻³.

Figure 16 Observed chlorophyll fluorescence (red dots) at 60 m depth at the Palm Passage site with a comparison to configurations vB3p0 (blue) and vB2p0 (pink).

Figure 17 Dissolved inorganic phosphorus (a), nitrate (b) and ammonium (c) concentrations at High Island, central GBR: observations (red dots), simulations vB2p0 (pink) and vB3p0 (blue).

11.2.3 Carbon chemistry

The model contains two state variables to represent the state of carbon chemistry, dissolved inorganic carbon and alkalinity, from which, at equilibrium and known temperature and salinity, other variables such as pH may be calculated. The biogeochemical model provides highly skilful predictions of pH and aragonite saturation (Fig. 18). This success is primarily due to the skill of the hydrodynamic model. The circulation produces a good representation of the alkalinity field, and the sea surface temperature sets an accurate difference in the partial pressure of CO₂ between the atmosphere and ocean. With these phenomena accurately calculated, the model needs only to correctly predict the time-averaged wind-speed-dependent air–sea flux. Further assessment of carbon chemistry properties along the entire length of the GBR (S5; Sect. 28) shows a bias ± RMSE of DIC of $-\mathrm{7.7}\pm \mathrm{34.2}$ mmol m⁻³. Further skill assessment from inshore sites is available in Mongin et al. (2016b).

Figure 18 Aragonite saturation state calculated from temperature, DIC and alkalinity at 20 m depth at the Integrated Marine Observing System (IMOS) Yongala mooring, central GBR: observations (red dots), simulations B2p0 (pink) and B3p0 (blue).

The outputs of all hindcasts in the eReefs project can be downloaded from http://dapds00.nci.org.au/thredds/catalogs/fx3/catalog.html (last access: 22 September 2020).

12 Relocatable Coast and Ocean Model (RECOM)

A web-based interface, RECOM, has been developed to automate the process of downscaling the EMS model using an existing hindcast as boundary conditions (https://research.csiro.au/ereefs/models/models-about/recom/, last access: 22 September 2020, including the RECOM User Manual). For the purposes of learning how to apply the EMS software available, RECOM provides the user with the ability to generate a complete test case of a domain situated along the northeast Australian coastline. Once a RECOM simulation has been generated using the web interface, the entire simulation including source code, forcing and initial condition files, model configuration files and the model output can be downloaded. This allows the user to repeat the model simulation on their own computing system and modify code, forcing and output frequency as required. The technical details of RECOM are detailed in Baird et al. (2018) and in the RECOM user manual.

13 Discussion

The EMS BGC model development has been a function of the historical applications of the model across a rage of ecosystems, so it is worth giving a brief history of the model development.

13.1 History of the development of the EMS biogeochemical model

The EMS biogeochemical model was first developed as a nitrogen-based model for determining the assimilative capacity for sewerage discharged into Port Philip Bay, the embayment of the city of Melbourne (Harris et al., 1996). This study saw a focus on sediment processes such as denitrification and demonstrated the ability of bay-wide denitrification to prevent change in the ecological state of the bay exposed to sewerage treatment plant loads (Murray and Parslow, 1997; Murray and Parlsow, 1999). The basic structure of the model and in particular the split of pelagic, epibenthic and sediment zones were in place for this project. This zonation generated the ability to resolve processes in shallow-water systems and in particular to consider benthic flora in detail.

The next major study involved simulating a range of estuarine morphologies (salt wedge, tidal, lagoon, residence times) and forcings (river flow seasonality, nutrient inputs, etc.) that were representative of Australia's 1000+ estuaries (Baird et al., 2003). At this point, carbon and phosphorus were included in the model, and the process of including physical limits to ecological processes begun (e.g. diffusion limitation of nutrient uptake and encounter-rate limitation of grazing).

Following studies in the phosphorus-limited Gippsland Lakes and macro-tidal Ord River system led to the refinement of the phosphorus absorption–desorption processes. Further studies of the biogeochemical–sediment interactions in the subtropical Fitzroy River (Robson et al., 2006) and investigation of the impacts of a tropical cyclone (Condie et al., 2009), saw a stronger link to remote observations. At this time, the use of offline transport schemes were also implemented (such as the Moreton Bay model), allowing for an order-of-magnitude faster model integration (Gillibrand and Herzfeld, 2016).

The next major change in the BGC model involved implementing variable $\mathrm{C}:\mathrm{N}:\mathrm{P}$ ratios of microalgae through the introduction of reserves of energy, nitrogen and phosphorus (Wild-Allen et al., 2010), allowing for more accurate prediction of the elemental budgets and impacts of natural and anthropogenic forcing of the Derwent River estuary, southeast Tasmania. This study was followed up by a number of studies developing scenarios to inform management strategies of the region (Wild-Allen et al., 2011, 2013; Skerratt et al., 2013; Hadley et al., 2015a, b).

From 2010 onwards, EMS has been applied to consider the impacts of catchment loads on the Great Barrier Reef. The focus on water clarity led to the development of a spectrally resolved optical model and the introduction of simulated true colour (Baird et al., 2016b). The eReefs project was the first EMS application to consider corals with the introduction of the host–symbiont coral system and equilibrium carbon chemistry (Mongin and Baird, 2014; Mongin et al., 2016b, a). Additionally, the calculation of model outputs that match remote-sensing observations allowed the model to be run in a data-assimilating system, where the observation–model mismatch was based on remote-sensing reflectance (Jones et al., 2016). The most recent application of the EMS BGC model has been for investigating the environmental impact of aquaculture in Los Lagos, Chile (Steven et al., 2019a). For the Los Lagos application, new processes for fish farms, dinoflagellates and benthic filter feeders were added, although these additions are not described in this document. As a demonstration of the ability to add and remove processes, the Los Lagos application was run with the same EMS C-executable file as the Great Barrier Reef application, just with the configuration files altered.

13.2 Comparison with other marine biogeochemical models

As introduced earlier, there are a number of complex marine biogeochemical models. The most similar model in scope and approach to EMS is the ERSEM (European Regional Seas Ecosystem Model) model (Butenschön et al., 2016). Both ERSEM and EMS consider in detail pelagic, benthic and sediment processes and could generally be described as functional group models. That is, the state variables and the processes that link them are chosen to represent groups of organisms that act in similar ways. This allows the complexity of real systems to be reduced to a tractable model. Many functional group style biogeochemical models exist and were in fact the earliest models developed (Riley, 1947; Fasham, 1993; Sarmiento et al., 1993). The most significant differences between EMS and ERSEM are (1) EMS concentrates more on benthic flora than ERSEM, while ERSEM considers lower-trophic-level ecosystem interactions such as fisheries that are not captured in EMS; and (2) while EMS and ERSEM have similar state variables and processes, EMS has a different set of governing equations that are based on geometric constraints of individuals while ERSEM, like most other functional biogeochemical models, has equations based on empirical relationships determined from population interactions.

The last two decades have seen addition modelling approaches emerge: trait-based models that consider changing processes rates as populations vary (Bruggeman and Kooijman, 2007); size-based models that determine rates based on organism size (Baird et al., 2007a); ecosystem-style models that consider a multiple “species” within a functional group, developing large food webs (Fulton et al., 2014); and models that consider a large number of functional groups that are refined through competition between groups (Follows et al., 2007). These new approaches are applied primarily in pelagic ecosystems, where the generic nature of pelagic interactions encourages overarching philosophies to model construction and with considerable success (Dutkiewicz et al., 2015). The awkwardness of the variety of benthic communities (corals, seagrass, kelp, etc.) and their prime role in shallow water, has meant that estuarine and coastal models have, like ERSEM and EMS, typically chosen the functional model approach (Madden and Kemp, 1996; Spillman et al., 2007).

13.3 Future developments in EMS

EMS has been developed to address specific scientific questions in Australia's coastal environment. As a result, the set of processes the EMS considers varies from those typically applied by other groups developing marine BGC models. Processes which have not been considered but often are considered in marine BGC models, include iron and silicate limitation (which are not common on the Australian continental shelf or estuaries), photo-inhibition of microalgae, explicit bacterial biomass. Each of these will be considered as the need arises.

A deliberate decision in the development of the EMS BGC model was made to avoid higher-trophic-level processes, such fish dynamics and reproduction of long-lived species. This decision was made because (1) including these longer timescale, often highly non-linear, processes reduces the ability of development to concentrate on BGC processes; and (2) it was recognised that CSIRO has developed a widely used ecosystem model (Atlantis, https://research.csiro.au/atlantis/, last access: 22 September 2020; Fulton et al., 2014), and that coupling the EMS with Atlantis takes advantage of complementary strengths of the two modelling systems.

A recent capacity introduced to EMS is the development of a relocatable capability (RECOM; Sect. 12), allowing model configurations (grid, river and meteorological forcing, ecological processes, boundary conditions) to be automatically generated. This capability will be a good test of the portability of the BGC model and in particular the use of geometric description of physical limits to ecological processes.

Future enhancements in the EMS BGC model for tropical systems are likely to continue to pursue those components at risk from human impacts, such as dissolution of marine carbonates affecting sediment substrate and herbicide interactions with photosystems. We also expect to continue to refine the optical model and in particular the relationship between particle size distribution and mass-specific scattering and absorption properties. In temperate systems, current and near-future deployments of EMS in Australia will be focused on coastal system characterisation for aquaculture, carbon sequestration and management decision support for the blue economy. Ongoing research includes improved methods for model validation against observations and translation of model outputs into knowledge that informs stakeholder decisions.

13.4 Concluding thoughts

The BGC model in the CSIRO EMS has developed unique parameterisation when compared to other marine biogeochemical models applied elsewhere due in part to a unique set of scientific challenges of the Australian coastline. It has proved to be useful in many applications, most notably the Great Barrier Reef where extensive observational datasets has allowed new process model development and detailed model skill assessment (Baird et al., 2016b, a; Mongin et al., 2016b; Skerratt et al., 2019 and http://ereefs.info/, last access: 22 September 2020). This document provides easy access to some of the novel process formulations that have been important in this success, as well as a complete scientific description of version B3p0.