Introduction

Terrigenous dissolved organic matter (DOM) plays a major role in the storage and cycling of carbon (C) at regional and global scales (Tranvik et al. 2009). The flux of C from soils to freshwater systems may be as high as 1.9–5.1 Pg C year−1 (Cole et al. 2007; Drake et al. 2018), of which a significant, although as yet not well constrained, proportion is as DOM (Raymond and Spencer 2015; Massicotte et al. 2017; Drake et al. 2018). Calculations of the global carbon budget (e.g., Ciais et al. 2013), including those derived from earth system models (ESMs), do not, however, explicitly represent the processes involved in the transfer of dissolved organic carbon (DOC) throughout the land-ocean-aquatic-continuum (LOAC), which spans freshwater, estuaries and the ocean. One reason is that sequestration of terrigenous DOC in the ocean was considered negligible for many years: isotopic signatures and absence of traditional terrestrial biomarkers (e.g., lignin) suggested that the large pool of oceanic DOM (~ 700 Pg C: Hansell and Carlson 1998) is produced mostly in situ. Recent results from molecular fingerprinting techniques have, however, demonstrated the presence of otherwise undetectable (in the pico-to nanomolar concentration range) terrestrially-derived compounds far into the open ocean (Medeiros et al. 2016, Riedel et al. 2016). A second reason why DOC processes have not been afforded an explicit representation in calculations of the global C budget is the sheer difficulty of developing mechanistic models of DOC turnover that span the entire LOAC, given the complexity of the system and associated uncertainties (Le Quéré et al. 2015). The requirement is for unified concepts of DOM turnover, based on mechanistic understanding of the underlying processes, in order to develop a unified model that has a single set of state variables and parameterisations that are applicable from freshwaters to the ocean. There are in existence models of DOC for individual parts of the LOAC (e.g., Futter et al. 2007; Bauer et al. 2013), but the different research communities have adopted differing concepts and nomenclature.

Here, we develop a new cross-system model, UniDOM (Unified model of Dissolved Organic Matter), that unifies process representations (i.e., provides a common set of concepts, state variables and parameterisations) across freshwater, estuarine and ocean systems. The development of UniDOM has two main aims. Our first and primary objective is to identify knowledge gaps in the turnover of terrigenous DOC across the LOAC, stimulating discussion and the evolution of new ideas, including associated modelling approaches. “By forcing one to produce formulas to define each process and put numbers to the coefficients, [a simulation of a natural ecosystem] reveals the lacunae in one’s knowledge…the main aim is to determine where the model breaks down and use it to suggest further field or experimental work” (Steele 1974). Our second aim is to present UniDOM as a fully functional mathematical model that represents the bulk DOC pool, and to use the model to provide a preliminary assessment of the extent to which terrigenous DOC traverses the LOAC and penetrates into the ocean. For this purpose, UniDOM is implemented in a simple physical framework with residence times for freshwater, estuaries and the ocean that are representative of United Kingdom waters.

UniDOM: conceptual basis

The main focus of UniDOM is DOC turnover across the entire LOAC. The requirement is to develop a model structure and associated parameterisations that capture the decline in concentration, along with associated changes in composition, as terrigenous DOC is acted on by different turnover processes, namely biological consumption, photooxidation and flocculation. The selection of state variables is of paramount importance in this regard and, after much debate (including two workshops), we chose to represent terrigenous DOC as two fractions, T1, and T2, with a third state variable for an autochthonous aquatic fraction, A (Fig. 1). The thinking behind this structure is described in this section, noting that a key consideration is that the different DOC fractions and associated rate processes should be amenable to straightforward and high-throughput measurement in the field and laboratory. A full mathematical description of UniDOM is presented later.

Fig. 1
figure 1

Flow diagram of the UniDOM model showing production and turnover of DOC within the LOAC (rivers and lakes, estuaries and ocean): terrigenous fractions T1 and T2, and autochthonous (A)

Full size image

Terrigenous DOM typically contains a large fraction of aromatic compounds, most notably lignin, but also a variety of secondary metabolites (Verhoeven and Liefveld 1997). These compounds are mostly of structural origin and thereby resistant to biological decomposition. Lignin, which is produced by terrestrial and aquatic vascular plants (Webster and Benfield 1986; Kirk and Farrell 1987; Kalbitz et al. 2003), and associated structural polysaccharides such as cellulose and hemicellulose (Webster and Benfield 1986), are good examples, although turnover rates of these compounds can nevertheless be significant (Ward et al. 2017). These chemical structures also contribute to the coloured or chromophoric DOM fraction (CDOM) and render the compounds prone to photooxidation by ultraviolet (UV) radiation (Cory et al. 2015; Koehler et al. 2016; Berggren et al. 2018). For example, Opsahl and Benner (1998) noted that 75% of lignin was photo-oxidised over 28 days in the Mississippi River, while Spencer et al. (2009) noted a 95% reduction when Congo River water was irradiated for 57 days.

Aromatic compounds are most readily precipitated by metal ions given their susceptibility to sorption and co-precipitation with reactive minerals (Vilge-Ritter et al. 1999; Riedel et al. 2012). Flocculation may be highest in estuaries for several reasons. High concentrations of particulate matter in estuaries promotes flocculation, especially in the turbidity maximum (Vinh et al. 2018). High ionic strength also facilitates flocculation by neutralising negative surface charges on DOM. Metal salts enhance flocculation, e.g., iron is precipitated as waters of terrestrial origin enter estuaries (Charette and Sholkovitz 2002). Flocculation requires molecules to collide and so, along with deflocculation, depends on hydrodynamic forcing, notably Brownian motion and fluid shear, i.e., turbulence (Kepkay 1994; Wang et al. 2013). The greater turbulence within estuaries relative to lakes and the ocean therefore further contributes to the potential for estuarine environments as hotspots for DOM flocculation (Geyer et al. 2008).

