Landscape dynamics and the Phanerozoic diversification of the biosphere – Nature

Global landscape model

Here we use goSPL^12,13, an open-source scalable parallel numerical model that simulates landscape and sedimentary basin evolution at the global scale (resolution about 5 km). It accounts for river incision and soil creep, both considered as the main drivers of long-term physiographic changes. goSPL also tracks eroded sediments from source to sink considering alluvial and marine deposition, sediment compaction and porosity change, and could be used to reconstruct basin stratigraphic records. To evaluate these processes on landscape dynamics, different forcing conditions could be imposed from spatially and temporally varying tectonics (both horizontal and vertical displacements) to multiple climatic histories (for example, precipitation patterns and sea-level fluctuations). The model’s main equation, the continuity of mass, has the following common form:

for which changes in surface elevation z with time t depend on tectonic forcing U (rock uplift rate, in metres per year), hillslope processes for which κ is the diffusion coefficient (set to 0.5 m² yr⁻¹)⁵⁷ and fluvial processes defined using the stream power law. m and n are dimensionless empirical constants (set to 0.5 and 1), ε is a precipitation-independent component of erodibility (set to 4.0 × 10⁻⁷ yr⁻¹ on the basis of the choice of m), and PA is the water flux combining upstream total area (A) and local runoff (P) obtained from palaeoclimate mean annual precipitation minus evapotranspiration⁵⁷. In our formulation, the weathering impact of runoff and its role on river incision enhancement is incorporated by scaling the erodibility coefficient with local mean annual runoff rate with a prefactor d (a positive exponent estimated from field-based relationships⁵⁸ and set to 0.42). It follows from equation (1) that deposition in flat plains or along gentle slopes is null. However, it simulates continental deposition in depressions and endorheic basins.

In goSPL, erosion occurring in upstream catchments is linked to basin sedimentation through a multiple-flow-direction algorithm that routes both water and sediment flux towards multiple downstream nodes, preventing the locking of erosion pathways along a single direction and helping the distribution of the corresponding flux in downstream regions. To solve the flow discharge globally (PA), we use a parallel implicit drainage area (IDA) method^59,60 in a Eulerian reference frame, expressed in the form of a sparse matrix composed of diagonal terms set to unity and off-diagonal terms corresponding to the immediate neighbours of each vertex composing the spherical mesh. The solution of the IDA algorithm is obtained using the Richardson solver with block Jacobian preconditioning⁵⁹, both available in PETSc⁶¹. Continental erosion and sediment transport solutions follow a similar approach⁶⁰. Some of the main advantages of goSPL lie in its design of implicit and parallel solutions of its constitutive equations⁶⁰, making it possible to increase the model stability even with large time steps, and to scale the simulation run time over hundreds of CPUs.

Palaeo-elevation and precipitation forcing

To reconstruct the past physiography, goSPL relies on time-evolving boundary conditions—climatic and palaeogeographic—that are used to compute the interplay between the solid Earth and the climate. To reconstruct high-resolution palaeo-elevations throughout the Phanerozoic, we use a series of 108 maps from the PALEOMAP palaeogeographic atlas¹⁴ ranging from the Holocene epoch to the Cambrian–Precambrian boundary (541 Ma). These palaeo-maps are defined at approximately 5-Myr intervals, and each of them is represented as a regular grid with a resolution of 0.1° × 0.1° (approximately 10-km cell width at the Equator). These palaeogeographic maps were initially based on information related to lithofacies and palaeoenvironmental datasets⁶² and supplemented and refined for more than 40 years with regional palaeogeographic atlases^14,63. We acknowledge that these maps bear some uncertainties and controversial aspects. For example, the very early Andean rise to their modern elevations consequently precociously lowers the predicted bulk sediment flux during the Neogene period while diversification continues to increase (Figs. 3 and 4). It is worth noting that even though the PALEOMAP dataset forms the basis of this study, from a methodology standpoint, other datasets^64,65,66 could be used.

To simulate riverine processes, the palaeo-precipitation dataset used was generated using a variant, HadCM3BL-M2.1aD (ref. ¹⁵), of the coupled atmosphere–ocean–vegetation Hadley Centre model. This climate model also uses the PALEOMAP Atlas⁶⁷ but at a lower resolution (3.75° × 2.5°). The reconstructed palaeo-precipitation regimes are obtained for each individual palaeo-elevation map after running the climate model for at least 5,000 model years to reach a dynamic equilibrium of the deep ocean¹⁵. In addition to palaeo-elevation grids, there are two additional time-dependent boundary conditions that were set in the climate model: the solar constant; and the atmospheric CO₂ concentrations. Regarding the last condition, two alternative CO₂ estimates are proposed¹⁵ and we chose the palaeo-precipitation and evapotranspiration maps generated from the set of HadCM3 climate simulations using the CO₂ local weighted regression curve from ref. ⁶⁸.

