# Emergence of the Southeast Asian islands as a driver for Neogene cooling

Yuem Park, Pierre Maffre, Yves Goddéris, Francis Macdonald, Eliel Anttila, Nicholas Swanson-Hysell

<figure><img src="image_1.png" /></figure>

## To cite this version:

Yuem Park, Pierre Maffre, Yves Goddéris, Francis Macdonald, Eliel Anttila, et al. Emergence of the Southeast Asian islands as a driver for Neogene cooling. Proceedings of the National Academy of Sciences of the United States of America, 2020, 117 (41), pp.25319-25326. 10.1073/pnas.2011033117. hal-03054120

HAL Id: hal-03054120

https://hal.science/hal-03054120v1

Submitted on 26 May 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>
Title

Emergence of the Southeast Asian islands as a driver for Neogene cooling

Permalink

https://escholarship.org/uc/item/1174m63s

Journal

Proceedings of the National Academy of Sciences of the United States of America, 117(41)

ISSN

0027-8424

Authors

Park, Yuem
Maffre, Pierre
Goddéris, Yves
et al.

Publication Date

2020-10-13

DOI

10.1073/pnas.2011033117

Copyright Information

This work is made available under the terms of a Creative Commons Attribution License,
available at https://creativecommons.org/licenses/by/4.0/

Peer reviewed
# Emergence of the Southeast Asian islands as a driver for Neogene cooling

Yuem Park<sup>a,1</sup>, Pierre Maffre<sup>a</sup>, Yves Goddéris<sup>b</sup>, Francis A. Macdonald<sup>c</sup>, Eliel S. C. Anttila<sup>c</sup>, and Nicholas L. Swanson-Hysell<sup>a</sup>

This manuscript is the author version of the paper generated using the PNAS LaTeX template.

Steep topography, a tropical climate, and mafic lithologies contribute to efficient chemical weathering and carbon sequestration in the Southeast Asian islands. Ongoing arc-continent collision between the Sunda-Banda arc system and Australia has increased the area of subaerially exposed land in the region since the mid-Miocene. Concurrently, Earth's climate has cooled since the Miocene Climatic Optimum, leading to growth of the Antarctic ice sheet and onset of Northern Hemisphere glaciation. We seek to evaluate the hypothesis that the emergence of the Southeast Asian islands played a significant role in driving this cooling trend through increasing global weatherability. To do so, we have compiled paleoshoreline data and incorporated them into GEOCLIM, which couples a global climate model to a silicate weathering model with spatially resolved lithology. We find that without the increase in area of the Southeast Asian islands over the Neogene, atmospheric $pCO_2$ would have been significantly higher than pre-industrial values, remaining above the levels necessary for initiating Northern Hemisphere ice sheets.

silicate weathering | weatherability | arc-continent collision | Neogene cooling | Southeast Asian islands

The Southeast Asian islands (SEAI) have an out-sized contribution to modern chemical weathering fluxes relative to its area. The confluence of steep topography, a warm and wet tropical climate, and the presence of mafic lithologies results in high fluxes of Ca and Mg cations in the dissolved load and associated $CO_2$ consumption (1–4). There has been a significant increase in the area of subaerially exposed land within the region since the mid-Miocene associated with ongoing arc-continent collision between Australia and the Sunda-Banda arc system (5–7). Concurrently, after the Miocene Climatic Optimum, a cooling trend began ca. 15 Ma and accelerated over the past 4 million years (m.y.) leading to the development of Northern Hemisphere ice sheets (8, 9). Many hypotheses have been proposed to explain this cooling trend including changes in ocean/atmosphere circulation (5, 10, 11), a decrease in volcanic degassing (12), or uplift in the Himalaya (13, 14). Here we seek to evaluate the hypothesis that emergence of the SEAI was a significant factor in driving long-term climatic cooling over the Neogene.

Over geologic time-scales, $CO_2$ enters Earth's ocean–atmosphere system primarily via volcanism and metamorphic degassing, and leaves primarily through the chemical weathering of silicate rocks and through organic carbon burial (15). Chemical weathering delivers alkalinity and cations to the ocean which drives carbon sequestration through carbonate precipitation. Steady-state $pCO_2$ is set at the $pCO_2$ level at which $CO_2$ sinks are equal to sources. As $CO_2$ sinks increase and $pCO_2$ falls, temperature decreases and the hydrological cycle is weakened, causing the efficiency of the silicate weathering sink to decrease until a new steady-state is achieved at lower $pCO_2$ (16).

Topography, climate, and lithology all effect chemical weathering. High-relief regions generally lack extensive regolith development, and thus tend to have reaction-limited weathering regimes that are more prone to adjust when climate changes (17–19). High physical erosion rates contribute to high chemical weathering fluxes in these high-relief regions (20). In warm and wet regions, mineral dissolution kinetics are faster leading to enhanced chemical weathering (18, 21). Mafic rocks have higher Ca and Mg concentrations and dissolution rates than felsic rocks, and thus have the potential to more efficiently sequester carbon through silicate weathering (22). These factors have led to the proposal that arc-continent collisions, which create steep landscapes that include mafic lithologies, within the tropical rain belt have been important in enhancing global weatherability, lowering atmospheric $pCO_2$, and initiating glacial climate over the past 520 m.y. (7, 23, 24) and perhaps in the Neoproterozoic as well (25).

Quantitatively estimating the magnitude of decrease in steady-state $pCO_2$ associated with the emergence of a region with a high carbon sequestration potential, such as the SEAI, requires constraints on changing tectonic context and accounting of associated feedbacks. As this region emerges, the total global silicate weathering flux will transiently exceed the volcanic degassing flux, causing $pCO_2$ to initially decline until a new steady-state is established where the total magnitude of the $CO_2$ sinks is the same as before the change. However,

## Significance Statement

The Southeast Asian islands are a modern-day hotspot of $CO_2$ consumption via silicate weathering. Since ~15 million years ago, these islands have been increasing in size at the same time that Earth's climate has been cooling. Here we test the hypothesis that this global cooling could have been driven by tectonic emergence of the Southeast Asian islands. Using a new compilation of paleoshorelines in conjunction with a coupled silicate weathering and climate model, we find that this emergence is associated with a large decrease in $pCO_2$. Without these changes in tropical island paleogeography, there would not have been large Northern Hemisphere ice sheets as a defining feature of Earth's climate over the past 3 million years.

Author contributions: Y.P., P.M., Y.G., and N.L.S.-H. designed and implemented the GEOCLIM model. F.A.M. and E.S.C.A. constructed the paleoshoreline database. Y.P. and P.M. executed the model and analyzed the data in consultation with N.L.S.-H and Y.G.. Y.P., P.M., Y.G., F.A.M., and N.L.S.-H. wrote the manuscript.

The authors declare no competing interest.

<footnote><sup>a</sup>Department of Earth and Planetary Science, University of California, Berkeley, CA 94720, USA; <sup>b</sup>Géosciences Environnement Toulouse, CNRS–Université Paul Sabatier - IRD, Toulouse 31400, France; <sup>c</sup>Department of Earth Science, University of California, Santa Barbara, CA 93106, USA</footnote>

<footnote n="1">¹To whom correspondence should be addressed. E-mail: yuempark@berkeley.edu</footnote>
the sensitivity of the silicate weathering flux in any particular
location to this change in pCO₂ is variable and dependent on
the specific topography, climate, and lithology at that location.
Furthermore, how regional climate responds to this change
in pCO₂ is itself spatially variable. Therefore, the magnitude
of pCO₂ change that is required to balance the total global
silicate weathering flux with the volcanic degassing flux will
depend on the specific spatial distribution of topography, cli-
mate, and lithology at the time of emergence. As a result, any
attempt to meaningfully estimate the decrease in steady-state
pCO₂ associated with emergence of the SEAI must model
spatially resolved climatology and silicate weathering fluxes in
tandem and account for the spatial distribution of the factors
that affect these inter-connected systems.

GEOCLIM Model

To estimate the decrease in steady-state pCO₂ associated with
the increase of subaerially exposed land area in the SEAI,
we use the global spatially resolved GEOCLIM model (26).
GEOCLIM estimates changes in steady-state pCO₂ associated
with coupled changes in erosion, chemical weathering, and
climatology by linking a silicate weathering model to climate
model runs at multiple pCO₂ levels.

<figure><img src="image_3.png" /><figcaption>**Fig. 1.** A schematic representation of the silicate weathering component of GEOCLIM in a single profile at steady-state. A rock particle leaves the unweathered bedrock with production rate *P*<sub>r</sub>, and transits vertically through a regolith of height *h*. Regolith production and physical erosion (*E*<sub>p</sub>) are equal at steady-state. As a particle transits upwards, some fraction of the primary phases (*x*) are chemically weathered (*W*), with the flux of dissolved Ca+Mg being *W* multiplied by the concentration of Ca+Mg in unweathered bedrock (*χ*<sub>CaMg</sub>). Details of the formulation for the silicate weathering component of GEOCLIM can be found in *Materials and Methods*.</figcaption></figure>

**Silicate Weathering Component.** The silicate weathering component of GEOCLIM calculates CO₂ consumption resulting from silicate weathering for subaerially exposed land. We assume that Ca and Mg are the only cations that consume CO₂ over geologic time-scales, such that each mole of Ca or Mg that is dissolved by silicate weathering consumes one mole of CO₂. While reverse weathering is another potential sink for Ca or Mg (27), its parameterization is unclear and it has been interpreted to be a relatively minor flux in the Cenozoic (28),

and we do not include it in our model. In previous versions of
the model, silicate weathering was a function of temperature
and runoff only, and all bedrock was assigned identical chem-
ical compositions (26). More recent versions of GEOCLIM
implement regolith development and soil shielding (Fig. 1),
which introduces a dependence on erosion rate (and therefore
topographic slope) (29). While this introduction of regolith
development into GEOCLIM is important for assessing the
impact of tropical arc-continent collisions on pCO₂, the rela-
tively high Ca+Mg concentration in arc rocks relative to other
lithologies must also be considered.

We therefore implement variable bedrock Ca+Mg concen-
tration into GEOCLIM (SI Appendix). The spatial distribu-
tion of lithologies is sourced from the Global Lithologic Map
(GLiM) (30) and is represented by 6 categories: metamorphic,
felsic, intermediate, mafic, carbonate, and siliciclastic sediment.
Each land pixel is assigned these lithologic categories at a res-
olution of 0.1° × 0.1°. The Ca+Mg concentrations of felsic,
intermediate, and mafic lithologies are assigned based on the
mean of data of these lithologic categories compiled in Earth-
Chem (www.earthchem.org/portal). Given that GLiM does
not distinguish ultramafic lithologies, such rocks are grouped
with mafic rocks. As a result, the Ca+Mg concentration is
likely an underestimate in regions of obducted ophiolites, such
that the estimated effect of these regions on changing steady-
state pCO₂ could be conservative (31). The weathering of
carbonate does not contribute to long-term CO₂ consumption
and its Ca+Mg concentration is ignored. The Ca+Mg concent-
trations of metamorphic and siliciclastic sediment lithologies
are more difficult to define, since their chemical composition
is strongly dependent on protolith composition and, in the
case of siliciclastic sediment, the degree of previous chemical
depletion. We explore a range of feasible Ca+Mg concentra-
tions for metamorphic rocks and siliciclastic sediment during
calibration of the silicate weathering component of GEOCLIM.

**Calibration.** The values of four parameters within the silicate weathering component that modify the dependence of silicate weathering on temperature, runoff, erosion, and regolith thickness are poorly constrained. Rather than prescribing single values, we select multiple values for each of these four parameters along with the Ca+Mg concentration of metamorphic and siliciclastic lithologies from within reasonable ranges (*SI Appendix*; Table S2). We then permute all possible combinations of these values for the six parameters, leading to 93,600 unique parameter combinations (i.e. permutations). For each combination, we compute spatially resolved long-term CO₂ consumption associated with Ca+Mg fluxes using present-day runoff, temperature, and slope. We sum computed CO₂ consumption over watersheds for which data-constrained estimates are available (1, 32), then calculate the coefficient of determination (*r*²) between computed and measured CO₂ consumption in each of these watersheds. After eliminating parameter combinations that result in low *r*², 573 parameter combinations remain (*SI Appendix*; Fig. S3). The resulting global CO₂ consumption of these filtered model runs all overlap with independently derived estimates of the global CO₂ degassing flux (33), as they should given that the long-term carbon cycle is in steady state (*SI Appendix*; Fig. S3).

**Climate Model Component.** Having calibrated the silicate weathering component of GEOCLIM, we use it to estimate the
<figure><img src="image_4.png" /><figcaption>**Fig. 2.** The emergence of the Southeast Asian islands (also referred to as the Maritime Continent in the climate science literature) from the mid-Miocene to present. Past shorelines at 5, 10, and 15 Ma are shown in A with associated land area summarized in B. A significant increase in area over the past 5 million years is coincident with cooling and the onset of Northern Hemisphere glaciation as reflected in the benthic oxygen isotope record (35) shown in C.</figcaption></figure>

decrease in steady-state $pCO_2$ associated with emergence of the SEAI. For the climate model component, we use temperature and runoff from a subset of the GFDL CM2.0 experiments (34) (*SI Appendix*). These experiments are well-suited for this analysis because all non-$CO_2$ forcings are held constant at values representative of pre-industrial conditions, allowing the effect of changing $pCO_2$ on climatology to be isolated. Furthermore, the experiments were run long enough for the final system to approximate steady-state.

## Paleoshorelines

To determine the position of paleoshorelines in the SEAI over the past 15 m.y., we use terrestrial and marine sedimentary deposits (Fig. 2; *SI Appendix*). The paleoshoreline data indicate that the Sunda-Banda Arc and New Guinea are primarily responsible for the increase in area since 15 Ma. Exhumation of the modern Sunda-Banda Arc is the result of ongoing arc-continent collision with the Australian Plate (37). Most of Sumatra and Java along with the non-volcanic islands of the Outer Banda Arc were elevated above sea level after

5 Ma (38). In New Guinea, emergence in the mid-Miocene is associated with collision between the Melanesian Arc and Australia's distal margin (39), which drove exhumation of the Irian-Marum-April Ophiolite Belt. Exhumation accelerated over the past 4 m.y. in the New Guinea Central Range due to slab-breakoff and buoyant uplift, and in eastern New Guinea due to jamming of the north-dipping subduction zone (39). We also include changes in areas of presently submerged continental shelves such as the Sunda Shelf that were previously exposed (*SI Appendix*; Fig. S7). These tectonic drivers and others throughout the region led to progressive emergence over the past 15 m.y. that accelerated following 5 Ma (Fig. 2B). This trend mirrors broad cooling over the Neogene that resulted in the initiation of Northern Hemisphere ice sheets (Fig. 2C).

We use GEOCLIM to estimate $pCO_2$ associated with the reconstructed subaerial extent of the SEAI at ca. 15, 10, and 5 Ma ("paleo-SEAI" scenarios; Fig. 3). Because we use a climate model forced with modern geography, the position of the tectonic blocks remain fixed. Although there has been motion of these tectonic blocks since 15 Ma, they have re-
<figure><img src="image_5.png" /><figcaption>**Fig. 3.** Steady-state *pCO*₂ estimates from GEOCLIM for the various scenarios discussed in the text. For each of the seven scenarios, each point represents an estimate from one of the 573 unique parameter combinations that most closely matched estimates of present-day CO₂ consumption in 80 watersheds around the world (SI Appendix). The box encloses the middle 50% of the *pCO*₂ estimates (i.e. the interquartile range), and the notch represents the median with its 95% confidence interval. The whiskers extend to the 2.5 and 97.5 percentile values. Glaciation thresholds (36) are shown on the x-axis.</figcaption></figure>

mained within tropical latitudes such that this fixed scenario is a good approximation of the paleogeography (SI Appendix; Fig. S10). We also test an end-member scenario, in which all islands associated with arc-continent collision in the region are removed (“removed SEAI” scenario; Fig. 3).

## pCO₂ Estimates

Using the 573 unique parameter combinations, the “paleo-SEAI” scenarios resulted in 526–678 ppm for 15 Ma, 457–516 ppm for 10 Ma, and 391–434 ppm *pCO*₂ for 5 Ma (Fig. 3). These results indicate a progressive decrease in *pCO*₂ over the Neogene associated with the emergence of the SEAI, and suggest that without this emergence, pre-industrial *pCO*₂ would have been ~526–678 ppm. These “paleo-SEAI” scenarios do not account for Neogene changes outside of the SEAI (e.g. changes in ocean/atmosphere circulation, volcanic degassing, and weathering fluxes elsewhere on Earth, discussed in *Alternative Mechanisms for Neogene Cooling*). Therefore, these results are not estimating *pCO*₂ at 15 Ma, but rather are quantifying *pCO*₂ change associated with emergence of the SEAI.

