CO₂ emissions from peat-draining rivers regulated by water pH

Southeast Asian peatlands represent a globally significant carbon store that is destabilized by land-use changes like deforestation and the conversion into plantations, causing high carbon dioxide (CO₂) emissions from peat soils and increased leaching of peat carbon into rivers. While this high carbon leaching and consequentially high DOC concentrations suggest that CO₂ emissions from peat-draining rivers would be high, estimates based on field data suggest they are only moderate. In this study, we offer an explanation for this phenomenon by showing that carbon decomposition is hampered by the low pH in peat-draining rivers. This limits CO₂ production in and emissions from these rivers. We find an exponential pH limitation that shows good agreement with laboratory measurements from high-latitude peat soils. Additionally, our results suggest that enhanced input of carbonate minerals increases CO₂ emissions from peat-draining rivers by counteracting the pH limitation. As such inputs of carbonate minerals can occur due to human activities like deforestation of river catchments, liming in plantations, and enhanced weathering application, our study points out an important feedback mechanism of those practices.

1 Introduction

Rivers and streams emit high amounts of carbon dioxide (CO₂) to the atmosphere (Cole et al., 2007), but estimates of these emissions (0.6–1.8 PgC yr⁻¹) are highly uncertain (Aufdenkampe et al., 2011; Raymond et al., 2013). Studies agree that more than three-quarters of global river CO₂ emissions occur in the tropics (Raymond et al., 2013; Lauerwald et al., 2015). River CO₂ emissions are controlled by the partial pressure difference between CO₂ in the atmosphere and the river water (Raymond et al., 2012). Riverine CO₂ is fed by the decomposition of organic matter leached from soils (Wit et al., 2015), by direct leaching of dissolved CO₂ from soil respiration (Abril and Borges, 2019; Lauerwald et al., 2020), and by photomineralization of dissolved organic carbon (DOC) (Nichols and Martin, 2021; Zhou et al., 2021). Studies suggest Southeast Asia as a potential hotspot for river CO₂ emissions (Raymond et al., 2013; Lauerwald et al., 2015) due to the presence and degradation of carbon-rich peat soils. However, measurements of river CO₂ emissions from this region are sparse.

More than half of the known tropical peatlands are located in Southeast Asia (Page et al., 2011; Dargie et al., 2017), whereby 84 % of these are Indonesian peatlands, mainly on the islands of Sumatra, Borneo, and Irian Jaya (Page et al., 2011). Already in 2010, land-use change affected 90 % of the peatlands located on Sumatra and Borneo (Miettinen and Liew, 2010) and turned them from CO₂ sinks to CO₂ sources (Hooijer et al., 2010; Miettinen et al., 2017; Hoyt et al., 2020). Enhanced decomposition in disturbed peatlands additionally increases the leaching of organic matter from soils into peat-draining rivers (Moore et al., 2013; Rixen et al., 2016; Cook et al., 2018). According to Regnier et al. (2013), land-use change remobilizes (1.0±0.5) Pg of soil organic carbon per year, of which 40 % is decomposed in rivers and emitted as CO₂ to the atmosphere. The resulting CO₂ emissions of 0.4 PgC yr⁻¹ represent 33 % of the total CO₂ emissions from rivers (Regnier et al., 2013).

Peat soils are rich in carbon, causing high concentrations of DOC in peat-draining rivers that increase with increasing peat coverage of the river catchments (Wit et al., 2015). However, despite high carbon leaching rates that cause DOC concentrations which can be more than 4 times higher than those in temperate regions (Butman and Raymond, 2011; Müller et al., 2015; Gandois et al., 2020), measured CO₂ fluxes from tropical rivers with high peat coverage (18–41 $\mathrm{gC}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$) hardly exceed those measured for rivers in temperate regions (18.5 $\mathrm{gC}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$, Müller et al., 2015; Butman and Raymond, 2011). Different reasons for this were suggested in literature. Müller et al. (2015) suggested short residence times of peat-derived DOC in rivers due to the location of peatlands near the coast as a possible cause. Other suggestions are the recalcitrant nature of DOC (Müller et al., 2016) and the lack of oxygen (O₂, Wit et al., 2015), which both lower the rate of DOC decomposition. Moreover, Borges et al. (2015) suggested a limitation of bacterial production and the resulting DOC decomposition in African peat-draining rivers as a consequence of low pH based on observations at rivers in the Congo basin.

The assumption of low O₂ concentrations and pH as a cause for moderate CO₂ emissions is supported by the regulating effect of these parameters on decomposition rates in peat soils. pH and O₂ are the key parameters that limit the activity of the decomposition-impelling enzyme phenol oxidase (Pind et al., 1994; Freeman et al., 2001). Phenol oxidase is needed to decompose phenolic compounds that are especially present in tropical peat soils (Hodgkins et al., 2018; Yule et al., 2018). Those phenolic compounds are more rapidly decomposed in the upper layer of peat soils than in deep peat (Gandois et al., 2014). Studies agree that the limiting effect of oxygen on decomposition is accurately represented by the Michaelis–Menten kinetics (Fang and Moncrieff, 1999; Pereira et al., 2017). This approach assumes that DOC decomposition is linearly limited for low O₂ concentrations but that there is no limitation for higher O₂ concentrations once they are sufficient to meet the decomposition demands (Keiluweit et al., 2016). Due to high rates of decomposition caused by the carbon-rich environment and low rates of photosynthesis caused by low nutrient concentrations and dark water colours that limit light availability to algae, peat-draining rivers are usually undersaturated with regard to atmospheric O₂ (Wit et al., 2015; Baum and Rixen, 2014). Still, their O₂ concentrations exceed those in peat soils due to gas exchange with the atmosphere (Müller et al., 2015; Rixen et al., 2008) and thus are assumed to limit decomposition less strongly than in peat soils (Pind et al., 1994). The same applies to the pH limitation, as pH in peat-draining rivers is usually higher than in peat soils (Pind et al., 1994). Other than for O₂ limitation, however, the form of the pH limitation is still subject to discussion. Linear (Sinsabaugh, 2010) and exponential (Williams et al., 2000; Kang et al., 2018) correlations have been stated in the literature.

This study aims at quantifying the limiting impact of pH and O₂ on the DOC decomposition in peat-draining rivers to explain the moderate CO₂ emissions observed from these rivers. We analysed data from 10 Southeast Asian peat-draining rivers with DOC concentrations between 200 and 3000 µmol L⁻¹ and pH and O₂ concentrations ranging from 3.8 to 7.1 and from 50 to 200 µmol L⁻¹, respectively.

2 Materials and methods

This study's methods were separated into two parts. The first part provides information on the study area, conducted measurement campaigns, and collected data that our analyses are based on. The second part describes the processes and equations used to quantify the decomposition dependency on O₂ and pH.

2.1 Measurement campaigns and study area

2.1.1 Study area

Southeast Asian peatlands store 42 Pg soil carbon across an area of 271 000 km² (Hooijer et al., 2010). More than 97 % of these peat soils are located in lowlands (Hooijer et al., 2006). The development of peatlands in Southeast Asia is favoured by its tropical climate with high precipitation rates that range between 120 mm in July and 310 mm in November with an annual mean of 2700 mm yr⁻¹ (Yatagai et al., 2020). Due to land-use changes like deforestation and the conversion into plantations, today less than one-third of those Southeast Asian peatlands remain covered by peat swamp forests, while in 1990 it was more than three-quarters (Miettinen et al., 2016). Southeast Asian rivers mostly originate in mountain regions and cut through coastal peatlands on their way to the ocean (Fig. 1). Measurement data included in this study were obtained in river parts that flow through peat soils to capture the influence of peatlands on the carbon dynamics in the rivers. The impact of sampling locations and seasonality is discussed in Appendix B.

Figure 1 Map of river catchments with the location of peat areas. Blue lines indicate the main rivers. Blue shaded areas outline the river basins, and brown areas indicate peatlands. Coloured data points indicate the sampling stations of the individual campaigns.