From the palaeo-elevation and palaeo-runoff maps, we generate the input files for goSPL by resampling the global temporal grids on an icosahedral mesh composed of more than 10 million nodes and 21 million cells (corresponding to an averaged resolution of about 5 km globally—about 0.05° resolution at the Equator). Inspired by techniques used in palaeoclimate modelling^15,67, we designed an approach to achieve a dynamic equilibrium (erosion rates balance tectonics; equation (1) and Extended Data Fig. 2) under steady boundary conditions (rainfall, tectonic uplift and erodibility). For each individual time slice, we run two sets of simulations over 168 CPUs to estimate their corresponding physiographic characteristics and associated water and sediment dynamics (Extended Data Fig. 1). A first simulation is carried out over an interval of 2 Myr under prescribed elevation and runoff conditions and simulates landscape evolution and associated water discharge and sediment flux assuming no other forcing. Under this setting, the role of surface processes is not counterbalanced by tectonics, and they excessively erode the reconstructed elevation, trimming part of the major long-lived orogenic belts and upland areas, and causing extensive floodplains. The resulting elevations are then corrected by assimilating the palaeo-elevation information¹³. Model predictions account for landscape evolution, and at this stage already contain a more detailed representation of terrestrial landforms (for example, canyons, valleys, incised channels and basins to cite a few) than the initial palaeo-elevation (Extended Data Fig. 2a). To preserve these morphological features during the correction step, high-amplitude and high-frequency structures are removed using a combination of moving average windows (ranging from 0.5° to 2°) that conserves the global distribution of the initial palaeo-elevation with minimal change of its hypsometry (≤0.5%; Extended Data Fig. 2c). We then derive a tectonic map (uplift and subsidence rates) by computing the local differences between the palaeo-elevation values and the adjusted ones.

A second simulation starts with previous palaeo-elevation and runoff conditions and additionally accounts for tectonic forcing. This simulation runs until dynamic equilibrium is reached (that is, erosion rates compensate tectonic uplift rates) within the first million years of landscape evolution (Extended Data Fig. 2c). The outputs of this second simulation are then used to evaluate water and sediment flux for the considered time slice, as well as the major catchment characteristics (river networks, drainage areas, and erosion and deposition rates).

The parametrization of equation (1) is obtained by calibrating its variables with modern estimates of average global erosion rates¹⁸ (mean value of 63 m Myr⁻¹ with a standard deviation of 15 m Myr⁻¹; Extended Data Fig. 6c) and of suspended sediment flux from the BQART model^19,20 (corresponding to 12.8 Gt yr⁻¹). Following calibration, we predict an average present-day global erosion rate of 71 m Myr⁻¹, and a sediment flux of 12.15 Gt yr⁻¹ (assuming an average density of 2.7 g cm⁻³). On the basis of the multiple-flow IDA approach used to integrate runoff over upstream catchments⁵⁹ (IDA algorithm), we also extract the spatial distribution of the largest water discharges and sediment flux (Extended Data Fig. 3a) and their respective basin characteristics.

Sr isotopic variations from mantle origin

We use the geochemical archive of oceanic sediments to test the validity of the model predictions. The ⁸⁷Sr/⁸⁶Sr isotopic ratio of seawater (Extended Data Fig. 4) reflects the balance between continental weathering and mantle dynamics (hydrothermalism at mid-ocean ridges and weathering of island arcs and oceanic islands)^69,70, making it a classic proxy to diagnose the relative importance of geodynamic and climatic forcings through time^71,72.

From the measured isotopic budget of the ocean (O), present-day low ⁸⁷Sr/⁸⁶Sr ratios (about 0.703) are produced from mantle sources (M), whereas high ratios (about 0.713) come from continental runoff (CR, measured from main rivers worldwide)⁶⁹. As a result, the strontium isotope oceanic mass balance has the following form:

in which ξ represents the mass fraction of the Sr coming from the mantle (Q_M/(Q_M + Q_S) with Q_S being the predicted net sediment flux to the ocean derived from our simulation) and the flux of mantle origin (Q_M). At the present day, Q⁰_M is given by the percentages defined above:

Assuming that weathering scales with erosion rates, our reconstructed global net sediment flux to the ocean (Q_S) offers an independent alternative to existing approaches evaluating Sr flux from tectonic origin and could be used to infer past tectonic activity^73,74,75,76. Under this assumption, differences (∆(⁸⁷Sr/⁸⁶Sr)) between our predicted Sr ratio and the one from the geological record³⁹ would reflect the dynamics of the Earth’s exogenic system, with positive ∆(⁸⁷Sr/⁸⁶Sr) values corresponding to periods of higher tectonic activities relative to the present day and negative ones coinciding with more quiescent periods. We then use the isotope oceanic mass balance to independently derive the mantle contribution to the ⁸⁷Sr/⁸⁶Sr ratio, relying on our estimates of terrigenous flux and on the observed ⁸⁷Sr/⁸⁶Sr ratio of seawater. We crudely use the Sr isotope oceanic mass balance to estimate the Sr flux of mantellic origin (Q_M; Extended Data Fig. 4) in the mass fraction ξ:

Note that we assume that the isotopic ratios \({\left[\frac{{}^{87}{\rm{Sr}}}{{}^{86}{\rm{Sr}}}\right]}_{{\rm{CR}}}\) and \({\left[\frac{{}^{87}{\rm{Sr}}}{{}^{86}{\rm{Sr}}}\right]}_{{\rm{M}}}\) remained equivalent to those of the present day. Whereas the Sr ratio from the mantle budget might change marginally (⁸⁷Sr/⁸⁶Sr about 0.703 for mid-ocean ridge hydrothermal and about 0.7035 for island arcs and oceanic islands⁶⁹), the contribution of the continental crust to the Sr ratio is highly dependent on the type of weathered materials⁷⁰ (⁸⁷Sr/⁸⁶Sr about 0.708 for limestones, compared to about 0.721 for silicates⁶⁹).

We find that the contributions of mantle and terrigenous sources relative to those of the present day covary at long wavelengths, with two periods of reinforced influx from both sources separated by a quieter period during Pangaea. This trend mirrors the Wilson cycle. It indicates that the periods of high erosion, coinciding with periods of continental aggregation and contraction, increased elevation and wetter climates, also match periods of reinforced tectonic activity. Seafloor expansion is faster during periods of continental dispersal, and the total length of ridges increases during breakup⁷⁷. As the mantle input to the ⁸⁷Sr/⁸⁶Sr ratio of seawater is partially driven by seafloor kinematics, we find several similarities between the predicted mantle Sr flux and subduction rates^72,74,76 that can be taken as a proxy for oceanic spreading rates; Extended Data Fig. 4). Notably, over the past 250 Ma, we deduce that mantle fluxes were low during Pangaea, and subsequently increased during Pangaea breakup; flux decreased during the late Palaeogene, mirroring the decrease in crustal production rates in the Atlantic and Pacific oceans⁷⁸ and the consumption of the East Pacific ridge.

Our model predictions of sediment flux compare at first order to the observed increase in the ⁸⁷Sr/⁸⁶Sr ratio of seawater over the past 150 Myr (Extended Data Fig. 4). Assuming that weathering scales with erosion rates, it corroborates the first-order impact of Cenozoic orogenesis⁷⁹ on the Sr ratio. Likewise, over the entire Phanerozoic, short-lived (20–40 Myr) increases in predicted erosion flux can explain the increase in the ⁸⁷Sr/⁸⁶Sr record during major orogenic phases⁸⁰, whereas the long-term variations of the record can be at first order explained by the coevolution of the terrigenous and mantle sources during the Wilson cycle. These results indirectly substantiate our model predictions of sediment flux to the oceans.

Limitations and sensitivity

In goSPL¹², erosion is defined using a first-order parametrization of the physics at play (equation (1)), which captures the long-term, large-scale landscape evolution^22,57,81. More sophisticated treatments directly linked to sediment transport theory and incorporating different erosional behaviours have been proposed (for example, by accounting for mobile sediment and bed resistance to erosion, or using different formulations of sediment transport based on transport-limited equations)^82,83,84.

In addition, the erodibility parameter does not consider spatiotemporal variations that could be induced by environmental (for example, temperature gradients and seasonal precipitation), geological (for example, soil composition and fault-induced bedrock weakening) or biological (for example, plant root growth and soil microbial activity) mechanisms^{6,85,86,87,88,89}. Instead, we assume uniform erodibility across all continents. Accounting for variable lithologies in model space could be achieved by assigning an erodibility prefactor depending on the surface rock composition in the stream power law term of equation (1) by an erodibility prefactor depending on the type of surficial lithology classes (typically with values ranging between 1.0 and 3.2; ref. ⁹⁰). However, this approach would require global palaeo-lithological surficial cover data that are difficult to obtain when looking into deep geological times. Although we do not account for seasonality variations, the weathering impact of precipitation and its role on river incision enhancement is incorporated by scaling the erodibility coefficient with the local mean net annual precipitation rate⁵⁸. One could also incorporate the effect of temperature on weathering of rocks according to rock composition by scaling the erodibility coefficient using the palaeoclimate temperature distribution. Such refinement possibilities are many, and although in principle desirable to better reproduce observations, adding those would necessarily add poorly controlled degrees of freedom in the parameterization, and lead to illusory enhanced predictive capabilities.

