How temporal patterns in rainfall determine the geomorphology and carbon fluxes of tropical peatlands

Soil respiration.

Soil respiration measurements were made with an automated dynamic chamber system (LI-8100 with LI-8150 multiplexer; LI-COR Biosciences), in which carbon dioxide concentrations were measured after periodic closure of four chambers connected to an infrared gas analyzer via a gas multiplexer. Each chamber was closed for 6 min and a flux computed every hour, using instrument software (LI-8100 software version 2.0.0). Chambers were supported 30 cm above the peat by four 2-inch PVC pipes, 1 m long, to avoid submergence of the chambers. Although temperature affected individual flux measurements around the diurnal cycle, the effect of temperature on daily mean fluxes was negligible. Water table height was simultaneously monitored using a logging pressure transducer as described below (piezometers).

Microtopography transect.

We used a total station [TC-GTS-105N (60536) 5“ Totalstation; B. L. Makepeace, Inc.] to survey peat microtopography and 12 piezometers along a 180-m transect at the field site in August 2012 (Fig. 2 D and E). We refer to the piezometers in this transect as “transect piezometers” to distinguish them from the 5 “trail piezometers” along the 2.5-km trail. We located the transect survey points in geographic space by rotating the points to match the compass orientation between the first two piezometers and then shifting the rotated points so that the position of the first piezometer matched its GPS coordinates (Fig. 2E).

To obtain the elevation of the survey points above mean sea level, we matched the survey data to the LiDAR data. The LiDAR digital terrain model (DTM) was defined by finding local minima in last-return elevations on a 20-m × 20-m grid. Therefore, we found local minima in surveyed peat surface points on the same grid and then aligned the (arbitrary) vertical coordinate of the minimal survey points and the DTM elevation in the corresponding grid cell (Fig. 3B).

We then defined a reference land surface p in the transect as follows. First, we fitted a line to survey points by orthogonal distance regression and then projected all points onto that line. To find the water table elevation in each piezometer, we subtracted the water table depth (Solinst Model 101 Water Level Meter; Solinst Canada) from the elevation of the piezometer casing. We then found the calibrated water table elevation vs. time in the piezometer by matching the measured water table elevation to the piezometer output at the time of measurement (Fig. 3 B and C). We defined the land surface p in the microtopography transect as a linear interpolant of the average water table in each piezometer, shifted up to align with the local minima in the surveyed peat surface points (Fig. 3B). By this definition, the land surface p is smooth because its shape comes from the water table, and it passes through the bottom of hollows because it is shifted to the local minima in the peat surface.

We next used the transect hydrologic and microtopographic data to define the land surface at each trail piezometer. Our definition of the land surface p as a smooth surface fitted through the bottom of hollows makes it easy to estimate the local land surface elevation by looking at a nearby hollow. However, because this reference land surface is smoothed, we cannot easily find the exact height of trail piezometers relative to the land surface. The trail piezometer nearest the transect piezometers (80–130 m away) had almost identical water table behavior to the transect piezometers during the time interval of the transect piezometer data (Figs. 3C and 5A ). Therefore, to find the land surface elevation for each trail piezometer, we first found the land surface for this nearest piezometer by matching the data on the overlapping time interval. Because the behavior of the water table was so similar in all trail piezometers except the bog plain piezometer, for plotting in Fig. 5 we simply matched the water table height ζ=H−p for these piezometers to that of the piezometer nearest the transect by subtracting the difference in mean. Although this may introduce some error for the piezometer outside the area of uniform water table behavior (the bog plain piezometer), we did not use that piezometer to determine the specific yield and transmissivity functions for simulations.

LiDAR data and topography.