The collective data for this study were derived from four rivers on Borneo (Sarawak, Malaysia) and six rivers on Sumatra (Indonesia). The investigated rivers on Borneo are the Rajang, Simunjan, Sebuyau, and Maludam, and the rivers surveyed on Sumatra are the Rokan, Kampar, Indragiri, Batang Hari, Musi, and Siak (Fig. 1). We additionally include data from the Siak's tributaries Tapung Kiri, Tapung Kanan, and Mandau. River peat coverages range from 4 % in the Musi catchment to 91 % in the Maludam catchment, whereby the bigger rivers that originate in the uplands generally have lower peat coverages than smaller coastal rivers.

2.1.2 River expeditions and measured parameters

Data were derived from a total of 16 campaigns in Sumatra and Sarawak (Fig. 1, Table A1 in the Appendix). For the Indonesian rivers, 10 measurement campaigns between 2004 and 2013 were conducted. We use published data from Baum et al. (2007) for the Mandau, Tapung Kanan, and Tapung Kiri rivers; from Wit et al. (2015) for the Siak, Indragiri, Batang Hari, and Musi rivers; and from Rixen et al. (2016) for the Rokan and Kampar rivers. CO₂ measurements are available for the campaigns performed after 2008.

For the Malaysian rivers, measurements were performed in six campaigns between 2014 and 2017. We use data published by Müller-Dum et al. (2019) and Martin et al. (2018) for the Rajang River and by Müller et al. (2015) for the Maludam campaigns in 2014 and 2015. Additional campaigns for this study were conducted in March 2015 at the Simunjan and Sebuyau rivers as well as in January 2016 and March and July 2017 at the Simunjan, Sebuyau, and Maludam rivers. Measurements of DOC, CO₂, and O₂ concentrations as well as water pH, water temperatures (T), and gas exchange coefficients (k₆₀₀) for these additional campaigns were performed in the same manner as during the 2014 Maludam campaign (Müller et al., 2015). However, due to instrumental problems, the CO₂, O₂, and pH data measured at the Simunjan River in 2016 were not available for our analysis. Table 3 lists the averaged river parameters, including the catchments' peat coverages and atmospheric CO₂ fluxes.

During the January 2016 and March and July 2017 campaigns, concentrations of particulate inorganic carbon (PIC) in the form of CaCO₃ were measured in addition to the other parameters. These data are not included in the before-mentioned studies. Therefore, we describe the measurement principle here. Discrete water samples, taken from approximately 1 m below the water surface, were filtered through pre-weighed and pre-combusted glass fibre filters (0.7 µm) to sample particulate material within the water volume. To determine the particulate carbon (organic and inorganic), the samples were catalytically combusted at 1050 ^∘C, and combustion products were measured by thermal conductivity using an Euro EA3000 elemental analyser. The PIC was determined from the difference between this total particulate carbon and particulate organic carbon that was measured after addition of 1 M hydrochloric acid in order to remove the inorganic carbon from the sample.

2.1.3 River catchment size and peat coverage

River catchment sizes were derived from Hydro-SHEDS (Lehner et al., 2008) at 15 s resolution in WGS 1984 Web Mercator Projection. Sub-basins belonging to the catchments were identified using the HydroSHEDS 15 s flow direction data set and added to the main basins. The estimated accuracy of final catchment lines is 0.4 %.

Catchment peat coverage was derived from peat maps downloaded from https://www.globalforestwatch.org/ (last access: 10 December 2018) for Indonesia and Malaysia. The Indonesian peatland map was published by the Indonesian Ministry of Agriculture (MoA, 2012). The Malaysian peat map was made available by (Wetlands International, 2004) and is based on a national inventory by the Land and Survey Department of Sarawak (1968). Both maps include peatlands in different conditions, from undisturbed peat swamp forest to disturbed peat under plantations, which is nowadays widespread in those countries. Peat coverage was determined from the areal extent of peatlands in the catchment divided by catchment size. Peat coverages derived using other peat maps are compared in Appendix C.

2.2 Analysis of decomposition dependency on pH and O₂

Decomposition dependencies on pH and O₂ were derived based on the assumption that concentrations of dissolved inorganic carbon (DIC) and O₂ in peat-draining rivers, as a first approximation, are derived from an equilibrium between gas exchange with the atmosphere and DOC decomposition in the river water. This approximation assumes photosynthetic CO₂ consumption, photomineralization, and direct CO₂ leaching from soils to be negligible. We discuss the impact of these processes later on. In this section, we introduce the calculation of atmospheric gas exchange fluxes and decomposition rates. Then we derive equations to quantify decomposition limitations by pH and O₂ based on an equilibrium between these two processes.

2.2.1 Gas exchange between rivers and the atmosphere

Atmospheric CO₂ fluxes from rivers were calculated from CO₂ gas exchange coefficients and river CO₂ concentrations according to

$$\begin{array}{}\text{(1)}& F_{{\text{CO}}_{\mathrm{2}}}=k_{{\text{CO}}_{\mathrm{2}}}\left(T\right)\cdot \left({\text{CO}}_{\mathrm{2}}-K_{{\text{CO}}_{\mathrm{2}}}\left(T\right)\cdot p{\text{CO}}_{\mathrm{2}}^{\text{a}}\right).\end{array}$$

Exchange coefficients for CO₂ ($k_{{\text{CO}}_{\mathrm{2}}}\left(T\right)$) were calculated from k₆₀₀ and water temperature according to Wanninkhof (1992) as

$$\begin{array}{}\text{(2)}& k_{{\text{CO}}_{\mathrm{2}}}=k_{\mathrm{600}}\cdot {\left({\displaystyle \frac{\mathrm{1911.1}-\mathrm{118.11}\cdot T+\mathrm{3.4527}\cdot T^{\mathrm{2}}-\mathrm{0.041320}\cdot T^{\mathrm{3}}}{\mathrm{600}}}\right)}^{-n}.\end{array}$$

An exponent of $n=\mathrm{1}/\mathrm{2}$ (valid for rough surfaces; Zappa et al., 2007) was used for the rivers. The temperature T is given in degrees Celsius. pCO${}_{\mathrm{2}}^{\text{a}}$ is the atmospheric partial pressure of CO₂ (≈400 µatm), and $K_{{\text{CO}}_{\mathrm{2}}}$ describes the temperature-dependent Henry coefficient for CO₂, which was calculated according to Weiss (1974) as

$$\begin{array}{}\text{(3)}& \ln K_{{\text{CO}}_{\mathrm{2}}}=-\mathrm{58.0931}+\mathrm{90.5069}\cdot {\displaystyle \frac{\mathrm{100}}{T}}+\mathrm{22.2940}\cdot \ln \left({\displaystyle \frac{T}{\mathrm{100}}}\right).\end{array}$$

$k_{{\text{CO}}_{\mathrm{2}}}$ and $K_{{\text{CO}}_{\mathrm{2}}}$, derived for the individual rivers, are listed in Table A2.

Atmospheric O₂ fluxes ($F_{{\text{O}}_{\mathrm{2}}}$) were derived analogously to $F_{{\text{CO}}_{\mathrm{2}}}$, with $k_{{\text{O}}_{\mathrm{2}}}\left(T\right)$ calculated according to Wanninkhof (1992) as

$$\begin{array}{}\text{(4)}& k_{{\text{O}}_{\mathrm{2}}}=k_{\mathrm{600}}\cdot {\left({\displaystyle \frac{\mathrm{1800.6}-\mathrm{120.10}\cdot T+\mathrm{3.7818}\cdot T^{\mathrm{2}}-\mathrm{0.047608}\cdot T^{\mathrm{3}}}{\mathrm{600}}}\right)}^{-n}\end{array}$$

and Henry coefficients for O₂ ($K_{{\text{O}}_{\mathrm{2}}}$) calculated according to Weiss (1970) as

$$\begin{array}{}\text{(5)}& \ln K_{{\text{O}}_{\mathrm{2}}}=-\mathrm{58.3877}+\mathrm{85.8079}\cdot {\displaystyle \frac{\mathrm{100}}{T}}+\mathrm{23.8439}\cdot \ln {\displaystyle \frac{T}{\mathrm{100}}}.\end{array}$$

$k_{{\text{O}}_{\mathrm{2}}}$ and $K_{{\text{O}}_{\mathrm{2}}}$ for the individual rivers in this study are also listed in Table A2.

