Rates of seafloor and continental weathering govern Phanerozoic marine phosphate levels

Phosphate is an essential macronutrient for all organisms with a key role in setting levels of marine primary productivity. Despite its importance for marine biogeochemical cycles and its role in shaping the evolution of marine organisms, the factors controlling phosphate bioavailability on geologic timescales remain poorly understood. Here we develop a statistical model of the coupled cycles of phosphate, carbon, oxygen and calcium to constrain the weathering-derived fluxes and seawater concentrations of phosphate through Phanerozoic time (541 million years ago to the present). Our model includes input parameters and time-dependent forcings derived from geologic and geochemical data. We find that the climate sensitivity of chemical weathering of the oceanic crust by low-temperature fluids exerts a first-order control on phosphate availability. Specifically, continental weathering is a source of the limiting nutrient phosphate, but seafloor weathering is considered to be a minor phosphate sink. Consequently, times in Earth history during which seafloor weathering constituted a large fraction of the total (seafloor + continental) weathering were also times during which phosphate influxes to and concentrations in the ocean were relatively low. Lower seawater phosphate levels during those times probably resulted in lower primary productivity and oceanic and atmospheric oxygen concentrations. Marine phosphate levels and biological productivity were lowest during the early Phanerozoic when seafloor weathering rates were high and continental weathering rates were muted, according to a statistical model of coupled elemental cycles.