LiDAR data were obtained from the Brunei Survey Department from a 2009 aerial mission. We constructed a DTM from the LiDAR data by taking the minimum last-return elevation within 20-m × 20-m grid cells and then smoothing the resulting surface, corresponding to our definition of a reference land surface p as a smooth surface fitted through hollows (local minima in p∼). We then created a contour vector map from the DTM (GRASS GIS 6.4.3; grass.osgeo.org/) and outlined a flow tube containing our piezometer and cores by constructing boundary flow lines from contours using the TAPES-C algorithm (35), tracing back to an estimated location of the groundwater divide on the Malaysian side of the border (Fig. 2C). Although the groundwater divide is on the Malaysian side of the border where LiDAR data are unavailable, the tapered shape of the flow tube ensures that the simulations are insensitive to the precise location of the divide (Fig. 6).

We sampled the magnitude of the gradient along each contour, using a spatial database (SpatiaLite 4.2.0; www.gaia-gis.it/gaia-sins/). By the divergence theorem, the average Laplacian of surface elevation within an enclosed region is equal to the normal surface gradient integrated around its boundary divided by the enclosed area. We estimated the Laplacian in the region around the piezometer by linear regression of the integrated normal gradient against the enclosed area in the four piezometers nearest the river, where the slope was approximately uniform (Fig. 6).

Peat cores.

Methods for the two peat cores and 13 samples taken for radiocarbon dating with a modified Livingstone corer are available in Dommain et al. (27) and those for the seven peat cores and 22 radiocarbon samples taken with a Russian peat sampler are in Gandois et al. (31). Median calendar ages were estimated from radiocarbon ages, using Calib 7.1.0 with the IntCal13 calibration curve (36, 37). Because cores were taken in hollows, the top of each core corresponds roughly to the reference land surface p, which is fitted through the bottoms of hollows.

Piezometers.

Water table height was recorded every 20 min at each of five piezometers along a 2.5-km transect, using logging pressure transducers (Solinst Levelogger Edge 3001; Solinst Canada) (Fig. 2C). Transducer measurements were corrected for barometric fluctuations, using a barometer (Solinst Barologger Gold 3001). Each transducer was suspended on a steel cable inside a 2-inch PVC pipe 1.5 m in length installed to 1.4 m depth and screened at 1.3–1.45 m below the top of the piezometer casing. Piezometer locations were determined with a GPS (GPSmap 60csx; Garmin Ltd.).

An additional 12 piezometers were installed along an 180-m transect 1 km from the river (Fig. 2 D and E and SI Methods). These piezometers were screened from 1 m below the peat to 5 cm above the peat and were attached to 2-inch PVC pipes anchored in the underlying clay. Data from these piezometers were recorded from August 17, 2012 through November 3, 2012.

Throughfall measurement.

Three siphoning tipping-bucket rain gauges (Texas Electronics TR-525S) were installed along the microtopography transect to record throughfall in the forest understory. Each gauge was mounted on an acrylic plate supported 50 cm above the peat surface on a 2-inch PVC tube pushed into the peat. Large fronds of understory plants (mostly P. andersonii) were cleared from above the gauges to reduce the spatial variability in measurements. Tip events were recorded by a switch closure event logger (UA-003-64; Onset Computer Corp.) enclosed in a small PVC case attached to each gauge.

Master recharge and recession curves.

The master recharge and recession curves were computed based on the uniform water table behavior among the four piezometers nearest the river. A related approach, without the assembly of master curves but using the division of water table time series into periods of heavy rain, days without precipitation, and nights without precipitation, has been used previously in other studies of wetland hydrology (28, 29). Under very heavy rain the groundwater divergence term ∇·(T∇H) in the groundwater flow equation (Eq. 1) becomes negligible,

and because specific yield Sy is a function of water table relative to the surface ζ only, an integral form of the simplified differential equation

corresponds to the master recharge curve, giving the cumulative rainfall depth P driving a change in water table height. During dry intervals, “net precipitation” consists entirely of evapotranspiration, and taking evapotranspiration as constant and equal to its average through the diurnal cycle ET¯ (nonnegative) yields a simplified differential equation without rain,

which, in its integral form,