Based on the above, we define T1 as the fraction that includes compounds such as lignin and other structural entities that are prone to photooxidation and flocculation, but relatively resistant to microbial decomposition. Conversely, we classify T2 compounds as those that are not readily photooxidised or flocculated, but which are more susceptible to microbial processing. The T1:T2 partition is further justified in context of carbon fluxes across the LOAC by the observation that it is aromatic fraction that dominates the signal of terrigenous DOM in the open ocean (Meyers-Schulte and Hedges 1986; Hedges et al. 1997; Opsahl and Benner 1997). It is tempting to use aromaticity as the specific criterion for defining the division of terrigenous DOM between T1 and T2, but we choose not to do so because photooxidation is not exclusively restricted to aromatic compounds. Carboxylic acids are also prone to photooxidation (Xie et al. 2004; Gonsior et al. 2013), as are some proteins (Janssen and McNeill 2014). Instead, we operationally define T1 as terrigenous DOM that strongly absorbs UV radiation, whereas T2 represents compounds that are only weakly or not absorbing for UV, such as many carbohydrates and proteins.

Autochthonous DOM is produced in situ throughout the LOAC via a range of processes including phytoplankton exudation, zooplankton grazing and detritus turnover (e.g., Anderson and Williams 1998). Algae typically do not contain lignin although lignin-like compounds, e.g., sporopollenin, are sometimes present (Gunnison and Alexander 1975). Recalcitrant compounds are nevertheless present in autochthonous DOM, such as detritus in marine systems that contains structural polysaccharides (Mann 1988) that require solubilisation by exoenzymes in order to be utilised by bacteria (Chróst 1990). CDOM is thought to be a by-product of microbial processing and metabolism (Rochelle-Newall and Fisher 2002; Nelson et al. 2004; Romera-Castillo et al. 2011; Kinsey et al. 2018), is present within autochthonous DOM and is ubiquitous in both freshwater systems (Kutser et al. 2005) and the world ocean (Nelson and Siegel 2013). A significant fraction of CDOM may be loosely classed as humic material (Nelson and Siegel 2013; Kellerman et al. 2018), although the aromatic content is usually lower than that of soil organic matter (McKnight et al. 2001). Autochthonous DOM also includes colloidal material such as gels and transparent exopolymer particles (TEP) that aggregates to form flocs of marine snow (Santschi 2018).

The inclusion of autochthonous DOC in UniDOM is required in order to meet our objective to provide a fully functional model that represents the bulk DOC pool. Noting that the main focus of our study is terrigenous DOC (T1 and T2), we adopt a simple approach whereby the autochthonous fraction is represented as a single state variable, A. As with DOM in general, autochthonous DOM comprises a spectrum of compounds of varying lability. A fraction is present as simple sugars that are rapidly consumed by bacteria on timescales of hours to days (Søndergaard and Middelboe 1995) and, as such, are unlikely to contribute significantly to lateral fluxes given that the transit time for DOM traversing the LOAC is usually at least a few days. Instead, and in common with other modelling studies of marine DOM (e.g., Schmittner et al. 2005; Salihoglu et al. 2008), our use of a single state variable for A represents semi-labile compounds that turn over on timescales of weeks to months. Compounds may be rendered semi-labile for a number of reasons including the requirement for exoenzyme hydrolysis, low concentrations of individual biomolecules, the presence of competing substrates, and bacterial community structure (Anderson et al. 2015a). Photooxidation is not usually considered as a loss process for labile or semi-labile DOM in marine biogeochemical models (Anderson et al. 2015a). A novel feature in UniDOM is that, given the contribution of CDOM, we assume that the A state variable is subject to photooxidation by UV. Rather than having two state variables, however, we instead assume that a fixed fraction of the total autochthonous pool is present as CDOM. In common with most ocean models, the formation of aggregates (TEP) is not included as a loss process for DOC in UniDOM (Anderson et al. 2015a).

It is important that UniDOM is amenable to linking with soil models in order to generate input fluxes of T1 and T2. Compounds that may be classified as T1 tend to be retained in the mineral soil layers through adsorption to mineral surfaces and co-precipitation (Kaiser and Kalbitz 2012). They are therefore typically released from topsoil during high rainfall events (Raymond and Saiers 2010) as relatively ‘young’ substrates, i.e., derived from photosynthetically fixed carbon during the last 20–30 years (Tipping et al. 2012). Older, but functionally similar, compounds may also be released due to ecosystem modifications such as drainage (Evans et al. 2014), cultivation and urbanisation (Butman et al. 2014), or permafrost melting (Neff et al. 2006). Fraction T2, on the other hand, is less susceptible to immobilisation and so a higher proportion percolates downwards through the soil profile. It therefore accrues from base flow emanating from the deep (mineral) soil layers, as well as being released during rapid flow events (Ward et al. 2012; Pereira et al. 2014). Man-made release of DOM into freshwaters and estuaries, notably from sewage, is an issue in many countries. For the sake of simplicity (avoiding the use of an additional state variable), we suggest that this DOM is assigned as T2 to fit with the strongly versus weakly UV-absorbing scheme for separating terrigenous organic matter.

Data

Fractions T1 and T2 should be amenable to measurement in order for UniDOM to link with observational programmes and field data. Here, we recommend the suitability of SUVA254 (mass-specific absorbance at 254 nm), due to its widespread measurement and use in the DOM literature. The approach involves a 2-component end-member calculation, noting that it can only be an approximation given the variability in the functional and optical composition of DOM (Her et al. 2003). Briefly, a two-component description of DOM is consistent with the monotonic relation observed for diverse sample of freshwater DOM between the ratio of absorbance at two UV wavelengths, and specific absorbance at the upper wavelength (for further details, see Supplementary Appendix 2). This feature of DOM absorbance spectra has been used to generate accurate predictions of DOC concentration for independent, globally diverse samples in which terrigenous DOM predominates (Thacker et al. 2008; Tipping et al. 2009; Carter et al. 2012). The SUVA254 end members employed are derived from a variety of empirical experimental and environmental data, a full description of which is given in Supplementary Appendix 2. They have values of 7.7 L mg C−1 m−1 for 100% T1, and 1.8 L mg C−1 m−1 for 100% T2, in which case the T1 proportion of DOC, fT1, is:

$$ {\text{f}}_{{{\text{T}}1}} \; = \;0.1695\;{\text{SUVA}}_{254} \;{-}\;0.3051 $$
(1)