2.2.2 Decomposition rates and their dependency on pH and O₂

The decomposition rate of DOC (R) is defined as molecules of CO₂ that are produced per available molecule of DOC during a specific time step and thus represents the proportionality factor between the CO₂ production rate and the DOC concentration:

$$\begin{array}{}\text{(6)}& R={\displaystyle \frac{\mathrm{\Delta }{\text{CO}}_{\mathrm{2}}}{\text{DOC}\cdot \mathrm{\Delta }t}}\Rightarrow {\displaystyle \frac{\partial {\text{CO}}_{\mathrm{2}}}{\partial t}}=R\cdot \text{DOC}.\end{array}$$

As discussed before, R can be limited by O₂ concentrations and by pH. We use an O₂ limitation factor that is based on the Michaelis–Menten equation ($L_{{\text{O}}_{\mathrm{2}}}=\frac{{\text{O}}_{\mathrm{2}}}{K_{\text{m}}+{\text{O}}_{\mathrm{2}}}$) as suggested by Pereira et al. (2017). For pH limitation, we consider an exponential limitation factor ($L_{\text{pH}}^{\text{exp}}=\exp \left(\mathit{\lambda }\cdot \left(\text{pH}-{\text{pH}}_{\mathrm{0}}\right)\right)$) as suggested by Williams et al. (2000) and a linear limitation factor ($L_{\text{pH}}^{\text{lin}}=\frac{\text{pH}}{{\text{pH}}_{\mathrm{0}}}$) as suggested by Sinsabaugh (2010). Considering the definition of pH as the negative decadic logarithm of H⁺ activity ({H⁺}), the exponential limitation factor is equivalent to a linear correlation with ${\left\{{\text{H}}^{+}\right\}}^{\frac{\mathit{\lambda }}{\ln \left(\mathrm{10}\right)}}$.

The CO₂ production rates due to DOC decomposition for the linear and the exponential pH limitation approach are thus defined as follows.

$$\begin{array}{}\text{(7)}& \begin{array}{rl}& {\left({\displaystyle \frac{\partial {\text{CO}}_{\mathrm{2}}}{\partial t}}\right)}_{\text{lin}}=R_{\text{max}}\cdot L_{{\text{O}}_{\mathrm{2}}}\cdot L_{\text{pH}}^{\text{lin}}\cdot \text{DOC}\\ & =R_{\text{max}}\cdot {\displaystyle \frac{{\text{O}}_{\mathrm{2}}}{K_{\text{m}}+{\text{O}}_{\mathrm{2}}}}\cdot {\displaystyle \frac{\text{pH}}{{\text{pH}}_{\mathrm{0}}}}\cdot \text{DOC}\\ & {\left({\displaystyle \frac{\partial {\text{CO}}_{\mathrm{2}}}{\partial t}}\right)}_{\text{exp}}=R_{\text{max}}\cdot L_{{\text{O}}_{\mathrm{2}}}\cdot L_{\text{pH}}^{\text{exp}}\cdot \text{DOC}\\ & =R_{\text{max}}\cdot {\displaystyle \frac{{\text{O}}_{\mathrm{2}}}{K_{\text{m}}+{\text{O}}_{\mathrm{2}}}}\cdot \exp \left(\mathit{\lambda }\cdot \left(\text{pH}-{\text{pH}}_{\mathrm{0}}\right)\right)\cdot \text{DOC}\end{array}\end{array}$$

R_max is the maximum decomposition rate. K_m is the Michaelis constant for O₂ inhibition. It is also called the half-saturation constant and gives the O₂ concentration at which O₂ limits decomposition by 50 % (Loucks and Beek, 2017). λ is the exponential pH inhibition constant, and pH₀ is a normalization constant that was set to 7.5 since this is reported to be the optimal pH for the activity of the decomposition-impelling enzyme phenol oxidase (Pind et al., 1994; Kocabas et al., 2008). Calculations of pH₀ based on our data and the exponential pH approach are described in Appendix D4. They yield an optimum pH of approximately 7.2 and thus agree well with the pH₀ of 7.5 used in this study.

$L_{{\text{O}}_{\mathrm{2}}}$ and L_pH can take values between 0 and 1. Thus, Eq. (7) is only valid for pH ⩽ pH₀. For higher water pH, a different approach would be needed. However, for the rivers in this study Eq. (7) is sufficient since their pH is <7.5 (Table 3). The limitation factors represent the fraction of decomposition that remains after the limitation by the parameter. Later on, we refer to the fraction by which decomposition is limited, which is (1−L_pH) for pH limitation and ($\mathrm{1}-L_{{\text{O}}_{\mathrm{2}}}$) for O₂ limitation. The total fraction by which pH and O₂ limit decomposition is given by ($\mathrm{1}-L_{\text{pH}}\cdot L_{{\text{O}}_{\mathrm{2}}}$). When O₂ concentrations and water pH are high enough not to limit the decomposition rate, Eq. (7) simplifies to Eq. (6) with R=R_max.

2.2.3 Least-squares optimization to quantify the pH and O₂ impact on decomposition rates

As mentioned before, we base our calculations on the assumption that DIC concentrations in peat-draining rivers result from an equilibrium between CO₂ emissions and CO₂ production by decomposition. Thus, we optimized the parameters in Eq. (7) such that the production of CO₂ in the water volume beneath a specific surface area equals the atmospheric CO₂ flux through this area. The CO₂ production is calculated by multiplication of Eq. (7) with the product of river depth d and surface area A, and the CO₂ emissions are calculated by multiplication of Eq. (1) with the surface area A:

$$\begin{array}{}\text{(8)}& \begin{array}{rl}& d\cdot A\cdot R_{\text{max}}\cdot L_{{\text{O}}_{\mathrm{2}}}\cdot L_{\text{pH}}\cdot \text{DOC}=A\cdot k_{{\text{CO}}_{\mathrm{2}}}\left(T\right)\\ & \cdot \left({\text{CO}}_{\mathrm{2}}-K_{{\text{CO}}_{\mathrm{2}}}\left(T\right)\cdot p{\text{CO}}_{\mathrm{2}}^{\text{a}}\right).\end{array}\end{array}$$

Analogously, river O₂ concentrations result from an equilibrium between the atmospheric O₂ flux and O₂ consumption due to decomposition. During decomposition, the O₂ consumption is proportional to the CO₂ production ($\mathrm{\Delta }{\text{O}}_{\mathrm{2}}=-b\cdot \mathrm{\Delta }{\text{CO}}_{\mathrm{2}}$). The proportionality factor b is usually <1 since a fraction of the O₂ used for decomposition is taken from the oxygen content in the dissolved organic matter (Rixen et al., 2008). Thus, the equilibrium between O₂ consumption within the water volume and O₂ flux through the surface area can be written as

$$\begin{array}{}\text{(9)}& \begin{array}{rl}& -b\cdot d\cdot A\cdot R_{\text{max}}\cdot L_{{\text{O}}_{\mathrm{2}}}\cdot L_{\text{pH}}\cdot \text{DOC}=A\cdot k_{{\text{O}}_{\mathrm{2}}}\left(T\right)\\ & \cdot \left({\text{O}}_{\mathrm{2}}-K_{{\text{O}}_{\mathrm{2}}}\left(T\right)\cdot p{\text{O}}_{\mathrm{2}}^{\text{a}}\right).\end{array}\end{array}$$

In order to compare these dependencies to measured data, Eq. (8) and Eq. (9) were analytically solved for CO₂ and for O₂, respectively. The resulting equations based on linear pH limitation ($L_{\text{pH}}^{\text{lin}}=\frac{\text{pH}}{{\text{pH}}_{\mathrm{0}}}$) are listed in Table 1. The analogously derived equations for CO₂ and O₂ based on the exponential pH approach $\left(L_{\text{pH}}^{\text{exp}}=\exp \left(\mathit{\lambda }\cdot \left(\text{pH}-{\text{pH}}_{\mathrm{0}}\right)\right)\right)$ are listed in Table 2.

Table 1 Equations to derive CO₂ and O₂ based on the linear pH approach.

