Impact of peatlands on carbon dioxide (CO₂) emissions from the Rajang River and Estuary, Malaysia

Tropical peat-draining rivers are known as potentially large sources of carbon dioxide (CO₂) to the atmosphere due to the high loads of carbon they receive from surrounding soils. However, not many seasonally resolved data are available, limiting our understanding of these systems. We report the first measurements of carbon dioxide partial pressure (pCO₂) in the Rajang River and Estuary, the longest river in Malaysia. The Rajang River catchment is characterized by extensive peat deposits found in the delta region, and by human impact such as logging, land use and river damming. pCO₂ averaged 2540±189 µatm during the wet season and 2350±301 µatm during the dry season. Using three different parameterizations for the gas transfer velocity, calculated CO₂ fluxes to the atmosphere were 1.5 (0.5–2.0) g C m⁻² d⁻¹ (mean, minimum – maximum) during the wet season and 1.7 (0.6–2.6) g C m⁻² d⁻¹ during the dry season. This is at the low end of reported values for Southeast Asian peat-draining rivers, but similar to values reported for Southeast Asian rivers that do not flow through peat deposits. In the Rajang River, peatlands probably do not contribute much to the CO₂ flux due to the proximity of the peatlands to the coast, which limits the opportunity for degradation of organic C during transport. Thus, we suggest that peat coverage is, by itself, insufficient as the sole predictor of CO₂ emissions from peat-draining rivers, and that other factors, like the spatial distribution of peat in the catchment and pH, also need to be considered.

1 Introduction

Tropical rivers transport large amounts of terrestrially derived carbon to the ocean (Dai et al., 2012) and the atmosphere (Aufdenkampe et al., 2011; Raymond et al., 2013). It has been estimated that 78 % of riverine carbon dioxide (CO₂) emissions occur in the tropics (Lauerwald et al., 2015). Tropical wetlands exert a particularly strong influence on the carbon budget of these rivers. Two regional studies independently showed that the partial pressure of CO₂ (pCO₂) in rivers increases with increasing wetland coverage in the catchment. Borges et al. (2015) established a relationship between wetland extent and pCO₂ for African rivers. Wit et al. (2015) presented an analog synthesis for Southeast Asian rivers, which flow through peatlands. Peatlands are a special type of wetland, where organic matter accumulates at rates that make them the most effective terrestrial carbon store on a millennial timescale (Dommain et al., 2011). Southeast Asian peatlands store 68.5 Gt carbon (Page et al., 2011). The highest riverine dissolved organic carbon (DOC) concentrations reported so far were found in Southeast Asian peat-draining rivers (Alkhatib et al., 2007; Moore et al., 2011; Müller et al., 2015), with an annual average of 68 mg L⁻¹ DOC found in an undisturbed peat-draining river (Moore et al., 2013). Because of these high DOC concentrations, Indonesian rivers may account for 75 % of the DOC flux into the South China Sea (SCS) while accounting for 39 % of the discharge (Huang et al., 2017). Surprisingly, CO₂ emissions from these rivers are not exceptionally high (Müller et al., 2015; Wit et al., 2015). This is attributed to a short residence time of the organic matter in the river, allowing little time for decomposition, and the resistance of peat-derived carbon to bacterial degradation. Nevertheless, the CO₂ flux from peat-draining rivers to the atmosphere increases with increasing peat coverage in the river basin (Wit et al., 2015), showing that these ecosystems exert an important influence on a river's carbon budget.

Most Southeast Asian peat-draining rivers are disturbed by human activities such as river damming, urbanization, deforestation (Milliman and Farnsworth, 2011) and discharge of untreated wastewater (Park et al., 2018). Anthropogenic change poses a new challenge to understanding carbon fluxes in Asian river systems, and more data are urgently needed to constrain the carbon budget for this important region (Park et al., 2018). In Malaysia, the country holding the second largest share of tropical peat (Page et al., 2011), river CO₂ emissions have only been studied in a small undisturbed peat-draining river (Müller et al., 2015), in estuaries (Chen et al., 2013; Müller et al., 2016) and in two river reaches which were not influenced by peat (Müller et al., 2016). In this study, the longest Malaysian river, the Rajang River on the island of Borneo, was investigated. This river flows through largely logged-over tropical rainforest (Gaveau et al., 2014), urban areas and disturbed peat swamps (Gaveau et al., 2016). The aim of this study was to assess the Rajang River and Estuary carbon load and to investigate the impact of peatlands on its CO₂ emissions. To this end, we surveyed longitudinal transects extending from river reaches that were not influenced by peat to the peat-covered delta.

2 Materials and methods

2.1 Study area