is equivalent to the master recession curve (Fig. 5F) turned on its side. The Laplacian ∇2p was obtained from topographic data as described above (SI Methods).

We assembled the master recharge and recession curves by using throughfall data (above) to identify intervals of the time series either with no rain or with heavy rain (rainfall intensity greater than 4 mm/h). Within each interval, we resampled a spline of the head vs. time data to yield time on a uniform grid of water table height. We then found an offset in time or rain depth for each interval to minimize the least-squares difference between overlapping intervals, creating the master recharge and recession curves (SI Master Recession and Recharge Curve Assembly). At low water tables, the water table was nearly constant at night and decreased linearly during the day (Fig. 5F), so we estimated evapotranspiration as 2.246 mm/d, using the slope of the daytime decrease in water table and specific yield obtained from the master recharge curve.

We estimated parametric functions for the specific yield (cubic spline) and transmissivity (piecewise-linear logarithm of hydraulic conductivity) by simultaneously fitting the recharge and recession curves to model results (Levenberg–Marquardt; www.pesthomepage.org). We simulated recharge [water table height vs. cumulative rainfall depth ζ(P)] or recession [water table height vs. time since rain ζ(t)] by solving the relevant differential equation (Eq. S5 for recharge, Eq. S7 for recession), using an integrating solver (SUNDIALS; computation.llnl.gov/projects/sundials) with the current parameter estimates. Although the data provide no information on how conductivity is distributed in the peat below the lowest excursion of the water table, the dynamics of the water table are always dominated by other factors, either the conductivity of peat higher in the soil profile or evapotranspiration.

Calibration for peat surface evolution parameters.

After estimation of hydrologic parameters T,Sy, we estimated the peat accumulation parameters fp,α by fitting the simulated shape of our peat dome to the LiDAR surface, using a nonlinear least-squares optimization routine (Levenberg–Marquardt; www.pesthomepage.org). The coastal peat swamp forests of Brunei Darussalam are underlain by a broad, very gently sloped mangrove clay (1, 27), so as the initial condition we used a flat, sloped surface corresponding to a gradient in the elevation of basal peat in our peat cores (Fig. 7A).

We drove our simulation of the peat dome using a precipitation time series derived from our own throughfall data, meteorological data, and nearby speleothem δ18O records, using an approach similar to that of Kurnianto et al. (16). Our goal was not to reconstruct the actual rainfall at our site over the last 3,000 y, but to reasonably approximate the fluctuation in rainfall at subannual, decadal, and millennial timescales. For subannual rainfall, we cycled through field throughfall intensity data from 2012 on a 20-min grid from our site. We then scaled throughfall intensities by a factor specific to each simulation year and month as follows. First, we followed the approach of Moerman et al. (40) and calculated a regression of 2-mo rainfall means at a nearby meteorological station [Brunei International Airport (BWN), 90 km northeast of our site] against δ18O at Gunung Mulu airport, 61 km southeast of our site, from August 2006 through April 2011 (P= 2.53006mm / d−δ18O× 1.00825mm / d, R2= 0.208513). We then used this regression to compute mean rainfall intensities for each decadal to centennial interval captured by the speleothems at Mulu, which are believed to characterize regional patterns of rainfall (40, 41). We then repeated 44 y of monthly precipitation totals from the meteorological station (1966–2009, mean annual rainfall 2.90 m) for 3,000 y, adjusting monthly means to match the linearly interpolated precipitation averages from the speleothem data regression. Finally, we scaled 20-min intensities in each month of the 3,000-y time series to match the month precipitation total determined by the meteorological data and speleothem record, scaled by a constant for all months and years so that precipitation in the measured year (2012) matched throughfall at our site. Evapotranspiration was modeled as zero at night and constant by day, equal to effective evapotranspiration estimated from piezometer data as described above. Because the details of changes in river stage over the last 3,000 y are unknown, we used dates in the core closest to the river to drive the gradually shifting surface elevation and assumed that average river stage remained in the same position relative to the surface at the boundary.