Estimates of T1 and T2 based on SUVA254 for four large UK river catchments—the Avon, Tamar, Conwy and Halladale—are presented in Table 1, in order to demonstrate the feasibility of using field data to derive these fractions, as well as to provide a means of initialising the model. Sampling took place on the freshwater edge of the tidal limits. It is unsurprising that these rivers have contrasting T1:T2 ratios because their catchments differ in the relative proportions of major soil types, i.e., from a dominance of peats through to arable land covers, the latter providing a proxy for mineral soil coverage. Sampling was undertaken monthly from January to December 2017. For comparison, the data are accompanied by estimates for two headwaters catchments within the Conwy River system with contrasting peat cover, sampled monthly in the 2011 and 2012 hydrologic years (October to September). Laboratory analysis was undertaken for DOC concentration and UV absorbance, permitting the calculation of SUVA254. For further details, see Supplementary Appendix 1. Results demonstrate the utility of SUVA254 for estimating the T1 fraction, fT1, and show that there is considerable variability between sites. The highest T1 fraction, fT1 = 0.65, is seen at Afon Ddu, where soils are almost entirely dominated by peat. In contrast, fT1 = 0.13, the lowest estimate, is seen for the Avon, which is largely bounded by mineral soils.

Table 1 Mean annual DOC concentration (mg L−1), SUVA254 (L mg C−1 m−1) and T1 proportion (fT1), along with catchment area (km2), peat% and arable%, for four large UK catchments at locations close to their tidal limits (Avon, 50.747° N, − 1.782°E; Tamar, 50.531° N, − 4.222° E; Conwy, 53.107° N, − 3.792° E; Halladale, 58.479° N, − 3.903° E) and two headwater catchments (Afon Ddu, 52.976° N, − 3.835° E; Maenan, 53.171° N, − 3.800° E)
Full size table

Mathematical realisation

A full mathematical version of UniDOM is now presented. Differential equations are described in this section for T1, T2 and A, along with lists of model variables and parameters in Tables 2 and 3. The very act of seeking out mathematical parameterisations of the various processes of DOM turnover is itself conducive to revealing conceptual insights. When we initially constructed UniDOM, constant rates were assigned to DOC turnover via microbes and photooxidation, as is common practice (at least for the former) in many models of DOM (Anderson et al. 2015a). Recent empirical evidence, however, suggests that average turnover rates should decline with time, i.e., with DOM age, because biologically and photochemically labile compounds are progressively removed from solution, leaving recalcitrant ones behind (Catalán et al. 2016; Evans et al. 2017). A novel feature of UniDOM is that the decline in DOC turnover rate with age, φL, is described by:

Table 2 Model variables
Full size table
Table 3 Model parameters
Full size table
$$ \varphi_{L} = \frac{1}{{L^{\alpha } }} = L^{ - \alpha } $$
(2)

where parameter α defines the rate of decline. Note that this relationship is invalid for L = 0 (division by zero). The solution is to truncate at L = L0 (the age at which φL starts to decline) giving:

$$ \varphi_{L} = f_{L} \left( {L_{0} ,\alpha } \right) = \left( {L - L_{0} + 1} \right)^{ - \alpha } ,\;{\text{L}}\; > \;{\text{L}}_{0} $$
(3a)
$$ \varphi_{L} \; = \;f_{L} \left( {L_{0} ,\alpha } \right)\; = \;1,\;{\text{L}}\; \le \;{\text{L}}_{0} $$
(3b)

The predicted decline in φL with age (Eqs. 3a, 3b) is shown in Fig. 2, for α = 0.38, showing φL = 1, 0.23, 0.17 and 0.09 at 1, 50, 100 and 500 days, respectively. Extending to longer timescales (not shown), φL = 0.11, 0.044 and 0.018 at 1, 10 and 100 years.

Fig. 2
figure 2

Predicted decline in normalised DOC turnover (φL) with age (L) for α = 0.38 and L0 = 1 day: a linear axes, b log axes (slope for log(age) > 0 is − 0.38)

Full size image

The value assigned to parameter α is based on two relationships shown graphically in Evans et al. (2017; Fig. 3a therein) where, using log–log axes, DOC turnover rate is plotted against water residence time (a proxy for DOM age). The first is based on laboratory dark-incubation data (Catalán et al. 2016) and has a slope of − 0.38, representing microbial turnover but with no contribution from photooxidation. The second relationship is for 82 predominantly European and North American lakes and reservoirs with water residence times ranging from a week to 700 years. We fitted a slope of -0.76 to the linear part of the relationship. The data in both cases are for net DOC removal and therefore, especially in the case of the lake data, may underestimate turnover, and overestimate α, because of in situ production of DOC (Köhler et al. 2013; Evans et al. 2017). We will investigate the sensitivity of UniDOM to the full range of α, from α = 0 (no effect of DOC age on lability), to α = 0.76 (assumes that the entire decrease in reactivity across the age gradient in the lake input–output dataset is due to declining reactivity, with no offset due to production), using an average value, α = 0.38, as standard.

The differential equation for T1 involves production via soil and loss terms for photo-oxidation, flocculation and microbial decomposition:

$$ \frac{{dT_{1} }}{dt} = \sigma_{T1} P_{soil} - f_{L} \left( {L_{0} ,\alpha_{pox} } \right)\emptyset \left( {\emptyset^{ref} ,z_{col} ,k_{UV} } \right)T_{1} - \gamma_{X} T_{1}^{2} - f_{L} \left( {L_{0} ,\alpha_{mic} } \right)\eta_{1} T_{1} + h\left( {z_{col} } \right)T_{1} $$
(4)

Terrigenous DOC (T1 + T2) is produced via soils and released into headwater catchments (flux Psoil; mmol C m−3 day−1). Fraction σT1 of this release is as T1 (remainder as T2) and will vary for different soil types. Parameter σT1 was set as a weighted mean of data in Table 1 for the Avon and Tamar (mineral soil), Conwy (organic-mineral) and Halladale (organic), based on fractional cover in the UK of 43.7, 22.3 and 29.7% for organic, organic-mineral and organic soils, respectively (Bradley et al. 2005), giving σT1 = 0.31.