Proxy-based estimates of the magnitude and trajectory of *pCO*₂ change from the Miocene to the Pliocene are variable between techniques and associated assumptions underlying their interpretation (SI Appendix; Fig. S11). The *pCO*₂ values from the 5 Ma “paleo-SEAI” scenario overlap with many proxy-based estimates (40) as well as values that emerge from approaches that assimilate climate and ice sheet model output with benthic δ¹⁸O data (41, 42). The modeled *pCO*₂ values for 15 Ma resemble the higher end of proxy-based *pCO*₂ estimates for the early to mid-Miocene, indicating that the increase in subaerially exposed land area and tectonic topography of the SEAI is sufficient to explain long-term cooling of Earth’s climate over the Neogene. The *pCO*₂ threshold for Antarctic glaciation is estimated to be ~750 ppm with that for Northern Hemisphere glaciation being significantly lower at ~280 ppm

(36). These modeled values of decreasing *pCO*₂ associated with emergence of the SEAI are therefore consistent with the record of Neogene climate with Miocene ice sheets on Antarctica (43) followed by Northern Hemisphere ice sheets developing in the Pliocene (44) as *pCO*₂ subsequently decreased.

The results of our “paleo-SEAI” scenarios highlight the importance of the combination of topography, runoff, and lithology in setting Earth’s climate state. To independently explore the effect of the modern-day surface exposure of lower-relief basaltic lavas on steady-state *pCO*₂ (45), we replace mafic volcanics associated with the Deccan Traps, Ethiopian Traps, and Columbia River Basalts with the Ca+Mg concentration of bulk continental crust in GEOCLIM (Fig. 3). The resulting *pCO*₂ is ~300–500 ppm, indicating that the presence of mafic rocks in these igneous provinces affects steady-state *pCO*₂ as has been suggested to be important for Paleogene cooling (45). However, the higher 526–678 ppm values for the 15 Ma “paleo-SEAI” scenario illustrate that higher relief and a wet tropical climate significantly increase the efficiency of CO₂ consumption, especially when paired with high Ca+Mg lithologies. As such, arc-continent collisions in the tropics are likely more important for driving long-term changes in *pCO*₂ than the eruption of flood basalts (7, 46).

Previous work has estimated that the decrease in *pCO*₂ since 5 Ma associated with the emergence of the SEAI and enhanced silicate weathering is ~19 ppm (5), in which case their emergence would be a relatively minor contributor to Neogene cooling. This 19 ppm estimate was obtained using an equation that assumes a direct linear relationship between mean global temperature and changes in weathering-rate-weighted land area, scaled by a factor that is intended to account for the influence of both runoff and temperature. *pCO*₂ was then estimated from the calculated temperature using a simple energy balance equation. However, the relationship between mean global temperature (or *pCO*₂) and weathering-rate-weighted land area is not linear. Furthermore, this simple linear relationship ignores spatial variability in topography and climatology,
and only crudely accounts for spatial variability in lithology. In fact, the 19 ppm estimate is closer in magnitude to the decrease in pCO₂ that we estimate if mafic volcanics associated with the Deccan Traps (a relatively flat area outside of the warm and wet tropics) are replaced with the Ca+Mg concentration of bulk continental crust (22–70 ppm; Fig. 3). The significant difference in steady-state pCO₂ estimated between the “removed Deccan Traps” scenario and the “paleo-SEAI” scenarios (Fig. 3) demonstrates that considering changes in the spatial distribution of lithologies alone is not adequate for estimating changes in steady-state pCO₂. Instead, spatially varying topography and climatology significantly modulates silicate weathering rates, and must be accounted for when estimating pCO₂ change associated with paleogeographic change.

An important caveat for these estimates of pCO₂ is that our modeling is determining the climatology in the GFDL CM2.0 model at which steady-state is achieved – a climatology that has an associated pCO₂ value in the model. However, climate models are variable in their response to changes in pCO₂. One way to summarize this variability is through the equilibrium climate sensitivity value – the steady-state change in global mean surface air temperature associated with a doubling of pCO₂. A range of 1.5 to 4.5°C per pCO₂ doubling was proposed in the landmark Charney report (47) and this range was considered to be the credible interval (>66% likelihood) in the last IPCC report (48). Integrating constraints both from understanding of climate feedback processes and the climate record, a recent comprehensive review estimates the 66% probability range of climate sensitivity to be 2.6 to 3.9°C per pCO₂ doubling with a 5 to 95% range of 2.3 to 4.7°C per pCO₂ doubling (49). The equilibrium climate sensitivity associated with the CM2.0 climate models is 2.9°C per pCO₂ doubling, which falls within these ranges although these ranges remain broad. An alternative way to consider the results from our analysis would be that an estimate of 572 ppm (2× pre-industrial pCO₂) for the 15 Ma “paleo-SEAI” scenario is implying that Earth would be ~2.9°C warmer. If Earth’s climate sensitivity is at the higher end of the probable range and higher than in the CM2.0 model, as it is in some climate models, this same amount of Neogene cooling resulting from the emergence of the SEAI could have been driven by a smaller change in pCO₂.

## Alternative Mechanisms for Neogene Cooling

**Ocean/Atmosphere Circulation.** Some hypotheses to explain ice sheet growth over the Neogene invoke changes in ocean/atmosphere circulation including: further climatic isolation of Antarctica due to strengthening of the circumpolar current (11); increased atmospheric moisture in the Northern Hemisphere due to intensified thermohaline circulation following Panama Isthmus emergence (10); and cooling of North America resulting from a strengthened Walker Circulation associated with emergence of the SEAI (5). Such changes in ocean/atmosphere circulation are likely to modulate pCO₂ thresholds for glacial initiation and ice sheet growth (36). However, the prolonged time-scale of the cooling trend since 15 Ma (Fig. 2C) is most readily attributable to decreasing pCO₂ associated with evolving geological sources and sinks of carbon, modulated by the silicate weathering feedback (16, 50–53).

**Volcanic Degassing.** A decrease in volcanic degassing (12) has also been proposed as a driver for Neogene cooling. However, proxy-based estimates of the evolution of volcanic degassing fluxes throughout the Neogene are inconsistent with each other, such that not even the sign of the change in volcanic degassing over the past ~15 m.y. is without ambiguity (26). For example, it has both been estimated that the volcanic degassing flux was ~25% lower (54) and ~10% higher (55) at 15 Ma relative to the present day.

Our model framework provides an opportunity to estimate the decrease in volcanic degassing flux necessary to achieve the same change in pCO₂ predicted for the increase in global weatherability associated with the emergence of the SEAI over the past 15 m.y. If we use the parameter combination that had the highest r² between computed and measured CO₂ consumption in watersheds around the world during calibration (Calibration and SI Appendix; Fig. S4), GEOCLIM estimates a pre-industrial volcanic degassing flux of 4.1×10¹² mol/yr to balance the silicate weathering flux at 286 ppm pCO₂. If we then assume that this volcanic degassing flux did not change over the past 15 m.y., then GEOCLIM estimates that the increase in global weatherability associated with the emergence of the SEAI led to a change in pCO₂ of ~280 ppm (“increase in weatherability only” scenario in Fig. 4). If we instead assume that global weatherability did not change over the past 15 m.y., then we estimate that the volcanic degassing flux needs to have been ~13% greater at 15 Ma relative to the pre-industrial to drive the same ~280 ppm change in pCO₂ (“decrease in degassing only” scenario in Fig. 4). This ~13% value is higher than 10%, the highest current estimate for the volcanic degassing flux at 15 Ma relative to the present day (55).

However, changes in the volcanic degassing flux would have modulated changes in pCO₂ associated with changes in global weatherability. For example, some proxy-based approaches as well as some model-data assimilation approaches estimate that mid-Miocene pCO₂ was lower than 568 ppm (SI Appendix; Fig. S11). Take a scenario in which pCO₂ was 400 ppm at 15 Ma. If we assume that the pCO₂ decrease to the pre-industrial value of 286 ppm was driven by the increase in global weatherability associated with emergence of the SEAI in conjunction with an increase in volcanic degassing which counteracts cooling by increasing the flux of CO₂ to the atmosphere (“increase in weatherability and degassing” scenario in Fig. 4), the volcanic degassing flux would have had to have been ~7% smaller than the pre-industrial. More robust constraints on pCO₂ (SI Appendix; Fig. S11) and/or volcanic degassing rates over the past 20 m.y. are needed to constrain which of the “increase in weatherability only” or “increase in weatherability and degassing” scenarios (Fig. 4) is more representative of the mechanisms driving Neogene cooling.

**Himalayan Uplift.** Marine $^{87}$Sr/$^{86}$Sr has overall been increasing since ca. 35 Ma (56). The traditional explanation for this trend is that it reflects increased weathering of radiogenic (i.e. high $^{87}$Sr/$^{86}$Sr) silicate rocks (13, 57). Associated with this explanation is the proposal that increasing weathering of radiogenic silicate rocks in the Himalayas was the primary driver of Neogene cooling (58). It could be argued that increasing marine $^{87}$Sr/$^{86}$Sr is inconsistent with the hypothesis that increasing weathering of juvenile (i.e. low $^{87}$Sr/$^{86}$Sr) silicate rocks in the SEAI was an important driver of Neogene cooling.
<figure><img src="image_6.png" /><figcaption>**Fig. 4.** Weatherability curves for the modern and "paleo-SEAI" scenarios shown in Figure 3. The lower panel expands the lower *pCO₂* range (x-axis) of the upper panel. Details on how these curves were generated are described in *Materials and Methods*. Each of the 4 curves represent a different tectonic boundary condition (i.e. the reconstructed paleoshorelines of the SEAI; Fig. 2A) and therefore a different global weatherability. The curves show the resulting *pCO₂* for a given volcanic degassing flux such that the input flux is balanced by the silicate weathering output flux. Point B represents the pre-industrial, in which *pCO₂* is 286 ppm. The arrow from Point A₁ to B represents the "increase in weatherability only" scenario, in which global weatherability increases as the SEAI emerge, but the volcanic degassing flux does not change over the past 15 m.y. In this scenario, the *pCO₂* decreases from the value dictated by the 15 Ma "paleo-SEAI" weatherability curve (568 ppm). The arrow from Point A₂ to B instead represents the "decrease in degassing only" scenario, in which global weatherability remains the same as the pre-industrial, but the same change in *pCO₂* as the "increase in weatherability only" scenario is achieved by decreasing the volcanic degassing flux from a value ~13% greater than the pre-industrial. The arrow from Point A₃ to B represents the "increase in weatherability and degassing" scenario, in which a change in *pCO₂* from 400 ppm to 286 ppm is achieved by increasing both global weatherability from our 15 Ma tectonic boundary condition as well as the volcanic degassing flux from a value ~7% smaller than the pre-industrial flux.</figcaption></figure>

However, the globally averaged ratio of silicate weathering fluxes from radiogenic cratonic rocks versus juvenile arc lithologies can be at least partially decoupled from marine $^{87}$Sr/$^{86}$Sr via the regional weathering of isotopically unique lithologies. For example, in addition to highly radiogenic granites and gneisses (57), unusually radiogenic carbonates are abundant in Himalayan strata, and it is estimated that ~75% of Sr coming from the Himalayas can be attributed to carbonate rather than silicate weathering (59–61). As such, there are challenges in interpreting the marine $^{87}$Sr/$^{86}$Sr record as a direct proxy for silicate weathering fluxes. Nevertheless, steadily increasing marine $^{87}$Sr/$^{86}$Sr is interrupted ca. 16 Ma when the slope of the $^{87}$Sr/$^{86}$Sr curve decreases (56). This decrease in slope has been attributed to coincident exhumation of relatively low $^{87}$Sr/$^{86}$Sr Outer Lesser Himalaya carbonates (62, 63), but

could also be at least partially driven by the emergence of low $^{87}$Sr/$^{86}$Sr lithologies in the SEAI during arc-continent collision. Increasing seawater Mg/Ca since ca. 15 Ma (64) is consistent with an increasing proportion of the global silicate weathering flux being derived from mafic and ultramafic sources.

Himalayan uplift would have affected geological carbon sinks, either via increased weathering of silicate rocks (58) or enhanced burial of organic matter in the Bengal Fan (14). Increased weathering of the emerging SEAI would have occurred in tandem with such changes in the Himalaya, such that the effects of these paleogeographic changes on geochemical proxy records, like marine $^{87}$Sr/$^{86}$Sr, become difficult to disentangle. In addition, given the large uncertainty associated with changes in regional climatology across Asia due to Himalayan orogeny, developing quantitative estimates of the evolution of global silicate weathering fluxes associated with Himalayan orogeny remains a major challenge.

## The Geologic Carbon Cycle

If geological carbon sources remain approximately constant, global alkalinity delivery from silicate weathering needs to be approximately constant as well to keep the long-term carbon cycle in steady-state (16). Enhanced silicate weathering in a region such as the SEAI is compensated by a decrease in silicate weathering elsewhere. Global alkalinity delivery from silicate weathering does not change, but occurs more efficiently and thereby at lower *pCO₂*. Given that carbonate weathering is disconnected from the long-term carbon-cycle mass balance, changes in carbonate accumulation through time (65) could be driven by changes in carbonate weathering.

The long-term carbon-cycle mass balance can be perturbed via mechanisms that are disconnected from changes in volcanic degassing and silicate weathering rates. For example, sulphide oxidation coupled to carbonate dissolution could act as a source of CO₂ on million year time-scales (66). Similarly, the weathering of sedimentary organic matter could serve as a source of CO₂ (67). On the other hand, enhanced burial of organic matter enabled by higher sediment and nutrient delivery could be an important sink of CO₂, as has been suggested in the Bengal Fan (14) and Taiwan (68). The fluxes of CO₂ represented by these processes are not accounted for in our model framework, and could have been affected by emergence of the SEAI and/or Himalayan orogeny. *pCO₂* changes that result from these processes would be superimposed on *pCO₂* changes associated with evolving silicate weathering fluxes. However, our coupled weathering-climate model indicates that the *pCO₂* change associated with increased global weatherability driven by emergence of the SEAI is sufficient to explain the majority of Neogene cooling (Fig. 3). Without this emergence, *pCO₂* would have remained above the level necessary for the growth of Northern Hemisphere ice sheets.

## Conclusions

Coupled geological constraints and modeling experiments demonstrate that the SEAI have been a growing hot spot for carbon sequestration due to silicate weathering from the Miocene to present. Changes in volcanic degassing and paleogeography elsewhere on Earth, particularly in the Himalaya and Central America, would have also affected geological carbon sources and sinks. Yet, not only does the history of
emergence of the SEAI coincide with Neogene cooling and the onset of Northern Hemisphere glaciation, but our coupled weathering-climate model also indicates that the associated steady-state pCO₂ change is sufficient to explain much of this cooling. These results highlight that the Earth's climate state is particularly sensitive to changes in tropical geography.

## Materials and Methods

The code for the GEOCLIM model used in this study can be found at: https://github.com/piermafrost/GEOCLIM-dynsoil-steady-state/releases/tag/v1.0. The code that generated the inputs and analyzed the output of the GEOCLIM model can be found at: https://github.com/Swanson-Hysell-Group/GEOCLIM_Modern.

**GEOCLIM Silicate Weathering Component.** The silicate weathering component of the GEOCLIM model has been modified from the previously published version (20). The new component implements the model of Gabet and Mudd (2009) (17) for the development of a chemically weathered profile. We refer to this chemically weathered profile as regolith where the base of the regolith is unweathered bedrock. In the model of Gabet and Mudd (2009), material enters the regolith and leaves either as a solute through chemical weathering of the material during its travel from the bedrock towards the surface, or as a physically weathered particle once it reaches the top. We use the DynSoil implementation of the Gabet and Mudd (2009) model, which integrates chemical weathering within the regolith (18). The transient time-varying version of this regolith model is described by three equations:

$$ \frac{dh}{dt} = P_r - E_p \quad [1] $$

$$ \frac{\partial x}{\partial t} = -P_r \frac{\partial x}{\partial z} - K \tau^\sigma x \quad [2] $$

$$ \frac{\partial \tau}{\partial t} = -P_r \frac{\partial \tau}{\partial z} + 1 \quad [3] $$

Equation 1 is a statement of material conservation, where h is the total height of the regolith (m), t is the model time (yr), P_r is the regolith production rate (m/yr), and E_p is the physical erosion rate (m/yr). Equation 2 describes how the residual fraction of weatherable phases (x, unitless) changes as a function of time (t, yr) and depth (z, m). Kτ^σ is the dissolution rate constant, which depends on the local climate (captured by K, yr⁻¹⁻^σ) and the time that a given rock particle has spent in the regolith (τ, yr) to some power σ (unitless) which implements a time-dependence. Equation 3 describes how the time that a given rock particle has spent in the regolith changes as time in the model progresses.