Local water balance is dominated by flows near the surface.

Eighteen months of data on water table height in five piezometers along a 2.5-km transect (Fig. 5) show two distinctive features of water table behavior in tropical peatlands. First, when the water table is high, it falls very rapidly; and second, the water table height relative to the land surface remains approximately uniform in all piezometers as the water table rises and falls, as observed elsewhere by Hooijer (28). In what follows, we use “water table height” ζ=H−p to refer to the water table height relative to the land surface, as distinct from the water table elevation H above mean sea level. Because the water table height ζ is approximately uniform, the water table behavior can be summarized by a pair of curves describing the uniform rise of the water table during heavy rain and the uniform decline of the water table during dry intervals between rains (Fig. 5 E and F). During heavy rain, the effects of evapotranspiration and outward flow are negligible, and the rainfall intensity vs. rate of increase in water table height gives the specific yield. Between rain events, the water table declines because of evapotranspiration and the divergence of groundwater flow.

Transmissivity T is a function of water table height ζ and controls the divergence of groundwater flow ∇·(T∇H). We determined the effect of water table height on transmissivity using our water table data. The water table declines during dry intervals because of a combination of evapotranspiration and the divergence of groundwater flow; however, the two are easily distinguished at low water tables because evapotranspiration ceases at night (Fig. 5D). Therefore, we can obtain the divergence of groundwater flow from the declining water table during dry intervals after accounting for evapotranspiration (refs. 28, 29; further details are in SI Methods). We find that transmissivity increases exponentially at high water tables, when water rises into hollows and flows through hummocks, but decreases dramatically at low water tables when water flows through fine pores in the peat matrix (Fig. 5C). Very high permeability near the peat surface is consistent with our observations of more void space higher in the peat profile and also with recent data from other tropical peatlands (30). The water table curves (Fig. 5 E and F) indicate that the near-surface permeability is so great that the total thickness of deeper peat is unimportant for groundwater flow. Therefore, transmissivity is approximately independent of peat depth and depends only on the water table height ζ, which is uniform in space (although highly variable in time).

Morphology of peat surface explains uniform water table behavior.

According to Boussinesq’s equation, uniform transmissivity is not, by itself, enough to explain the uniform fluctuation of the water table. Even in hydrologic systems where hydraulic properties are uniform, the water table can behave differently at different locations because of topography. For example, in most hydrologic systems a rainstorm drives a different water table response at a topographic divide than it does near where groundwater discharges to a river.

To understand the uniform water table behavior in peatlands, we refer back to Boussinesq’s equation (Eq. 1). If both the specific yield Sy and the transmissivity T depend only on the local water table height relative to the surface and not on position within the peatland, uniform water table movement occurs if the divergence of the peat surface gradient, or the peat surface Laplacian ∇2p, is uniform (Fig. S2 C–E). (The “Laplacian of the peat surface” ∇2p, or just “Laplacian,” is the scalar result of applying the Laplacian operator ∇2 to the land surface elevation p.) To see why a uniform land surface Laplacian explains uniform water table behavior, we rewrite Boussinesq’s equation (Eq. 1) in terms of the water table height relative to the land surface (ζ=H−p), instead of the water table elevation H:

We observe that water table height is uniform (∇ζ= 0). If transmissivity T is also spatially uniform, the groundwater divergence term simplifies to the transmissivity times the peat surface Laplacian (∇·[T∇(p+ζ)]=T∇2p). The time derivative ∂p/∂t of the land surface elevation is negligible because peat accumulation or loss is much slower than rise or fall of the water table, so the term p can be dropped from the time derivative. We observe that the fluctuations in water table height ∂ζ/∂t are uniform, as is net precipitation qn, so the groundwater divergence term T∇2p must also be spatially uniform. Thus, Boussinesq’s equation simplifies to an ordinary differential equation (ODE) describing the uniform fluctuation of the water table relative to the peat surface