Equations to derive CO₂ from measured temperature (T), DOC, pH, and O₂ as well as to derive O₂ from measured T, DOC, and pH. R_max, K_m, and b, derived via least-squares optimization using measured DOC, pH, T, O₂, and CO₂ data of the investigated rivers, are listed in Table 5.

Table 2 Equations to derive CO₂ and O₂ based on the exponential pH approach.

Equations to derive CO₂ from measured temperature (T), DOC, pH, and O₂ as well as to derive O₂ from measured T, DOC, and pH. R_max, K_m, λ, and b, derived via least-squares optimization using measured DOC, pH, T, O₂, and CO₂ data of the investigated rivers, are listed in Table 5.

Based on these equations, least-squares optimizations were performed to derive the decomposition parameters R_max, b, K_m, and λ such that CO₂(DOC, pH, O₂, T) and O₂(DOC, pH, T) are simultaneously optimized for the measured parameters of DOC, pH, T, CO₂, and O₂.

The equations in Tables 1 and 2 depend on the river gas exchange coefficients for CO₂ ($k_{{\text{CO}}_{\mathrm{2}}}$) and O₂ ($k_{{\text{O}}_{\mathrm{2}}}$), which both depend on k₆₀₀. Those exchange coefficients are poorly constrained and spatially as well as temporally extremely variable. The k₆₀₀ values we list in this study are based on a variety of techniques, including floating chamber measurements (Müller et al., 2015), calculations based on wind speed and catchment parameters (Müller-Dum et al., 2019), and balance models of water parameters (Rixen et al., 2008). Although all of those estimates remain highly uncertain, we find a fairly good agreement between k₆₀₀ and river depths (d, Fig. A1). We therefore use a fixed ratio of $k_{\mathrm{600}}/d=(\mathrm{7.0}\pm \mathrm{0.5})\times {\mathrm{10}}^{-\mathrm{6}}{\mathrm{s}}^{-\mathrm{1}}$ for the least-squares optimizations rather than individual exchange coefficients and depths of the rivers.

3 Results

3.1 Correlation with peat coverage

The data presented in Table 3 yield a linear increase in river DOC concentration with peat coverage (Fig. 2a) as well as a negative linear correlation between river pH and peat coverage (Fig. 2b). The river CO₂ concentration shows a strong increase for peat coverages <30 %. Despite further increase in DOC concentrations, CO₂ concentrations in rivers with peat coverage >30 % level off, resulting in a fairly constant CO₂ for peat coverages >50 % (Fig. 2c). The river O₂ shows an opposite behaviour to the CO₂. O₂ concentrations initially decrease with increasing peat coverage and show a decline in the regression rate for high peat coverages, resulting in a minimum O₂ concentration of approximately 65 µmol L⁻¹ (Fig. 2d).

Table 3 Measured data from the investigated rivers.

Values are means of river campaigns. Data variability is given by the standard deviation of the measurements. For rivers with campaigns exclusively before 2009 there are no data (n.d.) on river CO₂ concentrations or fluxes. ^a For the Simunjan, the March 2015 and July 2017 campaigns (Simunjan₁) were separated from the January 2016 and March 2017 campaigns (Simunjan₂) due to strong differences in the parameters. ^b Due to technical problems during the Simunjan campaign in January 2016, these values are only based on one measurement campaign.

Figure 2 Correlation of peat coverage (PC) with (a) DOC, (b) pH, (c) CO₂, and (d) O₂. Each data point represents one river. Variability is indicated by the error bars, which are given by standard deviation. For the Simunjan River, the January 2016 and March 2017 campaigns (Simunjan₂; see Tables 3 and 4), indicated by red data points, were separated from the other Simunjan campaigns (Simunjan₁) and excluded from the correlations due to strong deviations from the other campaigns that imply an additional process discussed in Sect. 4.4. Ordinary least-squares optimizations were used to calculate linear correlations with DOC and pH and exponential correlations with CO₂ and O₂. Rivers included in a previous study investigating these correlations (Wit et al., 2015) are indicated by squares.

Table 4 Data measured in the four Simunjan campaigns.

Values are means of measurements. Data variability is given by standard deviation of measurements. ^* Due to technical problems, the March 2017 pH, CO₂, and O₂ data need to be treated cautiously. ${}^{**}$ In March 2017 only a minimum CO₂ concentration could be derived.

The Simunjan River shows exceptions to these correlations. Although generally CO₂ concentrations stagnate for high peat coverages, extremely high CO₂ concentrations were measured during two campaigns in the Simunjan River (Fig. 2). In January 2016 and March 2017 DOC and CO₂ concentrations in the Simunjan River were significantly higher than in March 2015 and July 2017 (Simunjan₁, Table 4). O₂ concentrations during these campaigns were lower (≈50 µmol L⁻¹) than for the other Simunjan campaigns (≈107 µmol L⁻¹), while the water pH of 5.0 was only slightly lower than during the other campaigns (pH≈5.3). The Simunjan campaigns with high DOC and CO₂ concentrations were accompanied by high concentrations of particulate carbonate (CaCO₃, Table 4), while CaCO₃ concentrations in July 2017 were much lower.

3.2 Limitation of decomposition rates by pH and O₂

To gain a better understanding of the pH and O₂ impacts on decomposition rates, least-squares optimizations of the equations in Table 1 (linear pH limitation) and Table 2 (exponential pH limitation) were performed based on measured pH, DOC, CO₂, O₂, and temperature data. The resulting decomposition parameters for the two pH approaches are listed in Table 5. A quality assessment of the least-squares optimizations can be found in Appendix D.

Table 5 Decomposition parameters derived via least-squares optimization.

Data for linear (lin.) and exponential (exp.) approaches for the pH limitation of decomposition were derived via least-squares optimization of the equations in Table 1 and Table 2, respectively. R_max is the maximum decomposition rate, b is the fraction of O₂ consumption, K_m is the Michaelis constant for O₂ limitation, and λ is the exponential pH limitation constant.

For the linear pH limitation approach, the decomposition parameters result in a Michaelis constant for O₂ limitation of $K_{\text{m}}\approx \mathrm{400}\mathrm{\mu }\mathrm{mol}{\mathrm{L}}^{-\mathrm{1}}$, a maximum decomposition rate of $R_{\text{max}}\approx \mathrm{10}\mathrm{\mu }\mathrm{mol}{\mathrm{mol}}^{-\mathrm{1}}{\mathrm{s}}^{-\mathrm{1}}$, and a fraction of O₂ consumption of b≈90 % (Table 5). Thus, the fraction by which pH limits decomposition according to this linear approach ranges from 6 % in the Batang Hari to 49 % in the Maludam, and the limitation by O₂ ranges from 71 % in the Batang Hari to 88 % in the Maludam and Siak. In total, O₂ and pH would limit decomposition in a range between 71 % and 93 %. Limitation fractions for all rivers are listed in Table A3.

For the exponential pH limitation approach, the Michaelis constant for O₂ ($K_{\text{m}}\approx \mathrm{6}\mathrm{\mu }\mathrm{mol}{\mathrm{L}}^{-\mathrm{1}}$) is significantly smaller than the constant derived for linear pH limitation (Table 5). The maximum decomposition rate ($R_{\text{max}}\approx \mathrm{4}\mathrm{\mu }\mathrm{mol}{\mathrm{mol}}^{-\mathrm{1}}{\mathrm{s}}^{-\mathrm{1}}$) and the fraction of O₂ consumption (b≈81 %) are also smaller than for linear pH limitation, but on the same order of magnitude. The exponential pH limitation factor results in λ≈0.5. According to these parameters, O₂ limits decomposition in the investigated rivers by ≤10 %, while the fraction of pH limitation ranges from 20 % in the Batang Hari to 85 % in the Maludam. The total limitation by O₂ and pH ranges from 23 % to 87 % (Table A4).

To evaluate both decomposition approaches, CO₂ and O₂ concentrations calculated based on the equations in Tables 1 and 2 with the parameters in Tables 3 and 5 were compared to measured CO₂ and O₂ concentrations in the individual rivers. For the linear pH limitation approach, correlation coefficients for CO₂ and O₂ correlations are R²=0.80 and R²=0.88, respectively (Fig. 3). For the exponential pH limitation approach, the resulting correlation coefficients are similar, whereas the correlation for CO₂ (R²=0.89) is slightly stronger and the O₂ correlation (R²=0.85) is slightly weaker (Fig. 4).