The Rajang River is located in the Malaysian state of Sarawak in the northern part of the island of Borneo (Fig. 1a). Sarawak has a tropical climate with high temperatures (average 26.6 ^∘C, 1992–2016, in Sibu, DWD, 2018) and high precipitation (average 3578 mm yr⁻¹, 1992–2016 in Sibu, DWD, 2018). The region experiences two monsoonal periods: the northeastern monsoon, with enhanced rainfall and frequent floods, occurs between December and February (“wet season”; see Fig. 2a), while the southwestern monsoon from May until September is associated with relatively drier weather (“dry season”). However, despite the monsoon seasons, rainfall is high throughout the year (Sa'adi et al., 2017).

The Rajang River originates in the Iran mountains, a mountain range at the border between Malaysia and Indonesia (MacKinnon, 1996) with elevations of up to 1800 m (Milliman and Farnsworth, 2011). It drains an area of approximately 52 010 km² (Lehner et al., 2006; DID, 2017) whose geology is dominated by Cenozoic sedimentary and metamorphic rocks, consisting of siliciclastic rock with minor amounts of carbonates (Staub et al., 2000; Milliman and Farnsworth, 2011). The Rajang River flows approximately 530 km from east to west and discharges into the South China Sea (Milliman and Farnsworth, 2011). The main settlements along the river are the towns of Kapit, Kanowit and the city of Sibu (163 000 inhabitants) (see Fig. 1b). In addition, a large number of longhouses (traditional buildings inhabited by local communities) are located along the river and its tributaries (Ling et al., 2017). Hydroelectric power plants were built on two tributaries in the upper Rajang basin: The Bakun hydroelectric power plant commenced operation in 2011 and the Murum dam in 2015 (Sarawak Energy, 2013; see Fig. 1b). The construction of another hydroelectric power plant on a tributary in the southern Rajang basin is planned for the future (Sarawak Energy, 2013).

The Rajang delta system is comprehensively described in Staub and Gastaldo (2003). It is entirely surrounded by peatlands (Fig. 1b), which extend over an area that corresponds to approximately 11 % of the catchment size (Nachtergaele et al., 2009). Most of these peatlands have been converted to industrial oil palm plantations (Gaveau et al., 2016, Fig. 1b). The main distributary channels forming the delta (from north to south) are the Igan, Hulu Seredeng (which splits up into Lassa and Paloh), Belawai and Rajang, which have a maximum tidal range (spring tide) of 3–6 m (Staub and Gastaldo, 2003). Saltwater intrudes into the estuary approximately as far as the point where the Rajang River splits up into its four southernmost distributaries, a few kilometers downstream of Sibu (Fig. 1b), depending on the season. Tidal influence extends further inland approximately up to the town of Kanowit (Staub and Gastaldo, 2003).

Figure 1 Map of the Rajang basin on the island of Borneo (a) and a close-up of the basin (b) including the location of the peatlands (Nachtergaele et al., 2009), industrial oil palm plantations (Gaveau et al., 2016) and the stations during the cruise. The distributaries are marked with numbers: 1: Igan, 2: Lassa, 3: Paloh, 4: Belawai, 5: Rajang.

Figure 2 Monthly discharge calculated for the Rajang River (a) and average temperatures (b) for the city of Sibu (2^∘20^′ N, 111^∘50^′ E) for the years 2015 and 2016 and the long-term average from 1992 to 2016. Data from DWD (2018).

Monthly discharge (Fig. 2a) was estimated from monthly precipitation (1992–2016; DWD, 2018) and an evapotranspiration rate of 1545 mm yr⁻¹ (Kumagai et al., 2005) or 43.2 %. Annual average discharge from 1992 to 2016 was 3355 m³ s⁻¹, in good agreement with reported discharges of 3490 m³ s⁻¹ (Milliman and Farnsworth, 2011) and 3372 m³ s⁻¹ for the years 1991–2015 (Sa'adi et al., 2017).

2.2 Surveys

We sampled the Rajang River during two surveys, which were designed to get spatial coverage of both peat and nonpeat areas during the wettest and driest period of one year. The first survey took place at the peak of the monsoon season in January 2016 (“wet season”). The second one was performed during the dry season in August 2016 (“dry season”). In January 2016, we entered the Rajang River through the Rajang river mouth (distributary 5 in Fig. 1b), went upstream to the town of Kapit and back downstream to the town of Belawai at the Belawai river mouth (distributary 4 in Fig. 1b). In August 2016, we entered the Rajang River through the Rajang river mouth (5), went upstream to Kapit and back to Sibu. From there, we went out to the coast through the Lassa distributary (2), and back to Sibu through the Igan distributary (1). The last sampling stretch was from Sibu into the Paloh distributary (3) and back to Belawai (4). During this campaign, one stationary measurement was performed overnight in Sarikei in the Rajang distributary in order to assess tidal and diurnal variability.

2.3 CO₂ measurements

The setup on the boat was similar to the one described in Müller et al. (2016). Surface water was pumped through a shower-type equilibrator (Johnson, 1999) at a rate of approximately 15 L min⁻¹. In the beginning, the equilibrator headspace was connected to an FTIR analyzer (Griffith et al., 2012), which allows for the simultaneous measurement of CO₂, methane (CH₄), nitrous oxide (N₂O) and carbon monoxide (CO). During the cruise in January 2016, a failure of the FTIR analyzer occurred and measurements were continued (also in August 2016) using an Li-820 non-dispersive infrared (NDIR) analyzer for the measurement of CO₂ (Licor, USA). For calibration and inter-calibration of the two instruments, a set of gravimetrically prepared gas mixtures (DEUSTE Steininger) was measured, which were calibrated against the World Meteorological Organization (WMO) standard scale by the Max Planck Institute for Biogeochemistry in Jena, Germany. For the FTIR, spectra were averaged over 5 min and dry air mole fractions were retrieved using the software MALT5 (Griffith, 1996). Li-820 data were stored with a temporal resolution of 1 min. Gas partial pressure was determined using measurements of ambient pressure with a PTB110 barometer (Vaisala, Finland) and correction for removal of water according to Dickson et al. (2007). Water temperature was measured in the equilibrator and in the surface water and correction to water surface temperature was performed according to Dickson et al. (2007).

In August 2016, the internal pressure sensor of the Li-820 failed. Because the instrument performs an internal correction based on the cell pressure, this correction had to be reversed and recalculated with an assumed internal cell pressure. This procedure is described in the Supplement.

CO₂ fluxes (FCO₂, in g C m⁻² d⁻¹) across the water–air interface were computed using the gas transfer equation

$$\begin{array}{}\text{(1)}& {\text{FCO}}_{\mathrm{2}}=k{K}_{\mathrm{0}}(p{\text{CO}}_{\mathrm{2}}^{\text{water}}-p{\text{CO}}_{\mathrm{2}}^{\text{air}})\cdot f_{\mathrm{1}}{f}_{\mathrm{2}},\end{array}$$

where k is the gas transfer velocity (m s⁻¹), K₀ is the solubility (mol L⁻¹ atm⁻¹) calculated according to Weiss (1974), pCO${}_{\mathrm{2}}^{\text{water}}$ is the partial pressure of CO₂ in water, pCO${}_{\mathrm{2}}^{\text{air}}$ is the partial pressure of CO₂ in the overlying air (both in µatm), f₁ is a conversion factor from L⁻¹ to m⁻³, and f₂ is a conversion factor from µmol s⁻¹ to mg d⁻¹. The atmospheric mole fractions of CO₂ during the months of our measurements were derived from the NOAA ESRL Carbon Cycle Cooperative Global Air Sampling Network (Dlugokencky et al., 2018) for the closest station, which was Bukit Kototabang, Indonesia.

The gas transfer velocity (k) is a critical yet poorly constrained parameter. Many studies have attempted to relate k to the main drivers of turbulence, such as wind speed (e.g., Wanninkhof, 1992; Nightingale, 2000; Raymond and Cole, 2001) or, especially for rivers, catchment parameters like slope, water flow velocity and discharge (Raymond et al., 2012). We sampled mainly the downstream reaches of the Rajang River, which range in width from 271 m to several kilometers in the delta (Allen and Pavelsky, 2018). Therefore, k parameterizations that were developed for estuaries or big rivers were considered the most appropriate. We compared three parameterizations to constrain the CO₂ fluxes. The first parameterization is the one by Borges et al. (2004) for estuaries, which is driven both by wind speed and water flow velocity; the one by Alin et al. (2011), which was developed for rivers wider than 100 m and is driven by wind speed; and the one by Raymond and Cole (2001), which is driven by wind speed and was developed for big rivers and estuaries. Those parameterizations read as follows:

$$\begin{array}{}\text{(2)}& {\displaystyle }& {\displaystyle }{k}_{\mathrm{600},\text{B04}}=\mathrm{1.0}+\mathrm{1.719}{w}^{\mathrm{0.5}}{h}^{-\mathrm{0.5}}+\mathrm{2.58}{u}_{\mathrm{10}},\text{(3)}& {\displaystyle }& {\displaystyle }{k}_{\mathrm{600},\text{A11}}=\mathrm{4.46}+\mathrm{7.11}\cdot u_{\mathrm{10}},\text{(4)}& {\displaystyle }& {\displaystyle }{k}_{\mathrm{600},\text{R01}}=\mathrm{1.91}\cdot e^{\mathrm{0.35}{u}_{\mathrm{10}}},\end{array}$$

where w is the water flow velocity (cm s⁻¹), h is the depth (m) and u₁₀ is the wind speed at 10 m (m s⁻¹). The depth h was taken from the bottom sounder recordings of our boat, and w in the lower river reaches was measured by Staub and Esterle (1993) to be 0.7 m s⁻¹. A more recent study by Ling et al. (2017) reports w= 1.1 m s⁻¹ for the Rajang River upstream from Kapit. For the calculation of the gas exchange velocity k, we used the average of w=0.9 m s⁻¹. For u₁₀, on-site wind speed data were unfortunately not available. In such cases, other authors (e.g., Bouillon et al., 2012; some estuaries in Chen et al., 2013) have resorted to gridded wind data from the NOAA NCEP–NCAR reanalysis product (Kalnay et al., 1996). While we acknowledge the uncertainty introduced by using gridded data instead of in situ wind speed, we used this product as well, as it is the best one available for our study area. We retrieved daily wind speed at 10 m for the grid centered at 2.85^∘ N, 112.5^∘ E for the time of our measurements. In situ k is dependent on in situ salinity and temperature and was calculated from k₆₀₀, exploiting its relationship with the Schmidt number (Wanninkhof, 1992).

The calculation of fluxes using k parameterizations is associated with a large uncertainty, and it is difficult to determine the most suitable parameterization if none of them was developed in the study region. In addition, using input data from the literature (as for the water flow velocity) or gridded instead of measured wind data adds to this uncertainty (see Supplement). However, using three different parameterizations we are able to constrain the magnitude of CO₂ emissions from the Rajang River. We will report the average CO₂ fluxes from the three different parameterizations as well as minimum and maximum fluxes.

2.4 Ancillary measurements

In January 2016, individual water samples were taken at 15 stations between the river mouth and Kapit, including the distributary channels Rajang and Belawai. In August 2016, water samples were taken at 34 stations, with a higher sampling frequency and coverage in the delta (Rajang, Igan, Lassa, Paloh and Belawai, Fig. 1b). Water samples were taken from approximately 1 m below the surface using a Van Dorn water sampler. Particulate material was sampled on pre-weighed and pre-combusted glass fiber filters. From the net sample weight and the volume of filtered water, the amount of suspended particulate matter (SPM) was determined. For POC, 1 molar hydrochloric acid was added in order to remove inorganic carbon from the sample. For the determination of carbon, samples were catalytically combusted at 1050 ^∘C and combustion products were measured by thermal conductivity using a Euro EA3000 Elemental Analyzer. Repeatability for C content was 0.04 % (standard deviation).

In August 2016, water samples were also taken for the determination of dissolved inorganic carbon (DIC) and the isotopic composition (δ¹³C) of DIC, because the isotopic composition of DIC can help in identifying its sources (Das et al., 2005; Campeau et al., 2017, 2018). Samples were poisoned with 200 µL concentrated HgCl₂ and filtered through Whatman glass fiber filters (GF/F, pore size 0.7 µm). Samples were filled into 40 mL vials leaving no headspace, checked for the existence of bubbles, and stored refrigerated until analysis. Concentrations and δ¹³C of DIC were measured via continuous-flow wet-oxidation isotope ratio mass spectrometry (CF-WO-IRMS) using an Aurora 1030W TOC analyzer coupled to a Thermo Delta V Plus IRMS (Oakes et al., 2010). Sodium bicarbonate (DIC) of known isotope composition dissolved in helium-purged milli-Q was used for drift correction and to verify concentrations and δ¹³C values. Reproducibility for DIC was ±10 µmol L⁻¹ for concentrations and ±0.10 ‰ for δ¹³C (standard deviations).

During both surveys, dissolved oxygen and water temperature were continuously measured with a temporal resolution of 5 min using an FDO 925 oxygen sensor and a WTW 3430 data logger (Xylem Inc., USA). The oxygen sensor was calibrated by the manufacturer, and a routine function check was performed before the start of measurements using the check and calibration vessel (FDO^© Check) provided by the company. The reported accuracy of a dissolved oxygen measurement at 20 ^∘C in air-saturated water is 1.5 %, and the precision of the accompanying temperature measurement is 0.2 ^∘C (WTW, 2012). pH, salinity and temperature were measured at the stations, using a SenTix 940 pH sensor and a multiprobe (Aquaread AP-2000). The pH sensor was calibrated before the start of the measurements using NIST (National Institute of Standards and Technology) traceable buffers. Since salinity was only measured at the stations, we spatially interpolated salinity for the interpretation of pCO₂ data. This procedure is described in the Supplement.

Different geographical extents of the river were covered during the two campaigns, and the salt intrusion limits were different during the two seasons. In order to keep the results from the two surveys comparable, we report results for three categories: nonpeat (Kapit–Kanowit), peat (Kanowit–Sibu) and delta (downstream of Sibu). Their definition and properties are specified in Table 1. The nonpeat and peat areas are directly comparable between seasons, because the same spatial extent was covered during both surveys and they were non-saline during both seasons.

Table 1 Definition of the different areas considered in this study.

2.5 Data analysis and export calculations

Data analysis was performed with Python 2.7.15 and ArcMap 10.5. Averages of measured parameters are reported ±1 standard error unless stated otherwise. Errors for calculated parameters (e.g., river loads; see below) were determined with error propagation. For fluxes and derived quantities, we report the mean, minimum and maximum from the three k parameterizations. Seasonal differences were tested for significance using the Mann–Whitney U test from the Python Scipy Statistical Functions module. Data from the delta were excluded from the statistical tests due to the different geographical coverage achieved during the two surveys.

In order to calculate the total carbon export from the Rajang River for the months of our measurements, we derived DOC load, POC load, DIC load for the peat and nonpeat area, and CO₂ outgassing as follows:

The river loads of DOC, POC and DIC were calculated for the peat and nonpeat area combined, using

$$\begin{array}{}\text{(5)}& \text{RIVER LOAD}=C\cdot Q\cdot f_{\mathrm{3}},\end{array}$$

where C is the average concentration of DOC, POC and DIC (mg L⁻¹), Q is monthly discharge (m³ s⁻¹) and f₃ is a conversion factor from s⁻¹ to month⁻¹.

For DOC, we used the DOC concentrations reported by Martin et al. (2018). These data were acquired during 2017 downstream of Kanowit. Only freshwater values were considered (average for the wet and dry season: 2.0 and 2.1 mg L⁻¹). For POC, we used the area-weighted average concentration of peat and nonpeat river reaches determined during our surveys.

For DIC, we used an area-weighted average concentration as well, which was determined from our measurements during the dry season. For the wet season survey, DIC was calculated from pH and pCO₂ using the program CO₂sys (Lewis and Wallace, 1998). Note that pH measurements were only available at the stations, and sometimes we did not have parallel pCO₂ measurements. Therefore, the number of calculated DIC values for the peat and nonpeat area is 6.

For the calculation of total CO₂ emissions from FCO₂, the river surface area was required. River surface area was calculated from the GRWL (Global River Widths from Landsat) database (Allen and Pavelsky, 2018) using Esri's ArcMap 10.5. Missing segments in the delta were manually delineated using a Landsat satellite image and their surface area was determined. This procedure is described in the Supplement. With the surface area of individual river segments at hand, CO₂ emissions were calculated for the nonpeat area, peat area and the delta separately.

A Keeling plot (Keeling, 1958) was used to explore possible sources of DIC in the Rajang River. The Keeling plot method has been used to determine the isotopic signature in CO₂ from ecosystem respiration (Pataki et al., 2003; van Asperen et al., 2017) by plotting δ¹³C of CO₂ in an air sample vs. the inverse CO₂ concentration. From the y intercept of a linear regression, the isotopic signature of the source can be determined (Keeling, 1958; Pataki et al., 2003). The Keeling plot method assumes mixing of only two components: One background component (e.g., atmospheric background) and one additional source (e.g., respiration). While originally designed for atmospheric research, it has also been used in studies exploring possible sources of DIC in stream water (Campeau et al., 2017, 2018). However, it has to be interpreted with caution, because rivers are open systems where the basic assumption of two-component mixing is easily violated. When we interpret the Keeling plot in this study, we look at mixing between river DIC and added DIC from groundwater and peat-draining rivers. We thus have to assume that DIC in the main river does not change over the time of our measurements. Since the water is continuously replenished from upstream, this assumption is valid as long as the upstream (background) source's isotopic signature remains the same. This can be assumed during a relatively short sampling period (7 days). Nevertheless, we can only consider freshwater data, as mixing with seawater already constitutes a third component. Therefore, we show only freshwater data in our Keeling plot and caution that its interpretation is based on the relatively strong assumption that we can view the river's main stem δ¹³C of DIC as constant over the sampling period.

3 Results

3.1 General characterization of the Rajang River

Measured salinity ranged between 0 and 18.6 during the wet season and 0 and 32.1 during the dry season. Saltwater was detected further upstream during the dry season than during the wet season (Fig. 3a and b). Saltwater penetrated further inland in the Rajang and Belawai distributaries than in the Igan distributary (Fig. 3b), suggesting that most freshwater is discharged via the Igan distributary.

The Rajang River was slightly acidic (wet: 6.7; dry: 6.8), area-weighted mean for the peat and nonpeat area; see Table 2) and highly turbid, with area-weighted average SPM concentrations of 187.2±75.7 mg L⁻¹ (wet season) and 51.5±12.1 mg L⁻¹ (dry season; see Table 2). With higher SPM during the wet season (p<0.001), the organic carbon content of SPM was significantly decreased (1.5±0.4 % on average, p=0.01) compared to the dry season (2.1±0.6 % on average; see Table 2). POC ranged from 0.7 to 9.1 mg L⁻¹ during the wet season (average 2.6±0.6 mg L⁻¹) and from 0.3 to 1.9 mg L⁻¹ during the dry season (average 1.1±0.4 mg L⁻¹; see Table 2). The seasonal difference was significant (p=0.002).

Figure 3 Salinity (interpolated) and pCO₂ distribution in the Rajang River and delta during the wet season survey (a, c) and dry season survey (b, d).

Table 2 Average values of measured parameters for different river reaches (mean ± 1 SE) for all parameters except FCO₂, which is reported as mean (minimum − maximum).

^* Calculated, not measured values.

The river water was consistently undersaturated with dissolved oxygen (DO) with respect to the atmosphere. DO oversaturation was not observed. The area-weighted average DO was similar during the wet and dry seasons (wet: 81.1±5.5 %; dry: 79.6±3.5 %), with slightly lower DO in the delta (wet: 66.0±6.9 %; dry: 71.2±11.1 %).

Measured DIC in the dry season ranged from 153.7 µmol L⁻¹ in the nonpeat area to 2399.2 µmol L⁻¹ in the estuary and increased linearly with salinity (Fig. 4a, r=0.98, p<0.001). Concentrations averaged 177.9±22.4 µmol L⁻¹ in the peat and nonpeat area and 1302.3±749.2 µmol L⁻¹ in the delta (Table 2). Calculated DIC for the wet season averaged 289.8±32.1 µmol L⁻¹ (area-weighted mean for the nonpeat and peat area). δ¹³C-DIC ranged between −11.9 ‰ and −1.4 ‰ and averaged $-\mathrm{7.0}\pm \mathrm{1.5}$ ‰ in the peat and nonpeat area and −5.9 ‰ in the delta (Table 2). δ¹³C-DIC was positively correlated with DIC for the delta (r=0.81, p<0.001) and negatively correlated with DIC for the freshwater part (peat and nonpeat combined, $r=-\mathrm{0.98}$, p=0.004, Fig. 4b). The Keeling plot for freshwater samples revealed a linear relationship (Fig. 4c) with a y intercept (±SE) of $-\mathrm{18.6}\pm \mathrm{0.3}$ ‰.

3.2 Carbon dioxide

The Rajang River was oversaturated with CO₂ with respect to the atmosphere, with an average pCO₂ of 2531±188 µatm (wet season) and 2337±304 µatm (dry season) in the nonpeat area. pCO₂ was significantly higher in the peat area (p<0.001, both seasons), with 2990±239 µatm (wet season) and 2994±141 µatm (dry season). The area-weighted means for the peat and nonpeat areas were 2540±189 µatm (wet season) and 2340±301 µatm (dry season). In the delta, pCO₂ was more variable, and the average values of 3005±1039 µatm (wet season) and 2783±1437 µatm (dry season; see Table 2) were also significantly higher than in the nonpeat area (p<0.001, both seasons). pCO₂ values were strikingly similar between wet and dry seasons, and so were the spatial patterns in pCO₂ (Fig. 3c and d). Tidal variability of pCO₂ was observed at an overnight station in Sarikei in August 2016. During this time, pCO₂ increased from approximately 3000 to almost 6000 µatm during rising tide (see Supplement).

pCO₂ decreased with increasing salinity in the delta (Fig. 3). However, a big spread of data in both the high-salinity region (during the tidal measurement described above) and the freshwater region was observed. pCO₂ was correlated with DO (Fig. 5). An interesting pattern is consistently visible in both the wet and dry season data, by which main stem data can clearly be distinguished from those collected in the Belawai and Paloh distributaries. Unfortunately, our data did not allow identification of a diurnal signal for either pCO₂ or DO. In the tidal part of the river, we had only one stationary measurement overnight, when a diurnal signal could not be identified due to the strong tidal signal. In the non-tidal part of the river, we had insufficient nighttime data to make a statement about a day–night difference for pCO₂ and DO.

Wind speed in the grid centered at 2.85^∘ N, 112.5^∘ E averaged 0.57 m s⁻¹ during our campaign in January 2016 and 1.09 m s⁻¹ during our campaign in August 2016. The calculated gas exchange velocities for a Schmidt number of 600 using the B04 model (k_600,B04) were 8.2 and 9.6 cm h⁻¹, respectively. This compares to the A11 model with 8.5 and 12.2 cm h⁻¹ and to the R01 model with 2.3 and 2.8 cm h⁻¹ for the wet and dry season, respectively. The resultant CO₂ fluxes (FCO₂) to the atmosphere ranged between 0.5 and 2.4 g C m⁻² d⁻¹ in the wet season and between 0.6 and 3.5 g C m⁻² d⁻¹ in the dry season (per water surface unit area; see Table 2).

3.3 Carbon river load and CO₂ emissions

Discharge was above average during 2016 (Fig. 2a). Discharge during January 2016 (wet season) was in accordance with the long-term average, but discharge during August 2016 (dry season) was higher than usual. The Rajang River loads for the peat and nonpeat area were 104 Gg C (82–123 Gg C) in January 2016 and 65 Gg C (45–83 Gg C) in August 2016; another 31 Gg C (12–41 Gg C) (wet) and 34 Gg (12–51 Gg C) were emitted as CO₂ from the delta (Table 3). Of the river loads of carbon, 91 % (86–97 %) (wet) and 82 % (70–94 %) (dry) were exported laterally by discharge. Approximately half of the laterally transported carbon was in the organic form (wet: 57±12 %; dry: 60±18 %). River loads were similar during both seasons, except that POC export was 3-fold higher in the wet season. CO₂ emissions to the atmosphere accounted for 9 % (3–14 %) and 18 % (6–30 %) (wet and dry) of the total carbon load of the river.

Table 3 Average discharge and calculated carbon loads and emissions to the atmosphere, estimated for the months of the two surveys and the whole year 2016.

4 Discussion

4.1 Sediment yield and organic carbon load

The Asia-Pacific region is known for its high sediment yields, especially where rivers drain Cenozoic sedimentary and volcanic rock (Milliman and Farnsworth, 2011). Therefore, the high suspended sediment load in the Rajang River is not surprising. However, SPM concentrations during our expeditions were substantially lower (187.2±75.7 and 51.5±12.1 mg L⁻¹) than in July 1992 (613 mg L⁻¹; Staub and Esterle, 1993). This could be an effect of upstream dams (operational since 2011 and 2015), which trap sediment in their reservoirs, thereby reducing downstream sediment loads (Vörösmarty et al., 2003; Snoussi et al., 2002). In support of this, SPM concentrations were intermediate in the upper Rajang River in 2014–2015 (218.3 mg L⁻¹; Ling et al., 2017). These measurements were taken before, and the measurements in the current study were taken after, the Murum dam began full operation in the second quarter of 2015. Furthermore, SPM and POC concentrations (wet: 2.6; dry: 1.1 mg L⁻¹) in the Rajang River were similar to those in the Pearl River, China, (SPM: 70–247 mg L⁻¹, POC: 1.0–3.8 mg L⁻¹; Ni et al., 2008) and the Red River, Vietnam, (SPM: 294±569 mg L⁻¹ (wet) and 113±428 mg L⁻¹ (dry), POC: 3.7±2.0 mg L⁻¹ (wet) and 1.1±1.1 mg L⁻¹ (dry); Le et al., 2017), both of which are also affected by damming.

SPM was higher during the wet season than during the dry season, in agreement with observations at the Kinabatangan River, Malaysia (Harun et al., 2014). This can be attributed to enhanced erosion during the wet season. In logged-over forest, such as that found in most of the Rajang River basin, the energy impact of rain drops on the soil is higher than in densely vegetated areas, where rain is intercepted by the canopy before falling on the ground (Ling Lee et al., 2004). In agreement with this line of reasoning, Ling et al. (2016) showed that the amount of suspended solids in Malaysian streams draining areas with logging activities increased significantly after rain events. The decreased organic carbon content observed during the wet season further supports this, as it indicates a higher contribution of eroded mineral soil to SPM. This pattern is observed in many rivers in this region (Huang et al., 2017). Despite the changing carbon content, most POC was still exported during the wet season, as in other Southeast Asian rivers (Ni et al., 2008; Moore et al., 2011). The proportion of laterally transported carbon in the Rajang River that is in organic form (wet: 57±12 %; dry: 60±18 %) is similar to that reported for the carbon flux to the South China Sea (50±14 %, Huang et al., 2017).

4.2 Inorganic carbon load

4.2.1 DIC concentrations and sources

DIC concentrations measured during our dry season survey were comparable to those determined by Huang et al. (2017) for the Rajang River (201 and 487 µmol L⁻¹). Their values are based on seven measurements taken between 2005 and 2009 downstream of Sibu (Ting-Hsuan Huang, personal communication, 2018). DIC concentrations in the Rajang River thus appear to be substantially lower than in the Mekong River (1173–2027 µmol L⁻¹; Li et al., 2013) and the Pearl River (1740 µmol L⁻¹; Huang et al., 2017), but similar to those in the Musi River, Indonesia (250 µmol L⁻¹, Huang et al., 2017). This suggests that the Rajang River compares better to the equatorial Indonesian rivers than to rivers draining mainland Southeast Asia, probably because of the scarcity of carbonate rock, which has the highest weathering rate and is thus responsible for high DIC in rivers (Huang et al., 2012).

Figure 4 DIC vs. salinity (a), δ¹³C vs. DIC (b) and a Keeling plot of δ¹³C vs. 1∕DIC for nonpeat, peat and delta samples, respectively. Freshwater samples, including those from the delta region, are circled in panel (b), and only those are shown in the Keeling plot in panel (c). All data refer to dry-season samples.

Figure 5 Correlation of pCO₂ vs. dissolved oxygen (DO) in the wet and dry season for individual distributaries.

The source of DIC varied along the length of the Rajang River. In the estuarine part, the positive linear relationship between DIC and salinity (Fig. 4a) suggests that the main source of DIC in the estuary is marine. This is also supported by the relatively high δ¹³C-DIC of estuarine samples, as ocean DIC is more enriched, with δ¹³C-DIC between 0 ‰ and 2.5 ‰ (Rózanski et al., 2003). The positive correlation between DIC and δ¹³C in the delta thus implies an increasing contribution of marine DIC.

In the freshwater part of the Rajang River, δ¹³C values were more depleted ($-\mathrm{7.0}\pm \mathrm{1.5}$ ‰, Table 2) than downstream, but were more enriched than δ¹³C-DIC values reported for other rivers in the region (Lupar and Saribas Rivers in Sarawak, −15.7 ‰ to −11.4 ‰, Müller et al., 2016; Musi, Indragiri and Siak Rivers in Indonesia, −22.5 ‰ to −9.0 ‰, Wit, 2017).

Multiple sources and processes are likely to influence δ¹³C-DIC in the Rajang River. To start with, the y intercept of the Keeling plot for freshwater samples suggests that the DIC added to the main river has a δ¹³C of −18.6 ‰, which is consistent with δ¹³C values of bicarbonate from silicate weathering with soil CO₂ from C3 plants (−22.1 ‰ to −16.1 ‰, Das et al., 2005), i.e., groundwater input.

Another relevant source is rainwater DIC with a typical δ¹³C of −9.3 ‰ (Das et al., 2005). In a river with relatively low DIC (177.9 µmol L⁻¹) and large surface runoff due to heavy rain, this source term is presumably non-negligible, and a significant contribution would partially explain the relatively high δ¹³C-DIC values.

With regard to in-stream processes, photosynthesis increases and respiration decreases δ¹³C-DIC (Campeau et al., 2017). Due to the high turbidity, it can be assumed that photosynthesis in the Rajang River is negligible. In contrast, the correlation of DO and pCO₂ (Fig. 5) suggests that respiration is important. This assumption is supported by the negative correlation of δ¹³C and DIC for freshwater samples, because with increasing DIC, δ¹³C values get more depleted, suggesting that organic carbon (with a δ¹³C of around −26 ‰ for C3 plants; Rózanski et al., 2003) is respired to CO₂ within the river. However, we observed relatively high δ¹³C values overall. Two processes are likely to be responsible for the downstream increase in δ¹³C: (1) methanogenesis and (2) evasion of CO₂. (1) Campeau et al. (2018) observed a strong relationship between CH₄ concentration and δ¹³C-DIC in a boreal stream draining a nutrient-poor fen, suggesting that fractionation during methanogenesis leads to an increase in δ¹³C-DIC. As the peat soils in the Rajang delta are also anaerobic and nutrient-poor, it is likely that methanogenesis plays a role there as well. This is consistent with high reported soil CH₄ concentrations of up to 1465 ppm in a peat under an oil palm plantation in Sarawak (Melling et al., 2005). It would therefore be of great interest to investigate CH₄ concentrations in the Rajang River in the future. (2) CO₂ evasion is also known to lead to a gradual increase in δ¹³C-DIC values until in equilibrium with the atmosphere (with a slightly positive δ¹³C-DIC; Polsenaere and Abril, 2012; Venkiteswaran et al., 2014; Campeau et al., 2018). Due to intracarbonate equilibrium fractionation, dissolved CO₂ is more depleted in δ¹³C than the other carbonate species. Thus, if it is removed, δ¹³C of the remaining DIC increases. This effect depends on pH and is more pronounced in near-neutral waters and less strong in very acidic water (Campeau et al., 2018). We sampled the lower river reaches downstream of Kapit, which corresponds to approximately the last 200 km of the river. In addition, the terrain is much steeper upstream of Kapit than in the lower river reaches, so that CO₂ fluxes to the atmosphere are presumably much higher due to higher turbulence. This means that a large fraction of CO₂ had likely already been emitted from the river surface before reaching Kapit, leading to the observed high δ¹³C-DIC values.

Table 4 pCO₂, pH and CO₂ evasion of several Southeast Asian rivers flowing into the South China Sea.

^* The sampling points were located outside of the peat area, so the actual peat coverage at that point was zero.

4.2.2 pCO₂

pCO₂ in the Rajang River was higher than in most Southeast Asian rivers without peat influence, such as the Mekong River (pCO₂=1090 µatm; Li et al., 2013), the Red River (pCO₂=1589 µatm; Le et al., 2017) or the freshwater parts of the Lupar and Saribas Rivers (pCO₂=1274 µatm and 1159 µatm; Müller et al., 2016, Table 4, Fig. 6). This could be attributed to the peat influence. However, the Rajang River has a pCO₂ at the low end of values reported for peat-draining rivers (Fig. 6): Wit et al. (2015) report values between 2400 µatm in the Batang Hari (peat coverage = 5 %) and 8555 µatm in the Siak River (peat coverage = 22 %).

Figure 6 Comparison of average pCO₂ values in Southeast Asian rivers. Colors distinguish peat-draining rivers from nonpeat-draining rivers. The Rajang River (this study) is highlighted in red.

A meaningful comparison is also the one between the Rajang River and the Indragiri River, Indonesia, because they have a similar peat coverage (Rajang: 11 %, Indragiri: 12 %) and peat coverage has previously been considered as a good predictor of river CO₂ emissions (Wit et al., 2015). However, pCO₂ in the Indragiri (5777 µatm) was significantly higher than in the Rajang, which was associated with a lower pH (6.3, numbers from Wit et al., 2015). A simple exercise using CO₂Sys illustrates how important it is to consider both pH and pCO₂ when comparing different peat-draining rivers. At the given temperature, salinity and pH, the pCO₂ of 5777 µatm in the Indragiri corresponds to a DIC value of 327 µmol L⁻¹. At a hypothetical pH of 6.8, as measured in the Rajang River, this DIC value corresponds, under otherwise unchanged conditions, to a pCO₂ of 2814 µatm – which is very close to the average values measured in the peat area of the Rajang River. The close coupling of pH and pCO₂ implies that peat coverage in a river basin is insufficient as the sole predictor of CO₂ fluxes. Rather, the relationship between peat coverage and pCO₂ must be viewed in the context of the rivers' pH, and drivers of pH must be more carefully considered. Note also that peat coverage is usually reported for the entire catchment (e.g., Wit et al., 2015; Rixen et al., 2016) and does not reveal how much peat is found in estuarine vs. freshwater reaches, which complicates comparisons further.

δ¹³C-DIC in the Indragiri was lower (−16.8 ‰, Wit, 2017) than in the Rajang (−7.0 ‰). While this can be partially attributed to the different pH, it can also be interpreted as a greater contribution of respiratory CO₂ in the Indragiri, while the Rajang might be more strongly influenced by weathering. This would also help to explain the higher pH in the Rajang River.

While DOC and pCO₂ are positively related to discharge in most rivers (e.g., Bouillon et al., 2012), this pattern is sometimes reversed in peat-draining rivers. This is due to dilution, when rainfall exceeds the infiltration capacity of the wet soil and water runs off at the surface (Clark et al., 2008; Rixen et al., 2016). In the Rajang River, pCO₂ was slightly higher in the nonpeat area during the wet season, in agreement with many nonpeat-draining tropical rivers (Bouillon et al., 2012; Teodoru et al., 2015; Scofield et al., 2016). However, the seasonality of pCO₂ was very small, similar to other Malaysian rivers (Müller et al., 2016) and in line with the small seasonal variability of DOC concentrations in the Rajang River (Martin et al., 2018). Note that since our data were collected during two single surveys, they represent only a snapshot and do not allow strong claims about seasonality.

Due to an El Niño event, temperatures in Southeast Asia were unusually high in late 2015 and 2016, with a temperature extreme in April 2016 prevailing in most of Southeast Asia (Thirumalai et al., 2017) and unusually hot conditions also recorded in Sibu (Fig. 2b). Given that weathering rates increase at elevated temperatures, DIC from weathering could have been enhanced over other years, although this seems unlikely because DIC was relatively low and in agreement with previous studies. Another factor to be taken into account is that decomposition in the dry upper soil layer is more intense at higher temperatures, and with incipient rainfall, all the resultant DOC is flushed out to the river. Therefore, it is possible that during the year of our measurements, DOC concentrations were higher than usual, and respiratory CO₂ may therefore have been enhanced compared to other years.

4.2.3 Impact of the peatlands on the CO₂ emissions from the Rajang River

The fact that pCO₂ was significantly higher in the peat area than in the nonpeat area implies, at first glance, that the peat areas are a source of CO₂ to the river. However, this difference has to be interpreted with caution, as the entire peat area is under tidal influence (Fig. 1b). In the following, we will present several arguments that suggest that the peatlands exert only a small influence on the CO₂ emissions from the Rajang River.

• a.

  DOC. One indicator of peatland influence on a river's carbon budget is DOC. DOC concentrations in the Rajang River delta were reported to range between 1.4 and 3 mg L⁻¹ (Martin et al., 2018). This is at the low end of the range of DOC concentrations reported for peat-draining rivers in Indonesia: These range from 2.9 mg L⁻¹ in the Musi River (peat coverage = 3.5 %) up to 21.9 mg L⁻¹ in the Siak River (peat coverage = 22 %; Wit et al., 2015). Rivers whose catchment area is entirely covered by peat exhibit even higher DOC concentrations, with 52 mg L⁻¹ (wet season) and 44 mg L⁻¹ (dry season) in the Sebangau River, Indonesia (Moore et al., 2011), and 44 mg L⁻¹ in the Maludam River, Sarawak, Malaysia (Müller et al., 2015). The Rajang River compares rather to rivers like Lupar and Saribas, Malaysia, which exhibit DOC concentrations of 1.8 and 3.7 mg L⁻¹ in their freshwater parts (no peat influence, Müller et al., 2016). Consequently, DOC concentrations imply that the peatlands' influence on the Rajang's DOC is rather small.

• b.

  Mixing model. An alternative approach is to theoretically calculate the increase in Rajang River pCO₂ that would result from the influx from peatlands. For this, we created a simple model to simulate the mixing of two water masses (see Supplement), one with a pH of 6.8 and a pCO₂ of 2434 µatm (designed to resemble the Rajang River, nonpeat area) and the other with a pH of 3.8 and a pCO₂ of 8100 µatm (designed to resemble peat-draining tributaries, based on values for the peat-draining Maludam River in Sarawak, Müller et al., 2015). For simplicity, we assumed that mixing occurs at salinity = 0 and that the temperature in both water bodies is the same (28.4 ^∘C). From these values, DIC and total alkalinity (TA) of the two water bodies were calculated using CO₂Sys. Since they can be assumed to be conservative, we simply calculated DIC and TA of the mixture as

  $$\begin{array}{}\text{(6)}& {\displaystyle }& {\displaystyle }{\text{DIC}}_{S=\mathrm{0}}=\left(\mathrm{1}-\text{pc}\right)\cdot {\text{DIC}}_{\mathrm{1}}+\text{pc}\cdot {\text{DIC}}_{\mathrm{2}}\text{(7)}& {\displaystyle }& {\displaystyle }\text{and}\text{(8)}& {\displaystyle }& {\displaystyle }{\text{TA}}_{S=\mathrm{0}}=\left(\mathrm{1}-\text{pc}\right)\cdot {\text{TA}}_{\mathrm{1}}+\text{pc}\cdot {\text{TA}}_{\mathrm{2}},\end{array}$$

  where pc is the peat coverage in the basin (pc = 0.11). From DIC and TA, the pCO₂ of the mixture was computed (pCO₂=3058 µatm). This means that if all peat-draining tributaries in the Rajang River basin had a pCO₂ of 8100 µatm and a pH of 3.8, the pCO₂ in the peat area would be enhanced by around 600 µatm. However, this increase in pCO₂ is obviously gradual. For example, in the city of Sibu, peat coverage was estimated at around 2 %, for which the described mixing model yields a pCO₂ of 2548 µatm (Fig. S3). In most parts of the delta, dilution with seawater already plays a role. Therefore, the mixing model was extended, assuming that at pc = 3 %, salinity is still zero and then linearly increases until pc = 11 %, S=32, DIC${}_{S=\mathrm{32}}=\mathrm{2347}$ µmol L⁻¹ and TA${}_{S=\mathrm{32}}=\mathrm{2324}$ µmol L⁻¹ (two end-member mixing model; see Supplement). As a result, pCO₂ would theoretically not exceed 2605 µatm if the peat-draining tributaries were the only source of CO₂ downstream of Kanowit (Fig. S3).

  However, pCO₂ in the peat area was 2990 µatm (wet) and 2994 µatm (dry) and 3005 µatm (wet) and 2783 µatm (dry) in the delta, so there must be another source of CO₂. Since pCO₂ in Sarikei varied 2-fold with the tidal cycle, it seems likely that a large part of the difference in pCO₂ between the nonpeat area and peat area is actually a difference between river and tidal river. Tidal variability is often seen in estuaries (e.g., Bouillon et al., 2007; Oliveira et al., 2017), largely due to conservative mixing. However, during rising tide in Sarikei, we observed pCO₂ values of almost 6000 µatm, which is higher than the freshwater end-member, suggesting that other effects also play a role. Among those are decomposition of organic matter in intertidal sediments (Alongi et al., 1999; Cai et al., 1999) and subsequent transport of the produced CO₂ to the river, as well as groundwater input (Rosentreter et al., 2018).

• c.

  δ¹³C-DIC. If the peatlands acted as a significant source of CO₂ to the Rajang River, it would be expected that this had an impact on the δ¹³C-DIC values. In the peat-draining Maludam River, δ¹³C-DIC averaged −28.55 ‰ (Müller et al., 2015). Thus, the influx of peat-draining tributaries to the Rajang River (with an average δ¹³C-DIC of −7 ‰) would theoretically decrease δ¹³C-DIC. This was not observed. Instead, the main source identified by the Keeling plot was consistent with groundwater input, with a δ¹³C-DIC of −18.6 ‰. However, given the constraints on the applicability of the Keeling plot method discussed above, the effect of DIC inputs from peat on δ¹³C-DIC might simply not be detectable. Other processes which influence δ¹³C-DIC as discussed above might also prevent identification of a peat signal in the δ¹³C-DIC data.

In summary, we were unable to discern a large impact of peatlands on the DIC budget of the Rajang River. It is possible that, because the peatlands are located close to the coast in this system, mixing with seawater occurs before significant effects on the pCO₂ are theoretically possible. This means that not only the peat coverage in the catchment is relevant, but also how much of this peat is found in estuarine reaches. These findings support the arguments of Müller et al. (2015) and Wit et al. (2015) that material derived from coastal peatlands is swiftly transported to the ocean, explaining why peat-draining rivers may not necessarily be strong sources of CO₂ to the atmosphere.

5 Conclusions

The Rajang River is a typical Southeast Asian river, transporting large amounts of terrestrial material to the South China Sea. In contrast to other Southeast Asian rivers with similar peat coverage, the impact of the peatlands on the Rajang River's pCO₂ appeared to be rather small, probably due to the proximity of the peatlands to the coast. As a consequence, CO₂ emissions from the Rajang River were moderate compared to other Southeast Asian rivers and low if compared to Southeast Asian peat-draining rivers.