where the peat surface Laplacian ∇2p is uniform.

The peat surface Laplacian describes the curvature of the peat surface: It is equal to the sum of the second derivatives of the surface elevation in two perpendicular horizontal directions (∇2p=∂2p/∂x2+∂2p/∂y2). Thus, analysis of water table dynamics predicts uniform curvature of the peat surface where water table fluctuations are uniform. This uniformity of surface elevation curvature can be tested against elevation maps.

Maps of the peat surface Laplacian are highly sensitive to microtopographic noise in the surface elevation map because the Laplacian uses the second derivative of the surface elevation. However, by the divergence theorem, the average Laplacian within any closed contour is equal to the integral of the normal gradient along the contour divided by the enclosed area. Therefore, we can examine the uniformity of the surface Laplacian by studying the slope of a regression between the integrated normal gradient and the enclosed area (Fig. 6). Indeed, we find a linear relationship between the integrated normal gradient along each contour and the area enclosed by the contour in our LiDAR-derived peat surface elevation map, indicating a uniform surface Laplacian in the region of uniform water table behavior (Fig. 6). In contrast, outside the region of uniform Laplacian, the water table behaves differently (“bog plain piezometer” in Figs. 5 and 6).

Uniform surface Laplacian determines stable tropical peatland morphology.

The uniform peat surface Laplacian provides a remarkably simple way to compute a stable morphology for a tropical peat dome. By “stable morphology,” we mean a morphology in which the peat surface and water table continue to fluctuate with the vagaries of climate, but there is no long-term average change in the peat surface or water table elevation (they are stationary; ⟨∂p/∂t⟩= 0,⟨∂H/∂t⟩= 0). Uniform water table height is the simplest behavior that could make an entire peatland stable, because if the water table height is spatially uniform, the local rate of peat accumulation is also uniform. In a stable peatland, there is no long-term change in the water table height, so any water added by net precipitation must eventually be removed by groundwater flow

Thus, the Laplacian ∇2p∞ of the stable peatland surface p∞ is minus the average net precipitation divided by the average transmissivity

We can compute the stable topography of any tropical peatland by solving Poisson’s equation (Eq. 7) for the stable peat surface morphology p∞, using the appropriate Laplacian value for that climate. The average transmissivity ⟨T⟩ is a complicated function of the temporal pattern of rainfall and the hydrologic–biological system. However, for any rainfall regime, one can find the stable surface Laplacian ∇2p∞ by repeatedly simulating water table fluctuations (Eq. 5) with a trial Laplacian ∇2p and adjusting the Laplacian value until peat production balances decomposition (Eq. 3) everywhere in the peatland (SI Methods). In this way, one finds a shape parameter (∇2p∞) that describes stable peatland morphology under a given rainfall regime in any drainage network.

Climate and drainage network determine tropical peatland carbon storage capacity.

By specifying the stable peatland topography, the uniform-Laplacian principle gives the peat carbon storage capacity inside any drainage boundary and in any given climate. The volume under the surface satisfying Poisson’s equation times the mean carbon density of the peat gives the carbon storage capacity of the peatland. For example, the peat dome at our primary site currently has a mean peat depth of 3.88 m (max 4.92 m) and stores about 1,535 metric tons (t) C·ha−1; however, if the climate remains similar to the climate during its 2,300-y development, we predict that in about 2,500 y it will reach a stable shape with a mean peat depth of 4.54 m (max 7.10 m) and store 1,800 t C·ha−1 (Fig. S3; simulations of dynamics are described in the next section).