Figure 3 Correlation between measured and calculated concentrations of (a) CO₂ and (b) O₂. Grey lines indicate the 1:1 line. Calculations were performed based on the equations in Table 1 which represent linear pH limitation of decomposition rates. Each data point represents one river. Grey data points are excluded from the correlation since the data for these rivers are based on fewer than three campaigns within the same season. This includes the Simunjan campaigns with high carbon concentrations, which are excluded here due to figure scaling and further discussed in Appendix D5.

Figure 4 Correlation between measured and calculated concentrations of (a) CO₂ and (b) O₂. Grey lines indicate the 1:1 line. Calculations were performed based on the equations in Table 2 which represent exponential pH limitation of decomposition rates. Each data point represents one river. Grey data points are excluded from the correlation since the data for these rivers are based on fewer than three campaigns within the same season. This includes the Simunjan campaigns with high carbon concentrations, which are excluded here due to figure scaling and further discussed in Appendix D5.

4 Discussion

4.1 Carbon dynamics in peat-draining rivers and their dependencies on peat coverage

The linear correlations observed between peat coverage and DOC (Fig. 2a) as well as between peat coverage and pH (Fig. 2b) agree with results by Wit et al. (2015) and confirm the importance of peat soils as a major DOC source to these rivers, whereas the decomposition of DOC and leaching of organic acids lower the pH. The initial increase in CO₂ concentrations (Fig. 2c) and decrease in O₂ concentrations (Fig. 2d) with peat coverage can be explained by increased DOC decomposition due to higher DOC concentrations and also agrees with the results of Wit et al. (2015).

Previous studies included no data for rivers with peat coverage >25 % (Wit et al., 2015). In this study, we include additional campaigns at rivers with peat coverages up to 91 %. We observe that CO₂ concentrations in rivers of peat coverage >30 % level off to fairly constant values for peat coverage >50 % (Fig. 2c). This agrees with moderate CO₂ emissions that were stated for those rivers (Müller et al., 2015; Moore et al., 2013). We find that, according to Eq. (7), the stagnation can be explained by the pH and O₂ limitations of decomposition. A similar pattern of stagnating CO₂ concentrations has been observed in river sections of high DOC at the Congo River (Borges et al., 2015). The CO₂ and DOC concentrations measured in these rivers are comparable to those measured in our study, indicating that the underlying process is valid not only for Southeast Asian rivers but for tropical peat-draining rivers in general.

4.2 Decomposition in peat-draining rivers and its dependency on O₂ and pH

We were able to reproduce the stagnation in CO₂ and O₂ concentrations by introducing O₂ and pH limitations for decomposition rates in the rivers. Model approaches for both linear and exponential pH limitation factors reproduce the observed stagnation in CO₂ and O₂ concentrations and result in reasonably good correlations with the measured concentrations (Figs. 3 and 4). However, to evaluate the quality of the two approaches, the resulting parameters need to be further discussed.

The fractions of O₂ consumption by decomposition that we derived for both approaches, with $b=(\mathrm{81}\pm \mathrm{10})\mathit{\%}$ and $b=(\mathrm{90}\pm \mathrm{25})\mathit{\%}$, agree with the fraction of 0.8 that was calculated based on the oxygen-to-carbon ratio in Southeast Asian peat soils (Rixen et al., 2008).

The maximum decomposition rates of $\mathrm{4}\mathrm{\mu }\mathrm{mol}{\mathrm{mol}}^{-\mathrm{1}}{\mathrm{s}}^{-\mathrm{1}}$ for the exponential approach and $\mathrm{10}\mathrm{\mu }\mathrm{mol}{\mathrm{mol}}^{-\mathrm{1}}{\mathrm{s}}^{-\mathrm{1}}$ for the linear approach are higher than global soil phenol oxidase activity data published by Sinsabaugh et al. (2008) that stated global average soil phenol oxidase activity of 70.6 µmol h⁻¹ g⁻¹ organic matter. For a carbon content in organic matter of 38 mmol g⁻¹ (Sinsabaugh, 2010), this represents approximately $\mathrm{0.5}\mathrm{\mu }\mathrm{mol}{\mathrm{mol}}^{-\mathrm{1}}{\mathrm{s}}^{-\mathrm{1}}$, while sites of high phenol oxidase activity are listed with up to $\mathrm{3}\mathrm{\mu }\mathrm{mol}{\mathrm{mol}}^{-\mathrm{1}}{\mathrm{s}}^{-\mathrm{1}}$ (Sinsabaugh et al., 2008). Thus, the derived R_max values are slightly higher than measured decomposition rates and therewith of a realistic order of magnitude.

4.2.1 Functional dependency of decomposition on O₂

The two Michaelis constants for O₂ limitation of decomposition, derived for the linear and exponential pH limitation approaches, differ strongly (Table 5). As discussed before, the Michaelis constant represents the O₂ concentration at which O₂ availability limits decomposition by 50 %. In literature, Michaelis constants between 1 and 40 µmol L⁻¹ are suggested for the O₂ impact on phenol oxidase, depending on the phenolic species (Fenoll et al., 2002).

The linear pH limitation approach yields a Michaelis constant of $K_{\text{m}}\approx \mathrm{390}\mathrm{\mu }\mathrm{mol}{\mathrm{L}}^{-\mathrm{1}}$. This constant is higher than the O₂ concentration in atmospheric equilibrium ($\approx \mathrm{280}\mathrm{\mu }\mathrm{mol}{\mathrm{L}}^{-\mathrm{1}}$), which implies an oxygen deficit at atmospheric conditions that does not exist (Vaquer-Sunyer and Duarte, 2008). However, though the derived K_m value for this linear pH limitation is unrealistically high, this does not necessarily negate the linear pH approach. High parameter interdependence between K_m and R_max complicate the computation of these decomposition parameters (Sect. D1). To disentangle the impact of the intercorrelated parameters, additional least-squares optimizations at fixed K_m values ranging from 1 to 40 µmol L⁻¹ (Fenoll et al., 2002) were performed (Sect. D2). These optimizations result in maximum decomposition rates of R_max = (1.4–2.4) $\mathrm{\mu }\mathrm{mol}{\mathrm{mol}}^{-\mathrm{1}}{\mathrm{s}}^{-\mathrm{1}}$ and O₂ consumption factors of $b=(\mathrm{102}-\mathrm{109})\mathit{\%}$ and therewith do not agree with literature values of these parameters ($R_{\text{max}}\mathit{⩾}\mathrm{3}\mathrm{\mu }\mathrm{mol}{\mathrm{mol}}^{-\mathrm{1}}{\mathrm{s}}^{-\mathrm{1}}$ and b≈80 %; Sinsabaugh et al., 2008; Rixen et al., 2008).

The exponential pH limitation approach yields a Michaelis constant of $K_{\text{m}}\approx \mathrm{6}\mathrm{\mu }\mathrm{mol}{\mathrm{L}}^{-\mathrm{1}}$. This value is in good agreement with the literature data of 1 to 40 µmol L⁻¹ (Fenoll et al., 2002). Its large uncertainty (>400 %, Table 5) is mainly caused by relatively high concentrations of O₂ in the rivers. Due to exchange with atmospheric O₂, the concentrations in all rivers exceed the median O₂ threshold to lethal hypoxic conditions of 50 µmol L⁻¹ (Vaquer-Sunyer and Duarte, 2008). Thus, the O₂ limitation in peat-draining rivers is relatively small (between 3 % and 10 %, Table A4). Consequentially a majority of the decomposition limitation is caused by the low pH in peat-draining rivers that we found to limit the decomposition rates in rivers of high peat coverage (low pH) by up to 85 % (Table A4).

4.2.2 Functional dependency of decomposition on pH

Our results indicate the exponential pH limitation of decomposition to be more realistic than the linear pH limitation. The exponential limitation better represents river CO₂, especially for high CO₂ concentrations which are most strongly affected by the pH limitation. The exponential limitation is additionally supported by the unrealistically high O₂ limitation resulting from the linear pH approach. The strong collinearity between decomposition parameters in the linear pH limitation approach complicates the interpretation of the parameters mentioned above. Additional calculations of the parameters R_max and b for fixed K_m disagree with literature data and thus further disprove the linear approach (Sect. D2).