The net weathering rate in the regolith column (W, m/yr) can then be calculated with:

$$ W = \int_{0}^{h} K \tau^{\sigma} x \, dz \quad [4] $$

The regolith production rate can be expressed as the product of the optimal production rate (P₀) and a soil production function (f(h)):

$$ P_r = P_0 f(h) \quad [5] $$

$$ P_0 = k_{rp} q e^{-\frac{E_a}{R} \left(\frac{1}{T} - \frac{1}{T_0}\right)} \quad [6] $$

$$ f(h) = e^{-\frac{h}{d_0}} \quad [7] $$

P₀ is the 'optimal' regolith production rate (m/yr), which is defined to be the regolith production rate when there is no overlying regolith. In Equation 6, where k_{rp} is a proportionality constant (unitless), q is the runoff (m/yr), E_a is the activation energy (J/K/mol), R is the ideal gas constant (J/mol), T is the temperature (K), and T₀ is the reference temperature (K), we parameterize the 'optimal' regolith production rate (69). f(h) is the soil production function

(unitless), which describes how regolith production decreases as the thickness of the regolith increases. It takes an exponential form as is commonly applied in the literature (17). In Equation 7, d₀ is a reference regolith thickness (m) (70).

Our implementation of the erosion rate is parameterized based on runoff and slope (s; m/m):

$$ E_p = k_e q^m s^n \quad [8] $$

k_e is a proportionality constant ((m/yr)¹⁻²) and m and n are adjustable exponents that are kept as 0.5 and 1 (29). This formulation is directly inspired by the stream power law (71). This formulation and these exponent values are supported by compilations, but variability in the proportionality constant is difficult to capture at a global scale (72).

The K in the dissolution rate constant in Equation 2 describes the dependence of the chemical weathering on climate:

$$ K = k_d \left(1 - e^{-k_w q}\right) e^{-\frac{E_a}{R} \left(\frac{1}{T} - \frac{1}{T_0}\right)} \quad [9] $$

Equation 9 is an empirical simplification of mineral dissolution rates derived from kinetic theory and laboratory experiments (18), where k_d is a proportionality constant that modifies the dependence of dissolution rate on runoff and temperature (yr⁻¹⁻^σ), and k_w is a proportionality constant that modifies the dependence of dissolution rate on runoff (yr/m).

In this study, we are interested in obtaining the steady-state solution rather than the transient time-varying solution. The steady-state solution for DynSoil can be calculated analytically by setting the time derivatives equal to zero resulting in the following set of equations:

$$ h = \max \left( 0, d_0 \log \left( \frac{P_0}{E_p} \right) \right) \quad [10] $$

$$ x(z) = \exp \left( - \frac{K}{\sigma + 1} \left( \frac{z}{E_p} \right)^{\sigma + 1} \right) \quad [11] $$

$$ W = E_p(1 - x(h)) = E_p(1 - x_s) \quad [12] $$

x(z) is the abundance profile of primary phases inside the regolith, varying with height upward from the base of the regolith as shown in Figure 1. Setting z equal to the regolith thickness (h) gives x_s which is the proportion of primary phases remaining at the top of the regolith column.

**Weatherability Curves.** To create the 15 Ma “paleo-SEAI” curve shown in Figure 4, we use the reconstructed paleoshorelines of the SEAI at 15 Ma (Fig. 2A). We then select the parameter combination that had the highest r² between computed and measured CO₂ consumption in watersheds around the world during calibration (Calibration and SI Appendix; Fig. S4), and fix pCO₂ at the 3 pCO₂ levels at which the GFDL CM2.0 climate model experiments were computed (SI Appendix). We then run GEOCLIM at each of these pCO₂ levels until steady state is achieved (i.e. until the volcanic degassing flux is equal to the silicate weathering flux). We then repeat this process for the 10 and 5 Ma “paleo-SEAI” paleoshorelines and the present day shorelines to generate the 3 other weatherability curves. Each estimated pCO₂ in Figure 3 is the result of underlying weatherability curves that change with the different chemical weathering parameters.

**SI Appendix.** A detailed description of the implementation of lithology into the silicate weathering component of GEOCLIM, the calibration of the silicate weathering component of GEOCLIM, the GFDL CM2.0 climate model, and the paleoshoreline reconstructions can be found in the SI Appendix.

**ACKNOWLEDGMENTS.** Collaborative research between N.L.S.-H. and Y.G. was initially supported by a grant from the France-Berkeley Fund. N.L.S.-H. and F.A.M. gratefully acknowledge support through NSF FRES grants #1926001 and #1925990. We thank Alec Brenner, Sam Lo Bianco, Mariana Lin, and Judy Pu for their data compilation contributions to the paleoshoreline reconstructions.
1. J Gaillardet, B Dupré, P Louvat, CJ Allègre, Global silicate weathering and CO₂ consumption rates deduced from the chemistry of large rivers. *Chem. Geol.* **159**, 3–30 (1999).

2. J Hartmann, N Jansen, HH Dür, S Kempe, P Köhler, Global CO₂-consumption by chemical weathering: What is the contribution of highly active weathering regions? *Glob. Planet. Chang.* **69**, 185–194 (2009).

3. JD Millman, KL Farnsworth, *River Discharge to the Coastal Ocean: A Global Synthesis*. (Cambridge University Press), (2013).

4. J Hartmann, N Moosdorf, R Lauerwald, M Hinderer, AJ West, Global chemical weathering and associated P-release - the role of lithology, temperature and soil properties. *Chem. Geol.* **363**, 145–163 (2014).

5. P Molnar, TW Cronin, Growth of the Maritime Continent and its possible contribution to recurring ice ages. *Paleoceanography* **30**, 196–225 (2015).

6. R Hall, Southeast Asia: New views of the geology of the Malay Archipelago. *Annu. Rev. Earth Planet. Sci.* **45**, 331–358 (2017).

7. FA Macdonald, NL Swanson-Hysell, Y Park, L Lisiecki, O Jagoutz, Arc-continent collisions in the tropics set Earth's climate state. *Science* **364**, 181–184 (2019).

8. NJ Shackleton, et al., Oxygen isotope calibration of the onset of ice-rafting and history of glaciation in the North Atlantic region. *Nature* **307**, 620–623 (1984).

9. J Zachos, M Pagani, L Sloan, E Thomas, K Billups, Trends, rhythms, and aberrations in global climate 65 ma to present. *Science* **292**, 686–693 (2001).

10. GH Haug, R Tiedemann, Effect of the formation of the Isthmus of Panama on Atlantic Ocean thermohaline circulation. *Nature* **393**, 673–676 (1998).

11. AE Shevenell, Middle Miocene Southern Ocean cooling and Antarctic cryosphere expansion. *Science* **305**, 1766–1770 (2004).

12. RA Berner, AC Lasaga, RM Garrels, The carbonate-silicate geochemical cycle and its effect on atmospheric carbon dioxide over the past 100 million years. *Am. J. Sci.* **283**, 641–683 (1983).

13. ME Raymo, WF Ruddiman, PN Froelich, Influence of late Cenozoic mountain building on ocean geochemical cycles. *Geology* **16**, 649–653 (1988).

14. V Galy, et al., Efficient organic carbon burial in the Bengal fan sustained by the Himalayan erosional system. *Nature* **450**, 407–410 (2007).

15. LR Kump, SL Brantley, MA Arthur, Chemical weathering, atmospheric CO₂, and climate. *Annu. Rev. Earth Planet. Sci.* **28**, 611–667 (2000).

16. LR Kump, MA Arthur, *Global Chemical Erosion during the Cenozoic: Weatherability Balances the Budgets*, ed. WF Ruddiman. (Springer US, Boston, MA), pp. 399–426 (1997).

17. EJ Gabet, SM Mudd, A theoretical model coupling chemical weathering rates with denudation rates. *Geology* **37**, 151–154 (2009).

18. AJ West, Thickness of the chemical weathering zone and implications for erosional and climatic drivers of weathering and for carbon-cycle feedbacks. *Geology* **40**, 811–814 (2012).

19. K Maher, CP Chamberlain, Hydrologic regulation of chemical weathering and the geologic carbon cycle. *Science* **343**, 1502–1504 (2014).

20. Y Goddéris, et al., Onset and ending of the late Palaeozoic ice age triggered by tectonically paced rock weathering. *Nat. Geosci.* **10**, 382–386 (2017).

21. AC Lasaga, JM Soler, J Ganor, TE Burch, KL Nagy, Chemical weathering rate laws and global geochemical cycles. *Geochimica et Cosmochimica Acta* **58**, 2361–2386 (1994).

22. C Dessert, B Dupré, J Gaillardet, LM François, CJ Allègre, Basalt weathering laws and the impact of basalt weathering on the global carbon cycle. *Chem. Geol.* **202**, 257–273 (2003).

23. O Jagoutz, FA Macdonald, L Royden, Low-latitude arc-continent collision as a driver for global cooling. *Proc. Natl. Acad. Sci.* **113**, 4935–4940 (2016).

24. NL Swanson-Hysell, FA Macdonald, Tropical weathering of the Taconic orogeny as a driver for Ordovician cooling. *Geology* **45**, 719–722 (2017).

25. Y Park, et al., The lead-up to the Sturtian Snowball Earth: Neoproterozoic chemostratigraphy time-calibrated by the Tambien Group of Ethiopia. *GSA Bull.* **132**, 1119–1149 (2020).

26. Y Goddéris, Y Donnadieu, A sink- or a source-driven solid cycle at the geological timescale? Relative importance of palaeogeography versus solid Earth degassing rate in the Phanerozoic climatic evolution. *Geol. Mag.* **156**, 355–365 (2017).

27. P Michalopoulos, RC Aller, Rapid clay mineral formation in Amazon Delta sediments: Re-verse weathering and oceanic elemental cycles. *Science* **270**, 614–617 (1995).

28. TT Isson, NJ Planavsky, Reverse weathering as a long-term stabilizer of marine pH and planetary climate. *Nature* **560**, 471–475 (2018).

29. P Maffre, et al., Mountain ranges, climate and weathering, do orogens strengthen or weaken the silicate weathering carbon sink? *Earth Planet. Sci. Lett.* **493**, 174–185 (2018).

30. J Hartmann, N Moosdorf, The new global lithological DNA database GLIM: A representation of rock properties at the Earth surface. *Geochem. Geophys. Geosystems* **13**, 1–37 (2012).

31. H Schopka, L Derry, C Arcilla, Chemical weathering, river geochemistry and atmospheric carbon fluxes from volcanic and ultramafic regions on Luzon Island, the Philippines. *Geochimica et Cosmochimica Acta* **75**, 978–1002 (2011).

32. JS Moquet, et al., Temporal variability and annual budget of inorganic dissolved matter in Andean Pacific Rivers located along a climate gradient from northern Ecuador to southern Peru. *Comptes Rendus Geosci.* **350**, 76–87 (2018).

33. T Gerlach, Volcanic versus anthropogenic carbon dioxide. *Eos, Transactions Am. Geophys. Union* **92**, 201–202 (2011).

34. TL Delworth, et al., GFDL's CM2 global coupled climate models. part i: Formulation and simulation characteristics. *J. Clim.* **19**, 643–674 (2006).

35. JC Zachos, GR Dickens, RE Zeebe, An early Cenozoic perspective on greenhouse warming and carbon-cycle dynamics. *Nature* **451**, 279–283 (2008).

36. RM DeConto, et al., Thresholds for Cenozoic bipolar glaciation. *Nature* **455**, 652–656 (2008).

37. R Harris, Rise and fall of the Eastern Great Indonesian arc recorded by the assembly, dispersion and accretion of the Banda Terrane, Timor. *Gondwana Res.* **10**, 207–231 (2006).

38. R Hall, The palaeogeography of Sundaland and Wallacea since the Late Jurassic. *J. Limnol.* **72**, 1–17 (2013).

39. M Cloos, et al., Collisional delamination in New Guinea: The geotectonics of subducting slab breakoff. *Geol. Soc. Am. Special Pap.* **1–51** (2005).

40. GL Foster, DL Royer, DJ Lunt, Future climate forcing potentially without precedent in the last

420 million years. *Nat. Commun.* **8**, 1–8 (2017).

41. RSW van de Wal, B de Boer, LJ Lourens, P Köhler, R Bintanja, Reconstruction of a continuous high-resolution CO₂ record over the past 20 million years. *Clim. Past* **7**, 1459–1469 (2011).

42. CJ Berends, B de Boer, RSW van de Wal, Reconstructing the evolution of ice sheets, sea level and atmospheric CO₂ during the past 3.6 million years. (2020).

43. DE Sugden, et al., Preservation of Miocene glacier ice in East Antarctica. *Nature* **376**, 412–414 (1995).

44. GH Haug, et al., North Pacific seasonality and the glaciation of North America 2.7 million years ago. *Nature* **433**, 821–825 (2005).

45. DV Kent, G Muttoni, Modulation of Late Cretaceous and Cenozoic climate by variable drawdown of atmospheric pCO₂ from weathering of basaltic provinces on continents drifting through the equatorial humid belt. *Clim. Past* **9**, 525–546 (2013).

46. Y Park, NL Swanson-Hysell, FA Macdonald, L Lisiecki, Evaluating the relationship between the area and latitude of large igneous provinces and Earth's long-term climate state. *Earth-ArXiv* (2019).

47. JG Charney, et al., *Carbon dioxide and climate: a scientific assessment*. (National Academy of Sciences, Washington, DC), (1979).

48. T Stocker, et al., *Climate Change 2013: The Physical Science Basis. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change*. (Cambridge University Press), p. 1535 (2013).

49. S Sherwood, et al., An assessment of Earth's climate sensitivity using multiple lines of evidence. *Rev. Geophys.* (2020).

50. JCG Walker, PB Hays, JF Kasting, A negative feedback mechanism for the long-term stabilization of Earth's surface temperature. *J. Geophys. Res.* **86**, 9776–9782 (1981).

51. ME Raymo, GEochemical evidence supporting T. C. Chamberlin's theory of glaciation. *Geology* **19**, 344–347 (1991).

52. RA Berner, K Caldeira, The need for mass balance and feedback in the geochemical carbon cycle. *Geology* **25**, 955–956 (1997).

53. RA Berner, GEOCARB III: A revised model of atmospheric CO₂ over Phanerozoic time. *Am. J. Sci.* **301**, 182–204 (2001).

54. JP Cogné, E Humler, Trends and rhythms in global seafloor generation rate. *Geochem. Geophys. Geosystems* **7**, 1–17 (2006).

55. DG Van Der Meer, et al., Plate tectonic controls on atmospheric CO₂ levels since the Triassic. *Proc. Natl. Acad. Sci.* **111**, 4380–4385 (2014).

56. J McArthur, R Howarth, G Shields, Strontium isotope stratigraphy. *The Geol. Time Scale*, 127–144 (2012).

57. JM Edmond, Himalayan tectonics, weathering processes, and the strontium isotope record in marine limestones. *Science* **258**, 1594–1597 (1992).

58. ME Raymo, WF Ruddiman, Tectonic forcing of late Cenozoic climate. *Nature* **359**, 117–122 (1992).

59. AD Jacobson, JD Blum, LM Walter, Reconciling the elemental and Sr isotope composition of Himalayan weathering fluxes: insights from the carbonate geochemistry of stream waters. *Geochimica et Cosmochimica Acta* **66**, 3417–3429 (2002).

60. J Quade, N English, PG DeCelles, Silicate versus carbonate weathering in the Himalaya: a comparison of the Arun and Seti River watersheds. *Chem. Geol.* **202**, 275–296 (2003).

61. L Oliver, et al., Silicate weathering rates decoupled from the ⁸⁷Sr/⁸⁷Sr ratio of the dissolved load during Himalayan erosion. *Chem. Geol.* **201**, 119–139 (2003).

62. PM Myrow, et al., Neogene marine isotopic evolution and the erosion of Lesser Himalayan strata: Implications for Cenozoic tectonic history. *Earth Planet. Sci. Lett.* **417**, 142–150 (1995).

63. CL Colleps, et al., Zircon (U-Th)/He thermochronometric constraints on Himalayan thrust belt exhumation, bedrock weathering, and Cenozoic seawater chemistry. *Geochem. Geophys. Geosystems* **19**, 257–271 (2018).

64. J Higgins, D Schrag, Records of Neogene seawater chemistry and diagenesis in deep-sea carbonate sediments and pore fluids. *Earth Planet. Sci. Lett.* **357–358**, 386–396 (2012).