The uniformity of the stable peat surface Laplacian is an approximation that requires that (i) peat accumulation rate ∂p/∂t is a nondecreasing function of water table height, (ii) flow of water is proportional to water table gradient (Boussinesq’s equation), and (iii) transmissivity is independent of location because flow through deep peat is negligible compared with near-surface flow. In reality, groundwater flow through deeper peat will result in a deviation of the stable peat dome surface from the uniform-Laplacian shape in very large peat domes. Specifically, groundwater flow through deep, low-permeability peat will tend to flatten the dome center, because of slow infiltration of water into the deep peat, and steepen the dome margin, because of exfiltration of water back into the high-permeability near-surface peat near the boundary. Deep groundwater flow should be manifested as a downward (dome center) or upward (dome margin) trend in the water table during nights without rain when the water table is low; no such trend is apparent in our piezometer data (Fig. 5D), suggesting that deep groundwater flow is small. A small deep groundwater flow term is further supported by radiocarbon dating of porewater dissolved organic carbon at our site (31), which suggests a maximum downward velocity of water of about 1 m/y or at most a 1.4-mm water table decline during a single 12-h night, 1/16th of the 22-mm water table decline from evapotranspiration during the day (Fig. 5). (Evapotranspirative flux is about 1/10th of the rate of decline of the water table from evapotranspiration because about 1/10th of the deep peat cross-section is available for water flow; see specific yield curve in Fig. 5B.)

A shape parameter related to our stable peatland Laplacian (Eq. 7) appeared in Ingram’s model for temperate peatland morphology (10) assuming constant precipitation, uniform hydraulic conductivity, and simple river geometry (Ingram’s parameter is net precipitation divided by hydraulic conductivity, instead of average transmissivity). Our result is more general, because it handles varying rainfall and arbitrary landscapes, but is also mathematically simpler, because of our finding that transmissivity in tropical peatlands is approximately independent of peat depth.

Peat accumulation parameters regulate dome dynamics.

Our analysis shows how the rate of peat production fp and decomposition rate constant α affect both the stable morphology and the dynamics of tropical peat domes. These parameters of the peat accumulation function (Eq. 2) have an indirect but strong effect on the stable peat surface Laplacian and hence peatland carbon storage capacity via the mean transmissivity ⟨T⟩ (Eq. 7) because the mean water table depth must be equal to the ratio of the peat production rate to the decomposition rate constant (fp/α; Eq. 3). A higher decomposition rate constant implies a higher mean water table in stable peat domes, meaning a higher transmissivity, a smaller stable surface Laplacian, and less carbon storage. If both peat production fp and the decomposition rate constant α increase together, carbon storage capacity does not change, but peat dome dynamics are faster.

Fit parameters match literature data.

Peat accumulation parameters fitted to the topography of a peat dome at our Brunei field site agree with published data from other sites and also with our other field data (next section). We obtained peat accumulation parameters fp,α by simulating the evolution of the dome (Fig. 7) and minimizing the least-squares difference between the simulated peat surface and the modern peat surface measured by LiDAR. We then compared our calibrated peat accumulation function to literature data on subsidence in drained, vegetated peat swamps (4, 22). Our linear peat accumulation function was not calibrated to these subsidence data from the literature—only to the modern peat surface—but nonetheless matched the subsidence data almost exactly (Fig. 4A; fp= 1.46mm·y−1,α= 1.80d−1). Our soil CO₂ chamber measurements were also very similar to those from other sites, suggesting that the effect of water table on fluxes is similar at our site and in other tropical peatlands (Fig. 4B).

The uniform-Laplacian principle predicts a central bog plain and old peat near the surface at bog margins.

We find that a tropical peat dome reaches its stable shape first at its boundaries, because the stable dome surface is lowest there (Figs. 7 and 8 and Fig. S2). Meanwhile, the interior of the peat dome continues growing at an approximately uniform rate, forming a relatively flat (smaller-magnitude Laplacian) central bog plain. The vegetation of tropical bog plains may not be distinct (1), unlike the unforested bog plains of high-latitude peatlands (21); instead, we define the bog plain of a tropical peat dome as the central region that has not yet reached its stable Laplacian. Whereas the dome center continues to accumulate peat and sequester carbon, the margin has reached its stable shape and stopped growing, so peat near the surface is older there.