Phosphate is an essential macronutrient for all organisms with a key role in setting levels of marine primary productivity. Despite its importance for marine biogeochemical cycles and its role in shaping the evolution of marine organisms, the factors controlling phosphate bioavailability on geologic timescales remain poorly understood. Here we develop a statistical model of the coupled cycles of phosphate, carbon, oxygen and calcium to constrain the weathering-derived fluxes and seawater concentrations of phosphate through Phanerozoic time (541 million years ago to the present). Our model includes input parameters and time-dependent forcings derived from geologic and geochemical data. We find that the climate sensitivity of chemical weathering of the oceanic crust by low-temperature fluids exerts a first-order control on phosphate availability. Specifically, continental weathering is a source of the limiting nutrient phosphate, but seafloor weathering is considered to be a minor phosphate sink. Consequently, times in Earth history during which seafloor weathering constituted a large fraction of the total (seafloor + continental) weathering were also times during which phosphate influxes to and concentrations in the ocean were relatively low. Lower seawater phosphate levels during those times probably resulted in lower primary productivity and oceanic and atmospheric oxygen concentrations.
By limiting marine primary productivity and organic carbon (C) burial over geologic timescales, the long-term availability of phosphate (P) is a key determinant of atmospheric oxygen and carbon dioxide levels ( p O 2 and p CO 2 , respectively) 1-4 , climate 5 and marine faunal evolution and diversity patterns 6 . On geologic timescales (hundreds of thousands to millions of years), the size of the marine P reservoir is set by a balance between supply from continental weathering and loss mostly by burial in sedimentary rocks 7-9 . On timescales comparable to the mixing time of the ocean (thousands of years), the given marine P reservoir is distributed between the surface and deep ocean by a balance between P incorporation into newly fixed organic matter and P release during organic-matter remineralization 5 . The long-term source of P to the ocean is chemical weathering of continental silicate rocks [8][9][10] . Seafloor weathering, however, is considered to be a source of alkalinity only and a minor sink of P 11-13 . Away from seafloor-spreading axes, low-temperature (generally 5-20° C) water circulates in and reacts with the highly porous upper oceanic crust 14 . Pore-water and bulk-rock chemical analyses suggest that P is removed from seawater during these off-axis water-rock interactions, probably due to adsorption onto iron (oxy)hydroxides, secondary apatite precipitation and microbial activity in the crust 11 .
Rates of continental and seafloor weathering depend on climate. A higher surface temperature increases rates of continental weathering, both through temperature-dependent reaction kinetics 15-18 and through increased precipitation and run-off 19-21 . Off-axis weathering of the oceanic crust depends on deep-water temperatures 22-24 , which themselves depend on surface climate (Methods) 25,26 . The climate dependence of both continental and seafloor weathering has led to suggestions of negative climate feedbacks 15,22,24,[27][28][29][30][31] . According to these feedbacks, atmospheric CO 2 levels reach steady-state conditions when volcanic Article https://doi.org/10.1038/s41561-022-01075-1 weathering-derived alkalinity influx) stems from lower continental weatherability before the evolution of land plants. Land plants enhance continental weatherability in several ways, for example, by exuding organic acids and metal chelators, by increasing the surface area for weathering due to their extensive root systems and by recirculating water through evapotranspiration (Supplementary Information section 2.1) 33,34 . The lower continental weatherability results in a balance between the volcanic CO 2 source and the weathering CO 2 sink (comprising CaCO 3 and organic C burial) at a higher global-average temperature (Extended Data Fig. 3). If the climate sensitivity of seafloor weathering is not accounted for, as in some previous models, lower continental weatherability can be compensated only by changes in climate-dependent rates, notably the continental weathering rate. In a scenario of lower continental weatherability and climate-insensitive seafloor weathering, our model generates an early Palaeozoic steady-state climate sufficiently warmer that increased rates of continental weathering drive a combined carbonate and organic C burial flux equal to CO 2 outgassing (dark red lines in Fig. 2). However, accounting for the climate sensitivity of seafloor weathering, a warmer model climate accelerates both continental and seafloor weathering, and a balance between CO 2 sources and sinks is achieved at a lower global-average surface temperature than with climate-insensitive seafloor weathering. Furthermore, with low pre-land-plant continental weatherability, seafloor weathering becomes the dominant source of alkalinity to the ocean in most model simulations (colour contours on Fig. 2). As land plants expanded, ∼400-350 Ma, continental weatherability increased, and a balance between CO 2 inputs and outputs was reached at a lower surface temperature. This resulted in slower temperature-dependent kinetics of both seafloor and continental weathering, although net continental weathering rates increased due to the higher plant-related weatherability (Fig. 2).
During the assembly of the supercontinent Pangaea (~306 and 237 Ma), the average annual precipitation over land is thought to have decreased to a mid-Triassic minimum (Extended Data Fig. 2) due to the difficulty of delivering moisture to the continental interior (Supplementary Information section 2.2) 35 . Lower precipitation results in lower continental weathering rates at a given temperature, requiring a warmer climate for seafloor and continental weathering rates to match the CO 2 outgassing rates. As in the early Palaeozoic, the warmer climate accelerates both continental and seafloor weathering. Again, due to the low continental weatherability, which in this case is CO 2 outgassing is balanced by carbonate and organic C burial. These burial rates depend on weathering-derived alkalinity and nutrient fluxes to the oceans, respectively. Since continental silicate weathering is a source of both alkalinity and P, and seafloor weathering is only a source of alkalinity and a minor sink of P, the relative rates of seafloor and continental weathering influence the P cycle. For example, a gradual increase in p O 2 over 1,500-500 million years ago (Ma) has been suggested to reflect a gradual increase in the relative proportion of continental weathering, shifting the balance between the inorganic and organic C burial sinks 32 . The consequences of such dynamics for the Phanerozoic P cycle, to our knowledge, have not been fully explored.
In this Article, to study the evolution of the Phanerozoic P cycle, we developed a model for the coupled, long-term biogeochemical cycles of P, C, O 2 and Ca ( Fig. 1) and simulated the fluxes of these elements between the ocean-atmosphere and rock reservoirs in response to major geologic forcings (see Supplementary Information section 1 for flux parameterizations). The model outputs (Extended Data  Table 1) are the marine and sedimentary pools of these elements (for example, the P concentration) and their time-varying fluxes among the reservoirs (for example, organic C burial). Our model includes 35 input parameters (Extended Data Table 2 and Extended Data Fig. 1) and 11 time-dependent forcings (Extended Data Table 3 and Extended Data Fig. 2). The former include parameters such as volcanic CO 2 fluxes into the ocean-atmosphere and the C/P ratio of sediments in various settings. The latter include forcings related, for example, to palaeogeography, land cover and land-plant evolution (Supplementary Information section 2). To account for parameter and time-dependent forcing uncertainties, we randomly selected inputs across realistic parameter space and validated the results using geologic proxies for p CO 2 , the C content of oceanic crust and the C isotopic composition (δ 13 C) of carbonate minerals (Methods). The statistical distribution of valid model outputs constrains the likely marine P budget throughout the Phanerozoic.