65. W Si, Y Rosenthal, Reduced continental weathering and marine calcification linked to late Neogene decline in atmospheric CO₂. *Nat. Geosci.* **12**, 833–838 (2019).

66. MA Torres, AJ West, G Li, Sulphide oxidation and carbonate dissolution as a source of CO₂ over geological timescales. *Nature* **507**, 346–349 (2014).

67. RG Hilton, J Gaillardet, D Calmels, JL Birck, Geological respiration of a mountain belt revealed by the trace element rhodium. *Earth Planet. Sci. Lett.* **403**, 27–36 (2014).

68. SJ Kao, et al., Preservation of terrestrial organic carbon in marine sediments offshore Taiwan: mountain building and atmospheric carbon dioxide sequestration. *Earth Surf. Dyn.* **2**, 127–139 (2014).

69. S Carretier, Y Goddéris, T Delannoy, D Rouby, Mean bedrock-to-saprolite conversion and erosion rates during mountain growth and decline. *Geomorphology* **209**, 39–52 (2014).

70. AM Heimsath, WE Dietrich, K Nishiizumi, RC Finkel, The soil production function and landscape equilibrium. *Nature* **388**, 358–361 (1997).

71. P Davy, A Crave, Upscaling local-scale transport processes in large-scale relief dynamics. *Phys. Chem. Earth, Part A: Solid Earth* **25**, 533–541 (2000).

72. D Lague, The stream power river incision model: evidence, theory and beyond. *Earth Surf. Process. Landforms* **39**, 38–61 (2013).
# Supporting Information for “Emergence of the Southeast Asian islands as a driver for Neogene cooling”

Yuem Park<fnref n="1" />, Pierre Maffre<fnref n="1" />, Yves Goddéris<fnref n="3" />, Francis A. Macdonald<fnref n="2" />, Eliel S. C. Anttila<fnref n="2" />, Nicholas L. Swanson-Hysell<fnref n="1" />

1

2

3

4

5

6

These supplementary information materials provide details on the model framework used in this study. The code for the GEOCLIM model used in this study can be found at:

https://github.com/piermafrost/GEOCLIM-dynsoil-steady-state/releases/tag/v1.0. The code that generated the inputs and analyzed the output of the GEOCLIM model can be found at:

https://github.com/Swanson-Hysell-Group/2020_Southeast_Asian_Islands or

https://doi.org/10.5281/zenodo.4021653.

7

# Implementation of Lithology

8

9

10

11

12

13

14

15

We restricted the calculation of the weathering fluxes to the flux of dissolved Ca+Mg originating from continental silicate weathering, as they are essential for the long-term consumption of atmospheric CO₂ through silicate weathering. To calculate this flux, we need to assign the concentration of Ca and Mg ($\chi_{CaMg}$) within the unweathered bedrock (Fig. 1 in the main text).

Previous implementations of GEOCLIM have used a “diffuse lithology” where all exposed land is assigned the composition of bulk upper continental crust. These previous studies assumed that the weathering rates of all rocks was the same for each continental grid element, provided that the grid element was submitted to the same climatic conditions. Given that the hypothesis we seek to

<footnote n="1">¹ Department of Earth and Planetary Science, University of California, Berkeley, CA, USA</footnote>

<footnote n="2">² Department of Earth Science, University of California, Santa Barbara, CA, USA</footnote>

<footnote n="3">³ Géosciences Environnement Toulouse, CNRS-Université Paul Sabatier - IRD, Toulouse, France</footnote>
<table><caption>**Table S1.** Grouping of 16 lithologic categories in GLIM Hartmann and Moosdorf (2012) to 6 broader categories for GEOCLIM.</caption><thead><tr><th>GLIM ID</th><th>GLIM code</th><th>GLIM classification</th><th>GEOCLIM ID</th><th>GEOCLIM classification</th></tr></thead><tbody><tr><td>1</td><td>su</td><td>unconsolidated sediments</td><td>6</td><td>sediments</td></tr><tr><td>2</td><td>vb</td><td>basic volcanic rocks</td><td>4</td><td>mafics</td></tr><tr><td>3</td><td>ss</td><td>siliciclastic sedimentary rocks</td><td>6</td><td>sediments</td></tr><tr><td>4</td><td>pb</td><td>basic plutonic rocks</td><td>4</td><td>mafics</td></tr><tr><td>5</td><td>sm</td><td>mixed sedimentary rocks</td><td>6</td><td>sediments</td></tr><tr><td>6</td><td>sc</td><td>carbonate sedimentary rocks</td><td>5</td><td>carbonates</td></tr><tr><td>7</td><td>va</td><td>acid volcanic rocks</td><td>2</td><td>felsics</td></tr><tr><td>8</td><td>mt</td><td>metamorphics</td><td>1</td><td>metamorphics</td></tr><tr><td>9</td><td>pa</td><td>acid plutonic rocks</td><td>2</td><td>felsics</td></tr><tr><td>10</td><td>vi</td><td>intermediate volcanic rocks</td><td>3</td><td>intermediates</td></tr><tr><td>11</td><td>wb</td><td>water bodies</td><td>0</td><td>water/ice</td></tr><tr><td>12</td><td>py</td><td>pyroclastics</td><td>2</td><td>felsics</td></tr><tr><td>13</td><td>pi</td><td>intermediate plutonic rocks</td><td>3</td><td>intermediates</td></tr><tr><td>14</td><td>ev</td><td>evaporites</td><td>5</td><td>carbonates</td></tr><tr><td>15</td><td>nd</td><td>no data</td><td>1</td><td>metamorphics</td></tr><tr><td>16</td><td>ig</td><td>ice and glaciers</td><td>0</td><td>water/ice</td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>**Figure S1.** Distribution of lithologies at 0.1° × 0.1° resolution used in GEOCLIM modified from Hartmann and Moosdorf (2012).</figcaption></figure>

16 test involves the varying concentration of cations in different lithologies, we instead implement a

17 more realistic lithologically-resolved version of the model.

18 The spatial distribution of lithologies is sourced from the Global Lithologic Map (GLIM) of
Hartmann and Moosdorf (2012). The raw data takes the form of polygon vectors, where each
19
polygon is assigned one of 16 lithologic categories. We first group these 16 categories into 6
20
broader categories (metamorphic, felsic, intermediate, mafic, carbonate, and siliciclastic sediment;
21
Table S1). Note that the siliciclastic sediment lithologic category also includes sedimentary
22
sequences in which any carbonate is identified but is not the dominant lithology Hartmann and
23
Moosdorf (2012). We then rasterize the polygon vectors to 0.1° × 0.1° resolution, where each
24
pixel is assigned the lithologic category of the polygon that covers the greatest area in that pixel
25
(i.e. the ‘mode lithology’; Fig. S1). To improve the computing time of GEOCLIM, we decrease
26
the resolution of the raster to 0.5° × 0.5°. To do so, a 3-dimensional 720 × 360 × 7 matrix is
27
created, in which the fraction of each 0.5° × 0.5° pixel covered by each of the 6 lithologic
28
categories (or water/ice) is captured by the extra dimension. In this way, we calculate an
29
area-weighted mean Ca+Mg concentration of the surface in each 0.5° × 0.5° pixel.
30

The Ca+Mg concentrations of felsic (1,521 mol/m³), intermediate (4,759 mol/m³), and mafic
31
(10,317 mol/m³) lithologies are assigned based on the mean of data compiled from EarthChem
32
(www.earthchem.org/portal; calculations made in the code within the repository). For
33
metamorphic and siliciclastic lithologies, we explore a range of feasible Ca+Mg concentrations
34
during calibration of the silicate weathering component of GEOCLIM (GEOCLIM Calibration).
35
This approach makes the simplifying assumption that all pixels of a given lithologic category
36
share the same Ca+Mg concentration.
37

We also explore the sensitivity of the area-weighted mean Ca+Mg concentration of the surface in
38
each 0.5° × 0.5° pixel to the resolution of the underlying GLiM raster (Fig. S2). We find that the
39
difference in the area-weighted mean Ca+Mg concentration of the surface when using the 0.1° ×
40
0.1° GLiM raster versus the 0.05° × 0.05° GLiM raster can be well explained by a tight Gaussian
41
distribution (σ<200 mol/m³) about a mean of ~0 mol/m³. Therefore, using the lower-resolution
42
GLiM raster does not overall bias our results to higher or lower area-weighted mean Ca+Mg
43
concentrations of the surface, and only introduces some noise that has a magnitude that is
44
<figure><img src="image_8.png" /><figcaption>**Figure S2.** Histograms showing the difference in the area-weighted mean Ca+Mg concentration of the surface in each 0.5° × 0.5° pixel that results from starting with an 0.1° × 0.1° GLiM raster versus an 0.05° × 0.05° GLiM raster. The left panel is calculated over all land pixels, whereas the right panel is calculated over land pixels in the SEAIs only. Both panels use the parameter values that produced the highest r² during calibration: metamorphicCa+Mg = 2500 mol/m³ and sedimentCa+Mg = 2000 mol/m³.</figcaption></figure>

significantly smaller than the most Ca+Mg-poor lithologic category (i.e. felsics, at 1,521 mol/m³).

## GEOCLIM Calibration

Experimental determinations of the activation energy (E_a; Equations 6 and 9) associated with the weathering of silicate minerals are variable (Brantley, 2003). However, multiple efforts to invert for E_a in basaltic watersheds with varying temperature have yielded values (41.6 ± 3.2 kJ/mol in Li et al., 2016; 42.3 kJ/mol in Dessert et al., 2001) that are consistent with the lower end of activation energies of Ca+Mg bearing minerals in laboratory experiments such as that for diopside (40.5 ± 1.7 kJ/mol; Knauss et al., 1993) and for labradorite (42.1 kJ/mol; Carroll and Knauss, 2005). We use the value of 42 kJ/mol in our model runs. While we implement lithology-dependent Ca+Mg concentration, our implementation does not include lithology-dependent kinetics of mineral dissolution. Relative to felsic lithologies, mafic lithologies contain a higher concentration of minerals with faster dissolution kinetics (e.g. plagioclase).
<figure><img src="image_9.png" /><figcaption>**Figure S3.** A) Modeled global CO₂ consumption vs. the coefficient of determination (r²) between modeled and data-constrained CO₂ consumption in each of the watersheds. Each point represents model output using one of the 93,600 parameter combinations (Table S2). B) Same as A, but zoomed to the plotting space with positive coefficient of determination values. The black line represents the full range of estimates of the present-day non-anthropogenic global CO₂ emission rate compiled in Gerlach (2011). The black box represents the range of these estimates preferred by Gerlach (2011).</figcaption></figure>

However, Ca+Mg from both felsic and mafic lithologies is predominantly sourced from minerals with faster dissolution kinetics, and we therefore use the same chemical weathering formulation (including the same activation energy) across lithologies.

Within the equation that governs the physical erosion rate in GEOCLIM (Equation 8), the values of the exponents (m = 0.5 and n = 1) are supported by compilations (Lague, 2013).

However, there is uncertainty in the value of the proportionality constant (kₑ). Using the simplification of a uniform value, we therefore tune this value such that the total erosion flux in the model under present-day slope and runoff conditions (see below) matches the total erosion flux estimated in Milliman and Farnsworth (2013) (20 × 10¹² kg/yr). Assuming that the density of eroded materials is 2500 kg/m³, our tuned kₑ value is 0.0029110 m¹⁻m/yr¹⁻m.

Within the equations that govern the chemical weathering rate in GEOCLIM, we identify the less constrained parameters: the proportionality constant that modifies the dependence of dissolution rate on runoff and temperature (k₉; Equation 9), the proportionality constant that modifies the dependence of dissolution rate on runoff only (kₛ; Equation 9), the power constant
<figure><img src="image_10.png" /><figcaption>**Figure S4.** Modeled vs. data-constrained CO₂ consumption in watersheds around the world. **A)** Each point represents a single watershed, and the y-value of the point shows the mean value of the 573 parameter combinations that produce individual watershed CO₂ consumption fluxes that approximate those estimated in the literature for the present-day (i.e. the orange points in Figure S3). Whiskers extend to the minimum and maximum modeled watershed CO₂ consumption fluxes for the 573 parameter combinations. Watersheds within the compilation of Gaillardet et al. (1999) that are in the SEAIs (Fly, Kikori, Purari, and Sepik watersheds, all in New Guinea) are indicated in green. **B)** Same as A, but only showing the parameter combination that produced that highest r². The chemical weathering and regolith thickness maps of this parameter combination are shown in Figure S5.</figcaption></figure>

that modifies the dependence of dissolution rate on the time that a rock particle has spent in the regolith ($\sigma$; Equation 4), and the proportionality constant that modifies the dependence of regolith production on runoff and temperature ($k_{rp}$; Equation 6). Furthermore, the Ca+Mg concentrations of metamorphic and siliciclastic sediment grid cells are difficult to define. We allow these parameters to vary within reasonable bounds during the calibration stage of GEOCLIM. In total, we test 93,600 unique parameter combinations (Fig. S4; Table S2).

We compute spatially-resolved long-term CO₂ consumption (i.e. Ca+Mg fluxes) using present-day runoff (UNH/GRDC Composite Runoff Fields V1.0; Fekete et al., 1999), temperature (CRU TS v.4.03; Harris et al., 2013), and slope (Shuttle Radar Topography Mission; Farr et al.,
<figure><img src="image_11.png" /><figcaption>**Figure S5.** Chemical weathering and regolith thickness maps associated with the GEOCLIM model run using the parameter combination that produced the highest $r^2$ (Fig. S4B). Parameter values used for this model run are: $k_d = 5 \times 10^{-4}$, $k_w = 1$, $\sigma = -0.4$, $k_{rp} = 1 \times 10^{-2}$, metamorphic$Ca+Mg = 2500$ mol/m$^3$, sediment$Ca+Mg = 2000$ mol/m$^3$.</figcaption></figure>

2007) fields. As described in the main text, we sum the computed CO$_2$ consumption over

large-scale watersheds that appear in the global compilation of Gaillardet et al. (1999), as well as

smaller-scale watersheds of the Amazon Basin (HYBAM network) in the compilation of Moquet
<table><thead><tr><th>k<sub>d</sub><br>unitless</th><th>k<sub>w</sub><br>unitless</th><th>σ<br>unitless</th><th>k<sub>rp</sub><br>unitless</th><th>metamorphic<sub>Ca+Mg</sub><br>mol/m<sup>3</sup></th><th>sediment<sub>Ca+Mg</sub><br>mol/m<sup>3</sup></th></tr></thead><tbody><tr><td>1×10<sup>-5</sup></td><td>1×10<sup>-3</sup></td><td>-0.4</td><td>1.2×10<sup>-3</sup></td><td>1500</td><td>500</td></tr><tr><td>2×10<sup>-5</sup></td><td>2×10<sup>-3</sup></td><td>-0.2</td><td>2×10<sup>-3</sup></td><td>2000</td><td>1000</td></tr><tr><td>5×10<sup>-5</sup></td><td>5×10<sup>-3</sup></td><td>-0.1</td><td>3×10<sup>-3</sup></td><td>2500</td><td>1500</td></tr><tr><td>1×10<sup>-4</sup></td><td>1×10<sup>-2</sup></td><td>0</td><td>5×10<sup>-3</sup></td><td>3000</td><td>2000</td></tr><tr><td>2×10<sup>-4</sup></td><td>2×10<sup>-2</sup></td><td>0.1</td><td>1×10<sup>-2</sup></td><td>3500</td><td>2500</td></tr><tr><td>5×10<sup>-4</sup></td><td>5×10<sup>-2</sup></td><td>0.3</td><td>1.5×10<sup>-2</sup></td><td>4000</td><td>3000</td></tr><tr><td>1×10<sup>-3</sup></td><td>1×10<sup>-1</sup></td><td></td><td></td><td></td><td></td></tr><tr><td>2×10<sup>-3</sup></td><td>2×10<sup>-1</sup></td><td></td><td></td><td></td><td></td></tr><tr><td>5×10<sup>-3</sup></td><td>5×10<sup>-1</sup></td><td></td><td></td><td></td><td></td></tr><tr><td>1×10<sup>-2</sup></td><td>1</td><td></td><td></td><td></td><td></td></tr></tbody></table>

**Table S2.** Values tested for poorly constrained parameters in the silicate weathering component of GEOCLIM. Every permutation of the listed values were tested (except those permutations where the Ca+Mg concentration of the sediments are higher than that of the metamorphics), resulting in 93,600 unique parameter combinations.