Older peat near dome margins has not been predicted before, so we collected 22 additional radiocarbon dates from basal and near-surface peat samples to test this prediction. These radiocarbon dates confirmed that near-surface peat was older near dome margins than at the same depths toward the interior of the same domes (Fig. 7C). We also compared radiocarbon dates in deeper peat to simulated ages at the same locations and depths, excluding basal samples from the mangrove peat before the establishment of the peat swamp forest (Fig. 7 and SI Methods) (1, 27). Radiocarbon dates and simulated ages at the same locations and depths matched well (Fig. 7B). We did not expect radiocarbon dates from cores to match simulated peat ages exactly because (i) the drainage network may have shifted during the 2,300 y of dome growth, (ii) tree root growth may inject young carbon into peat below the surface, and (iii) tree falls in peat swamp forests remove older peat to form tip-up pools that then fill with younger peat. In an earlier study, we estimated that replacement of older peat by younger peat in tip-up pools would bias radiocarbon dates of deep peat to about 500 y later than when material was first deposited in that stratum (figure 11 in ref. 27), consistent with the offset between measured radiocarbon dates and ages simulated by our model (Fig. 7B).

Carbon sequestration rate is proportional to bog plain area.

The centripetal pattern of dome development makes the rate of carbon sequestration roughly proportional to the area of the central bog plain (Fig. 8). Under a given climate, the rate of sequestration decreases as the dome approaches its stable shape and the central region of peat accumulation—the bog plain—shrinks in area. For example, our simulations imply that the current rate of CO₂ sequestration at our site (0.80 t · ha⁻¹ · y⁻¹, 100-y average) is less than one-quarter of its initial rate about 2,300 y ago (3.81 t · ha⁻¹ · y⁻¹), and CO₂ sequestration is more than five times faster at the dome interior (1.89 t · ha⁻¹ · y⁻¹, 6.37 km from the river) than at its edge (0.36 t · ha⁻¹ · y⁻¹, 1 km from the river; Fig. S3). The mechanism of tropical peat dome development that we describe therefore creates landscape-scale patterns in local carbon fluxes and radiocarbon date profiles. Local measurements of carbon fluxes or radiocarbon dates cannot be upscaled to regional fluxes without considering dome morphology because the flatter interior of each peat dome sequesters carbon whereas the margins do not (Fig. 8). Old peat near the peatland surface (2), although in some cases caused by local climate change or disturbance, also can be expected at the margin of any peat dome.

Subdivision of a peatland by drainage canals reduces carbon storage.

The average surface elevation of a stable peat dome is proportional to the area of the dome because of the uniform-Laplacian principle. If we scale the area of a peat dome by some factor k by multiplying both x and y coordinates by k, the surface elevation p must increase by the same factor k to keep the same Laplacian. Therefore, the carbon storage capacity of a peat dome scales with its area. For example, a peat dome that is cut into halves of approximately the same shape as the original dome will have one-half the carbon storage capacity (half the mean stable peat depth) of the original dome. This provides a straightforward way to estimate the long-term impacts of artificial drainage networks that are now affecting over 50% of the peatlands of Southeast Asia (32) and from which a robust quantification of carbon emissions is urgently needed (6).

The dynamic response of a peat dome to changes in rainfall and sea level also depends on its area because of the centripetal pattern of dome development (Fig. 8). Because of their higher stable mean depth, larger-area domes reach their stable shape more slowly than smaller-area domes.

Relative effects of climate change on carbon storage capacity are independent of drainage network.

Although peatland drainage networks play a central role in determining absolute carbon storage and dynamics, we can calculate the proportional effect of climate change on long-term carbon storage of a tropical peatland independent of the drainage network. Poisson’s equation (Eq. 7) must be solved in each drainage boundary to obtain the topography of the stable peat surface. However, we can then predict the effects of changes in climate independent of the drainage network because of the linearity of the Laplacian operator. By the definition of linearity for a mathematical operator, a peat surface Laplacian that is larger by some factor k corresponds to a peat surface that is vertically stretched by the same factor (k∇2p=∇2kp) and therefore has a mean peat depth that is larger by the same factor. Thus, carbon storage capacity per area p¯∞ is proportional to the stable peat surface Laplacian p¯∞∝∇2p∞.