The exponential pH coefficient results in $\mathit{\lambda }=\mathrm{0.5}\pm \mathrm{0.1}$. Thus, in terms of H⁺ activity the correlation is given by ${\left\{{\text{H}}^{+}\right\}}^{\frac{\mathrm{0.5}}{\ln \left(\mathrm{10}\right)}}$, which roughly equals the fifth root of H⁺ activity. The derived limitation coefficient is similar to coefficients reported for high-latitude peat soils (λ=0.65 and λ=0.77) that were determined via laboratory measurements of phenol oxidase activity (Williams et al., 2000). The fact that the exponential inhibition by pH can be found in those high-latitude peat soils, as well as in tropical peat-draining rivers, suggests that the investigated correlations and processes are also relevant in other regions and that soil pH and water pH are important regulators of global carbon emissions.

4.3 Impact of additional processes

Our results neglect the direct leaching of CO₂ from soils as well as the photo-mineralization of DOC and the consumption of CO₂ by autotrophic production within rivers. CO₂ leaching rates are likely higher for peat soils than for mineral soils (Kang et al., 2018; Abril and Borges, 2019), and autotrophic production is limited in peat-draining rivers (Wit et al., 2015). Thus, both of these processes would work against the observed stagnation in CO₂ concentrations, and the exclusion of these processes could cause underestimation of the limitation factors rather than overestimation.

A recent study by Nichols and Martin (2021) found low phenol oxidase activity in Southeast Asian peat-draining rivers and low degradation of DOC from those rivers in an additional incubation experiment. They concluded that the remineralization of peat-derived DOC in Southeast Asian aquatic systems is likely dependent on photodegradation rather than microbial respiration (Nichols and Martin, 2021). This is supported by photolability of DOC from those regions (Martin et al., 2018). However, photomineralization rates would not be impacted by river pH or O₂. Thus, with photomineralization as the main cause of DOC degradation, no stagnation in CO₂ is expected. Accordingly, photomineralization of DOC, like the before-mentioned processes, would work against the observed CO₂ stagnation and could cause underestimation of the limitation parameters.

4.4 Disruption of the pH limitation by carbonates

Typically, concentrations of particulate carbonate in peat-draining rivers are low (Wit et al., 2018). However, we observed high CaCO₃ concentrations for two of the four campaigns in the Simunjan River. These two campaigns also show high DOC and CO₂ concentrations (Table 4). Possible causes for high carbonate concentrations during these campaigns could be increased erosion of mineral soils due to deforestation in mountain regions upstream or liming practices in plantations along the river. In either case, high carbonate concentrations at such a low pH indicate high dissolution of carbonates, which might have counteracted a more prominent decrease in pH due to decomposition of DOC. At the same time, the low river pH causes transformation of dissolved carbonates to CO₂ and thus additionally increases CO₂ concentrations. These processes seem to have suspended the natural pH limitation of CO₂ production in peat-draining rivers. This explains the high CO₂ concentrations observed during those two Simunjan campaigns (Table 4). Calculation based on the derived decomposition dependencies would indicate even higher CO₂ concentrations than measured. This indicates that the river carbon parameters were not in thermodynamic equilibrium during these campaigns as is further discussed in Sect. D5.

5 Conclusions

Our study shows that CO₂ concentrations in and emissions from Southeast Asian rivers stagnate for high peat coverages of the river catchments. Despite further increase in river DOC concentrations, CO₂ concentrations are fairly constant for peat coverages >50 %. We find that this stagnation is caused by a natural limitation of DOC decomposition in these rivers. This process provides an explanation of moderate CO₂ emissions measured from rivers of high carbon content

Correlation to measured data indicates that the limitation in decomposition is mainly caused by low river pH. Data reveal an exponential limitation of DOC decomposition by pH as the most realistic scenario. This reduces the CO₂ production in rivers of high peat coverage by up to 85 %. The limiting impact of O₂ on decomposition in the rivers is comparatively small with <10 %. This derived limitation of decomposition should be included to improve model studies and accurately capture river CO₂ emissions from tropical peat areas.

Campaigns with high carbon loads in the Simunjan River indicate that the natural CO₂ limitation can be suspended by high input of DOC and carbonates. Data from campaigns with enhanced concentrations of DOC and suspended carbonates reveal CO₂ emissions that were increased by almost 100 %. Here, the high DOC concentrations enhance decomposition, and the input of carbonates counteracts the pH decrease associated with large inputs of CO₂. Possible sources for enhanced carbonate concentrations can be soil erosion upstream of coastal peatlands, or liming practices in plantations along the rivers, which are common practice to improve plant growth on acidic soils. This carbonate impact should be considered for anthropogenic activities like liming and enhanced weathering.

Our study is based on measurements in Southeast Asian peat-draining rivers. However, comparison to data from African rivers and laboratory studies of decomposition in temperate peat soils suggests that the investigated correlations and processes are also relevant in other regions and that soil pH and water pH are important regulators of global carbon emissions.

Appendix A

Figure A1 Correlation of river depth with (a) atmospheric exchange coefficients (k₆₀₀) and (b) catchment peat coverage. A linear correlation between river depth and exchange coefficient reveals a slope of $k_{\mathrm{600}}/d=(\mathrm{2.5}\pm \mathrm{0.2})\mathrm{cm}{\mathrm{h}}^{-\mathrm{1}}{\mathrm{m}}^{-\mathrm{1}}=(\mathrm{7.0}\pm \mathrm{0.5})\times {\mathrm{10}}^{-\mathrm{6}}{\mathrm{s}}^{-\mathrm{1}}$.

Figure A2 Functional dependencies of CO₂ on (a) DOC, (b) pH, (c) temperature (T), and (d) O₂ according to the equation in Table 1 (linear pH approach) and in Table 2 (exponential pH approach).

Figure A3 Functional dependencies of O₂ on (a) DOC, (b) pH, and (c) temperature (T) according to the equation in Table 1 (linear pH approach) and in Table 2 (exponential pH approach).

Figure A4 Sensitivity of calculated CO₂ to (a) the maximum decomposition rate (R_max) and (b) the Michaelis constant for O₂ concentration (K_m) according to the equation in Table 1 for the linear pH approach and the parameter ranges given in Table 5.

Figure A5 Sensitivity of calculated CO₂ to (a) the maximum decomposition rate (R_max), (b) the Michaelis constant for O₂ concentration (K_m), and (c) the exponential pH limitation constant (λ) according to the equation in Table 2 for the exponential pH approach and the parameter ranges given in Table 5.

Figure A6 Sensitivity of calculated O₂ to (a) the maximum decomposition rate (R_max), (b) the Michaelis constant for O₂ concentration (K_m), and (c) the fraction of O₂ consumption by decomposition (b) according to the equation in Table 1 for the linear pH approach and the parameter ranges given in Table 5.

Figure A7 Sensitivity of calculated O₂ to (a) the maximum decomposition rate (R_max), (b) the Michaelis constant for O₂ concentration (K_m), (c) the exponential pH limitation constant (λ), and (d) the fraction of O₂ consumption by decomposition according to the equation in Table 2 for the exponential pH approach and the parameter ranges given in Table 5.

Table A1 List of river campaigns.

Table A2 Temperature-dependent exchange coefficients k and Henry coefficients K of CO₂ and O₂ for the individual rivers.

Exchange coefficients were derived from measured k₆₀₀ (Table 3) and water temperature (T) according to $k_{\text{X}}=k_{\mathrm{600}}\cdot (S{c}_{\text{X}}/\mathrm{600}{)}^{-n}$ with Schmidt numbers $S{c}_{{\text{CO}}_{\mathrm{2}}}$ and $S{c}_{{\text{O}}_{\mathrm{2}}}$ derived according to the equations in Wanninkhof (1992). An exponent of $n=\mathrm{1}/\mathrm{2}$ (valid for rough surfaces, Zappa et al., 2007) was used for the rivers. Henry coefficients were derived based on water temperature from the equations stated in Weiss (1974) for CO₂ and the equations stated in Weiss (1970) for O₂.

