Carbon storage capacity of tropical peatlands in natural and artificial drainage networks

Tropical peatlands store over 75 gigatons of carbon as organic matter that is protected from decomposition and fire by waterlogging if left undrained. Over millennia, this organic matter builds up between channels or rivers into gently mounded shapes called peat domes. Measurements of peat accumulation and water flow suggest that tropical peat domes approach a steady state in which the peat surface morphology is described by a uniform curvature, setting a limit on the carbon that a peatland can store. We explored the maximum amount of carbon that can accumulate as water-saturated peat in natural and artificial drainage networks of northwest and southern Borneo. We find that the maximum volume of peat accumulation in a channel-bounded parcel is proportional to the square of the parcel area times a scale-independent factor describing the shape of the parcel boundary. Thus, carbon capacity per area scales roughly with mean parcel area in the peatland. Our analysis provides a tool that can be used to predict the long-term impacts of artificial drainage, and to devise optimal strategies for arresting fires and greenhouse gas emissions in tropical peatlands.


Introduction
Tropical peatlands are estimated to store over 75 Gt carbon [1][2][3][4], preserved as peat where decomposition of woody organic matter is slowed by waterlogging [5]. Analogous ecosystems of the Carboniferous and Miocene created some of the world's largest deposits of coal [6,7], and many extant peatlands have continuously sequestered carbon for 5 000 years or more [8,9], thereby exerting a cooling effect on the global climate system [10,11]. However, in the last 30 years, artificial drainage of tropical peatlands for agriculture has been rapidly releasing this carbon via organic decomposition [12][13][14] and fire [15,16]. Decomposition and fires in drained peatlands of Southeast Asia are now a major source of greenhouse gas emissions 6 Mott MacDonald Singapore Pte Ltd 189721 Singapore, Singapore. [13] and transboundary smoke haze [15,16], having released about 2.5 Gt C since 1990 [13] and exposed millions of people to unhealthy levels of airborne particulates [17,18]. Tropical peatlands of Africa [3] and the Americas [2] could become similarly vulnerable if drained [19]. Because drainage has led to large emissions, restoration of tropical peatlands by blocking artificial drainage networks has been proposed as a cost-effective way to help meet climate targets [14,20,21]. Planning the restoration of tropical peatlands and understanding their role in Earth's carbon cycle require an understanding of how much carbon peatlands can preserve by waterlogging, and how it is affected by the geometry of natural and artificial drainage networks.
In tropical peat domes, peat can remain waterlogged in deposits that rise over 10 m above the water level in their bounding channels [22] because water that drains to dome margins is replaced by rainfall [5,[23][24][25]. Nonetheless, a peat dome cannot grow in height indefinitely, because if peat is piled too steeply, the surface peat cannot remain waterlogged and will decompose. Thus, the maximum height of a peat dome is limited by the size and shape of the channel-bounded parcel on which it accumulates. In the 20th century, this constraint was explored in raised bogs, which are analogous to tropical peat domes, through analyses of peatland hydrology and peat accumulation [26][27][28][29]. These analyses predict the profile of equilibrium peat elevation in a crosssection through a bog, assuming that the boundaries of the bog are long, straight lines formed by two parallel rivers [28,29], as may occur, for example, in post-glacial landscapes [30,31]. These assumptions do not apply where boundaries created by drainage networks are complex, and thus these profiles cannot directly predict overall limits on peat accumulation in most drainage networks.
In earlier work [25], we derived a model for the limiting morphology of tropical peat domes with boundaries of any size or shape. The model is based on field measurements showing that flow in tropical peatlands occurs predominantly in the loosestructured matrix of peat near the surface, where conductivity is very high [32,33], so that the efficiency of flow through the peatland, or hydraulic transmissivity, varies greatly in time but is approximately independent of peat thickness. Using this approximation, we showed that the limiting morphology of a tropical peat dome satisfies Poisson's equation, in which the curvature κ of the peat surface is uniform throughout the peatland (figure 1). This prediction was validated by collecting radiocarbon dates from a peat dome in Brunei Darussalam Darussalam, showing that the dome was at or near carbon equilibrium where it conformed to this limiting morphology. However, the analysis examined part of just one of many domes in the drainage network, and did not examine the overall limits on carbon storage imposed by the geometric shape and size of parcels in the peatland.
Over half of the tropical peat landscapes in Peninsular Malaysia, Sumatra and Borneo have been further partitioned by canals and ditches [34] primarily because crops such as oil palm [35] and acacia require that surface peat is dry. Ditch spacing in plantations is designed to ensure that the maximum water level in the peat is close to the water level in ditches [35,36]. This criterion implies that the maximum height of waterlogging above the drainage network in plantations, and therefore the limiting peat height, is small. However, drainage in tropical peatlands has generally been evaluated not by its effect on the capacity of the peatland to store carbon, but by rates of peat thickness loss or carbon emissions. Data on peat thickness loss per time, peat carbon density, and typical water levels are used to derive emission factors for different land cover classes, which are then multiplied by land cover area to derive regional estimates of emissions [13,20,[37][38][39][40]. The emission factor approach provides a snapshot of overall rates of peat loss in drained peatlands, but it does not attempt to predict changes in peatland morphology, which control the long-term evolution of carbon emissions.
In this paper, we apply the analysis developed by Cobb et al [25] to examine how morphology and the long-term capacity for carbon storage are constrained in peatland landscapes. We describe the boundaries of parcels created by natural and artificial drainage networks at sites in northwest and southern Borneo, and solve Poisson's equation in each parcel to analyze how distributions of parcel shape and size control the potential for peatland carbon storage. To illustrate how our approach can be used to evaluate the effects of land management decisions, we also calculate the long-term impact of a canal on the morphology and carbon storage capacity of a peat dome in northwest Borneo. Our analysis shows that fires and greenhouse gas emissions from Southeast Asian peatlands can be understood as a side effect of their re-equilibration after dissection by canals and ditches, and demonstrates how drainage network geometry can be taken into account when evaluating land-use impacts and mitigation strategies.

Carbon storage capacity
To evaluate the capacity of tropical peatland drainage networks to store carbon, we used an equation derived in [25] that relates equilibrium peat dome morphology (figure 1) to the efficiency of lateral flow through peat, or transmissivity, T; time-averaged precipitation P; and evapotranspiration ET. This analysis showed that the curvature κ of the peat surface, defined using the Laplacian of the peat surface elevation (κ =−∇ 2 p), reaches an ecosystem-and climatespecific equilibrium curvature κ * given by κ * = ⟨P − ET⟩/⟨T⟩ where ⟨⟩ indicates time averaging (a brief derivation is provided in the appendix). The equation is based on the following approximations, supported by literature data: (1) peat oxidation increases with water table depth [40]; (2) there is a sharp increase in hydraulic transmissivity as the water table rises [32,33], so that to a good approximation, transmissivity of the peat is a function only of the water level relative to the surface; and (3) the water level relative to the surface remains spatially uniform across the domed peatland even as it fluctuates with time [24].
Perhaps surprisingly, the form of the equilibrium peat morphology equation makes it possible to write a parcel's maximum carbon storage C as a product of independent factors describing the ecosystem and climate on one hand, and the geometry of the parcel boundary on the other  The maximum volume of peat that can be stored per area in a parcel can be broken down into the area of the parcel (horizontal axis) and a scale-independent parcel 'shape factor' (slope of line). Because of the smaller shape factor of a 4:1 rectangle than a square, the 1 4 km 2 rectangular parcels in (g) cannot store as much carbon as the 1 4 km 2 squares in (f), and creation of canals as d → e → g results in the loss of more potential for carbon storage than creation of canals as d → e → f.
where the peat carbon density ρ c and equilibrium curvature κ * are determined by the ecosystem and climate [25]. The other terms SA 2 are purely geometric, and provide a climate-and ecosystem-independent measure of the capacity for carbon storage in a parcel relative to a 1 km 2 disk-shaped parcel in the same ecosystem and climate (appendix). The parcel shape factor S is a scale-independent factor determined by the geometric shape of the parcel and A is the parcel area ('parcel' refers to any channel-bounded area where peat accumulates). The parcel shape factor S must be computed numerically for realistic shapes, but is independent of the parcel size; it is at most 1, for a disk (the shape that can support the highest water table in a fixed area), and approaches 0 in a parcel that is infinitesimally narrow (figure 2).

Analysis of parcel geometry in Borneo peatlands
To probe the distributions of peatland parcel areas A and shape factors S in real landscapes, we examined parcels created by artificial and natural drainage networks at sites in northwest and southern Borneo: in northwest Borneo, oil palm plantations near the Baram River ('Baram'), and natural peat domes near the Belait River ('Belait'); and in southern Borneo, an area drained for the ex-Mega Rice Project area ('EMRP'), and natural peat domes near the Katingan and Sebangau Rivers ('Sebangau'; figure 3). At all four sites, peat domes were built up during the Holocene by floristically similar hardwood tropical peat swamp forests [8]. At the two sites in northwest Borneo (Baram and Belait), peat domes were initiated between about 5 000 and 1 000 cal BP on mangrove and marine clay deposits [8,23]. The Baram site was drained for oil palm plantation in 2012-2013 [20], while the drainage network at the Belait site has not been altered. At the two sites in southern Borneo (EMRP and Sebangau), peat domes were initiated on podzolic terraces between about 14 000 and 9 000 cal BP, during a period of increasing rainfall and rising sea levels [8]. The EMRP site was drained by large primary canals excavated beginning in 1995 as part of a large rice cultivation project. The project was subsequently abandoned [20] without construction of the network of field drains typical of plantation agriculture.
For each site, we digitized parcel boundaries using Google Earth (www.google.com/earth) at a resolution of about 100 m cm −1 , except for areas in southern Borneo without high-resolution images, where we used 200 m cm −1 . We demarcated parcels using visible streams and rivers where possible, excluding flooded forest that can occur at confluences and the edges of large rivers. In the Belait and Sebangau regions, we looked for parcels minimally affected by drainage, but used extant drainage boundaries, even if artificial. In some areas of the Sebangau region, artificial channels were indistinct because of a combination of lower-resolution imagery and land use and we instead drew the parcel boundary to exclude the general area that appeared drained. We used a geospatial data translator library (http://gdal.org) to convert each polygon to a suitable projected coordinate system (EPSG:5247 GDBD2009 for northwest Borneo, EPSG:32 749 WGS 84 / UTM zone 49S for southern Borneo) and load it into a spatial database (http://postgis.net), which we used to measure the area A of each polygon. We then generated a quadrilateral mesh inside each parcel [41] and solved Poisson's equation on the mesh with a solver written in Cython (http://cython.org) using the deal.II finiteelement library [42]. We then sampled the numerical solution (using VTK Python, http://vtk.org/) on a square grid chosen to give at least 2000 points within each polygon (grid spacing of 1 m, 8 m, 20 m, 70 m, and 100 m for Baram, EMRP, Belait, Badas, and Sebangau regions, respectively) to determine the shape factor S (figure 4(a)).

Carbon capacity in a peat dome cut by a canal
As an example of the effect of a canal on peatland carbon storage capacity, we examined a site in northwest Borneo where a large peat dome, the Badas peat dome, was divided by a canal before 1972. To calculate the carbon storage capacity of the Badas peat dome before and after its drainage by the canal, we drew a boundary around the entire dome and then around the two smaller parcels created by the canal using Google Earth, as described above. We then solved Poisson's equation on the pre-canal parcel and the two post-canal parcels, and evaluated the limiting dome volume as done for other parcels. We used the equilibrium peat surface curvature κ * = 2.36 × 10 −3 km −1 from an earlier topographic analysis of a peat dome 7 km to the south [33], which hosts the same monodominant forest types and has the same developmental history as the other Belait domes [23]. The total limiting volume of each dome above the drainage network was calculated as the area of the boundary polygon multiplied by the mean limiting peat elevation from solving Poisson's equation, sampled every 70 m on an east-west, north-south grid as described above. Finally, we estimated the total carbon storage capacity of the pre-and post-canal drainage boundaries at this site by multiplying the pre-and post-canal volumes by the carbon density ρ c = 35 kg m −3 [43] measured in cores taken from the same dome used for estimation of equilibrium peat surface curvature.

Results and discussion
Our analysis shows how the partitioning of a peatland into parcels by a drainage network constrains its capacity for peat accumulation (figures 1 and 2). In practice, plantation parcels are much smaller than natural peat domes (figure 3): in the peatlands we examined in northwest Borneo, plantation parcels are smaller in mean area than natural parcels in nearby undrained peatlands by a factor of more than 7180 (0.005 65 vs. 40.6 km 2 ). The plantation parcels also have smaller shape factors because they are long and narrow (mean S = 0.277 vs. 0.667). Thus the mean equilibrium peat height in natural peat dome parcels of the Belait district is 3.19 m and their carbon storage capacity per area averages 125 kg C m −2 , whereas the mean equilibrium peat height above drainage boundaries in the Baram plantation parcels is less than 1 mm, and their capacity to store carbon above the drainage network is negligible (less than 6 g C m −2 ). The negligible capacity of these plantation parcels for carbon storage above the drainage network reflects the almost complete oxidation of carbon stocks that will occur in the long term if these drainage networks are maintained.
The studied parcels in the EMRP area in southern Borneo are exceptionally large (figure 3) because the area was never developed for large-scale plantation agriculture [20], and are more isodiametric in shape (mean S = 0.686) than the natural parcels of Sebangau (mean S = 0.475). Nonetheless, on a geometric basis alone, the carbon storage capacity per area in the Sebangau domes is larger on average by a factor of 76.4 (522 vs. 6.83) because of their vastly greater size (mean parcel area 1370 km 2 vs. 6.86 km 2 ).
Analysis of parcels in the four sample peatlands shows that while a parcel's shape affects its capacity for carbon storage, parcel area ranges over more than 6 orders of magnitude (0.002 0-4087 km 2 ), and therefore has a dominant effect on parcel storage capacity (figure 4). By definition, the parcel shape factor cannot exceed a value of 1, but can be arbitrarily small. Thus, in principle, a parcel the size of our largest (4087 km 2 ) could have a smaller capacity than our smallest (0.002 km 2 ), if the large parcel were sufficiently long and thin (shape factor less than 5 × 10 −7 ). In practice, the largest parcel shape factor measured (0.966) was 6.14 times the smallest (0.157, figure 4(c)). Thus, although the precise shape factor for a complex parcel can only be computed numerically, the relative capacity of two networks in the same climate and ecosystem can be roughly estimated based on just parcel area: maximum peat volume scales with total peatland area times mean parcel area. If a dome is cut into one hundred parcels of similar shape, the maximum volume of each parcel will be roughly 0.01% of the dome's original limiting volume, and the total carbon storage capacity of the parcels will be roughly 1% that of the original dome.
Our analysis of carbon storage capacity separates the effects of drainage network geometry from ecosystem and climate (equation (1)). In different regions and ecosystems, peat carbon density may vary; for example, published data suggest that typical carbon densities in the southern Borneo sites are about 40% higher than in coastal peatlands of Borneo and Sumatra (61 vs. 44 kg m −3 [8]). Limiting peat surface curvatures κ * are likely to differ too because of contrasts in rainfall patterns, peat hydraulic conductivity, peatland productivity and decomposition rates [25]. Thus, to compare natural carbon storage capacities across regions and ecosystems, an estimate of the limiting peat surface curvature using the methods of Cobb et al [25] is required as well, preferably from topographic data with a high vertical resolution. For example, although the parcels of the Sebangau region have a maximum volume per area 15.4 times that of  3 and 4), each relative to a 1 km 2 disk in their respective climates, we do not expect that the maximum carbon storage per area in Sebangau is larger by the same factor because stronger seasonality of rainfall [44,45] implies a smaller limiting surface curvature κ * in Sebangau [8,20].
The parcels of Sebangau are not only much (33.7 times) larger than those of the Belait peatlands, but also have been accumulating peat for about twice as long [8], suggesting that some process of parcel coalescence may occur over long periods of time whereby a river that previously divided two parcels becomes a site of peat accumulation, as described for higher-latitude bogs [28]. Such processes are examples of how dynamics in peatland river networks [46] could, like shifts in climate, cause very long-term changes in peatland carbon storage capacity.
Computing limiting peat morphologies can also explain spatial patterns in peat emissions and flammability driven by changes in drainage networks. Our calculation for the Badas peat dome shows that the canal through the dome reduces its capacity for carbon storage by about 257 million cubic meters of peat, or about 10.1 megatons of carbon, by removing a volume shaped like an inverted saddle from its original limiting morphology (figure 5). The largest differences between the pre-canal and post-canal equilibrium morphologies occur in the vicinity of the canal, but the post-canal equilibrium morphology is lower across the entire dome, both in the west (average 1.09 m lower) and in the east part (average 2.74 m lower). Though a detailed topography of the Badas dome is not available, it is clear that the peat surface has been approaching this post-canal equilibrium because visible subsidence and recurring fires have created a V-shaped depression around the canal. The Badas canal was recently blocked as part of a restoration project. If the canal had been maintained, we would anticipate continued subsidence and fire to a distance of about 5 km from the canal (figure 5(f)) until the peatland reached a new limiting morphology with two smaller peat domes.
Although our approach does not predict the details of the transition of a peatland towards a different morphology, it yields a rough lower bound on the time that the transition will take. This lower bound can be computed by dividing the maximum difference in elevation between the morphologies by a typical rate of peat accumulation or loss. By our analysis, the maximum peat height in a parcel (like the maximum peat volume per area) scales with parcel area, so a parcel with twice the area will take roughly twice the time to equilibrate. The EMRP parcels are on average 0.08% the area of natural parcels in the Sebangau area; if the original morphology of the drainage network in the EMRP resembled that of Sebangau, the EMRP parcels would approach their equilibrium morphologies in very roughly 0.08% of the time that would be required to completely destroy the equilibrium relief of the original domes. The faster equilibration of the small parcels created by canals in the EMRP is apparent in high-resolution topographic data from the area [47], which show a domed morphology in each EMRP parcel superimposed on a large natural dome that will gradually disappear from peat loss to decomposition and fire.
As parcels become very large, processes in peat domes may begin to deviate from the approximations underlying the theory we apply. Flow through deeper peat will account for an increasing proportion of the water budget in larger domes, and could begin to limit carbon storage because water is effectively drained from the dome interior to its edges through the deep peat [48]. Anaerobic decomposition in peat below the water table could also begin to limit height in sufficiently large domes [49], though available evidence suggests that anaerobic decomposition is extremely slow in tropical peats [8]. These refinements would predict a 'leveling-off ' of the storage capacity vs. area relationship in very large peat domes due to a more flattened shape than predicted by the theory we apply here.
A similar method could be developed to evaluate long-term drainage effects in raised bogs of the Northern Hemisphere. The outline of the procedure would be identical: (1) Validate a model for the equilibrium shapes of raised bogs in arbitrary boundaries; (2) Use this model to compute limiting morphologies of bogs with existing boundaries; (3) Compute altered limiting morphologies given altered drainage; and then (4) Examine the (spatially distributed) difference between the two. The model for equilibrium bog shape (1) might be simply the analysis of Ivanov [28] and Ingram [29] extended to three dimensions [36] and solved numerically, although it might need to be refined regionally or modified because bog types are diverse [50] and raised bog cross-sections have not always matched modeled profiles [31]. The Ivanov / Ingram model assumes that the equilibrium shape of the bog is controlled by flow in permanently saturated peat below the drought year water table, and thus implies a sublinear (square-root) storage capacity vs. area curve due to the leveling-off mechanism we have just described. Whatever the model for equilibrium bog shape, a storage capacity analysis would bring the same strengths and weaknesses in northern peatlands as in the tropics: it would not predict instantaneous fluxes or detailed dynamics, but would provide spatially explicit predictions of long-term impact, and complement an emission factor approach by describing how much carbon can be preserved in the long term under different management scenarios.

Conclusions
Tropical peat domes progressively sequester carbon as waterlogged organic matter over thousands of years. Classic studies explored how the hydrologicecosystem feedback that drives peat accumulation also limits the height and shape of raised bog cross-sections. Our analysis applies an analogous approach to landscapes, and shows the limiting shapes of tropical peat domes in three dimensions. The maps produced by this analysis convey intuition about not only how drainage networks constrain the morphology of tropical peatlands, but also how re-configuration of these networks will transform peatland morphology over time. Although we do not analyze dynamics, our equilibrium model can be combined with typical rates of peat accumulation and loss to yield a rough lower bound on the time required for re-equilibration of a peatland after alteration of the drainage network. The approach can predict spatial patterns in decomposition and flammability, as well as the long-term effects of interventions such as blocking of peatland ditches and canals. Our analysis shows how the critical factors of parcel size and shape can be taken into account when designing mitigation strategies to arrest or reverse tropical peatland greenhouse gas emissions.

A.1. Equilibrium peat morphology
In [25], we derived a peat equilibrium morphology equation from Boussinesq's equation for water table elevation H subject to precipitation P and evapotranspiration ET where both the specific yield S y (ζ) and the transmissivity T(ζ) are functions of the local water level ζ, defined as the elevation of the water table H relative to the surface p (ζ = H − p). Because net peat loss shifts to net accumulation as the water level rises, there is a mean water level at which the rate of peat production by the ecosystem is balanced by the rate of decomposition in surface peat (figure 1). The peat surface morphology p that maintains a stationary and uniform water level ζ is derived by substituting ζ into the groundwater flow equation (A1) and setting the timeaveraged change in local water storage ⟨S y ∂H/∂t⟩ to zero 0 = ⟨P − ET⟩ + ∇ · ⟨T ∇ (p + ζ)⟩ .
The water level ζ is approximately spatially uniform, so its gradient ∇ζ is zero. The transmissivity T is also approximately uniform and is independent of p, so it can be moved outside the divergence operator ∇·, yielding the equilibrium peat morphology equation −∇ 2 p = κ * where κ * = ⟨P − ET⟩/⟨T⟩.

A.2. Carbon capacity
We first describe the carbon capacity of a diskshaped parcel, which we use as a reference for other calculations. The shape of a surface with uniform Laplacian ∇ 2 p = −κ * above a flat, circular boundary with radius R is a parabolic dome with equation p = κ * ( R 2 − r 2 ) /4. The limiting volume V o of this reference dome can be calculated by integration as V o = πκ * R 4 /8. Substituting in the area of the parcel A = πR 2 , the limiting volume of the dome inside can equivalently be written V o = κ * A 2 /(8π). The carbon capacity C o of this disk-shaped parcel is the peat carbon density ρ c times the limiting peat volume, We define the parcel shape factor S as the ratio of a parcel's carbon capacity to that of a disk of the same area (S = C/C o ), yielding (1).

A.3. Scaling of carbon capacity with curvature and area
We have just shown (A2) that the carbon capacity of a disk-shaped parcel is proportional to the equilibrium peat surface curvature κ * and the parcel area squared A 2 . The same is true for a parcel of arbitrary shape; we illustrate this fact informally here.
Consider a parcel of arbitrary shape with area A. The boundary of the parcel is flat and the height of