Continental and seafloor weathering through the Phanerozoic
Our analysis suggests that during the early Palaeozoic era, continental weathering rates were low and seafloor weathering rates were high relative to their present-day values for most successful model runs (colour contours on Fig. 2). The lower contribution of continental weathering-derived alkalinity influx to the ocean (out of the total  Article https://doi.org/10.1038/s41561-022-01075-1 related to lower precipitation over land, seafloor weathering becomes the dominant source of alkalinity to the ocean in most model runs (Fig. 2b). By contrast, during the break-up of Pangaea, which started ~237 Ma, continental weathering rates increased and seafloor weathering decreased to present-day levels (Fig. 2). We emphasize that, as in the early Palaeozoic, seafloor weathering can become the dominant source of alkalinity to the ocean only when its climate dependence is accounted for in the model. Our proposed geologic history of seafloor weathering rates is consistent with observed changes in the CO 2 content of the oceanic crust 14,22,36,37 . During seafloor weathering reactions, Ca is leached from the basalts, increasing the Ca concentration in crustal pore fluids and promoting CaCO 3 precipitation in veins and void fillings 14,15,22,36 . The rate of C addition to altered oceanic crust, as inferred from the measured CO 2 content of the crust 22 , was higher in the Mesozoic than over Cenozoic time (Fig. 2b markers and Extended Data Table 4). Measurements of 87 Sr/ 86 Sr in the crustal carbonate minerals suggest that the ~fivefold C enrichment in Mesozoic drill cores results not from longer accumulation times but from faster temperature-dependent kinetics of seafloor weathering 22 . Our model trajectory with climate-sensitive seafloor weathering rates matches the available Mesozoic-Cenozoic observations, providing confidence in this trajectory over Phanerozoic time.

The Phanerozoic phosphate budget
Because continental weathering is a source of both alkalinity and P, whereas seafloor weathering is a source of alkalinity and a minor sink of P, the variation in the relative importance of continental and seafloor weathering influences the marine P budget. Our results (of ∼10 6 random parameter draws; Fig. 3) indicate that the high relative importance of seafloor weathering in the early Palaeozoic (Fig. 2b) led to weathering-derived P fluxes and deep-ocean P concentrations that were lower than present by a factor of ~2 in most model runs. When land plants evolved (~400-350 Ma, Supplementary Information section 2.1) and continental weatherability increased, weathering-derived P influxes and seawater concentrations increased in response to the higher relative importance of continental weathering.  in lower weathering-derived P influxes and seawater concentrations. By contrast, since the onset of the break-up of Pangaea (~237 Ma), P delivery rates and concentrations in seawater increased towards their present-day values. Weathering-derived P influxes and the seawater P concentrations depend on assumptions about the climate sensitivity of seafloor weathering rates. Neglecting the climate sensitivity of seafloor weathering rates, only alkalinity influxes from continental weathering may respond to changes in CO 2 outgassing or continental weatherability through a change in the global-average surface temperature (Fig. 4 right). As continental weathering is also a source of P, in this case, higher CO 2 influxes translate directly to higher P influxes, and vice versa. For this reason, with climate-insensitive seafloor weathering rates, the Palaeozoic weathering-derived P influxes and seawater concentrations are higher than with climate-sensitive seafloor weathering, and they closely follow the temporal evolution of the CO 2 outgassing rate (dark red envelope on Fig. 3). By contrast, with climate-sensitive seafloor weathering, the climate response to a change in CO 2 outgassing or continental weatherability is muted as a smaller change in the global-average surface temperature is necessary to restore mass balance in the C and alkalinity cycles through the silicate weathering feedback. In this case, higher CO 2 influxes translate into an increase in both continental and seafloor weathering rates, an increase in overall P delivery to the ocean, but a decrease in the ratio of P to alkalinity delivery ( Fig. 4 lower left). Lower continental weatherability translates into an increase in the relative contribution of seafloor weathering to the total alkalinity flux and an actual reduction in the P influx ( Fig. 4 upper left).
Next, we compare the trajectories predicted by our model with available estimates of P weathering fluxes and seawater P concentrations. Our predicted history of a ~twofold increase in P weathering fluxes differs from a recent prediction based on the 87 Sr/ 86 Sr record in marine carbonates, which has been interpreted to suggest a Phanerozoic decrease in P influxes by a factor of ~1.5 (ref. 38 ). We note that this contrasting prediction was made under an assumption that seawater 87 Sr/ 86 Sr reflects a simple mass balance between radiogenic and non-radiogenic strontium inputs from the continents and mid-ocean ridges, respectively. However, it has been shown that to derive weathering fluxes from the carbonate 87 Sr/ 86 Sr record, knowledge of continental weatherability, lithology, palaeogeography and off-axial seafloor weathering fluxes is necessary 29,39 , and these factors were unaccounted for in the aforementioned prediction. To our knowledge, there is a limited number of proxies for the marine P concentration. An increase in P/Fe ratios in iron-oxide-rich marine sediments, from a pre-Cretaceous average of 0.38 to a Cretaceous-Cenozoic average of 2.55, has been interpreted to reflect a change in the silica cycle rather than an increase in the seawater P concentration 4 . Specifically, a decrease in the seawater silica concentration due to the expansion of silicifying phytoplankton (diatoms) was inferred to have reduced the competition with P for adsorption sites on the surface of iron oxides, resulting in enhanced P adsorption and higher P/Fe without a major change in the seawater P concentration. However, the Phanerozoic silica cycle is poorly constrained 40   The heights of the bars represent the relative magnitude of the following steady-state fluxes: CO 2 outgassing (grey), alkalinity (Alk) input from continental weathering (red), alkalinity input from seafloor weathering (blue) and P input from continental weathering (pink). From an initial steady state (bars in the centre of the figure), a change in continental weatherability or in CO 2 outgassing rates ultimately leads to different steady-state fluxes through the action of climate feedbacks. The lighter-toned bars represent the difference between the initial and final steady-state fluxes. At lower continental weatherability (for example, before the evolution of land plants or during the formation of the supercontinent Pangaea; upper panels), or at higher CO 2 outgassing (lower panels), atmospheric CO 2 levels and global temperature increase until steadystate alkalinity fluxes to the ocean from weathering match the CO 2 input. If both seafloor and continental weathering rates are sensitive to climate (as in reality), then both rates increase (lighter-toned bars in the left panels). If the climate sensitivity of seafloor weathering is not accounted for, as in some previous models, only continental weathering rates increase (blue bar does not change in the right panels). Since the P flux is proportional to continental weathering rates, the change in the P influx is smaller when seafloor weathering is climate sensitive relative to the case of climate-insensitive seafloor weathering (compare the heights of the pink bars in the left and right panels).
Article https://doi.org/10.1038/s41561-022-01075-1 P/Fe. Instead, at least part of the P/Fe increase in sedimentary rocks may reflect the Triassic-Cretaceous increase in seawater P concentrations predicted by our model in response to the break-up of Pangaea (Fig. 3). Indirect evidence for an increase in P availability through the Phanerozoic is suggested by the temporal succession of dominant phytoplankton taxa. Green algae and cyanobacteria dominated the early Palaeozoic, followed by dinoflagellates and coccolithophorids in the Mesozoic and Palaeogene and by an explosive diversification of diatoms during the Neogene 41,42 . In cultures of modern representatives of these major groups, taxa from groups that emerged later in Earth history tend to prefer progressively more P-rich conditions 44 and tend to produce biomass richer in P (refs. [45][46][47] ). If these nutrient requirements and utilization patterns are conserved group-specific traits, they suggest that P availability increased over the Phanerozoic 6 .

Implications for the carbon and oxygen cycles
Predicted marine primary productivity (TmolC yr -1 ) tracks the dissolved surface P concentration (Fig. 5a, b), as primary productivity is calculated from the surface-ocean P concentration and the elemental composition of the primary producers (here assumed to be constant at the Redfield ratio; Supplementary Information section 1.7). Accordingly, early Palaeozoic primary productivity is ~2 times lower than today (Fig. 5b). Although the burial of marine-derived organic C is proportional to marine primary productivity, early Palaeozoic organic C burial is lower than today by only ~30% (Fig. 5c). This muted response to lower-than-present primary productivity is related to lower early Palaeozoic p O 2 (Fig. 5e, f), which led to a higher marine organic C burial efficiency 48 . Lower early Palaeozoic p O 2 stems from the absence of land plants and the associated terrestrial organic C burial flux (Fig. 5d), which has C/P approximately an order of magnitude higher than marine organic C (Supplementary Information section 1.12).
The δ 13 C values of marine carbonate rocks are widely used to constrain the relative proportions of organic C and CaCO 3 burial 35,49,50 . With the emergence of land plants, our model predicts a ~threefold increase in organic C burial, driven mostly by the onset of terrestrial organic C burial. Over the same period, inorganic C burial is predicted to increase by only a factor of ~2 (Extended Data Fig. 3). Nevertheless, model carbonate δ 13 C values increase only modestly over this interval, consistent with the marine carbonate δ 13 C record (Extended Data Fig. 3). The modest increase is due to the stabilizing effect of oxidative weathering of organic C-bearing sedimentary rocks (for example, ref. 51 ), which delivers isotopically depleted C to the oceans more rapidly upon the land-plant-related p O 2 increase. Our trajectory for deep-ocean O 2 concentrations is consistent with molybdenum isotope and iron-based proxies for deep-ocean O 2 concentrations, as well as with previous model estimates of atmospheric p O 2 , which suggest a late Palaeozoic increase in deep-ocean O 2 concentrations from a few tens µM to near-modern values 52-55 .
We identify the relative contributions of continental and seafloor weathering to the oceanic alkalinity influx as a key factor controlling the long-term P budget. During periods of low continental weatherability, the relative contribution of seafloor weathering increases, resulting in relatively low rates of P delivery to the ocean from continental weathering. We suggest that this mechanism has driven changes in P delivery rates to the ocean and seawater P concentrations over Phanerozoic time, with a notable role for two evolutionary and tectonic events. The first was the Devonian evolution of land plants, which increased the extraction efficiency of nutrients from rocks, thereby enhancing P delivery to the ocean and increasing seawater P concentrations. The second was the late Palaeozoic assembly and subsequent Mesozoic break-up of the supercontinent Pangaea, which resulted in a decrease and then an increase in rainfall over the continents. We suggest that the associated low continental weatherability during the tenure of Pangaea (~306-237 Ma) led to relatively low seawater P concentrations, which then increased to near-modern levels upon the supercontinent's break-up. Our biogeochemical model suggests that the general increase in the supply of P to the ocean over Phanerozoic time led to an increase in marine primary productivity, with implications for the C and O 2 cycles and for evolutionary trajectories. Overall, our study highlights the crucial role played by the rock-CO 2 cycle in nourishing life on Earth.

Online content
Any methods, additional references, Nature Research reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/s41561-022-01075-1.   Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

References
Springer Nature or its licensor (e.g. a society or other partner) holds exclusive rights to this article under a publishing agreement with the author(s) or other rightsholder(s); author self-archiving of the accepted manuscript version of this article is solely governed by the terms of such publishing agreement and applicable law.

Model development
To constrain the Phanerozoic evolution of P fluxes and concentrations in seawater in response to evolutionary, tectonic and other processes, we developed a model for the coupled, long-term biogeochemical cycles of P, C, O 2 and Ca, as well as total pools of sedimentary organic matter and carbonate rocks (Fig. 1). The model simulates the fluxes of these elements between the ocean-atmosphere and rock reservoirs and computes the speciation of dissolved inorganic carbon in seawater ( Supplementary Information section 1). Our model includes 35 input parameters (Extended Data Table 2 and Extended Data Fig. 1) and 11 time-dependent forcings (Extended Data Table 3 and Extended Data Fig. 2). Using the model, we simulated the evolution of the marine and sedimentary pools of these elements through the Phanerozoic by adopting some established forcings and parameterizations and introducing new parameterizations ( Supplementary Information  section 1). The model outputs are the estimates of the long-term sources and sinks, as well as the concentrations and reservoir sizes of the state variables (for example, seawater P concentration, p CO 2 , p O 2 , carbonate rocks and organic C in sedimentary rocks; Extended Data Table 1) through time. We present these outputs as colour contours in Figs. 2, 3 and 5 and Extended Data Fig. 3. To account for parameter and time-dependent forcing uncertainties, we randomly selected inputs across realistic parameter space and validated the results using geologic proxies for p CO 2 , the C content of oceanic crust and the δ 13 C of carbonate minerals (see Numerical solution for more information about the statistical approach). The mass balances of P and Ca in the ocean, C and O 2 in the oceanatmosphere and the sedimentary pools of organic matter (C org ) and carbonate rocks are: Reactive P enters the ocean from continental weathering (F wP,ocean ) and is removed via organic P burial (F borg,P ), apatite burial (F bap ), adsorption onto iron oxides in sediments (F bFeP ) and in hydrothermal plumes (F plume,P ) and removal during off-axis seafloor weathering reactions (F off,P ). Carbon enters the ocean-atmosphere by volcanic outgassing, which includes mid-ocean-ridge fluxes (F vmor ), mantle plumes (F vplume ) and metamorphic flux of organic-matter-bearing and CaCO 3 -bearing rocks (F varc,org and F varc,carb , respectively), and weathering of sedimentary rocks containing CaCO 3 and organic matter (F wcarb and F worg , respectively). The major removal processes of C are inorganic C burial (F bcarb ), marine organic C burial (F borg ) and land-plant-derived organic C burial (F borg,land ). Oxygen is produced by the burial of organic C and consumed by oxidation of organic C during weathering of organic-mattercontaining sedimentary rocks (F worg ) and by oxidation of reduced volcanic gases such as H 2 , CO and CH 4 (F ored ). The main source of Ca is continental weathering of carbonate and silicate rocks (F wcarb and F wsil , respectively). Furthermore, a substantial amount of Ca enters the ocean from off-axis seafloor weathering reactions (F wsf ). The main Ca removal mechanism is carbonate mineral burial. Organic-matter-containing sedimentary rocks are formed by the burial of marine and terrestrial organic C and are lost by weathering and metamorphism (F morg ). Calcium-carbonate-containing sedimentary rocks are formed by carbonate deposition (chemically or biologically driven) and lost by weathering and metamorphism (F mcarb ). To simulate the evolution of the model variables through time, we adopted major Phanerozoic evolutionary, tectonic, palaeogeographic and other forcings. The full set of forcings is shown in Extended Data Table 3 and Extended Data Fig. 2 and described in Supplementary Information section 2. The full parameter ranges and initial conditions are presented in Extended Data Table 2,  Supplementary Tables 6 and 7 and Extended Data Fig. 1. Off-axial seafloor weathering is a critical process that highly affects our results. In the following, therefore, we describe the related model development.
The rest of the model parameterizations and scaling relationships are fully described and discussed in Supplementary Information section 1.

Off-axial seafloor weathering
Seafloor weathering occurs when seawater circulates through the upper ~600 m of the highly permeable oceanic basalts in off-axial (low-temperature) systems. Like continental weathering, seafloor weathering is a net source of alkalinity, which releases Ca 2+ ions, a fraction of which precipitate as CaCO 3 in veins and void fillings 22,56 and a fraction of which are injected into the ocean. Global circulation models relate the temperature of deep-ocean water to the surface temperature at high latitudes, where deep water forms 25 . Given several indications of a linear relationship between deep-water temperature and pore-water (aquifer) temperature 23,24 , an approximately linear relationship between aquifer temperature and surface climate is expected 25 .
Adopting aquifer temperatures that are linearly proportional to the deep-ocean temperature, which is itself linearly proportional to surface temperature, we scale seafloor weathering rates with temperature in an Arrhenius-type dependence. The total rate of Ca 2+ leaching due to seafloor weathering (F wsf ) is further scaled with seafloor-spreading rate (f SF ; Supplementary Information section 2.3), assuming that faster seafloor spreading enhances the seafloor weathering rate by exposing fresh crust to weathering: Here F wsf is a unitless parameter representing the enhancement of seafloor weathering rates due to the temperature relative to the present day, R is the gas constant, T 0 pore is the present-day, average pore-water temperature (K), aT is a parameter that describes the linear relationship between the surface and pore-water temperatures extracted from proxy data and an ensemble of coupled global climate models 26 and ΔE bas is the effective activation energy of basalt dissolution (kJ mol -1 ). F 0 wsf is present-day flux of calcium into the ocean from seafloor weathering ( Supplementary Information section 3.1) and f SF is the seafloor-spreading rate ( Supplementary Information section 2.3).
For parameter values, see Extended Data Table 2.
The activation energy of basalt dissolution controls the strength of the seafloor-climate feedback. A weaker feedback could be a result of a greater role of geothermal heating and/or sediment cover in governing pore-water temperature 57 . To account for the uncertainties in the strength of this feedback, we draw the activation energy for basalt dissolution from a uniform distribution between 40 and 120 kJ mol -1 https://doi.org/10.1038/s41561-022-01075-1 (Extended Data Table 2 and Supplementary Table 2). Sensitivity analysis indicates that when this feedback is stronger (at higher activation energy for basalt dissolution), P concentrations in the early Palaeozoic are lower than the case of a weaker feedback ( Supplementary Fig. 1, panel 'ΔE bas ').

Numerical solution
Our model contains many parameters such as present-day fluxes and rate constants (Extended Data Tables 1-3). Some of these parameters represent a global average and have been estimated using a range of methods in past studies. Therefore, there are great uncertainties associated with some of the parameter values. To account for the uncertainties associated with the parameters and time-dependent forcings, we performed ~10 6 simulations. In these simulations, we randomly drew parameters and time-dependent forcings from distributions that represent the uncertainty in their values (Extended Data Tables 2 and  3 and Extended Data Figs. 1 and 2). To construct these distributions, we conducted an extensive literature review of existing estimates of the model parameters 7,[9][10][11][12][14][15][16][17][18][19][20][21][23][24][25]28,31,36,50,56, . We explain the procedure for estimation of these parameters and their values in Supplementary Information sections 2 and 3 and in Extended Data Tables 2 and 3. For most of the parameters, there exist only a limited number of estimations (for example, Supplementary Tables 1 and 2). We therefore draw these parameters from uniform distributions with lower and upper bounds constrained by the minimum and the maximum values in the literature. For example, we found four different estimates of the present-day burial of organic C on land (F borg,land = 4.3, 3.6, 7.9, 6.6 TmolC yr -1 ; Supplementary Table 1). Therefore, in each simulation, we randomly draw this parameter from a uniform distribution between 3.6 and 7.9 TmolC yr -1 . Estimates of other parameters are more numerous, and we fit normal distributions to these parameters. For example, in each simulation, we draw the present-day mid-ocean-ridge volcanic CO 2 source (F 0 vmor ) from a normal distribution fitted to 21 available estimates (Supplementary Table 3). Our model also includes time-dependent geologic forcings, such as the seafloor-spreading rate and weathering enhancement due to land-plant evolution 28,29,34,55,85,131-143 (Extended Data Table 3 and Extended Data Fig. 2). To account for uncertainties in these time-dependent forcings, in each simulation, we draw them also from distributions. We present these distributions in Extended Data Fig. 2, and we explain how we chose these in Supplementary Information section 2.
We exclude simulations that yielded results out of allowable ranges of several present-day observables. The allowable ranges include ±50% pre-industrial p CO 2 (µatm), ±15% of modern p O 2 (atm), ±15% of modern deep-ocean P concentrations (µM), ±15% of modern oceanic Ca concentration (mM), ±4 K of modern global surface temperature and 2σ around the mean of the C isotopic composition of carbonate minerals in modern marine sediments (‰) (Supplementary Table 8). We note that the exclusion method did not change the randomly sampled parameter and forcing space, only excluded some combinations of parameters (Extended Data Fig. 1). Since our results do not narrow the uncertainty on most model parameters and forcings, they faithfully represent the explored parameter space.
To test the sensitivity of the model prognostics to the different parameters, we varied each parameter (total 35) within its range and recorded the deep-ocean P concentration at two time points, 0 and 500 Ma. The results are presented in Supplementary Fig. 1. We further tested the model prognostics' sensitivity to the different geologic forcings by applying each of the forcings separately, together with the Phanerozoic evolution of solar luminosity, which was included in all model sensitivity simulations. The results are presented in Supplementary Figs. 2-4. The marine P concentration, which is the main prognostic of interest in this study, is sensitive to the parameters and forcings that affect the delivery of alkalinity and CO 2 to the oceanatmosphere. The reason for this sensitivity is that the influxes of P to the ocean depend on the absolute rates and relative contributions of continental and seafloor weathering (sources of alkalinity), which are required to balance the outgassing source of CO 2 by a combination of carbonate and organic C burial. Thus, the parameters and forcings that affect the delivery of alkalinity and CO 2 to the ocean-atmosphere ultimately control the sources and major sinks of P.

Model validation
In Extended Data Fig. 3, we compared our model output with a compilation of δ 13 C values in marine carbonate rocks 50 , which contains 43,972 samples. Data before the mid-Jurassic are derived mostly from shallow-water (platform), near-shore carbonates. Data from the Jurassic through the modern are derived mostly from CaCO 3 secreted by planktonic foraminifera and calcareous nanoplankton, which was extracted mainly from oceanic drill cores. Overall, the envelope covering 90% of the model δ 13 C results is mostly within one standard deviation of a 10-million-year moving average of δ 13 C values measured in marine carbonate rocks.
In Extended Data Fig. 3, we compared our model p CO 2 predictions with ∼1,240 newly calibrated proxy p CO 2 estimates derived from 112 studies of stomatal density in fossilized leaves, the δ 13 C values in alkenones, liverworts and carbonate minerals in palaeosols, and marine boron isotopes 144 . The uncertainty envelope on the proxy p CO 2 estimates (the 5th to 95th percentiles) includes the error associated with each proxy and age determination 144 . The overall pattern in the proxy records is of high p CO 2 (maximum probability between 1,000 and 2,000 ppm) during the Early Devonian, a decline to present-day levels by ~320 Ma, a rise to high p CO 2 (maximum probability between 1,000 and 2,000 ppm) during the early Mesozoic and then a decline to present-day levels 144 . This pattern is roughly captured by our model envelope over the period covered by the proxy records. The model produces elevated p CO 2 (~1000-4,000 ppm) during the Early Devonian, low p CO 2 (~300-400 ppm) between ~350 and 252 Ma, higher early Mesozoic p CO 2 (~1,000-3,000 ppm) and a Cenozoic decline to present-day p CO 2 (Extended Data Fig. 3).
We note that the evolution of the model prognostics is a result of tectonic processes such as continental drift, subduction and metamorphism and of major evolutionary events, such as the expansion of land plants. These processes and events operate on 100-million-year timescales, and our model is, therefore, capable of reproducing the long-term evolution of atmospheric CO 2 and the δ 13 C of carbonate minerals. Shorter-timescale forcings and events, such as orbitally driven changes in sea level or ~million-year-long ocean anoxia events, are not represented. Thus, even though such events probably affected atmospheric CO 2 and carbonate δ 13 C values over Earth history 35 , they are not accounted for in the model and cannot produce short-timescale variation in model CO 2 and δ 13 C values.
In our model, variations in surface temperature are derived from solar luminosity and p CO 2 ( Supplementary Information section 1.1). Our model predicts an overall decline in surface temperature of ~15° C over the Phanerozoic (Extended Data Fig. 3). Furthermore, our model suggests a difference of ~8° between the Mesozoic and Cenozoic surface temperatures, which agrees with variations in deep-water and pore-water temperatures based on a compilation of the oxygen isotope composition of carbonate minerals deposited in the upper oceanic crust 23 . We adopt observation-supported suggestions that seafloor temperatures are offset from surface temperatures by an approximately constant amount 23,25,26 . Therefore, variations in surface temperatures are expected to manifest as similar variations in seafloor temperatures.
In Fig. 2b we compared the Ca flux from the altered oceanic crust, estimated from the bulk CO 2 content of oceanic crust at different ages and an assumption of a 1/1 molar ratio of Ca/C in the crustal carbonates that yielded the CO 2 , with the fluxes of Ca leaching from the oceanic crust obtained in our model. The bulk CO 2 content of the altered crust was obtained from table 1 in ref. 22 . We note that the parameterization Article https://doi.org/10.1038/s41561-022-01075-1 of the Ca leaching fluxes in our model (Off-axial seafloor weathering) is unrelated to the observations of ref. 22 , making this comparison an independent test of the predicted Ca fluxes. To convert from CO 2 (wt%) to C flux from the ocean into the crust (TmolC yr -1 ), which is approximately equivalent to the Ca flux from alteration of the oceanic crust, we used the following equation: F CO 2 = ΔCO 2 MW CO 2 × d × ρ × R 0 SFS × f SF × 10 × 10 −12 , where F CO 2 is the flux of CO 2 from the ocean into the crust (TmolC yr -1 ), ΔCO 2 is the inorganic C uptake by the oceanic crust (as measured, not including sediment) (wt% = g 100 g -1 ), d is the crust thickness (600 m), ρ is the density of the crust (2,900 kg m −3 ), R 0 SFS is the present-day rate of crust production (3.45 × 10 6 m 2 yr -1 ), MW CO 2 is the molecular weight of CO 2 (44 g mol -1 ), f SF is a unitless, time-dependent parameter that describes the enhancement in the rate of seafloor production relative to the present day (Supplementary Information section 2.3 and Extended Data Fig. 2), 10 is a multiplier to convert wt% units to 100 g kg -1 and 10 -12 is to convert from mol to Tmol.

Data availability
All study data are included in this article and/or in the Supplementary Information. Source data are provided with this paper.

Code availability
The code used in this study is available via Zenodo (https://doi. org/10.5281/zenodo.6874786).    Table 3 for description of these arrays. (a) Weathering enhancement due to landplant evolution (f E ). f E was drawn from uniform distributions at any given time, where the boundaries of the uniform distribution change over time to account for the evolution of land plants described in Section 2.1, SI. (b) The effect of paleogeography on runoff throughout the Phanerozoic (f PG , Section 2.2, SI). (c) The seafloor spreading rate (f SF , Section 2.3, SI). (d) A time-dependent forcing that accounts for the enhancement of the climate sensitivity during cold periods (f glac , Section 2.4, SI). f glac is drawn from a uniform distribution between unity and two when there is evidence for long-lived glaciations, and is set to unity over the rest of the Phanerozoic. (e) The effect of continental configuration, latitude and vegetation cover on surface albedo, and consequently, on average continental temperature (T geog , Section 2.5, SI). T geog is drawn from a normal distribution, where the mean was adopted from Goddéris et al. (2012) 85 , and the standard deviation was adopted from the range of climate predictions associated with the CMIP5 models 145 . (f) Uplift (f U , Section 2.6, SI). (g) Land area (f A , Section 2.7, SI). (h) The fraction of land area covered by carbonates (f L , Section 2.7, SI). (i) The effect of volcanic land cover on continental weatherability (f volc , Section 2.8, SI). ( j) A time-dependent forcing that represents the colonization and expansion of terrestrial biomass (f cland , Section 2.9, SI). f cland was drawn from uniform distributions at any given time, where the boundaries of the uniform distribution change over time to account for the evolution of land plants described in the SI. (k) The effect of the evolution of pelagic calcifiers on CaCO 3 subduction (f C , Section 2.10, SI). All the time-dependent variables are unitless, and normalized to the present day. Unless mentioned otherwise, all the timedependent variables were drawn from normal distributions, where the mean was adopted from different sources (black lines, SI), and 2σ is taken to be 25% of the mean (the upper and lower boundaries of the gray shaded area).