et al. (2011), Moquet et al. (2016), and Moquet et al. (2018). The latter are nested watersheds, which requires upstream weathering fluxes to be subtracted from downstream fluxes. Watersheds for which this subtraction yields an aberrant value are not considered. 80 watersheds in total are used in this study. We calculate the coefficient of determination ($r^2$) between computed and measured CO$_2$ consumption in each of these basins:

$$r^2 = 1 - \frac{\sum [\log_{10}(M_i) - \log_{10}(O_i)]^2}{\sum [\log_{10}(O_i) - \overline{\log_{10}(O)}]^2}$$

$M_i$ is the modeled CO$_2$ consumption over watershed $i$, $O_i$ is the observed CO$_2$ consumption over watershed $i$, and $\overline{\log_{10}(O)}$ is the mean of the log of observed CO$_2$ consumption over all watersheds.

The majority of the original 93,600 parameter combinations produce CO$_2$ consumption maps with poor fits to the measured watershed data (Fig. S3). Given that the coefficient of determination ($r^2$) is calculated using Equation rather than fitting a linear model, many of the combinations associated with particularly poor fits result in negative $r^2$ values (Fig. S3).

However, given the right permutation of parameters, GEOCLIM produces CO$_2$ consumption maps that fit the measured watershed data reasonably well. We eliminate all parameter combinations
that produce a $r^2 \le 0.5$, which leaves 573 unique parameter combinations (Figs. S3 and S4).

Global CO₂ consumption calculated from these 573 parameter combinations ranges from 3.7×10¹² mol/yr to 6.3×10¹² mol/yr, with a mean of 5.2×10¹² mol/yr. These estimates fall within the range of independently estimated outgassing rates for the present-day: the full range of estimates of the present-day non-anthropogenic global CO₂ emission rate compiled in Gerlach (2011) is 3.0–10.0×10¹² mol/yr, with the preferred range of those estimates being 3.4–5.9×10¹² mol/yr (Fig. S3B).

Each parameter combination predicts a different total CO₂ consumption that should match the CO₂ degassing at steady-state. Hence, in the non-calibration GEOCLIM experiments, each parameter combination is used with its corresponding steady-state CO₂ degassing. In these estimates, we are not including the effects of reverse weathering (Michalopoulos and Aller, 1995) as it is not clear how it should be parameterized and it is interpreted to be a relatively minor flux in the Cenozoic (Isson and Planavsky, 2018). Formation of authigenic clays that scales with Ca and Mg concentration of riverine waters could decrease the CO₂ consumption associated with silicate weathering although the overall trend seen in this study would remain.

<figure><img src="image_12.png" /><figcaption>**Figure S6.** Ca, Mg, Na, and K concentrations within mafic, intermediate, and felsic lithologies, as calculated based on the mean of data compiled from EarthChem.</figcaption></figure>

A related issue is that our analysis focuses on the weathering of Ca and Mg silicates. Sources of
alkalinity to the ocean include the four major cations of Ca, Mg, Na and K. A common approach for the long-term carbon cycle follows that of Berner et al. (1983) where an emphasis is placed on Ca and Mg with the weathering of Na and K silicates considered to not be a long-term sink of CO2. The reasoning behind this assumption is that while weathering of these silicates can lead to transient uptake of CO2, this uptake is balanced by release. The ultimate mechanism through which silicate weathering sequesters carbon is the incorporation of Ca+Mg carbonates into the lithosphere. Given that Na and K are not precipitating carbonate, when simulating long term steady-state solutions, as we are here, it is valid to follow the approach of Berner et al. (1983) and not consider the Na and K cycles. Complexity arises given the long residence time of Na and K in seawater (~45–80 m.y. for Na, and ~8–10 m.y. for K; Emerson and Hedges, 2008; Lécuyer, 2016) and the lag between the generation of alkalinity and its consumption that could modulate the carbonate system. Dynamic modeling of the effects of the Na and K cycles on the long term carbon cycle is an area ripe for future research although uncertainty in their sinks and the timescale of the processes that control them make them difficult to parameterize. An alternative approach for estimating carbon sequestration associated with silicate weathering that sought to incorporate all four cations was taken by France-Lanord and Derry (1997) who used a formulation of ΔCO2 = ΔMg + ΔCa + 0.15ΔNa + 0.1ΔK. This relationship, also applied by Schopka et al. (2011), assumes that all Ca2+ and Mg2+ eventually forms carbonates and applies an estimate that 20% of K+ and 30% of Na+ exchanges for Ca and Mg and ultimately produces carbonate (and considers charge balance). Different parameterizations such as that of France-Lanord and Derry (1997) could modulate our estimates, but the importance of the emergence of the SEAIs on Neogene pCO2 would remain unchanged. The relative importance of mafic lithologies also remains given that with significantly higher Ca+Mg concentrations and subequal Na+K concentrations to intermediate/felsic lithologies they have the highest Ca+Mg+Na+K concentrations (Fig. S6).
136 Climate Model