By design, our simulation is sensitive to both the climatic and palaeo-elevation reconstructions. Although other palaeogeography reconstructions exist^62,64,65,66, many are restricted to specific geological intervals and, to our knowledge, are not tied to a series of palaeo-precipitation maps for the entire Phanerozoic. Consequently, we chose to evaluate model sensitivity on palaeo-elevation using a single time slice (Aptian period about 120 Ma) by comparing our predicted sediment flux with a different set of palaeogeography and palaeoclimate reconstructions⁹¹. The results highlight several differences at the regional scale (Extended Data Fig. 5). For instance, a more humid equatorial climatic zone in the second reconstruction⁹² induces higher erosion rates on the northern and central part of Gondwana. The palaeo-elevation differences also redefine the drainage network and the volumes of sediment transported to the oceans or stored in continental basins. This is the case in Antarctica and on the eastern part of Eurasia where we note higher erosion rates or an increase in terrestrial sediment accumulation depending on the considered palaeogeography. Despite these local variations, those disparities become more tenuous when evaluating the global response. As an example, the percentage of endorheic basins varies from 24 to 26.5% between the two simulations and the net sediment flux to the ocean changes from 2.72 to 2.26 km³ yr⁻¹ (Extended Data Fig. 5). This suggests that although regional differences exist and if the imposed forcing conditions are not too dissimilar (both in terms of palaeoclimates and palaeogeographies), the reported global evolution and global trends that are presented in the study should remain relatively unchanged between reconstructions.

To evaluate the model sensitivity to palaeo-runoff, we ran a full series of simulations throughout the Phanerozoic. The palaeoclimate reconstructions from ref. ¹⁵ have been carried out with two different CO₂ conditions (the atmospheric concentrations from ref. ⁶⁸ and a smoothed curve), but the mean global continental runoff remains very similar in both cases; Extended Data Fig. 6a) and should only very marginally change our results. Instead, we opted for the recent Phanerozoic palaeo-precipitation reconstruction from ref. ⁹³ that was run at 10-Myr intervals using the PALEOMAP palaeogeographic atlas¹⁴ with a much higher resolution (1°) than those of ref. ¹⁵. The release from ref. ⁹³ does not contain the evapotranspiration time slices, and we could therefore use the total palaeo-precipitation only as a proxy for runoff. Global mean continental runoff from ref. ⁹³ exhibits a similar temporal trend to the ones from ref. ¹⁵ but, because evapotranspiration could not be accounted for, with higher values (about 0.3 m yr⁻¹ on average over the Phanerozoic; Extended Data Fig. 6a). As erosion scales with runoff (equation (1)), this inflated runoff directly affects the global net sediment flux to the ocean (about 0.84 km³ yr⁻¹ on average) and to a lesser extent the continental deposition flux (<0.1 km³ yr⁻¹ on average; Extended Data Fig. 6b). The spatial distributions of these two runoff scenarios and their relative impact on denudation rates show substantial spatial differences over time (Extended Data Fig. 7). At the continental scale, the higher resolution in ref. ⁹³ should better account for the control of topography on the spatial variability in precipitation. For example, we note at 40 Ma (Extended Data Fig. 7) the orographic control of the Himalayas on the regional rainfall regime with its implications on erosion and deposition.

The sensitivity analysis provides two important pieces of information. First, both simulations show similar responses in terms of global sediment flux and denudation rates (Pearson correlations of 0.94 and 0.92, respectively). This suggests that irrespectively of the chosen palaeoclimatic reconstruction, our interpretation of the relationships between predicted sedimentary flux and biodiversity still holds. Second, the runoff has a stronger effect on the amplitudes in net sediment flux to the ocean and denudation rates when considering similar palaeo-elevation reconstruction (Extended Data Fig. 6b). Instead, differences in palaeo-elevations affect continental sediment cover and sediment accumulation in endorheic basins (Extended Data Fig. 5).