The second term in Eq. (4) is loss due to photooxidation that depends on a default maximum rate (ϕref) that decays with time (Eqs. 3a, 3b) and the rate of UV attenuation with depth (zcol, m) according to the UV extinction coefficient (kUV, m−1). Parameter ϕref can be estimated from measurements of DOC decay when exposed to UV irradiation: 14% decrease over 8 h, ϕ = 0.22 day−1 (Bertilsson and Tranvik 1998; Swedish Lakes), 40% over 21 days, ϕ = 0.025 day−1 (Shiller et al. 2006; Mississippi River), 84% over 67 days, ϕ = 0.043 day−1 (Helms et al. 2013; North Pacific water), 81% over 110 days, ϕ = 0.028 day−1 (Helms et al. 2014; Chesapeake, USA) and 50% over 1 day, ϕ = 0.35 day−1 (Gonsior et al. 2014; oligotrophic ocean waters). These ϕ values were calculated by integrating Eqs. (3a, 3b), with α = 0.38, assuming that half the loss is due to photooxidation and calculating the rate at t = 1 day (= L0). We use an average of these values, ϕref = 0.13 day−1, noting that all are from low latitude sites except Bertilsson and Tranvik (1998).

Attenuation of UV-B, the most photo-active part of the spectrum, is strong in natural waters (Graneli et al. 1996) and can be described by Beer’s Law. The extinction coefficient kUV (m−1) was derived using a typical z10% (depth at which 10% of UV remains) for open ocean waters of 20 m (Tedetti and Sempéré 2006), which gives a rate of extinction of 0.12 m−1 (parameter kUVw). The rate of UV extinction increase with DOC concentration can be described from SUVA350, which is the specific UV absorption at 350 nm per unit DOC (m2 (mmol C)−1; Massicotte et al. 2017), in which case kUV is:

$$ k_{UV} = k_{UVw} + SUVA_{350} \left( {T_{1} + \xi A} \right) $$
(5)

Note that a fraction of autochthonous DOM, ξ, is also assumed to absorb UV (see below). A typical value of SUVA350 for freshwater systems is 1.0 m2 (g C)−1 (Massicotte et al. 2017), which is equivalent to 0.012 m2 (mmol C)−1). This value is for total DOC and so, if instead normalised to T1, then a typical value could be 0.039 m2 (mmol C)−1, assuming that T1 is a 31% of the total (parameter σT1). For example, this value of SUVA350, combined with kUVw = 0.12 m−1, gives an extinction, kUV, of 7.92 m−1 for a DOC concentration of 200 mmol m−3, meaning that 55% of UV is attenuated at 0.1 m and 98% at 0.5 m. The average rate of photo-oxidation within a water column of depth zcol (m), ϕ(ϕref, zcol,kUV) is then:

$$ \emptyset \left( {\phi^{ref} ,z_{col} ,k_{UV} } \right) = \frac{{\emptyset^{ref} }}{{z_{col} }}\mathop \smallint \limits_{0}^{{z_{col} }} e^{{ - k_{UV} z}} dz $$
(6)
$$ \mathop \smallint \limits_{0}^{{z_{col} }} e^{{ - k_{UV} z}} dz = \frac{1}{{k_{UV} }} - \frac{{e^{{ - k_{UV} z_{col} }} }}{{k_{UV} }} $$
(7)

Flocculation is assumed only to apply to fraction T1, given that aromatic compounds are prone to removal from solution via sorption and co-precipitation with reactive minerals (Vilge-Ritter et al. 1999; Riedel et al. 2012). After due consideration, we concluded that a mechanistic model of flocculation involving electrostatic chemistry, molecular interactions, turbulent mixing, etc., is not practically feasible, especially across the entire LOAC, because of the difficulties of reliable parameterisation. Unlike microbial turnover and photooxidation, we resorted to the application of different flocculation rates in freshwater (F), estuaries (E) and the ocean (O), parameters γF, γE and γO, calculated as a square function of T1 (encounter rates between molecules increase with density). A flocculation rate of 0.05 day−1 for estuaries was observed by Asmala et al. (2014) at a salinity of 1 (DOC declined from 15.0 to 14.3 mg L−1 over 24 h). This therefore equates to 0.05/1250 (15 mg L−1 = 1250 mmol C m−3), giving γE ≈ 4.0 × 10−5 day−1 (mmol C m−3)−1. The extent of flocculation may, however, be highly variable from site to site. In some cases, DOC shows conservative behaviour in estuaries, with little turnover (e.g., Mantoura and Woodward 1983). We tentatively set γE ≈ 2.0 × 10−5 day−1 (mmol C m−3)−1, a value half way between zero and that derived from Asmala et al. (2014), recognising that there is considerable uncertainty in this parameter. The flocculation rate for lakes was assigned a value of 0.0002 day−1 (g m−3)−1 by Tipping et al. (2016); converting units gives γF = 2 × 10−6 day−1 (mmol C m−3)−1. We tentatively set the flocculation parameter for the ocean to the same value, i.e., γO = 2 × 10−6 day−1, while noting that realised rates will be much lower because of the lower concentrations involved. The formation of transparent exopolymer particles (TEP) by particle aggregation is not included in UniDOM, as is the case for most marine ecosystem models.

The third loss term in Eq. (4) is DOC utilisation by bacteria to fuel growth and metabolism. The rate of microbial degradation at time 1 day (L0 = 1) was estimated from the linear relationship for dark incubation experiments shown in Evans et al. (2017; their Fig. 3a, dashed line), giving 0.03 day−1 (note that Catalán et al. (2016) show data with the same relationship for residence time as low as 1 day). Separate decay rates are specified for T1 and T2. We assume that microbial degradation of fraction T1 is slow because of its high lignin content and thus use a value 1/3 for that of fraction T2 based on a comparison of bacteria respiration rates in peat-influenced water versus clear mountain waters (Berggren and del Giorgio 2015). If the ratio of T2:T1 is 69:31 (parameter σT1 = 0.31) and the average rate is 0.03 day−1, then η1 = 0.013 day−1 (for T1) and η2 = 0.038 day−1 (for T2). Finally, h(zcol) calculates changes in concentration as T1 is diluted when the depth of the water column increases as water moves from rivers to the estuary and the ocean (see next section).