Table A3 pH and O₂ limitations in the individual rivers based on linear pH approach.

Fraction by which the decomposition is lowered due to the impact of pH and O₂, calculated based on the limitation factors in Eq. (7) and the parameters in Table 5 according to pH $\text{lim.}=(\mathrm{1}-L_{\text{pH}})$, O${}_{\mathrm{2}}\text{ lim.}=(\mathrm{1}-L_{{\text{O}}_{\mathrm{2}}})$, and total $\text{lim.}=(\mathrm{1}-L_{\text{pH}}\cdot L_{{\text{O}}_{\mathrm{2}}})$.

Table A4 pH and O₂ limitations in the individual rivers based on exponential pH approach.

Fraction by which the decomposition is lowered due to the impact of pH and O₂, calculated based on the limitation factors in Eq. (7) and the parameters in Table 5 according to pH $\text{lim.}=(\mathrm{1}-L_{\text{pH}})$, O${}_{\mathrm{2}}\text{ lim.}=(\mathrm{1}-L_{{\text{O}}_{\mathrm{2}}})$, and total $\text{lim.}=(\mathrm{1}-L_{\text{pH}}\cdot L_{{\text{O}}_{\mathrm{2}}})$.

Appendix B: Impact of data limitation on study results

B1 Impact of sampling location

The data for this study were collected from samples taken in river sections that flow through peat soil. This ensures that the impact of peat soils on the river parameters is captured.

Concentrations measured in the small Malaysian rivers (Maludam and Sebuyau and Simunjan), with the exception of the Simunjan campaigns in January 2016 and March 2017 (Table 4, Fig. B3), show little variation over the river path and between campaigns (Figs. B1, B2 B3). However, the larger rivers drain mineral soils for the majority of their path and only reach peat regions close to the coast. Those rivers exhibit stronger differences in carbon concentrations along the length of the river. Rixen et al. (2010) found that DOC concentrations in the Siak River are higher by a factor of up to 4 in coastal peat regions than in the upstream river. CO₂ concentrations in the large Sumatran rivers were not measured outside of the coastal peat regions. Due to the lower pH in river parts that cut through peat and the related pH limitation of DOC decomposition, the difference in CO₂ concentrations along the river is likely lower than the difference in DOC concentrations. This is also indicated by CO₂ measurements in the Rajang River that reveal CO₂ concentrations in the peat-draining river sections to be only 15 %–20 % higher than CO₂ concentrations upstream of the peat regions (Müller-Dum et al., 2019).

B2 Impact of seasonality

The Southeast Asian study area is impacted by the Malaysian–Australian monsoon that causes the presence of moisture-loaded air with high precipitation rates from October to April while dry air dominates from May to September. To catch the impact of these rain and dry seasons on river carbon dynamics, campaigns in different months of the year were performed (Table A1). However, the seasonal data coverage is not dense enough to clearly identify or disprove a seasonal pattern in the measured data (Figs. B4, B5, and B6).

Figure B1 Individual DOC, CO₂, O₂, and pH measurements in the Maludam River versus the distance from the river mouth. Different colours stand for the individual river campaigns, and the grey shaded area indicates regions of S>1 that were excluded from the data in this study.

Figure B2 Individual DOC, CO₂, O₂, and pH measurements in the Sebuyau River versus the distance from the river mouth. Different colours stand for the individual river campaigns, and the grey shaded area indicates regions of S>1 that were excluded from the data in this study.

Figure B3 Individual DOC, CO₂, O₂, and pH measurements in the Simunjan River versus the distance from the river mouth. Different colours stand for the individual river campaigns, and the grey shaded area indicates regions of S>1 that were excluded from the data in this study.

Figure B4 Average CO₂ concentrations for individual campaigns compared to monthly temperature (Natawa, 2018a; Yasutomi et al., 2011) and precipitation (Natawa, 2018b; Yatagai et al., 2020) data (2005–2015 average) at the location of the respective river. Each panel represents one river.

Figure B5 Average DOC concentrations for individual campaigns compared to monthly temperature (Natawa, 2018a; Yasutomi et al., 2011) and precipitation (Natawa, 2018b; Yatagai et al., 2020) data (2005–2015 average) at the location of the respective river. Each panel represents one river.

Figure B6 Average O₂ concentrations for individual campaigns compared to monthly temperature (Natawa, 2018a; Yasutomi et al., 2011) and precipitation (Natawa, 2018b; Yatagai et al., 2020) data (2005–2015 average) at the location of the respective river. Each panel represents one river.

Appendix C: Comparison of different peat coverage estimates

Different peat maps are available for Southeast Asia, and the approaches to determine peat coverage of river catchments were inconsistent among different studies cited in our paper. We want to show here that the choice of a data product is crucial for the determination of peat coverage. We compare three different products (Table C1): the FAO Digital Soil Map of the World, country products downloaded at Global Forest Watch, and the Center for International Forestry Research (CIFOR) Wetlands distribution.

Those three products lead to highly different results (Table C2). We observed a tendency that CIFOR leads to smaller peat coverage than FAO and GFW. This is because CIFOR misses some, but not all peat areas that are known to be under industrial plantations. Gumbricht et al. (2017a) already pointed out that their model underestimates peatland area in Sumatra because peats are largely drained, which the model does not capture. However, in the Musi and Batang Hari catchment, CIFOR sees larger peat areas than FAO and GFW, which means that some peatlands might be missing in those maps.

We decided to use the GFW maps for several reasons: (1) CIFOR seems to miss peat under industrial plantations, which is still relevant for river carbon dynamics. Therefore, we chose not to use the CIFOR maps. (2) Between GFW and FAO, GFW is more recent than FAO for Indonesia. For Sarawak (Malaysia), both are based on the 1968 soil map by the Land Survey Department, but FAO uses a 10-fold coarser scale than the 1968 soil map (1:5 000 000 compared to 1:5 000 000). Thus, the GFW product was used. (3) GFW maps are based on official information, and we believe that the local authorities would know best about the peatland distribution in their country.

Similar to the peat coverage, the publications from which we use data in our study all had different approaches to determining catchment size – either including (Müller-Dum et al., 2019) or excluding (Wit et al., 2015) smaller sub-catchments. In our study, we aimed to unify those different approaches. Therefore, we recalculated catchment areas from one single data product (HydroSHEDS, Lehner et al., 2008), including sub-catchments that were identified using HydroSHEDS flow directions. The Simunjan catchment is included in the bigger Sadong catchment in HydroSHEDS. Therefore, it was manually delineated using HydroSHEDS flow directions.

(FAO-UN, 2007)(MoA, 2012)(Wetlands International, 2004)Gumbricht et al. (2017b)Gumbricht et al. (2017a)

Table C1 Different data products used to assess peatland extent in the catchments.

Table C2 Results for peat coverage (PC) in the different catchments using the three different data products.

Table C3 Correlations between the derived parameters.

Positive ($\mathit{\rho }\ge \mathrm{0.5}:+$), negative ($\mathit{\rho }\le -\mathrm{0.5}:-$), and non-significant ($-\mathrm{0.5}\le \mathit{\rho }\le \mathrm{0.5}:{}^{\circ }$) correlations between the parameters are indicated in the bottom left. The top right bold numbers represent the numerical Pearson correlation coefficients (ρ) between the parameters. The correlations are derived from the least-squares optimization of the linear pH approach (left table) and the exponential pH approach (right table). R_max is the maximum decomposition rate, b is the fraction of O₂ consumption by decomposition, K_m is the Michaelis constant for O₂ concentrations, and λ is the exponential pH limitation constant.

Appendix D: Quality assessment of least-squares optimizations

Uncertainty sources in the least-squares optimizations are interdependencies between the fitted parameters and noise in the measured data. We try to minimize the impact of measurement noise by including relative uncertainties (σ) of measured CO₂ and O₂ in the least-squares optimization. Thus, data from rivers with higher variation in measured parameters are constrained less rigidly in the optimization. The parameter interdependence is a more important source of uncertainties for our optimization as they cause interdependencies between the fitted parameters as well. This is especially relevant for the linear approach, where the functional dependencies of CO₂ and O₂ on the different river parameters are more similar than for the exponential approach (Fig. A2).