Another limitation of our approach is to propose a hypothesis by comparing time series of mean model outputs with independent variables, but similar trends could possibly be expected for any model that accounts for plate tectonics—with essentially two cycles of continental aggregation and dispersal over the Phanerozoic eon—and subsequent climate reconstructions. The highly relevant correlations that we found can be however advocated to hierarchize these studies. As plate tectonics—and the breakup of Pangaea in particular—form the cornerstone of such studies^11,40,41, we thus compared continental fragmentation⁴¹ with Phanerozoic marine diversity, which yields only moderate correlations (Pearson coefficients of 0.46 with number of marine families²; 0.54 to 0.58 with the number of genera^2,5; Extended Data Fig. 8).

Physiographic diversity and variety

To evaluate the relationships between physiography and the Phanerozoic climate and tectonics, we define a unique continuous variable (named the physiographic diversity index) that encapsulates several of the reconstructed morphometric attributes. Simulation outputs are first interpolated from the icosahedral mesh on a regular 0.05° grid (open-access online release, HydroShare¹⁶).

First, we quantify palaeo-landforms by calculating the topographic position index on each cell i (TPI_i) that measures the relative relief⁹⁴:

in which z_i is the considered elevation at cell i and z_k is the mean of its surrounding cells (z_k), with n being the number of cells contained inside an annulus neighbourhood. Topographic position is an inherently scale-dependent calculation⁹⁵. To circumvent this problem, TPI is computed considering two scales of observation, a fine one ranging between 0.05° and 0.15° and a coarser one from 0.25° to 0.5°. As elevation is generally spatially autocorrelated, TPI values increases with scale, making it difficult to compare both scales of observation directly. To overcome this issue, we calculate a standardized TPI_S (equation (5)), in which \(\overline{{\rm{TPI}}}\) is the mean over the entire grid and σ_TPI is its standard deviation⁹⁶ (top map in Extended Data Fig. 9a).

We also extract the slope and water discharge for each time slice (bottom maps in Extended Data Fig. 9a). Note that we selected these three variables—TPI, slope and discharge—as morphometric indicators of the physiographic diversity but other ones such as aspect (that is, the direction of maximum gradient), and erosion and deposition rates could equally be added⁹⁷. From these continuous variables, we then derive categorical variables by defining seven categories from the slope, five from the water flux and ten from the TPI_S (Extended Data Fig. 9a,b and Supplementary Table 1; ref. ⁹⁸).

From the categorical variables, we quantify the physiographic diversity index P_DIV (Supplementary Fig. 1a,b) using Shannon’s equitability (continuous variable [0,1]), which is calculated by normalizing the Shannon–Weaver diversity index (d_SW):

with p_k being the proportion of observations of type k in each neighbourhood and \({\mathcal{C}}\) being the number of categorical variables (here \({\mathcal{C}}\) = 3). The physiographic diversity index is calculated at each location for the 108 reconstructed time slices spanning the Phanerozoic (Supplementary Fig. 1c). The variations of P_DIV can be described either spatially for a given period (Supplementary Fig. 1c) or across time (Supplementary Fig. 2). As it characterizes how landscape complexity evolves over the geological timescale⁹⁹, this index can be used to infer the role physiography, and changes thereof, might have played on species migration, dispersal or isolation at both global and continental scales. Biodiversity richness emerges from many abiotic and biotic interactions; however, the role physiography might have played has often been overlooked¹⁰⁰. The high-resolution global maps generated from our simulation could be used as an additional palaeoenvironmental layer in mechanistic eco-evolutionary models^101,102.

To compare the temporal evolution of physiographic diversity with biotic^2,5 or abiotic geochemical and geophysical proxies, the mean value of P_DIV is insufficient, as it ignores the variety of landforms, which is better accounted for by the probability density function: the probability density function can be usefully reduced to a unique scalar—physiographic variety—given by the interquartile range P_VAR = Q₃ − Q₁ (that is, the range between the first quartile (Q₁; 25%) and the third quartile (Q₃;75%); Supplementary Fig. 2c).

As an example, we observe an increase in physiographic diversity of South America from the Upper Cretaceous to the Palaeocene related to the Andean mountain building period (Supplementary Fig. 2b,c) and concomitant with microflora diversity patterns in northern South America¹⁰³ and plant diversification in Patagonia¹⁰⁴. We also note periods of low P_VAR values around 50 Ma and 28 Ma related to the intermittent development of internal seas or lakes in the central part of the Amazon basin. Periods of increasing P_VAR follow each of these episodes, suggesting that the Pebas basin could have acted as a permeable biogeographical system favouring biotic exchange and adaptation in the region¹⁰⁵.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.