The differential equation for state variable T2, is:

$$ \frac{{dT_{2} }}{dt} = (1 - \sigma_{T1} )P_{soil} + f_{L} \left( {L_{0} ,\alpha_{pox} } \right)\beta \emptyset \left( {\phi^{ref} ,z_{col} ,k_{UV} } \right)T_{1} - f_{L} \left( {L_{0} ,\alpha_{mic} } \right)\eta_{2} T_{2} + h\left( {z_{col} } \right)T_{2} $$
(8)

A fraction, β, of photo-oxidised T1 is allocated to T2 (the remainder is lost as CO2), with an assigned value of 0.24 (e.g., Medeiros et al. 2015; Mostovaya et al. 2016). There is assumed to be no flocculation of T2.

Our main focus is the turnover of terrigenous DOC (fractions T1, T2) and so we adopt a simple description of the autochthonous fraction, A, which completes the representation of bulk DOC in UniDOM. It is possible to describe autochthonous DOC production in terms of a range of processes including phytoplankton exudation, zooplankton “messy feeding”, detritus turnover, etc., and to describe its turnover using separate pools such as labile, semi-labile and refractory (Anderson et al. 2015a). We use a single state variable, A, with production specified as a simple fraction, σA, of primary production, and which is subject to microbial turnover and photooxidation as loss terms. The inclusion of photooxidation is a novel feature of UniDOM, introduced to provide consistency with the parameterisation of T1, given that a fraction of autochthonous DOC is coloured. Note that it is entirely possible for users of UniDOM to substitute the simple parameterisation of A presented herein with more complex representations of the autochthonous pool of DOC.

The equation for autochthonous DOC (A) is:

$$ \frac{dA}{dt} = \sigma_{A} PP - \phi \left( {\phi^{ref} ,z_{col} ,k_{UV} } \right)\xi \left( {1 - \beta } \right)A - \eta_{A} A $$
(9)

where the terms are production and losses due to photooxidation and microbial consumption (note that A is unaffected by dilution: see next section). The ultimate source is primary production (PP, mmol C m−3 day−1) that occurs throughout the LOAC and is calculated from a surface rate, PP0, that is integrated over depth according to the attenuation coefficient, kPAR:

$$ PP = \frac{{PP_{0} }}{z}\mathop \smallint \limits_{0}^{{z_{col} }} e^{{ - k_{PAR} z}} {\text{dz}} $$
(10)

The calculation does not include PP from macrophytes and benthic algae in rivers. We use PP0 = 1.4 mmol C m−3 day−1 based on primary production of 200 g C m−2 year−1 for the North Sea (Capuzzo et al. 2018) and kPAR of 0.046 m−1. State variable A is assumed to have properties equivalent to those of the semi-labile pool in the ocean in terms of its production and biological turnover. Parameter σA is then calculated as semi-labile DOC normalised to oceanic PP, giving a value of 0.40 when derived from the model of Anderson and Ducklow (2001), assuming that 90% of non-phytoplankton release processes are semi-labile (Anderson and Williams 1998).

The first loss term in Eq. (9) represents photooxidation of CDOM in which, as for photo-oxidation of T1, fraction β is returned to the DOC pool, hence only (1 − β) is lost to CO2. A novel aspect of the representation of autochthonous DOM in UniDOM is that a fixed fraction of A, ξ, is assumed to be subject to photo-oxidation (Coble 2007); we use 20%, i.e., ξ = 0.2 (note that this fraction influences UV attenuation in surface waters and thereby influences the rate of degradation of T1). Semi-labile DOM turns over on seasonal timescales (e.g., Anderson and Williams 1999; Hansell 2013). A microbial degradation rate of 0.006 day−1 at 0 °C was used by Tipping et al. (2016) when modelling macronutrient processing within temperate lakes, giving a mean residence time of 1/0.006 = 166 days. We use ηA = 0.012 day−1, assuming a typical temperature of 10 °C and a Q10 for temperature dependence of 2. Finally, the turnover rate of A is assumed not to decline with increasing DOM age, unlike fractions T1 and T2. It is not necessary to model the ageing fraction of autochthonous DOC in UniDOM, given our focus on fate of long-lived terrigenous DOM.

Model setup

The UniDOM model is tested in a simple physical framework that involves a three-phase representation of water movement from freshwater to estuary to coastal ocean (Fig. 3). Residence times are specified for each domain, RF1, RF2, RE and RO days for freshwater (rivers and lakes), estuaries and ocean, respectively. Water enters the freshwater domain at time t = 0, with initial concentrations for T1 and T2, and travels for RF days until it reaches the estuary. It is assumed that the major terrigenous inputs of DOC are from wet, uncultivated organic-rich soil types including peats. In our case study of UK waters, these occur predominantly in the upper reaches of catchments and so the current model configuration assumes that terrigenous DOC only enters the LOAC at time t = 0, without further additions en route. This simplifying assumption will need to be modified if the model were applied to other areas of the world, e.g. when DOC-rich water is produced in lowland wetlands (e.g., Wiegner and Seitzinger 2004), or where non-peat DOC sources are proportionally more important. The water then passes sequentially through the estuarine and ocean systems.

Fig. 3
figure 3

Physical setup. T1 and T2 are input at time t = 0 and water then travels through freshwater, estuary and ocean with residence times RF1 + RF2, RE and RO. Autochthnous DOC, A, is produced throughout the LOAC (not shown)

Full size image