Figure D1 Correlation between measured and calculated concentrations of (a) CO₂ and (b) O₂. Grey lines indicate the 1:1 line. Calculations were performed for exponential pH limitation without O₂ limitation. Each data point represents one river. Grey data points are excluded from the correlation since the data for these rivers are based on fewer than three campaigns within the same season.

Figure D2 Correlation between measured and calculated concentrations of (a) CO₂ and (b) O₂. Grey lines indicate the 1:1 line. Calculations were performed based on the equations in Table 1, which represent linear pH limitation of decomposition rates. Each data point represents one river. Grey data points are excluded from the correlation since the data for these rivers are based on fewer than three campaigns within the same season. The Simunjan campaigns with high carbon concentrations (Simunjan₂, Table 4) are indicated by red circles.

Figure D3 Correlation between measured and calculated concentrations of (a) CO₂ and (b) O₂. Grey lines indicate the 1:1 line. Calculations were performed based on the equations in Table 2, which represent exponential pH limitation of decomposition rates. Each data point represents one river. Grey data points are excluded from the correlation since the data for these rivers are based on fewer than three campaigns within the same season. The Simunjan campaigns with high carbon concentrations are excluded from these figures and further discussed in the Appendix D5. The Simunjan campaigns with high carbon concentrations (Simunjan₂, Table 4) are indicated by red circles.

D1 Parameter collinearity

The functional CO₂ dependencies on pH, O₂, and DOC are more similar to each other for the linear than for the exponential pH approach (Fig. A2). This is also reflected in higher parameter uncertainties derived from the linear pH approach (Table 5). However, investigation of the correlation coefficients between the individual parameters reveals a strong positive correlation between the maximum decomposition rate (R_max) and the Michaelis constant for O₂ (K_m) in both the linear and the exponential pH approaches (Table C3). Additionally, there is a significant negative correlation between the exponential pH limitation constant (λ) and K_m (Table C3).

For the linear pH approach, the extremely high correlation between R_max and K_m (ρ=0.99) makes it impossible to meaningfully disentangle the individual impacts of these parameters. To test the possibility of a linear pH limitation in decomposition, least-squares correlations with fixed K_m parameters within literature values (1–40 µmol L⁻¹, Fenoll et al., 2002) were performed (Appendix D2).

For the exponential approach, while the parameters show a strong correlation (ρ=0.82 for R_max and K_m and $\mathit{\rho }=-\mathrm{0.86}$ for K_m and λ; Table C3), the functional dependencies are distinct enough to disentangle the parameter's impacts comparatively well, and the comparison to literature values supports the exponential pH limitation. The high uncertainty in the K_m parameter for this approach is only of small relevance as the O₂ limitation is comparatively weak. In fact, the pH limitation alone is able to reproduce the measured parameters quite well (Appendix D3).

D2 Least-squares optimizations of linear pH approach with fixed K_m parameters

To test the possibility of a linear pH limitation in decomposition, least-squares correlations for fixed O₂ Michaelis constants within literature values of 1–40 µmol L⁻¹ (Fenoll et al., 2002) were performed (Table D1). The good agreement for all K_m values is caused by the strong collinearity between K_m and R_max that enables a change in R_max to compensate for changes in K_m. R_max values for the fixed K_m values range between 1.4 and $\mathrm{2.4}\mathrm{\mu }\mathrm{mol}{\mathrm{mol}}^{-\mathrm{1}}{\mathrm{s}}^{-\mathrm{1}}$. These maximum decomposition rates are lower than high decomposition rates derived based on global phenol oxidase activity by Sinsabaugh et al. (2008). The fraction of O₂ consumption ranges between 102 % and 109 %, indicating that for this approach to be true, more O₂ would need to be consumed than CO₂ is produced. In reality, additional O₂ can be taken from the organic matter, reducing the needed O₂ compared to CO₂ production. Thus, despite the good correlation to measured data (Table D1) the derived parameters for the linear approach do not agree well with literature data, which makes this approach unlikely.

D3 Least-squares optimizations of exponential pH approach without O₂ limitation

Our study revealed that low pH is the main decomposition-impelling parameter in peat-draining rivers. To study whether this parameter alone can explain the observed stagnation on CO₂ and O₂ for high peat coverage (Fig. 2), a least-squares optimization without O₂ limitation was performed. This optimization yields decomposition parameters that differ only insignificantly from the parameters derived for exponential pH limitation with additional limitation by O₂ (Table D2). The correlation of measured CO₂ and O₂ to concentrations derived based on these parameters and the equations in Table 2 reveals a good agreement (Fig. D1). Only for the Kampar River does neglect of the O₂ limitation yield negative river O₂ concentrations (Fig. D1). This indicates that for O₂ concentrations in the examined rivers (O${}_{\mathrm{2}}>\mathrm{50}\mathrm{\mu }\mathrm{mol}{\mathrm{L}}^{-\mathrm{1}}$), the pH limitation alone is sufficient to explain the majority of the observed stagnation.

D4 Validation of optimal pH for phenol oxidase activity

To validate the optimal pH for decomposition (pH₀) in our study area, a least-squares optimization of the exponential pH approach (Table 2) including the parameter pH₀ was performed. The resulting pH value of pH₀≈7.2 agrees well with the literature value of 7.5 used in our study (Table D3). However, it reveals a high collinearity to R_max that causes high parameter uncertainties.

D5 Decomposition approach for abnormal Simunjan campaigns

In the correlation figures Figs. 3 and 4, the Simunjan campaigns of January 2016 and March 2017 (Table 4) were excluded due to scaling of the figures. Here we show the correlation figures with inclusion of those campaigns (Figs. D2 and D3). Since the data are based on only one campaign, they were excluded from the least-squares optimization. Calculated CO₂ concentrations based on both limitation approaches result in significantly higher concentrations than measured during the campaign (Figs. D2 and D3). At the same time, calculated O₂ concentrations are lower than measured concentrations in the rivers.

This indicates that the parameters in these campaigns are not in equilibrium based on the processes of atmospheric gas exchange and decomposition. This could be caused by additional processes of CO₂ sources and sinks during these anomalous campaigns. However, since the observed events are temporal, we consider it likely that the river parameters simply had not reached a state of equilibrium yet. With such high carbon yields it is also possible that the river cannot reach a state of equilibrium before the water discharges into the ocean. However, as mentioned before, the data are mainly based on one campaign. To validate this assumption, further studies would be needed.

Table D1 Decomposition parameters for fixed K_m values.

Decomposition parameters derived for the linear pH limitation approach via least-squares optimization of the equations in Table 1 with fixed Michaelis constant for O₂ concentrations given in µmol L⁻¹. R_max is the maximum decomposition rate stated in $\mathrm{\mu }\mathrm{mol}{\mathrm{mol}}^{-\mathrm{1}}{\mathrm{s}}^{-\mathrm{1}}$, and b is the fraction of O₂ consumption by decomposition stated in percent. Literature values for R_max and b are taken from Sinsabaugh et al. (2008) and Rixen et al. (2008), respectively. Additionally, coefficients of determination (R²) and root-mean-square errors (RMSEs) of the correlation between measured CO₂ and CO₂ derived using these parameters and the equations in Table 1 are listed to indicate the quality of the fit.

Table D2 Decomposition parameters for exponential pH limitation without O₂ limitation.

Decomposition parameters derived for exponential pH limitation with and without additional O₂ limitation. Parameters were derived via least-squares optimization of the equations in Table 2 to measured data. R_max is the maximum decomposition rate, b is the fraction of O₂ consumption by decomposition, K_m is the Michaelis constant for O₂ limitation, and λ is the exponential pH limitation constant.

Table D3 Decomposition parameters for exponential pH limitation without O₂ limitation.

Decomposition parameters derived for exponential pH limitation with freely set pH₀ and with a fixed pH₀ of 7.5. Parameters were derived via least-squares optimization of the equations in Table 2 to measured data. R_max is the maximum decomposition rate, b is the fraction of O₂ consumption by decomposition, K_m is the Michaelis constant for O₂ limitation, and λ is the exponential pH limitation constant.