137 For the climate model component of GEOCLIM, we use temperature and runoff fields from a 138 subset of the GFDL CM2.0 experiments (Delworth, 2006; Delworth et al., 2006; available for 139 download at https://nomads.gfdl.noaa.gov/dods-data/gfdl_cm2_0/) These experiments were 140 performed in order to explore the effect of various changes in forcing agents on climate since ca. 141 1860 at 2.0° × 2.5° resolution. In the “1860 control” experiment, forcing agents representative of 142 conditions ca. 1860 (including CO2, CH4, N2O, O3, sulfates, carbon, dust, sea salt, solar 143 irradiance, and the distribution of land cover types) are held constant for 500 years after reaching 144 equilibrium. pCO2 in 1860 is assumed to be 286 ppm. In the “+1%/yr to 2×” experiment, initial 145 conditions are taken from the “1860 control” experiment, then pCO2 is prescribed to increase 146 from 286 ppm at a compounded rate of +1% per year for 70 years, when pCO2 reaches double 147 (572 ppm) of the initial value. pCO2 is then held constant until the end of the 280 year 148 experiment. All non-CO2 forcing agents are held constant. The “+1%/yr to 4×” experiment is 149 identical to the “+1%/yr to 2×” experiment, except that pCO2 is prescribed to increase for 140 150 years, when pCO2 reaches quadruple (1144 ppm) of the initial value. pCO2 is then held constant 151 for 160 years. We take the mean of the last 100 years of each of these three experiments (when 152 pCO2 is being held constant at its final level) to obtain temperature and runoff fields for 286, 572, 153 and 1144 ppm pCO2 respectively. Temperature and runoff fields associated with pCO2 values 154 between the three modeled pCO2 levels (286, 572, and 1144 ppm) are obtained through linear 155 interpolation between the model results. 156

156 Given that precipitation in the SEAIs is primarily driven by large-scale convective upwelling 157 leading to high precipitation over both land and water (Donohoe and Voigt, 2017), the current 158 approach of using a climate model forced by modern paleogeography for all SEAIs scenarios is a 159 reasonable approximation for this region. However, for paleogeographic change that would have 160 very significantly modified precipitation over broad swathes of land such as growth of the 161 Himalaya, this approach of using modern climate model results is not viable. The climate model
also assumes the present-day extent of ice-sheets for all $pCO_2$ levels. However, without ice-sheets, the global climate would be warmer for a given $pCO_2$ (hysteresis loop; Pollard and DeConto, 2005), but the amount of emergent land would be lower. We expect this bias to be minimal, as all of our modeled $pCO_2$ estimates fall below the glaciation threshold for Antarctica.

## Southeast Asian Islands Scenarios

We use geological data to quantify the changes in the area of the SEAIs over the past 15 m.y. (described in detail in Paleoshoreline Reconstruction and Geological Synthesis). To generate the lithologic map used in the “removed SEAIs” scenario (Fig. S7), we remove all land associated with arc-continent collision in the SEAIs. To generate the lithologic map used in the “paleo-SEAIs” scenarios (Fig. S7), we take the estimated paleoshorelines of the SEAIs for 15, 10, and 5 Ma, then remove all land that falls outside of these bounds. The slope and lithologic classification of each pixel is identical to that of the present day, provided that these pixels fall within the estimated paleoshorelines of the SEAIs for their respective time slice. This means that while we are changing the amount of emergent land, we are not changing the latitudinal or longitudinal position and are keeping the slope of emerged pixels the same. Given that similar distributions of lithology and slope are reasonable in the past and that latitudinal translation has been relatively minor, we consider this simplification to be reasonable to first-order (Fig. S10). We then add land in the region that is not exposed today. The assignment of lithology, runoff, and slope for pixels on the Sunda Shelf is not trivial, given that the shelf is not currently exposed. As discussed in Paleoshoreline Reconstruction and Geological Synthesis, the Sunda Shelf is a flat-lying, relatively stable platform of continental crust. Islands between Malaysia/Sumatra and Borneo are composed of granite and are surrounded by shallow marine siliciclastic sediments (Darmadi et al., 2007; Hall, 2009, 2013b). However, basement highs of granite do not appear further to the northwest in the Gulf of Thailand, and seismic and drill core data suggest siliciclastic fluvial systems draining out toward the South China Sea when the Sunda Shelf was
<figure><img src="image_13.png" /><figcaption>**Figure S7.** 0.5° × 0.5° lithologic maps of the Southeast Asian islands (SEAIs) used to force GEOCLIM. Only the 15 Ma scenarios are shown here, but the 10 and 5 Ma scenarios would be similar, using the 10 and 5 Ma shorelines shown in Figure 1 of the main text instead. Each column represents a tested scenario. Each row represents a lithologic category. Solid black lines show present-day shorelines. Total areas of each lithologic category within the SEAIs for each of these tested scenarios is shown in Figure S8.</figcaption></figure>

exposed (Darmadi et al., 2007). We therefore assign pixels of the Sunda Shelf between
Malaysia/Sumatra and Borneo to be 50% felsic and 50% sediment, and pixels further to the north
to be 100% sediment. For the slope, we take the mean slope value (0.0043 m/m) of relatively flat
land on eastern Sumatra, and assign this value to all pixels of the Sunda Shelf. For the runoff, we
linearly interpolate along latitude bands between present-day land pixels to obtain runoff values
<figure><img src="image_14.png" /><figcaption>**Figure S8.** Total area of the lithologic categories within the SEAIs for each tested scenario. Note that the grey bars in the middle panel indicate the fractional area of each lithologic category for the entire Earth.</figcaption></figure>

for the presently-submerged Sunda Shelf. This approach preserves large N-S variation in runoff associated with the large-scale Hadley circulation and therefore should be reasonable to first-order, although it does not accurately capture the smaller-scale interactions between land and ocean/atmosphere circulation. For other minor islands in the region that were exposed 15, 10, or 5 m.y. ago, but are currently submerged, we take the same approach as that described for the Sunda Shelf to generate the runoff field. For the slope and lithology, we take the mean slope and Ca+Mg concentrations from neighbouring pixels. We also remove the Greenland ice sheet in the “paleo-SEAIs” scenarios, to be consistent with Northern Hemisphere ice sheets developing after 5 Ma (Haug et al., 2005). The runoff field over Greenland is generated using the same
<figure><img src="image_15.png" /></figure>

**Figure S9.** Steady-state *pCO*₂ estimates from GEOCLIM for the various scenarios discussed in the text. For each of the 13 scenarios, each point represents an estimate from one of the 573 unique parameter combinations that resulted in reasonable total global CO₂ consumption and most closely matched estimates of present-day CO₂ consumption in 80 watersheds around the world. The box encloses the middle 50% of the *pCO*₂ estimates (i.e. the interquartile range), and the notch represents the median with its 95% confidence interval. The whiskers extend to the 2.5 and 97.5 percentile values. This figure is identical to Figure 2 in the main text, except for the addition of the scenarios used to test the sensitivity of the results to the inclusion of the Sunda Shelf and the removal of the Greenland ice sheet.

method as that for the Sunda Shelf. For the slope, we take the mean slope value (0.0529 m/m) of
land for which slope data exists around the edges of Greenland, and assign this value to all pixels
of Greenland covered by the ice sheet. We also assign the Ca+Mg concentration of bulk
continental crust to these pixels.

To explore the sensitivity of our results to the inclusion of the Sunda Shelf and the removal of
<figure><img src="image_16.png" /></figure>

**Figure S10.** Paleogeographically-reconstructed paleoshorelines for the SEAIs at 5, 10, and 15 Ma, using the paleoshoreline reconstruction of this study coupled to the paleogeographic model of Matthews et al. (2016).

the Greenland ice sheet, we also tested “paleo-SEAIs - Greenland” scenarios for 15, 10, and 5 Ma, which are identical to that of the “paleo-SEAIs” scenarios without the removal of the Greenland ice sheet. We also tested “paleo-SEAIs - Greenland - Sunda Shelf” scenarios for 15, 10, and 5 Ma,

206

207

208
which are identical to that of the “paleo-SEAIs” scenarios without the inclusion of the Sunda
Shelf and the removal of the Greenland ice sheet. We find that the estimated steady state $pCO_2$’s
for the “paleo-SEAIs - Greenland” and “paleo-SEAIs - Greenland - Sunda Shelf” scenarios are
only marginally higher than those of the “paleo-SEAIs” scenarios (Fig. S9), suggesting that our
results are relatively insensitive to the inclusion of the Sunda Shelf and the removal of the
Greenland ice sheet. This insensitivity can be attributed to the low Ca+Mg concentrations of
felsic and sediment lithologies in conjunction with the low relief across the Sunda Shelf, and the
low runoff values at high latitudes.

<figure><img src="image_17.png" /><figcaption>**Figure S11.** Proxy-based $pCO_2$ estimates for the past 20 m.y. For data that were compiled within Foster et al. (2017), error bars indicate standardized uncertainties. For the Ji et al. (2018) data, error bars represent the 16-84^th^ percentile of resampled estimates. Individual data points from Cui et al. (2020) were not published with uncertainties. The Bereiter et al. (2015) data come from the EPICA Dome C ice core.</figcaption></figure>

Proxy-based $pCO_2$ estimates for the past 20 m.y. are shown in Figure S11 (Bereiter et al., 2015; Foster et al., 2017; Ji et al., 2018; Cui et al., 2020). Even when stringent quality criteria and the latest understanding of each of the $pCO_2$ proxies have been applied to available $pCO_2$ records, both significant uncertainty in the estimated $pCO_2$ for any given data point as well as
disagreement between techniques remain. Statistical methods that are often utilized in an
attempt to extract trends from this data (e.g. locally-weighted scatter plot smoothing, LOWESS)
evaluate the mean (with some weighting) of the highly scattered data in a given time interval.

However, such methods are only appropriate for estimating the “true” value at any given time
interval in datasets in which individual data points are being drawn from a single distribution
about a mean value. The pCO₂ proxy compilation shown in Figure S11 consists of several distinct
techniques and samples, each with their own set of assumptions and associated probability
distributions that influence their reported values. Therefore, running a regression through these
values neglects the fact that some of these techniques/samples are likely more reliable than the
others. Constraining how robust each pCO₂ proxy technique is remains an important challenge,
and as such the “true” pCO₂ could plausibly lie anywhere within the full range of proxy-based
pCO₂ estimates for any given time interval. Nevertheless, our modeled pCO₂ values for 15 Ma
(Fig. S9) resemble the higher end of proxy-based pCO₂ estimates for the early-mid-Miocene.

Given this scatter in the proxy-based pCO₂ estimates, we instead infer the Neogene cooling trend
from the Miocene benthic foramin oxygen isotope record shown in Figure 1 in the main text.

236 Paleoshoreline Reconstruction and Geological Synthesis

Geological maps, stratigraphic data, and previous paleoshoreline compilations were used to
calculate the changes in different types of subaerially exposed rocks in the SEAIs. Following
Molnar and Cronin (2015), we analyzed the area changes of islands that are larger than ~200 km²
(Fig. 1 in main text). We also included changes in areas of submerged continental shelves that
were previously exposed, like the Sunda Shelf. Larger islands were further divided into regions
based upon their position on microplates as defined by Matthews et al. (2016). Miocene to
present stratigraphic columns were compiled from each region to develop an age model for
paleoenvironmental indicators. In general, carbonate strata and thin-bedded siliciclastic strata
were assumed to be deposited in a subaqueous marine environment, whereas evidence for
exposure and terrestrial siliciclastic deposits were used as indicators of a subaerial
paleoenvironment. Following previous paleoshoreline compilations and these environmental
indicators, paleoshorelines were outlined in QGIS to calculate areas for the Early Pliocene (5 Ma),
Late Miocene (10 Ma), and Middle Miocene (15 Ma). The paleoshoreline reconstructions within
the paleogeographic models broadly coincide with sub-Epoch boundaries (i.e. Early Miocene =
23-16 Ma, Middle Miocene = 15.99-11.61 Ma, Late Miocene = 11.6-5.4 Ma, Pliocene = 5.39-3.8
Ma), but we recognize that the biostratigraphic resolution is an additional source of uncertainty.

Several new sources for paleoshoreline data have become available since Molnar and Cronin (2015). Particularly, for Sulawesi, we follow the recent stratigraphic compilation and paleoshorelines delineated by Nugraha and Hall (2018), and in New Guinea, we follow Gold et al. (2017) and Harrington et al. (2017). Although the outlines of our paleoshorelines are significantly different from those used for 5 Ma by Molnar and Cronin (2015) our calculated area is comparable.

Malay Peninsula and Sunda Shelf

Paleogeographic reconstructions of Peninsular Malaysia suggest that it has been largely exposed over the last 20 Ma (Hall and Nichols, 2002; Hall, 2013b). Although the majority of the Sunda Shelf is currently submerged, large portions of this flat-lying, relatively stable platform of continental crust were also emergent throughout the Neogene, and were repeatedly subaerially exposed during the Pleistocene as recently as the Last Glacial Maximum (Hall and Nichols, 2002).

Similarities between the terrestrial biotas of Borneo and mainland Southeast Asia confirm the
existence of land bridges on the Sunda Shelf throughout the Miocene (Moss and Wilson, 1998).

Eocene through Miocene rifting of the South China Sea resulted in subsidence throughout the Sunda Shelf (Morley and Morley, 2013). Grabens associated with this subsidence became freshwater rift lakes that later transitioned to partially enclosed inland seas and extensive brackish or saline wetland environments. Palynological analysis suggests widespread swamp
environments persisting through the Late Miocene (Morley and Morley, 2013). Basement highs (e.g. the Natuna Arch, currently subaerially exposed as the largely granitic island of Natuna off SW Borneo, as well as the granitic Tin Islands south of the Malay Peninsula) are typically bounded by Paleogene to Neogene basins dominated by shallow marine clastic fill (Darmadi et al., 2007). Most paleogeographic reconstructions of the region incorporate some degree of exposure of the Sunda Shelf from 20 Ma onwards (Hall, 2013b; Madon et al., 2013), although it is omitted in the 5 Ma paleogeographic model of Molnar and Cronin (2015).

## Borneo and Palawan

Our reconstructions of the paleoshorelines of Borneo and Palawan are largely informed by extant paleogeographic constructions detailed in Hall (2001, 2013a,b), as well as geologic maps and local shoreline reconstructions found in van de Weerd and Armin (1992); Witts et al. (2012); Madon et al. (2013); Kessler and Jong (2015). Comprised of Paleozoic to Mesozoic crustal components that were largely accreted to Sundaland by the late Cretaceous (Metcalfe, 2013), the southern and western portions of Borneo (SW Kalimantan) have been subaerially exposed throughout the Neogene (Hall, 2013a,b). A collision along the northern margin of Borneo associated with the initiation of rifting in the South China Sea resulted in the late Eocene uplift of the Central Range mountains (Hutchison, 1996), which provide sediment to basins along the Southern (Kalimantan) and northern (Sabah and Sarawak) coasts of the island. In the east, the Barito, Kutei, and Tarakan basins developed as a single area of subsidence (associated with the opening of the Makassar Straits in the Eocene) before the basins were isolated by Oligocene faulting and Miocene uplift (Witts et al., 2012). The Kutei Basin is characterized by eastward-prograding deltaic and shallow shelf deposits that have been steadily supplied with sediment from the Schwaner Mountains of SW Kalimantan and the Central Range (van de Weerd and Armin, 1992). North of the Kutei Basin, the deltaic deposits found in the Tarakan Basin similarly prograde eastward from the mid Miocene onward, fed by the northern drainages of the Central Range
(Satyana et al., 1999). In SE Kalimantan, the Barito Basin is bounded to the east by the ophiolite-bearing Meratus Mountains. Sedimentological data suggests that the Meratus Mountains were not emergent until the Late Miocene (Witts et al., 2012).

Separated from the Paleozoic continental core of SW Borneo by the Lupar Line suture zone, the northern portion of Borneo (Sabah) is underlain by ophiolitic basement that extends to Palawan (Hall et al., 2008; Ilao et al., 2018). The Sarawak Basin and NW Borneo trough offshore Sabah host >10 km-thick Neogene sedimentary sequences, indicative of the extent and duration of exhumation in northern Borneo (Hall et al., 2008). The Late Oligocene/Early Miocene Sabah Orogeny resulted in the uplift of both Sabah and southern Palawan (Hall, 2013a), as well as the obduction of the Palawan ophiolite and Telupid ophiolite of northern Sabah. Under-thrusting of thinned passive-margin continental crust beneath these suprasubduction ophiolites resulted in Early to Middle Miocene exhumation and offshore unconformities (Hall et al., 2008). Early Miocene sediments in northern Sabah have a provenance from Palawan (Suggate et al., 2014), indicating that Palawan, like northern Borneo, experienced uplift during the Sabah orogeny. However, initiation of back-arc extension in the Sulu Sea around 19 Ma (Hall, 2013a) resulted in the subsidence of the eastern Palawan Mountains, eliminating or reducing the role of Palawan as a sediment source for the Borneo trough. The 14 Ma Capoas Granite on Palawan and the 7.5 Ma Kinabalu Granite in Sabah are interpreted to be the result of crustal thinning related to pulses of regional backarc extension (Hall, 2013a). Thermochronological data from the Sabah highlands suggest extremely rapid uplift (7 km/m.y.) and exhumation during the latest Miocene and early Pliocene, resulting in the recent formation of Mt. Kinabalu, the tallest peak in Borneo (Cottam et al., 2013).

## Sulawesi

For the paleoshorelines of Sulawesi, we largely followed the recent stratigraphic compilation and paleoshorelines compilation of Nugraha and Hall (2018). Geographically, Sulawesi can be divided
into a central highland region flanked by North, South, Southeast, and East Arms. These arms have high relief (>3 km) that are separated by deep basins and broadly define Sulawesi's seven tectonic provinces: 1) the West Sulawesi magmatic arc of the South Arm, 2) the Central Sulawesi metamorphic Belt, 3) the Sangihe arc of the North Arm, 4) the East Sulawesi Ophiolite of the East Arm, 5) the southeast metamorphic belt of the Southeast Arm, and the microcontinental blocks of 6) Banggai-Sula and 7) Buton-Tukang Besi (Hamilton, 1979; Katili, 1978). Like SW Borneo, the North and South Arms, and much of Central Sulawesi, are underlain by continental blocks that rifted off of the Australian-Birds Head margin in the Jurassic and collided in the Cretaceous with Eurasian basement of Sundaland as part of the Woyla Arc system (Parkinson, 1998b; Hennig et al., 2016; Hall, 2017; Hennig et al., 2017a). After collision, subduction polarity reversed and a Cretaceous to Miocene SE-facing volcanic arc developed (Polvé et al., 1997; Elburg and Foden, 1998) that was connected to the paleo-Sunda Arc (Hall, 2002), also referred to as the Great Indonesian arc (Harris, 2006). The ophiolites of Sulawesi were generated in this arc system in a back-arc to intra-arc setting (Monnier et al., 1995) (but see Kadarusman et al., 2004 for an alternative view) and the largest fragments of the East Sulawesi ophiolite were detached from their metamorphic sole in the late Oligocene (32–28 Ma; Parkinson, 1998a) and thrust to the east above a west dipping slab (East Sulawesi block; Villeneuve et al., 2001). During the Early to Middle Miocene, the Tukang Besi-Buton block began to collide obliquely with the southeastern end of the Sunda Arc, causing ophiolite emplacement in West and SE Sulawesi, including Buton (Smith and Silver, 1991; Bergman et al., 1996).

The North Arm consists of Late Pliocene and younger volcanic rocks built upon Eocene to Early Miocene oceanic basalt, basaltic andesite, pelagic sediments, and metamorphic rocks (Elburg and Foden, 1998). Eocene to Early Miocene volcanic rocks formed above a NW dipping slab in an arc system that extended to West and South Sulawesi (van Leeuwen and Muhardjo, 2005). Gorontalo Bay is an extensional basin that formed in the Pliocene, and prior to that time, the East Sulawesi was attached to the North Arm of Sulawesi. After the Oligocene to Early Miocene accretion of the East Sulawesi ophiolite, the arc system stepped out to the SE forming a Neogene (23–16 Ma)
volcanic belt on the East arm (Kadarusman et al., 2004). In the Late Miocene, SE-dipping subduction was initiated below the North Arm, which continues today at the North Sulawesi Trough. North Sulawesi was not substantially emergent until the Pliocene (van Leeuwen and Muhardjo, 2005). Active volcanoes extend from the North Arm through the Sangihe Arc into the southern Philippines.

The East and Southeast Arms consist predominantly of mafic and ultramafic rocks of the East Sulawesi Ophoilite that are exposed for more than 10,000 km² (Monnier et al., 1995). The East Arm preserves a complete ophiolitic sequence underlain by a metamorphic sole, mélange, imbricate continental margin and crystalline basement with a blueschist metamorphic overprint (Silver et al., 1983; Monnier et al., 1995; Parkinson, 1998a). Locally, the structural thickness of the ophiolite exceeds 15 km with surface relief over 3 km (Kadarusman et al., 2004). Seventeen igneous K-Ar dates from the East Sulawesi ophiolite range from 93–32 Ma, clustering between 60 Ma and 40 Ma (Parkinson, 1998b). K-Ar dates on hornblende from the metamorphic sole yielded cooling ages between 36–23 Ma (Parkinson, 1998a; Villeneuve et al., 2001), which dates the initial emplacement of the East Sulawesi ophiolite as the north Sulawesi volcanic arc was underthrusted by the Sula spur (Australian crust) (Silver et al., 1983; Parkinson, 1998a); however, despite these relatively old ages of emplacement, it appears that the ophiolites on the East Arm were not substantially subaerially exposed and exhumed until the Miocene.

Collision of the Sula-Banggai block with the East margin of Sulawesi began in the latest Miocene and uplift accelerated in the early Pliocene (5.2 Ma to 3.8 Ma) and is associated with a major pulse of sedimentation in adjacent basins (Davies, 1990; Villeneuve et al., 2000). Off the northeast margin of the east arm of Sulawesi, Miocene platform carbonates on the Sula-Banggai block are overlain by Late Miocene to Early Pliocene ophiolite detritus in the Celebes mélange (Davies, 1990). Thus, although the East Sulawesi ophiolite was trapped/emplaced onto the composite Sundaland margin (i.e. a fragment of previously accreted Australian crust) between 36 Ma and 23 Ma, mafic and ultramafic rocks on the south and west arms appear to have not
been substantially subaerially exposed in the south and west arms until after 15 Ma, and on the east arm until after 5 Ma.

West Sulawesi rifted from Borneo during the Eocene forming the Makassar Straits back arc basin behind a southwest-facing arc (Polvé et al., 1997). Eroded fragments of ophiolite and extensive belts of volcanic rocks are preserved on the West and South Arms of Sulawesi (Bergman et al., 1996; van Leeuwen et al., 2010). During the Middle Miocene (ca. 15–13 Ma), extensional faults in the Bone basin reversed, which was accompanied by uplift and erosion of the Bone Mountain ophiolite and Lamasi complex in West Sulawesi (Bergman et al., 1996; van Leeuwen et al., 2010). Shortening was likely due to the collision of the leading edge of the Buton-Tukang Besi block, which collided with Buton and the SE Arm of Sulawesi (Smith and Silver, 1991). Uplift and erosion is recorded by the presence of a major Middle Miocene unconformity and sedimentary breccias in marginal basins (Bergman et al., 1996; van Leeuwen et al., 2010). Uplift was diachronous, not effecting units below and to the west of the Lamasi ophiolite until the Middle Miocene (ca. 13 Ma) (van Leeuwen et al., 2010). Alkali volcanism ensued at ca. 11 Ma and is associated with a second phase of extension and exhumation from the Late Miocene to Pliocene with fission track ages implying deep exhumation (Smith and Silver, 1991; Bergman et al., 1996; van Leeuwen et al., 2010).

Fission track ages from granitoids in central Sulawesi indicate rapid uplift of Central Sulawesi (200–700 m/m.y.) starting at about 5 Ma associated with movement on the Palo-Koro fault (Bellier et al., 2006). The fault system also shows a normal component with rapid exhumation of rock west of the fault in western Sulawesi – all fission track dates are younger than 5 Ma (Bellier et al., 2006). Just east of the Palo-Koro fault, the Palu Metamorphic Complex was exhumed in the Late Miocene to early Pliocene in the north (ca. 5.3 Ma) and later Pliocene in the south (ca. 3.1–2.7 Ma) at rates of up to 400 m/m.y. (Hennig et al., 2017b).
New Guinea and Halmahera

Paleoshoreline maps were georeferenced from several previous tectonic and paleogeographic analyses, most notably those of Nichols and Hall (1991); Cloos et al. (2005); Gold et al. (2017); Harrington et al. (2017). These studies were all based to varying degrees on lithological distributions, biostratigraphic, borehole data, and tectonic models. We complimented these data with our stratigraphic compilations for the region to provide further paleoenvironmental context.

Northern New Guinea, including the Melanesian Arc, was emplaced above the Australian plate during the Miocene (Hamilton, 1979; Cloos et al., 2005; van Ufford and Cloos, 2005; Baldwin et al., 2012). Two major ophiolite belts marking the suture – the Irian-Marum ophiolite belt (including the April ultramafics), and the Papuan Ultramafic Belt (PUB) – are preserved along the Central Range and Peninsular Range, respectively. South of the Irian-Marum ophiolites, the Ruffaer Metamorphic Complex constructed from the accretionary wedge, forms the spine of the Central Range (up to ~5 km elevation today).

The Middle Miocene (16–14 Ma) basal Makats Formation contains siliciclastic sediment that was transported from the south into the forearc basin associated with the Irian-Marum ophiolite belt (Visser and Hermes, 1962; Cloos et al., 2005). The beginning of widespread synorogenic sedimentation to the south and on the Australian continental basement was later at ca. 12 Ma (van Ufford and Cloos, 2005). Mountain building began in Late Miocene time, ca. 8–7 Ma (van Ufford and Cloos, 2005; Baldwin et al., 2012), but major relief was not generated until the Pliocene (Weiland and Cloos, 1996). Similarly, the Marum ophiolite was uplifted between 8–5 Ma with 3–4 km of denudation (Crowhurst et al., 1996). An estimated 80–100 km of shortening has been accommodated by deformation on the south side of the Central Range (Hill and Gleadow, 1989; Cloos et al., 2005).

Along the Papuan Peninsula of Eastern New Guinea, the PUB was obducted above a north dipping slab during Oligocene arc-continent collision between Australian continental fragments
and the Melanesian arc, but remained largely subaqueous until the Miocene (Davies and Smite, 1971; Davies and Jaques, 1984; van Ufford and Cloos, 2005). Miocene to present arc-continent collision in New Guinea has progressed from west to east. Exhumation of the Central Range accelerated over the past 4 m.y. which is interpreted to be the result of slab-breakoff and buoyant uplift (Cloos et al., 2005). A change to left-lateral lateral motion on the northern coast of New Guinea at this time, caused the exhumation of ophiolites both along the coast (Monnier et al., 1999) and on the islands Obi (Ali et al., 2001) and Halmahera (Hall et al., 1988; Ballantyne, 1992). In eastern New Guinea, progressive jamming of the north-dipping subduction zone has caused major uplift over the past 4 Ma (van Ufford and Cloos, 2005), which is well dated with a change in provenance from continental to volcanic detritus (Abbott et al., 1994) and thermochronology (Hill and Gleadow, 1989).

In the Bird's Neck of western New Guinea, the Lengguru fold belt formed during Middle to Late Miocene clockwise rotation and obduction of the Weyland Terrane (Bailly et al., 2009). Recent counter-clockwise rotation has exhumed several core complexes and has produced Pliocene compressional deformation to the southwest in the Misool-Onin-Kumawa Ridge (Sapin et al., 2009).

## Philippines

Our reconstructions of the paleoshorelines of the Phillipines were primarily informed by the paleogeographic reconstructions of Hall (1997), coupled with geologic maps of the archipelago (Geologic Survey Division, 1963). The Philippine ophiolites have been grouped in four distinct belts (Yumul, 2007). We combine belts 1 and 2 in an eastern suture zone and belts 3 and 4 in a western suture zone. Belts 1 and 2 were juxtaposed in the eastern suture zone during sinistral transpression, which resulted in Early Miocene uplift and deposition of coarse clastic sediments (Pubellier et al., 1991). The Palawan microcontinental block and Philippine mobile belt collided during late Early Miocene to early Middle Miocene resulting into the emplacement of the western
ophiolites, with exhumation extending from the Late Miocene to the present (Yumul et al., 2013).
In addition to collision, arc magmatism contributed significantly to crustal growth in the
Philippines (Dimalanta and Yumul, 2006).

Sunda-Banda Arc

Paleoshoreline reconstructions of the Sunda-Banda Arc system relied largely on the work of Hall (2013a). The Sunda-Banda Arc system is composed of thick sequences of mixed sedimentary rocks and basement intruded and overlain by igneous rocks derived from continual arc magmatism through much of the Cenozoic (Hall, 2017). Sumatra is the largest island within the Sunda-Banda Arc system that stretches from the Andaman Sea in the northwest to the Banda Sea in the east. The western portion of the Sumatra is underlain by the Wolya ophiolite, which was obducted in the Cretaceous and exhumed in the Eocene to Oligocene (Allen et al., 2008), whereas the southeast corner on the islands of Bangka and Belitung is composed of granite basement of the Sunda Shelf (Hall, 2009, 2013b). Java is composed of a complex of E-W-striking deformational and magmatic belts (Audley-Charles, 2004). The belts have a volcanic arc in the south and a continental shelf to the north with an adjacent sedimentary basin (Audley-Charles, 2004). Bali and Flores are mostly volcanic with some sedimentary cover (Hall, 2009).

Exhumation of the modern Sunda-Banda Arc is the result of ongoing arc-continent collision with the subducting Australian Plate (Harris, 2006). From ca. 20–10 Ma, the majority of Sumatra and nearly all of Java were submerged, although Bangka and Belitung were exposed in addition to a portion of the greater Sunda Shelf (Hall, 2009, 2013b). Six major islands were subaerially exposed before 5 Ma: Sumatra, Belitung, Bangka, Java, Bali, and Flores. However, most of Sumatra and Java were elevated above sea level and emerged to their present exposures only since 5 Ma (Hall, 2009, 2013b). Most of the non-volcanic islands of the Outer Banda Arc emerged after 5 Ma, associated with slab roll-back and collision with the Australian continental margin (Audley-Charles, 2004; Harris, 2006; Hall, 2013b). In Timor and Sumba, arc-continent
473 collision resulted in rapid uplift of deep marine sedimentary rocks to elevations >1 km above sea

474 level with estimated average uplift rates of 1.5 km/m.y. (Audley-Charles, 1986).
475

References

476

477

478

479

Abbott, L. D., Silver, E. A., Thompson, P. R., Filewicz, M. V., Schneider, C., and Abdoerrias, 1994, Stratigraphic constraints on the development and timing of arc-continent collision in northern Papua New Guinea: Journal of Sedimentary Research, vol. 64, pp. 169–183, doi:10.1306/D4267F82-2B26-11D7-8648000102C1865D.

480

481

482

Ali, J., Hall, R., and Baker, S., 2001, Palaeomagnetic data from a Mesozoic Philippine Sea Plate ophiolite on Obi Island, Eastern Indonesia: Journal of Asian Earth Sciences, vol. 19, pp. 535–546, doi:10.1016/s1367-9120(00)00053-5.

483

484

485

486

487

Allen, R., Carter, A., Najman, Y., Bandopadhyay, P., Chapman, H., Bickle, M., Garzanti, E., Vezzoli, G., Andò, S., Foster, G., and et al., 2008, New constraints on the sedimentation and uplift history of the Andaman-Nicobar accretionary prism, South Andaman Island: Special Paper 436: Formation and Applications of the Sedimentary Record in Arc Collision Zones, pp. 223–255, doi:10.1130/2008.2436(11).

488

489

490

Audley-Charles, M., 2004, Ocean trench blocked and obliterated by Banda forearc collision with Australian proximal continental slope: Tectonophysics, vol. 389, pp. 65–79, doi:10.1016/j.tecto.2004.07.048.

491

492

493

Audley-Charles, M. G., 1986, Rates of Neogene and Quaternary tectonic movements in the Southern Banda Arc based on micropalaeontology: Journal of the Geological Society, vol. 143, pp. 161–175, doi:10.1144/gsjgs.143.1.0161.

494

495

496

Bailly, V., Pubellier, M., Ringenbach, J.-C., de Sigoyer, J., and Sapin, F., 2009, Deformation zone “jumps” in a young convergent setting; the Lengguru fold-and-thrust belt, New Guinea Island: Lithos, vol. 113, pp. 306–317, doi:10.1016/j.lithos.2009.08.013.

497

498

499

Baldwin, S. L., Fitzgerald, P. G., and Webb, L. E., 2012, Tectonics of the New Guinea region: Annual Review of Earth and Planetary Sciences, vol. 40, pp. 495–520, doi:10.1146/annurev-earth-040809-152540.

500

501

502

Ballantyne, P., 1992, Petrology and geochemistry of the plutonic rocks of the Halmahera ophiolite, eastern Indonesia, an analogue of modern oceanic forearc: Geological Society, London, Special Publications, vol. 60, pp. 179–202, doi:10.1144/gsl.sp.1992.060.01.11.

503

504

505

506

Bellier, O., Sébrier, M., Seward, D., Beaudouin, T., Villeneuve, M., and Putranto, E., 2006, Fission track and fault kinematics analyses for new insight into the Late Cenozoic tectonic regime changes in West-Central Sulawesi (Indonesia): Tectonophysics, vol. 413, pp. 201–220, doi:10.1016/j.tecto.2005.10.036.

507

508

509

Bereiter, B., Eggleston, S., Schmitt, J., Nehrbass-Ahles, C., Stocker, T. F., Fischer, H., Kipfstuhl, S., and Chappellaz, J., 2015, Revision of the EPICA Dome C CO₂ record from 800 to 600 kyr before present: Geophysical Research Letters, vol. 42, pp. 542–549, doi:10.1002/2014gl061957.

510

511

Bergman, S. C., Coffield, D. Q., Talbot, J. P., and Garrard, R. A., 1996, Tertiary tectonic and magmatic evolution of western Sulawesi and the Makassar Strait, Indonesia: evidence for a
512
513

Miocene continent-continent collision: Geological Society, London, Special Publications, vol. 106, pp. 391–429, doi:10.1144/gsl.sp.1996.106.01.25.

514
515
516

Berner, R. A., Lasaga, A. C., and Garrels, R. M., 1983, The carbonate-silicate geochemical cycle and its effect on atmospheric carbon dioxide over the past 100 million years: American Journal of Science, vol. 283, pp. 641–683.

517
518

Brantley, S. L., 2003, Reaction kinetics of primary rock-forming minerals under ambient conditions: Treatise on Geochemistry, pp. 73–117, doi:10.1016/b0-08-043751-6/05075-1.

519
520
521

Carroll, S. A. and Knauss, K. G., 2005, Dependence of labradorite dissolution kinetics on CO₂(aq), Al(aq), and temperature: Chemical Geology, vol. 217, pp. 213–225, doi:10.1016/j.chemgeo.2004.12.008.

522
523
524

Cloos, M., Sapiie, B., Quarles van Ufford, A., Weiland, R. J., Warren, P. Q., and McMahon, T. P., 2005, Collisional delamination in New Guinea: The geotectonics of subducting slab breakoff: Geological Society of America Special Papers, pp. 1–51, doi:10.1130/2005.2400.

525
526
527

Cottam, M. A., Hall, R., Sperber, C., Kohn, B. P., Forster, M. A., and Batt, G. E., 2013, Neogene rock uplift and erosion in northern Borneo: evidence from the Kinabalu granite, Mount Kinabalu: Journal of the Geological Society, vol. 170, pp. 805–816, doi:10.1144/jgs2011-130.

528
529
530

Crowhurst, P. V., Hill, K. C., Foster, D. A., and Bennett, A. P., 1996, Thermochronological and geochemical constraints on the tectonic evolution of northern Papua New Guinea: Geological Society, London, Special Publications, vol. 106, pp. 525–537, doi:10.1144/gsl.sp.1996.106.01.33.

531
532

Cui, Y., Schubert, B. A., and Jahren, A. H., 2020, A 23 m.y. record of low atmospheric CO₂: Geology, doi:10.1130/g47681.1.

533
534
535

Darmadi, Y., Willis, B., and Dorobek, S., 2007, Three-dimensional seismic architecture of fluvial sequences on the low-gradient Sunda Shelf, offshore Indonesia: Journal of Sedimentary Research, vol. 77, pp. 225–238, doi:10.2110/jsr.2007.024.

536
537
538

Davies, H. L. and Jaques, A. L., 1984, Emplacement of ophiolite in Papua New Guinea: Geological Society, London, Special Publications, vol. 13, pp. 341–349, doi:10.1144/gsl.sp.1984.013.01.27.

539
540

Davies, H. L. and Smite, I. E., 1971, Geology of eastern Papua: Geological Society of America Bulletin, vol. 82, pp. 3299–3312, doi:10.1130/0016-7606(1971)82[3299:goep]2.0.co;2.

541
542
543

Davies, I. C., 1990, Geological and exploration review of the Tomori PSC, eastern Indonesia: Indonesian Petroleum Association Proceedings, 19th Annual Convention, doi:10.29118/ipa.1845.41.67.

544
545

Delworth, T., 2006, Preface to the special section: Journal of Climate, vol. 19, pp. 641–641, doi:10.1175/jcli9016.1.
546

547

548

549

Delworth, T. L., Broccoli, A. J., Rosati, A., Stouffer, R. J., Balaji, V., Beesley, J. A., Cooke, W. F., Dixon, K. W., Dunne, J., Dunne, K. A., and et al., 2006, GFDL's CM2 global coupled climate models. part i: Formulation and simulation characteristics: Journal of Climate, vol. 19, pp. 643-674, doi:10.1175/jcli3629.1.

550

551

552

553

Dessert, C., Dupré, B., François, L. M., Schott, J., Gaillardet, J., Chakrapani, G., and Bajpai, S., 2001, Erosion of Deccan Traps determined by river geochemistry: impact on the global climate and the ⁸⁷Sr/⁸⁶Sr ratio of seawater: Earth and Planetary Science Letters, vol. 188, pp. 459-474, doi:10.1016/S0012-821X(01)00317-X.

554

555

556

Dimalanta, C. B. and Yumul, G. P., 2006, Magmatic and amagmatic contributions to crustal growth in the Philippine island arc system: Comparison of the Cretaceous and post-Cretaceous periods: Geosciences Journal, vol. 10, pp. 321-329, doi:10.1007/bf02910373.

557

558

Donohoe, A. and Voigt, A., 2017, Why future shifts in tropical precipitation will likely be small: Geophysical Monograph Series, pp. 115-137, doi:10.1002/9781119068020.ch8.

559

560

561

Elburg, M. and Foden, J., 1998, Temporal changes in arc magma geochemistry, northern Sulawesi, Indonesia: Earth and Planetary Science Letters, vol. 163, pp. 381-398, doi:10.1016/s0012-821x(98)00143-5.

562

563

Emerson, S. and Hedges, J., 2008, Chemical oceanography and the marine carbon cycle: Cambridge University Press.

564

565

566

Farr, T. G., Rosen, P. A., Caro, E., Crippen, R., Duren, R., Hensley, S., Kobrick, M., Paller, M., Rodriguez, E., Roth, L., and et al., 2007, The shuttle radar topography mission: Reviews of Geophysics, vol. 45, pp. 1-33, doi:10.1029/2005rg000183.

567

568

569

Fekete, B. M., Vörösmarty, C. J., and Grabs, W., 1999, Global, composite runoff fields based on observed river discharge and simulated water balances: Global Runoff Data Centre Koblenz, Germany.

570

571

572

Foster, G. L., Royer, D. L., and Lunt, D. J., 2017, Future climate forcing potentially without precedent in the last 420 million years: Nature Communications, vol. 8, pp. 1-8, doi:10.1038/ncomms14845.

573

574

France-Lanord, C. and Derry, L. A., 1997, Organic carbon burial forcing of the carbon cycle from Himalayan erosion: Nature, vol. 390, pp. 65-67, doi:10.1038/36324.

575

576

577

Gaillardet, J., Dupré, B., Louvat, P., and Allègre, C. J., 1999, Global silicate weathering and CO₂ consumption rates deduced from the chemistry of large rivers: Chemical Geology, vol. 159, pp. 3-30, doi:10.1016/S0009-2541(99)00031-5.

578

579

Geologic Survey Division, 1963, Geological map of the Philippines: Tech. rep., Bureau of Mines, in coordination with the Board of Technical Surveys and Maps.

580

581

Gerlach, T., 2011, Volcanic versus anthropogenic carbon dioxide: Eos, Transactions American Geophysical Union, vol. 92, pp. 201-202, doi:10.1029/2011eo240001.
582

583

584

Gold, D., White, L., Gunawan, I., and BouDagher-Fadel, M., 2017, Relative sea-level change in western New Guinea recorded by regional biostratigraphic data: Marine and Petroleum Geology, vol. 86, pp. 1133–1158, doi:10.1016/j.marpetgeo.2017.07.016.

585

586

Hall, R., 1997, Cenozoic plate tectonic reconstructions of SE Asia: Geological Society, London, Special Publications, vol. 126, pp. 11–23, doi:10.1144/gsl.sp.1997.126.01.03.

587

588

589

Hall, R., 2001, Cenozoic reconstructions of SE Asia and the SW Pacific: changing pattern of land and sea: *In* Metcalfe, I., Smith, J., Morwood, M., and Davidson, I., eds., Faunal and Floral Migrations and Evolution in SE Asia-Australasia, Swets and Zeitlinger, pp. 35–56.

590

591

592

Hall, R., 2002, Cenozoic geological and plate tectonic evolution of SE Asia and the SW Pacific: computer-based reconstructions, model and animations: Journal of Asian Earth Sciences, vol. 20, pp. 353–431, doi:10.1016/s1367-9120(01)00069-4.

593

594

Hall, R., 2009, The Eurasian SE Asian margin as a modern example of an accretionary orogen: Geological Society, London, Special Publications, vol. 318, pp. 351–372, doi:10.1144/sp318.13.

595

596

Hall, R., 2013a, Contraction and extension in northern Borneo driven by subduction rollback: Journal of Asian Earth Sciences, vol. 76, pp. 399–411, doi:10.1016/j.jseaes.2013.04.010.

597

598

Hall, R., 2013b, The palaeogeography of Sundaland and Wallacea since the Late Jurassic: Journal of Limnology, vol. 72, pp. 1–17, doi:10.4081/jlimnol.2013.s2.e1.

599

600

601

Hall, R., 2017, Southeast Asia: New views of the geology of the Malay Archipelago: Annual Review of Earth and Planetary Sciences, vol. 45, pp. 331–358, doi:10.1146/annurev-earth-063016-020633.

602

603

604

Hall, R., Audley-Charles, M. G., Banner, F. T., Hidayat, S., and Tobing, S. L., 1988, Basement rocks of the Halmahera region, eastern Indonesia: a Late Cretaceous–early Tertiary arc and fore-arc: Journal of the Geological Society, vol. 145, pp. 65–84, doi:10.1144/gsjgs.145.1.0065.

605

606

607

Hall, R. and Nichols, G., 2002, Cenozoic sedimentation and tectonics in Borneo: climatic influences on orogenesis: Geological Society, London, Special Publications, vol. 191, pp. 5–22, doi:10.1144/gsl.sp.2002.191.01.02.

608

609

Hall, R., van Hattum, M. W., and Spakman, W., 2008, Impact of India–Asia collision on SE Asia: The record in Borneo: Tectonophysics, vol. 451, pp. 366–389, doi:10.1016/j.tecto.2007.11.058.

610

611

Hamilton, W. B., 1979, Tectonics of the Indonesian region: Professional Paper, p. 345, doi:10.3133/pp1078.

612

613

614

Harrington, L., Zahirovic, S., Flament, N., and Müller, R. D., 2017, The role of deep Earth dynamics in driving the flooding and emergence of New Guinea since the Jurassic: Earth and Planetary Science Letters, vol. 479, pp. 273–283, doi:10.1016/j.epsl.2017.09.039.

615

616

617

Harris, I., Jones, P., Osborn, T., and Lister, D., 2013, Updated high-resolution grids of monthly climatic observations - the CRU TS3.10 dataset: International Journal of Climatology, vol. 34, pp. 623–642, doi:10.1002/joc.3711.
Harris, R., 2006, Rise and fall of the Eastern Great Indonesian arc recorded by the assembly, dispersion and accretion of the Banda Terrane, Timor: Gondwana Research, vol. 10, pp. 207–231, doi:10.1016/j.gr.2006.05.010.

Hartmann, J. and Moosdorf, N., 2012, The new global lithological map database GLiM: A representation of rock properties at the Earth surface: Geochemistry, Geophysics, Geosystems, vol. 13, pp. 1–37, doi:10.1029/2012GC004370.

Haug, G. H., Ganopolski, A., Sigman, D. M., Rosell-Mele, A., Swann, G. E. A., Tiedemann, R., Jaccard, S. L., Bollmann, J., Maslin, M. A., Leng, M. J., and Eglinton, G., 2005, North Pacific seasonality and the glaciation of North America 2.7 million years ago: Nature, vol. 433, pp. 821–825, doi:10.1038/nature03332.

Hennig, J., Breitfeld, H. T., Hall, R., and Nugraha, A. S., 2017a, The Mesozoic tectono-magmatic evolution at the Paleo-Pacific subduction zone in West Borneo: Gondwana Research, vol. 48, pp. 292–310, doi:10.1016/j.gr.2017.05.001.

Hennig, J., Hall, R., and Armstrong, R. A., 2016, U-Pb zircon geochronology of rocks from west Central Sulawesi, Indonesia: Extension-related metamorphism and magmatism during the early stages of mountain building: Gondwana Research, vol. 32, pp. 41–63, doi:10.1016/j.gr.2014.12.012.

Hennig, J., Hall, R., Forster, M. A., Kohn, B. P., and Lister, G. S., 2017b, Rapid cooling and exhumation as a consequence of extension and crustal thinning: Inferences from the Late Miocene to Pliocene Palu Metamorphic Complex, Sulawesi, Indonesia: Tectonophysics, vol. 712-713, pp. 600–622, doi:10.1016/j.tecto.2017.06.025.

Hill, K. C. and Gleadow, A. J. W., 1989, Uplift and thermal history of the Papuan Fold Belt, Papua New Guinea: Apatite fission track analysis: Australian Journal of Earth Sciences, vol. 36, pp. 515–539, doi:10.1080/08120098908729507.

Hutchison, C. S., 1996, The “Rajang accretionary prism” and “Lupar Line” problem of Borneo: Geological Society, London, Special Publications, vol. 106, pp. 247–261, doi:10.1144/gsl.sp.1996.106.01.16.

Ilao, K. A., Morley, C. K., and Aurelio, M. A., 2018, 3D seismic investigation of the structural and stratigraphic characteristics of the Pagasa Wedge, Southwest Palawan Basin, Philippines, and their tectonic implications: Journal of Asian Earth Sciences, vol. 154, pp. 213–237, doi:10.1016/j.jseaes.2017.12.017.

Isson, T. T. and Planavsky, N. J., 2018, Reverse weathering as a long-term stabilizer of marine pH and planetary climate: Nature, vol. 560, pp. 471–475, doi:10.1038/s41586-018-0408-4.

Ji, S., Nie, J., Lechler, A., Huntington, K. W., Heitmann, E. O., and Breecker, D. O., 2018, A symmetrical CO₂ peak and asymmetrical climate change during the middle Miocene: Earth and Planetary Science Letters, vol. 499, pp. 134–144, doi:10.1016/j.epsl.2018.07.011.
Kadarusman, A., Miyashita, S., Maruyama, S., Parkinson, C. D., and Ishikawa, A., 2004, Petrology, geochemistry and paleogeographic reconstruction of the East Sulawesi Ophiolite, Indonesia: Tectonophysics, vol. 392, pp. 55–83, doi:10.1016/j.tecto.2004.04.008.

Katili, J. A., 1978, Past and present geotectonic position of Sulawesi, Indonesia: Tectonophysics, vol. 45, pp. 289–322, doi:10.1016/0040-1951(78)90166-x.

Kessler, F.-L. and Jong, J., 2015, Tertiary uplift and the Miocene evolution of the NW Borneo Shelf margin: Berita Sedimentologi - Indonesian Journal of Sedimentary Geology, vol. 33, pp. 21–46.

Knauss, K. G., Nguyen, S. N., and Weed, H. C., 1993, Diopside dissolution kinetics as a function of pH, CO₂, temperature, and time: Geochimica et Cosmochimica Acta, vol. 57, pp. 285–294, doi:10.1016/0016-7037(93)90431-u.

Lague, D., 2013, The stream power river incision model: evidence, theory and beyond: Earth Surface Processes and Landforms, vol. 39, pp. 38–61, doi:10.1002/esp.3462.

Lécuyer, C., 2016, Seawater residence times of some elements of geochemical interest and the salinity of the oceans: Bulletin de la Société Géologique de France, vol. 187, pp. 245–260, doi:10.2113/gssgbull.187.6.245.

Li, G., Hartmann, J., Derry, L. A., West, A. J., You, C.-F., Long, X., Zhan, T., Li, L., Li, G., Qiu, W., Li, T., Liu, L., Chen, Y., Ji, J., Zhao, L., and Chen, J., 2016, Temperature dependence of basalt weathering: Earth and Planetary Science Letters, vol. 443, pp. 59–69, doi:10.1016/j.epsl.2016.03.015.

Madon, M., Kim, C. L., and Wong, R., 2013, The structure and stratigraphy of deepwater Sarawak, Malaysia: Implications for tectonic evolution: Journal of Asian Earth Sciences, vol. 76, pp. 312–333, doi:10.1016/j.jseaes.2013.04.040.

Matthews, K. J., Maloney, K. T., Zahirovic, S., Williams, S. E., Seton, M., and Müller, R. D., 2016, Global plate boundary evolution and kinematics since the late Paleozoic: Global and Planetary Change, vol. 146, pp. 226–250, doi:10.1016/j.gloplacha.2016.10.002.

Metcalfe, I., 2013, Gondwana dispersion and Asian accretion: Tectonic and palaeogeographic evolution of eastern Tethys: Journal of Asian Earth Sciences, vol. 66, pp. 1–33, doi:10.1016/j.jseaes.2012.12.020.

Michalopoulos, P. and Aller, R. C., 1995, Rapid clay mineral formation in Amazon Delta sediments: Reverse weathering and oceanic elemental cycles: Science, vol. 270, pp. 614–617, doi:10.1126/science.270.5236.614.

Milliman, J. D. and Farnsworth, K. L., 2013, River Discharge to the Coastal Ocean: A Global Synthesis: Cambridge University Press, doi:10.1017/cbo9780511781247, URL http://dx.doi.org/10.1017/CBO9780511781247.
Molnar, P. and Cronin, T. W., 2015, Growth of the Maritime Continent and its possible contribution to recurring ice ages: Paleoceanography, vol. 30, pp. 196–225, doi:10.1002/2014PA002752.

Monnier, C., Girardeau, J., Maury, R. C., and Cotten, J., 1995, Back-arc basin origin for the East Sulawesi ophiolite (eastern Indonesia): Geology, vol. 23, pp. 851–854, doi:10.1130/0091-7613(1995)023(0851: baboft)2.3.co;2.

Monnier, C., Girardeau, J., Pubellier, M., Polvé, M., Permana, H., and Bellon, H., 1999, Petrology and geochemistry of the Cyclops ophiolites (Irian Jaya, East Indonesia): Consequences for the Cenozoic evolution of the north Australian margin: Mineralogy and Petrology, vol. 65, pp. 1–28, doi:10.1007/bf01161574.

Moquet, J.-S., Crave, A., Viers, J., Seyler, P., Armijos, E., Bourrel, L., Chavarri, E., Lagane, C., Laraque, A., Casimiro, W. S. L., and et al., 2011, Chemical weathering and atmospheric/soil CO₂ uptake in the Andean and Foreland Amazon basins: Chemical Geology, vol. 287, pp. 1–26, doi:10.1016/j.chemgeo.2011.01.005.

Moquet, J.-S., Guyot, J.-L., Crave, A., Viers, J., Filizola, N., Martinez, J.-M., Oliveira, T. C., Sánchez, L. S. H., Lagane, C., Casimiro, W. S. L., and et al., 2016, Amazon River dissolved load: temporal dynamics and annual budget from the Andes to the ocean: Environmental Science and Pollution Research, vol. 23, pp. 11,405–11,429, doi:10.1007/s11356-015-5503-6.

Moquet, J.-S., Guyot, J.-L., Morera, S., Crave, A., Rau, P., Vauchel, P., Lagane, C., Sondag, F., Lavado, C. W., Pombosa, R., and et al., 2018, Temporal variability and annual budget of inorganic dissolved matter in Andean Pacific Rivers located along a climate gradient from northern Ecuador to southern Peru: Comptes Rendus Geoscience, vol. 350, pp. 76–87, doi:10.1016/j.crte.2017.11.002.

Morley, R. J. and Morley, H. P., 2013, Mid Cenozoic freshwater wetlands of the Sunda region: Journal of Limnology, vol. 72, pp. 18–35, doi:10.4081/jlimnol.2013.s2.e3.

Moss, S. J. and Wilson, M. E. J., 1998, Biogeographic implications of the Tertiary palaeogeographic evolution of Sulawesi and Borneo: In Hall, R. and Holloway, J. D., eds., Biogeography and Geological Evolution of SE Asia, Backbuys Publishers Leiden, The Netherlands, vol. 133, pp. 133–163.

Nichols, G. J. and Hall, R., 1991, Basin formation and Neogene sedimentation in a backarc setting, Halmahera, eastern Indonesia: Marine and Petroleum Geology, vol. 8, pp. 50–61, doi:10.1016/0264-8172(91)90044-2.

Nugraha, A. M. S. and Hall, R., 2018, Late Cenozoic palaeogeography of Sulawesi, Indonesia: Palaeogeography, Palaeoclimatology, Palaeoecology, vol. 490, pp. 191–209, doi:10.1016/j.palaeo.2017.10.033.

Parkinson, C., 1998a, Emplacement of the East Sulawesi Ophiolite: evidence from subophiolite metamorphic rocks: Journal of Asian Earth Sciences, vol. 16, pp. 13–28, doi:10.1016/s0743-9547(97)00039-1.
727
728
729

Parkinson, C., 1998b, An outline of the petrology, structure and age of the Pompangeo Schist Complex of central Sulawesi, Indonesia: The Island Arc, vol. 7, pp. 231–245, doi:10.1046/j.1440-1738.1998.00171.x.

730
731

Pollard, D. and DeConto, R. M., 2005, Hysteresis in Cenozoic Antarctic ice-sheet variations: Global and Planetary Change, vol. 45, pp. 9–21, doi:10.1016/j.gloplacha.2004.09.011.

732
733
734
735

Polvé, M., Maury, R., Bellon, H., Rangin, C., Priadi, B., Yuwono, S., Joron, J., and Atmadja, R., 1997, Magmatic evolution of Sulawesi (Indonesia): constraints on the Cenozoic geodynamic history of the Sundaland active margin: Tectonophysics, vol. 272, pp. 69–92, doi:10.1016/s0040-1951(96)00276-4.

736
737
738
739

Pubellier, M., Quebral, R., Rangin, C., Deffontaines, B., Muller, C., Butterlin, J., and Manzano, J., 1991, The Mindanao collision zone: a soft collision event within a continuous Neogene strike-slip setting: Journal of Southeast Asian Earth Sciences, vol. 6, pp. 239–248, doi:10.1016/0743-9547(91)90070-e.

740
741
742
743

Sapin, F., Pubellier, M., Ringenbach, J.-C., and Bailly, V., 2009, Alternating thin versus thick-skinned decollements, example in a fast tectonic setting: The Misool-Onin-Kumawa Ridge (West Papua): Journal of Structural Geology, vol. 31, pp. 444–459, doi:10.1016/j.jsg.2009.01.010.

744
745
746
747

Satyana, A. H., Nugroho, D., and Surantoko, I., 1999, Tectonic controls on the hydrocarbon habitats of the Barito, Kutei, and Tarakan Basins, Eastern Kalimantan, Indonesia: major dissimilarities in adjoining basins: Journal of Asian Earth Sciences, vol. 17, pp. 99–122, doi:10.1016/s0743-9547(98)00059-2.

748
749
750
751

Schopka, H., Derry, L., and Arcilla, C., 2011, Chemical weathering, river geochemistry and atmospheric carbon fluxes from volcanic and ultramafic regions on Luzon Island, the Philippines: Geochimica et Cosmochimica Acta, vol. 75, pp. 978–1002, doi:10.1016/j.gca.2010.11.014.

752
753
754

Silver, E. A., McCaffrey, R., Joyodiwiryo, Y., and Stevens, S., 1983, Ophiolite emplacement by collision between the Sula Platform and the Sulawesi Island Arc, Indonesia: Journal of Geophysical Research: Solid Earth, vol. 88, pp. 9419–9435, doi:10.1029/jb088ib11p09419.

755
756
757

Smith, R. B. and Silver, E. A., 1991, Geology of a Miocene collision complex, Buton, eastern Indonesia: Geological Society of America Bulletin, vol. 103, pp. 660–678, doi:10.1130/0016-7606(1991)103(0660:goamcc)2.3.co;2.

758
759
760
761

Suggate, S. M., Cottam, M. A., Hall, R., Sevastjanova, I., Forster, M. A., White, L. T., Armstrong, R. A., Carter, A., and Mojares, E., 2014, South China continental margin signature for sandstones and granites from Palawan, Philippines: Gondwana Research, vol. 26, pp. 699–718, doi:10.1016/j.gr.2013.07.006.

762
763
764

van de Weerd, A. A. and Armin, R. A., 1992, Origin and evolution of the Tertiary hydrocarbon-bearing basin in Kalimantan (Borneo), Indonesia: AAPG Bulletin, vol. 76, pp. 1778–1803.
765 van Leeuwen, T. M. and Muhardjo, 2005, Stratigraphy and tectonic setting of the Cretaceous and Paleogene volcanic-sedimentary successions in northwest Sulawesi, Indonesia: implications for the Cenozoic evolution of Western and Northern Sulawesi: Journal of Asian Earth Sciences, vol. 25, pp. 481–511, doi:10.1016/j.jseaes.2004.05.004.

769 van Leeuwen, T. M., Susanto, E. S., Maryanto, S., Hadiwisastra, S., Sudijono, Muhardjo, and Prihardjo, 2010, Tectonostratigraphic evolution of Cenozoic marginal basin and continental margin successions in the Bone Mountains, Southwest Sulawesi, Indonesia: Journal of Asian Earth Sciences, vol. 38, pp. 233–254, doi:10.1016/j.jseaes.2009.11.005.

773 van Ufford, A. Q. and Cloos, M., 2005, Cenozoic tectonics of New Guinea: AAPG Bulletin, vol. 89, pp. 119–140, doi:10.1306/08300403073.

775 Villeneuve, M., Cornée, J.-J., Gunawan, W., Janin, M.-C., Butterlin, J., Saint-Marc, P., and Samodra, H., 2000, Continental block collision in the eastern arm of Sulawesi (Indonesia). 776 Structure and geodynamic interpretation: Comptes Rendus de l'Académie des Sciences - Series IIA - Earth and Planetary Science, vol. 330, pp. 371–378, doi:10.1016/s1251-8050(00)00163-4.
778

779 Villeneuve, M., Gunawan, W., Cornee, J.-J., and Vidal, O., 2001, Geology of the central Sulawesi belt (eastern Indonesia): constraints for geodynamic models: International Journal of Earth Sciences, vol. 91, pp. 524–537, doi:10.1007/s005310100228.

782 Visser, W. A. and Hermes, J. J., 1962, Geological results of the exploration for oil in Netherlands New Guinea, carried out by the Netherlandsche Nieuw Guinee Petroleum Maatschappij, 783 1935-1960: Staatsdrukkerij-en Uitgeverijbedrijf.
784

785 Weiland, R. J. and Cloos, M., 1996, Pliocene-Pleistocene asymmetric unroofing of the Irian fold belt, Irian Jaya, Indonesia: Apatite fission-track thermochronology: Geological Society of America Bulletin, vol. 108, pp. 1438–1449, doi:10.1130/0016-7606(1996)108(1438:ppauot)2.3.co;2.

789 Witts, D., Hall, R., Nichols, G., and Morley, R., 2012, A new depositional and provenance model for the Tanjung Formation, Barito Basin, SE Kalimantan, Indonesia: Journal of Asian Earth Sciences, vol. 56, pp. 77–104, doi:10.1016/j.jseaes.2012.04.022.

792 Yumul, G. P., 2007, Westward younging disposition of Philippine ophiolites and its implication for arc evolution: Island Arc, vol. 16, pp. 306–317, doi:10.1111/j.1440-1738.2007.00573.x.

794 Yumul, G. P., Dimalanta, C. B., Tamayo, R. A., and Faustino-Eslava, D. V., 2013, Geological features of a collision zone marker: The Antique Ophiolite Complex (Western Panay, Philippines): Journal of Asian Earth Sciences, vol. 65, pp. 53–63, doi:10.1016/j.jseaes.2012.08.017.