We use default residence times for UK waters of 1, 0, 3 and 730 days for RF1, RF2, RE and RO, respectively (i.e., with no lakes; RF2 = 109 days will be examined separately). Residence times for UK rivers, without lakes or impoundments, range from hours to several days, e.g. Worrall et al. (2014) give a mean of 26.7 h. They are highly variable for estuaries, often < 3 days, although sometimes considerably longer (Uncles et al. 2002; Yuan et al. 2007). We use RE = 3 days, which is the median value of 34 UK estuaries tabulated in Uncles et al. (2002). The 730 day residence time for the (coastal) ocean is, for our illustrative case study, representative of the North Sea, which directly or indirectly receives much of the DOM exported from the UK land mass (Painter et al. 2018) and where the age distribution of water varies between 0 and 4 years (Prandle 1984; Blaas et al. 2001). In a second simulation, we include an average lake residence time (RF2) of 109 days, based on average properties of UK lakes obtained from the Centre for Ecology and Hydrology (CEH) database (https://eip.ceh.ac.uk/apps/lakes/), which gives a median retention time and depth of 108.8 days and 4.42 m, respectively.

Depths are assigned for the freshwater system (zF1 and zF2 for rivers and lakes) and the estuary (zE). We use a river depth, zF1, of 1 m (e.g., see National River Flow Archive: http://nrfa.ceh.ac.uk/), zF2 of 4 m (see above) and an estuary depth of 10 m (e.g., Prandle and Lane 2015). Whenever depth increases, e.g., as water enters the estuary, T1 and T2 are diluted assuming that the diluting water has zero concentrations of these variables. Fraction A is assumed to be present in all waters and so there is no dilution. Maximum ocean depth, zO, is set to 100 m, which is a characteristic value for the North Sea (Blaas et al. 2001), with an assumed linear taper from zE to zO throughout period RO (Fig. 3).

The complications of modelling DOC release from soils and its entry into freshwater are avoided by setting initial concentrations of T1 and T2 at time t = 0. If a typical DOC concentration in catchment water is 674 mmol m−3 (based on a weighted average of the values in Table 1, as for σT1) and σT1 = 0.31, then T1 and T2 are 209 and 465 mmol C m−3 respectively. Fraction T2 is also generated throughout the LOAC by photooxidation of T1 (the second term in Eq. 8). This newly produced T2 is assumed to have a “reactivity age”, L, of zero. A cohort approach is therefore employed whereby a new state variable for T2 is created for each model day. The model is coded in R, using the EMPOWER ecosystem model testbed (Anderson et al. 2015b; see this reference for information about the structure of the code, input and output files, and instructions for use). The UniDOM code, including supporting files, is listed in Supplementary Appendix 3.

Results

Predicted quantities of T1, T2 and A (mmol C m−2) across the LOAC, with residence times of 1, 3 and 730 days for freshwater, estuary and ocean, respectively, are shown in Fig. 4. Note that results are shown on a depth-integrated basis (m−2) because the decrease in terrigenous DOC over time is then solely determined by DOC loss processes, whereas concentration units (m−3) are also affected by dilution as the water column deepens. Relatively little degradation of DOC is predicted over the first 4 days (i.e., by the time water exits the estuary), with T1 declining from 209 to 191 mmol C m−2, and T2 from 465 to 414 mmol C m−2, after 4 days (Fig. 4a). Photooxidation, flocculation and microbial consumption accounted for 52.6, 2.0 and 45.4% of turnover of T1, respectively, while microbes were the sole loss term for fraction T2. The longer residence time in the ocean led to greater removal, with T1 decreasing to 40.5 and 20.2 mmol C m−2, and T2 to 46.7 and 13.3 mmol C m−2, after 1 and 2 years of ocean transit, respectively. Thus, after 2 years in the ocean, predicted terrigenous DOC is 33.5 mmol C m−2, i.e., equivalent to 5.0% of the DOC released into catchment waters. Photooxidation again accounted for the majority of the loss of T1 (50.2%), with 49.7% from microbial turnover and only 0.1% from flocculation. Predicted concentrations of T1 and T2 after 2 years in the ocean are 0.20 and 0.13 mmol C m−3, within a water column 100 m deep. Autochthonous DOM (A) was produced throughout the LOAC, leading to a concentration of 1.4 mmol C m−3 exiting the estuary (after 4 days) and 9.0 mmol C m−3 after 2 years in the ocean. Microbial consumption accounted for 86.9% of turnover of autochthonous DOM, with the remaining 13.1% due to photooxidation.

Fig. 4
figure 4

Predicted quantities of T1, T2 and A (mmol C m−2) across the LOAC, with relative contributions of photo-oxidation, flocculation and microbial turnover shown in pie charts: a residence times RF1 (river) 1 day, RF2 = 0, RE (estuary) 3 days, RO (ocean) 2 years, b RF1 = 1 day, RF2 = 109 days (lake), RE = 3 days, RO = 2 years

Full size image

The inclusion of a lake with a 109-day residence time led to a much greater decrease in predicted DOC within the freshwater domain, with predicted T1 and T2 at day 110 of 75.1 and 162.7 mmol C m−2, respectively (Fig. 4b). The percentage losses of T1 and T2 are almost equal: 64% versus 65%. Predicted quantities of T1 and T2 further declined to 10.8 and 9.7 mmol C m−2, respectively, after 2 years of oceanic transit.

Sensitivity analysis was carried out by varying model parameters ± 10%, focusing on predicted terrigenous DOC (T1 + T2) after 734 days (residence times as in Fig. 4a; the standard output, i.e., with parameters as in Table 2, is 20.2 + 13.3 = 33.5 mmol C m−2). Only six parameters showed sensitivity > 5% (i.e., T1 + T2 at day 734 changed by ± 5% or more): σT1 (initial T1 as fraction of T1 + T2), photooxidation parameters ϕref (default rate) and SUVA350 (attenuation of UV due to coloured DOM), the microbial decay rates of T1 and T2 (η1, η2) and α (age-related decline in turnover). The sensitivity of predicted terrigenous DOC (T1 + T2) after 2 years in the ocean (day 734, with residence times as in Fig. 4a) to parameters ϕref, SUVA350, η1 and η2 is shown in Fig. 5. Decreasing the rate of photooxidation and/or increasing SUVA350 (leading to faster attenuation of UV in the water column and so less photooxidation; note that there is no link between parameters SUVA350 and ϕref in the model) results in predicted terrigenous DOC increasing from the standard value of 33.5 mmol C m−2 (default parameters) to values over 50 mmol C m−2 (Fig. 5a). Even larger quantities were predicted in response to lowering rates of microbial consumption, with DOC > 120 mmol C m−2 (Fig. 5b). Sensitivity is high for parameter η2 because microbial consumption is the sole loss term for fraction T2.

Fig. 5
figure 5

Predicted T1 + T2 at day 734 as influenced by a photooxidation parameters ϕref (default rate) and SUVA350 (UV attenuation due to DOM) and b microbial degradation rates (parameters η1 and η2 for T1 and T2 respectively; at L = 0)

Full size image

A key feature of UniDOM is that decay rates due to photooxidation and microbial degradation decline over time because T1 and T2 are assumed to become progressively more recalcitrant with increasing age. Model sensitivity to the parameter that controls this decline, α, led to dramatic changes in predicted fluxes of DOC across the LOAC (Fig. 6). Setting α = 0 means that parameters ϕref, η1 and η2 remain fixed, independent of DOM age, at their default values of 0.13, 0.013 and 0.038 day−1, respectively. In this case, predicted terrigenous DOC declined to < 1% of the initial 674 mmol C m−2 by day 129, to < 0.01% by day 266, with only a trace (≈ 10−6 mmol C m−2) left by day 734. In contrast, α = 0.76 (i.e., a rapid decline in degradation rates with DOM age) gives rise to relatively little degradation over the entire LOAC, with predicted T1 and T2 of 137 and 246 mmol C m−2 remaining after 734 days.

Fig. 6
figure 6

Model sensitivity of predicted T1 and T2 to rates of decline of photooxidation and microbial degradation (parameter α; default is α = 0.38). Residence times are as for Fig. 4a (i.e., without the inclusion of lakes)

Full size image

Discussion

Our main aim of developing UniDOM was to highlight knowledge gaps regarding the turnover of DOC across the LOAC, from which to present new ideas to progress this field of research. Understanding DOC turnover is an ongoing challenge for scientists because DOC consists of a plethora of compounds that differ in their reactivity and residence times, and which undergo extensive transformations during transport within rivers, estuaries and the ocean (Hopkinson et al. 1998; Geeraert et al. 2016). Furthermore, the concept of turnover (lability) is a complex one, involving biochemical structure, substrate concentration, presence/absence of competing substrates, microbial requirements (and hence microbial community composition) and the presence or absence of light (Anderson et al. 2015a).

Terrigenous DOC is subject to several turnover processes as it transits the LOAC—microbial consumption, photooxidation and flocculation—and the relative influence of these processes may vary considerably across freshwaters, estuaries and the ocean. Many models use a single state variable to represent terrigenous DOM (e.g., Futter et al. 2007; Rowe et al. 2014; Tipping et al. 2016), whereas we suggest that at least two are required to capture the complexity of DOC cycling across the LOAC. A novel feature of UniDOM is that terrigenous DOC is divided into two pools, T1 and T2, based on UV-absorbance (essentially, coloured and non-coloured). Fraction T1 is strongly UV-absorbing and so prone to photooxidation as well as flocculation, but relatively resistant to microbial consumption, whereas T2 is weakly or not UV-absorbing, is not prone to photooxidation and flocculation and is relatively amenable to biological degradation. We propose that the T1:T2 ratio can be quantified by measuring SUVA254 and demonstrated feasibility by presenting preliminary data for UK rivers. The measurement is simple and pragmatic, suitable for use across a range of contrasting and complex aquatic matrices. The calculation of T1:T2 involves a 2-component end-member approach based on UV absorbance (for details, see Supplementary Appendix 2). Nominally, T1 represents terrigenous DOM that includes a large fraction of aromatic compounds, notably lignin, that are prone to photooxidation by UV (Cory et al. 2015; Koehler et al. 2016; Berggren et al. 2018). Photooxidation is not, however, restricted to these compounds and so the definitions of T1 and T2 based on UV absorbance do not necessarily map precisely with specific molecular structures. The relationship between T1 and SUVA254 cannot therefore be precise and universal.

The development of UniDOM highlights a perennial difficulty when it comes to understanding and quantifying DOM turnover, namely that DOM consists of a heterogeneous mix of substrates, and that the average composition tends to become relatively more recalcitrant with age as labile substrates are stripped out. This phenomenon takes on particular significance when attempting to quantify a progressively declining terrigenous DOC flux across the LOAC, including its circulation within the ocean. One solution is to divide DOM into pools of differing lability (Worrall and Moody 2014) such as the labile, semi-labile, refractory, terminology that is often used for ocean DOC (Kirchman et al. 1993; Carlson and Ducklow 1995; Hansell 2013). Any such division is, however, arbitrary and will have difficulty in capturing a long-lived pool that becomes progressively recalcitrant. Reactivity continuum modelling is an alternative approach where degradation rates are expressed as a continuous probability distribution of reactivity (Vähätalo et al. 2010; Koehler and Tranvik 2015; Vachon et al. 2017). It has not, however, been adopted in biogeochemical models because of the mathematical complexity involved. We developed a novel approach in UniDOM, where degradation rates of terrigenous DOC, both microbial and by photooxidation, are inversely proportional to DOM age as parameterised from empirical estimates of turnover versus water residence time for a wide variety of lakes (Catalán et al. 2016; Evans et al. 2017). Using α = 0.38 (the parameter that defines the rate of decline with DOC age), the turnover rate normalised to that of a freshly-produced substrate (age = 0) is 1.0, 0.l1, 0.044 and 0.018 day−1 after 0, 1, 10 and 100 years, respectively. Model predictions for DOC transfer across the LOAC showed high sensitivity to the value of α, with α = 0.38 (our standard value) resulting in 5% of the DOC released in catchments remaining after 2 years in the ocean, whereas 57% remains if α is doubled to 0.76. If, on the other hand, the decline in lability with age is ignored (α = 0), terrigenous DOC is rapidly depleted, with only a trace quantity (≈ 10−6 mmol C m−2) remaining after 2 years in the ocean. The whole issue of how multiple processes, notably microbial turnover and photooxidation, interact to influence lability is complex. Photooxidation may make otherwise recalcitrant biomolecules available for use by microbes (e.g. Moran and Zepp 1997; Anesio et al. 2005; Cory et al. 2014), although the process may also consume some molecules that would otherwise have been available for biodegradation (Bowen et al. 2019; Bittar et al. 2015). The two processes may alternatively operate more or less independently on different components of the DOM pool (Benner and Kaiser 2011). Our work highlights the value of the empirical syntheses of DOC turnover versus water residence times by Catalán et al. (2016) and Evans et al. (2017) for understanding the fate of long-lived DOC, while emphasising the need for further work to understand the mechanistic underpinning of these relationships in terms of underlying processes and their parameterisation in models.

Model results were sensitive not only to the age-related parameterisation, but also to base turnover rates for freshly produced DOC (age = 0): we used η1, η2 = 0.013, 0.038 day−1 (biological turnover, T1, T2) and ϕref = 0.13 day−1 (photooxidation, fraction T1). It is unsurprising that a wide range of values for these rates is present in the literature, given variability between soil types, temperature, irradiance, etc., between sites. Photooxidation was predicted to account for 50% of turnover of the T1 fraction across the entire LOAC, including the coastal ocean. Our work thus emphasises the need to include this process in models of the coastal ocean. Further, a novel feature of UniDOM is that a fraction (20%) of the autochthonous DOC pool was assumed to be coloured and subject to photooxidation, accounting for 13% of the resulting loss. Biological lability depends on multiple factors including low concentrations of individual biomolecules, presence of competing substrates and the structure and physiology of the bacterial community (Anderson et al. 2015a). The last of these could be of particular importance for DOC fluxes across the LOAC, although it is by no means easy to incorporate into UniDOM or other models. Priming is an example of the physiological flexibility of microbes, depending on circumstances (Bengtsson et al. 2018). The addition of simple sugars, such as those produced by algal exudation, can facilitate the breakdown of more refractory DOM (Bianchi et al. 2015; Ward et al. 2016).

A perhaps surprising outcome of our study is that flocculation was predicted to make only a minor contribution to the turnover of terrigenous DOC. Maximum rates were predicted for estuaries, but even there flocculation accounted for only 2.5% of DOC losses. The literature is equivocal on the quantitative significance of this process (e.g., Mantoura and Woodward 1983; Asmala et al. 2014) and our work serves to emphasise the need for further studies to resolve this uncertainty. We were unable to develop a mechanistic parameterisation of flocculation in UniDOM given that this process depends on a combination of chemical, physical and biological processes (Wang et al. 2013). Instead, simple quadratic functions were applied separately to freshwaters, estuaries and the ocean, where the values assigned to the rate parameters are only tentative.

The division of terrigenous DOC into fractions T1 and T2, and the associated parameterisations for turnover, were conceptualised on a universal basis. UniDOM is parameterised using literature data, substantially from the boreal northern hemisphere, and the model is investigated for a United Kingdom setting. Application of UniDOM at the global scale is desirable to develop a more complete understanding of DOM dynamics and land to ocean transfers. We currently lack suitable data sets to rigorously test UniDOM across the globe in a fully representative manner. A strength of the T1:T2 approach is that SUVA measurements correlate strongly with DOM aromaticity (Weishaar et al. 2003) and can therefore be used across both temperate and tropical environments (Graeber et al. 2015). Future work could extend SUVA measurements to a variety of contrasting systems in terms of soil type, residence time, etc., enabling model verification across a range of geographical settings. Relationships between model parameters and geographical characteristics could be sought in order to aid extrapolation to river systems that currently lack measurements. For example, microbial and photooxidation controls on DOM turnover could be related to temperature and irradiance.

The results presented herein provide a preliminary assessment, using a simple physical framework, of DOC transfer across the LOAC and the extent to which terrigenous DOC penetrates into the ocean. The predicted turnover of DOC was minimal in river systems, consistent with recent observations (Huntington et al. 2019). After 2 years of ocean transit, 5.0% of terrigenous DOC was predicted to remain, following 1 and 3 day transit times for freshwater and estuary, respectively. This DOC flux decreased significantly, to 3.0%, when an “average” UK lake, with residence time of 109 days, was included in the simulation, highlighting the importance of freshwaters in turnover of DOM and release of CO2 to the atmosphere. The residence time of DOM in freshwaters is hugely heterogeneous, ranging from a few days in short rivers draining high rainfall areas (such as some of those in UK) that lack significant lakes or impoundments, increasing in larger, dryer and/or impounded catchments, and rising to years, decades or even centuries in some lake-dominated areas (e.g., Evans et al. 2017 and references therein). This variation may give rise to order-of-magnitude variability in the predicted export of terrigenous DOM from land to the ocean, when combined with the uncertainty in the reactivity versus age relationships discussed above. Despite this uncertainty, the results presented indicate that the contribution of terrigenous DOC to the ocean pool may by no means be negligible. This finding is consistent with measureable concentrations of dissolved lignin in the open ocean (Hernes and Benner 2002; Medeiros et al. 2016) that indicate some terrigenous DOC traverses the LOAC without remineralisation, although this flux nevertheless remains poorly quantified and understood (Fichot et al. 2014; Medeiros et al. 2017).

In conclusion, UniDOM is the first model that is designed to explore the removal of terrigenous DOC across the entire LOAC, with preliminary results indicating that ~ 5% of DOC produced in soils may penetrate the open ocean. Its unique structure incorporates separate state variables to represent terrigenous DOM, T1 (strongly-UV-absorbing) and T2 (non- or weakly-UV-absorbing), which show different susceptibilities to microbial consumption, photooxidation and flocculation. A novel parameterisation was derived for microbial turnover and photooxidation whereby rates of these processes were inversely related to DOC age. Model results were sensitive to this parameterisation, as well as values assigned to parameters for rates of these processes for freshly-produced DOC. Predicted rates of flocculation were surprisingly low, although the sheer complexity of this process makes a mechanistic representation in models a nigh impossible task. Our work highlights the need for ongoing research into how different turnover processes impact DOC lability and how lability changes with age, all in context of the different environments of the LOAC. We anticipate that our findings will guide future attempts to refine understanding of DOC transport, remineralisation and transformation as it journeys from land to ocean.