Dynamic simulations converge to new stable morphologies after changes in conditions.

Our simulations of peat dome dynamics demonstrate the convergence of initially stable domes to new, stable, uniform-Laplacian morphologies after perturbations (Fig. 9). The simulations show the effect of increased total rainfall (Fig. 9 A and E), which is a recognized climate feedback for tropical peatlands (12), and also show that artificial drainage for agriculture (Fig. 9D) can dominate all natural feedbacks if not curtailed (4, 16). In addition, our simulations demonstrate a third feedback: The increase in rainfall variability from warming climates (33) can cause peat loss if not compensated by an increase in total rainfall (Fig. 9 C and F). For these simulations, we generated new rainfall time series as similar to current rainfall as possible but with larger annual and El Niño–Southern Oscillation (ENSO) fluctuations (Fig. S4 A and B and SI Setting Annual and ENSO Amplitudes of Rainfall). Either greater seasonality or a stronger ENSO decreased peatland carbon storage capacity, but an increase in seasonality had a larger maximum effect, partly because the magnitude of the ENSO fluctuation is smaller. In contrast, sea level rise could drive peat accumulation in the long term by elevating the tidal rivers draining most peat domes (Fig. 9 B and E). In general, losses can be much more rapid than accumulation (Fig. 9E), because subsidence of drained peatlands can be far faster than typical accumulation rates (4). For example, the estimated area-averaged current CO₂ sequestration rate at our site is 0.80 t · ha⁻¹ · y⁻¹, whereas Hooijer et al. (5) estimated CO₂ emissions of at least 73 t · ha⁻¹ · y⁻¹ from tropical peatlands under plantation agriculture.

Intermittency of rainfall reduces tropical peatland carbon storage.

We find that fluctuations in net precipitation on timescales from hours to years can reduce long-term peat accumulation. We further explored the effects of variability in rainfall seen in our dynamic simulations (Fig. 9) by computing the effect of interstorm arrival time and annual and ENSO fluctuations on peatland carbon storage capacity (Fig. 10). The simulations demonstrate that long-term peat accumulation is controlled by variation in rainfall, not only by mean rainfall, because fluctuations in the water table cause exponential changes in groundwater flow. The high outward flow during peak water tables is not compensated by low flow rates after the water table declines. For example, a steady drizzle at the same average intensity as the intermittent rainfall actually observed at our site would sustain more than 10 times more long-term carbon storage (19.5 kt·ha−1 vs. 1.80 kt·ha−1; Fig. S4 D and E). The intermittency of tropical convective storms significantly affects long-term carbon storage: Carbon storage capacity can decrease by one-third depending on whether convective storms arrive every 14 h on average, as at our site, or every 24 h, with the same mean rainfall (Fig. 10A).

Our simulations with smoothed rainfall intensity and evapotranspiration show that models must consider the effects of subdiurnal fluctuations in rainfall to correctly predict the long-term evolution and carbon storage of tropical peatlands. The exact details of the fluctuations in rainfall are not important, in the sense that many distinct rainfall time series can give the same stable surface Laplacian and the same carbon storage capacity. However, carbon storage capacity can be severely overestimated by simulations that entirely ignore the effects of fluctuations in rainfall. We explored the effects of neglecting fluctuations in rainfall by computing the stable surface Laplacian after averaging net precipitation on hourly and longer intervals. Treating rainfall intensity and evapotranspiration as constant each hour, instead of every 20 min, increased the simulated stable surface Laplacian by a few percent, but averaging over 1 d led to an overestimate by 20%, over 1 wk by 100%, over 1 mo by 400%, and over 1 y by more than 1,000% (Fig. S4 D and E).