Phanerozoic seawater Mg/Ca variations driven by supercontinent cycles

Domain: biology
Published: 2026-03-20
Editorial outcome: major_revision

Abstract
Seawater Mg/Ca has oscillated markedly over the Phanerozoic, closely associated with secular variations in sedimentation, icehouse-greenhouse fluctuations, and the evolution of calcifying invertebrates. However, the drivers of these long-term oscillations remain unresolved. Here, we develop an inverse model that couples seawater Mg concentration and Mg isotope composition variations. By exploiting the contrasting Mg isotopic fractionation associated with Mg-silicate formation and dolomitization, this model quantifies the major marine Mg sinks. Model results show that increases in seawater Mg/Ca during supercontinent assembly were explained by reduced rates of silicate alteration and dolomitization. In contrast, decreases in seawater Mg/Ca during periods of supercontinent stasis and early breakup were driven by intensified Mg-silicate formation and dolomitization. Variations in marine Mg-silicate and dolomite formation were further regulated by seafloor spreading, continental configuration, and climate throughout the supercontinent cycle. Together, our results identify the supercontinent cycle as a first-order regulator of Phanerozoic seawater Mg/Ca through its control on marine Mg-silicate and dolomite fluxes.

Introduction
Seawater chemistry is regulated by a dynamic balance between elemental inputs and outputs via river discharge, hydrothermal fluids, seawater-rock reactions, and the formation of authigenic and biogenic minerals in marine sediments1,2. Magnesium and calcium are the second- and third-most abundant cations in seawater, and proxy data indicate their concentrations have varied considerably over the Phanerozoic, resulting in a maximum seawater Mg/Ca ratio of approximately 5 at present and minima of <1 during the Jurassic–Cretaceous and Cambrian–Carboniferous intervals (Fig. 1a). These periodic changes in seawater Mg/Ca have driven oscillations in the dominant mineralogy of marine carbonates, late-stage evaporites, and invertebrate skeletons throughout the Phanerozoic, and closely coincide with transitions between icehouse and greenhouse climate states (Fig. 1)1. It is well established that riverine input is a major process introducing Mg and Ca into the oceans, while both dolomitization and seafloor alteration remove Mg and add Ca, and the marine Ca budget is additionally regulated by calcium carbonate burial3–7, yet the relative importance of these processes in driving Phanerozoic seawater Mg/Ca variations remains difficult to quantify1,8. The prevailing view attributes seawater Mg/Ca variations to fluctuations in mid-ocean ridge hydrothermal fluxes9 or the rate of dolomitization3,4, since both processes remove Mg and release Ca to seawater, whereas some studies emphasize the role of authigenic clay formation in removing Mg from seawater during the Cenozoic10,11. Uncertainty in the relative contributions of these processes to seawater Mg/Ca arises from the incomplete nature of the geological record, which limits accurate reconstructions of past elemental fluxes1,4,6. Additionally, attributing changes in seawater Mg/Ca to specific processes is challenged by the difficulty in resolving Mg and Ca flux variations associated with individual sources and sinks1.

Integrating the seawater Mg/Ca record with Mg isotopic compositions (δ²⁶Mg, defined as [(²⁶Mg/²⁴Mg)sample / (²⁶Mg/²⁴Mg)standard – 1]×1000, where the standard is DSM3) enables identification of the processes controlling the marine Mg budget, because the formation of Mg-silicates and dolomite in the ocean overall fractionates Mg isotopes in opposite directions12. Consequently, variations in seawater δ²⁶Mg can distinguish between changes in seawater [Mg] driven primarily by Mg-silicate versus dolomite formation and associated geological processes10. This approach has been successfully applied to elucidate mechanisms underpinning seawater Mg/Ca variations during several key geological periods10,11,13,14.

Here, we develop an inverse model of the marine Mg cycle to quantitatively constrain the major processes that have driven seawater Mg/Ca variations throughout the Phanerozoic. By integrating the latest seawater [Mg] and δ²⁶Mg records into the model and comparing the results with Phanerozoic tectonic and climate reconstructions, we find that rising Mg/Ca ratios during supercontinent assembly result from reduced rates of silicate alteration and dolomitization, whereas declining ratios during breakup stem from intensified removal into these sinks. These results demonstrate that the long-term evolution of seawater Mg/Ca is fundamentally regulated by the supercontinent cycle and the associated changes in tectonics and climate.

Results
Inverse model development
Magnesium is a conservative element in seawater with a modern residence time of ca. 12 Myr15. Variations in seawater Mg/Ca are driven by long-term changes in continental weathering-derived riverine input, and by changes in the rates of Mg removal through i) silicate formation via high-temperature (>90 °C) hydrothermal reactions at the mid-ocean ridges (HT), ii) silicate formation via low-temperature (<90 °C) processes including seafloor alteration on ridge flanks and clay authigenesis (LT), and iii) dolomitization (dol)1,8. According to mass balance, the oceanic Mg inventory (N) and Mg isotopic composition (δSW) at a given time (t+Δt) are determined by the seawater Mg budget and δ²⁶Mg at time t, and the cumulative effects of sources and sinks over this period:

N(t+Δt) = N(t) + Δt [Friv(t) - FHT(t) - FLT(t) - Fdol(t)] (1)

N(t+Δt) δSW(t+Δt) = N(t) δSW(t) + Δt [Friv(t) δriv(t) - FHT(t) (δSW(t) + ΔHT-SW) - FLT(t) (δSW(t) + ΔLT-SW) - Fdol(t) (δSW(t) + Δdol-SW)] (2)

By applying these mass-balance equations, we can quantify the fluxes of Mg sources and sinks over the Phanerozoic, since key parameters are well constrained and their uncertainties can be robustly propagated. Seawater [Mg] and δ²⁶Mg are known throughout the Phanerozoic (Fig. 1), and N can be calculated using seawater [Mg] and assuming a constant Phanerozoic seawater mass of 1.4 × 10²¹ kg15. The Mg sinks associated with the formation of Mg-silicate and dolomite have known fractionation factors (Supplementary Table 1). For Mg-silicate formation at high temperature, ΔHT-SW is set as 0‰ due to the quantitative removal of Mg from fluid to crust6,16. For low-temperature Mg-silicate formation, including both seafloor alteration and clay authigenesis, we assign a ΔLT-SW of 0.8 ± 0.2‰ to cover a wider range of previous constraints based on marine clays, altered crust, and hydrothermal fluids5,17–19. For dolomitization, Δdol-SW has been experimentally constrained20: Δdol-SW = -0.1554(±0.0096) × 10⁶/T², with a calculated Δdol-SW of -1.8 ± 0.1‰ at 20 °C. We allow Δdol-SW to vary with Phanerozoic temperature21 and propagate the associated uncertainty.

The modern weighted average global riverine δ²⁶Mg (δriv) is -1.09 ± 0.05‰, primarily controlled by the mixture between carbonate and silicate weathering-derived Mg22. In our model, δriv is taken as increasing from -0.9 ± 0.1‰ in the Cambrian to the modern value of -1.09 ± 0.1‰ (Supplementary Fig. 1), following a modest Phanerozoic increase in carbonate weathering13.

To determine the unknown fluxes, riverine Mg input (Friv) and Mg removal via Mg-silicate formation (FHT and FLT) and dolomitization (Fdol), we performed Monte Carlo simulations, sampling across their plausible ranges and retaining combinations that reproduce the observed seawater [Mg] and δ²⁶Mg records. For each 1-Ma time step backwards from 0 to 550 Ma, 1,000,000 random parameter sets were generated within ranges of 4.5–6 Tmol/yr for Friv, 1–3 Tmol/yr for FHT, 0.6–2.5 Tmol/yr for FLT, and 1–4 Tmol/yr for Fdol that cover previous estimates (see Methods and Supplementary Table 2 for reference). At each time step, we computed the predicted N and δSW using Eqs. (1–2), and retained only those simulated parameter sets whose outputs matched the observed records within a relative 1% tolerance (accepted number > 20,000 per step). All Monte Carlo simulations and data analyses were performed in R (v4.3.1), and the code is provided in the Code Availability section. Uncertainties in accepted solutions for each flux are shown as 68% and 95% confidence intervals in Fig. 2. Since modeled FHT and FLT are nearly identical at each step and highly correlated over the Phanerozoic (R² = 0.84; Supplementary Figs. 2−3), they are combined as a total Mg-silicate sink (Fsil; Fig. 2c) in subsequent discussions.

Our inverse model avoids reliance on fragmentary sedimentary records or uncertain tectonic reconstructions and enables the simultaneous quantification of the major Mg fluxes into and out of ocean water over the Phanerozoic. Moreover, by including the variation of seawater [Mg] in Eq. (1), our model explicitly incorporates temporal variations in Mg residence time. To test the robustness of our results, we performed local sensitivity tests in which each input variable was varied within its data-driven uncertainty and the resulting ensembles were evaluated using Partial Rank Correlation Coefficients (PRCC) and local Δrange effect sizes to quantify short-term, monotonic controls on each flux. These local tests capture how small, realistic deviations in seawater [Mg], δ²⁶MgSW, Δdol-SW, ΔLT-SW, δ²⁶Mgriv, and prior input ranges propagate into Friv, Fdol, and Fsil at each time step (see Methods for details). For the secular analysis, we first calculated time-step mean fluxes and then aggregated them into 20-Myr bins, relating the binned mean fluxes to the slowly varying boundary conditions, including seawater [Mg], δ²⁶MgSW, Δdol-SW, and δ²⁶Mgriv, using a Deming regression that accounts for uncertainties in both variables.

Our sensitivity tests show that, for the integrated time step flux variations, the directions of flux changes are primarily governed by δ²⁶MgSW whereas their magnitudes are mainly controlled by seawater [Mg], with the prescribed Friv range additionally modulating its own variability (Fig. 3; Supplementary Figs. 4–9). Over the whole Phanerozoic, δ²⁶MgSW continues to dictate the direction of dolomite and riverine flux, while seawater [Mg] mainly determines the direction of Fsil. Variations in dolomite and riverine fluxes are mainly controlled by δ²⁶MgSW, while Fsil variation is dominated by seawater [Mg] (Fig. 4). We additionally explored scenarios in which the FLT is higher than the FHT (FLT = 1–3 Tmol/yr and FHT = 0.6–2.5 Tmol/yr). The resulting solutions are indistinguishable from our primary outputs (Supplementary Fig. 10), demonstrating that the assumption of a higher FHT relative to FLT does not influence our results. Together, these analyses demonstrate that the inferred controls on all Mg fluxes arise predominantly from the intrinsic mass-balance framework and from seawater [Mg] and δ²⁶MgSW variations, rather than from prior ranges of input-flux or isotope-fractionation factors.

Major drivers of phanerozoic seawater Mg/Ca variation
The Phanerozoic evolution of seawater [Mg] can be divided into six stages according to our fitted curve, comprising two full cycles of decreasing, fluctuating (minor variation), and increasing seawater [Mg] (Fig. 2a). Seawater [Mg] decreased from 550 to 500 Ma, remained relatively stable with minor fluctuations between 500 and 390 Ma, and increased from 390 to 310 Ma. The second cycle began with a decline in seawater [Mg] from 310 to 200 Ma, followed by a fluctuating interval between 200 and 90 Ma, and a final increase from 90 to 0 Ma. These intervals closely parallel secular variations in seawater Mg/Ca and the timing of aragonite–calcite transitions (Fig. 2a). In our model, we force source-sink net Mg imbalance (Fluxsource-sink) to co-vary with seawater [Mg] to maintain mass balance, ensuring positive Fluxsource-sink values during seawater [Mg] increase periods (Fig. 2e). Over the Phanerozoic, average Mg flux ranges are 5.2–5.4 Tmol/yr for Friv, 1.4–2.8 Tmol/yr for Fdol, and 2.9–3.9 Tmol/yr for Fsil (Fig. 2b–e).

Seawater [Mg] varies in concert with seawater Mg/Ca throughout the Phanerozoic, and the same geological drivers that regulate seawater [Mg] also control seawater Mg/Ca. Seawater [Mg] and [Ca] are independently reconstructed by measuring Mg and Ca concentrations of fluid inclusions in early-stage primary marine halite23, and their inverse evolution, together with the co-variation of seawater [Mg] and Mg/Ca (Fig. 1a–b), suggests that these secular variations are driven by the same geological processes. During intervals of decreasing and increasing seawater [Mg], seawater [Mg] and Mg/Ca show strong positive correlations (R² = 0.93–0.99; Supplementary Fig. 11), indicating that processes responsible for [Mg] changes also regulate Ca and the Mg/Ca ratio. This is reasonable since two major geological processes, seafloor alteration and dolomitization, involve 1:1 Mg-Ca exchange, such that seawater [Mg] increases correspond to [Ca] decreases and vice versa. In contrast, during intervals of minor seawater [Mg] fluctuations, correlations between [Mg] and Mg/Ca weaken (R² = 0.01–0.53; Supplementary Fig. 11) while seawater [Ca] varied significantly (Fig. 1b). This indicates that independent Ca-related processes, such as calcium carbonate burial, also considerably regulate seawater Mg/Ca during these periods. For instance, the development of pelagic calcification since the Mesozoic may have significantly contributed to declining seawater [Ca]24. In the following sections, we focus on the primary controls on major seawater Mg/Ca transitions between aragonite and calcite seas, while secondary fluctuations within prolonged calcite-sea intervals (ca. 500–390 and 200–90 Ma) may be better constrained through future studies of the marine Ca cycle.

The limited Friv variation in our results does not suggest invariant Friv in the geological record; instead, it reflects the relative insensitivity of Mg mass balance to Friv within plausible ranges. The modest Friv increases during the 420–370 Ma and 50–0 Ma periods coincide with the assembly of the supercontinent Pangea and Alpine-Himalayan collision (Fig. 2b), suggesting enhanced weathering flux during supercontinent assembly. The increase in Friv during 550–500 Ma may reflect intensified weathering and a reorganized clay factory associated with widespread felsic terrains, high continental topography, and elevated CO₂ outgassing14,25–27. However, our model cannot resolve the precise magnitude of these variations. This insensitivity does not substantially affect our interpretation because weathering generally alters Mg and Ca inputs in the same direction unless global lithology changes substantially, and the Phanerozoic continental lithologic composition appears to have remained broadly stable28. Nevertheless, short-lived shifts in lithology can influence riverine Mg/Ca, because calcium carbonate (Mg/Ca ≈ 0.1)4, dolomite (Mg/Ca = 1.0)29, and silicate rocks (Mg/Ca ≈ 0.4 to 0.7)30 contribute distinct Mg/Ca signatures. For instance, the extensive exposure of basalt-mafic rocks during the formation of Pangea at the late Permian31, the enhanced development of carbonate platforms during the breakup of Pangea at the Jurassic-Cretaceous32, and the Alpine-Himalayan Collision during the Cenozoic mountain uplift33 may involve significant changes in global continental lithology. Improved constraints on continental weathering fluxes and lithologic evolution across these intervals would further refine estimates of riverine Mg delivery and Mg/Ca modulation.

Finally, Phanerozoic seawater Mg/Ca variations can be primarily attributed to substantial changes in dolomite and Mg-silicate fluxes. Between 550 and 500 Ma, intense Mg-silicate removal coincided with net Mg loss (e.g., negative Fluxsource-sink), leading to a decline in seawater Mg/Ca, while subsequent weakening of this sink mitigated this decline; the dolomite sink remained nearly constant (Fig. 2). From 500 to 390 Ma, a slight decrease in Fdol combined with an initial rise and subsequent fall in Fsil yielded slightly positive and later negative Fluxsource-sink, accounting for fluctuating seawater Mg/Ca within this prolonged calcite sea period (Fig. 2a–d). Between 390 and 310 Ma, rising seawater Mg/Ca was first driven by a sharp reduction in the dolomite sink (390–340 Ma), followed by weakening of the Mg-silicate sink at 340–310 Ma (Fig. 2a–d). From 310 to 200 Ma, concurrent peaks in both dolomite and Mg-silicate fluxes caused renewed Mg removal and a further decline in Mg/Ca. Despite decreases in both sinks after 240 Ma, Fluxsource-sink remained negative, sustaining the downward trend (Fig. 2a–d). During 200–90 Ma, reduced dolomitization was offset by enhanced Mg-silicate removal, maintaining fluctuating seawater Mg/Ca (Fig. 2a–d). Since 90 Ma, rising seawater [Mg] was driven by a persistent decline in the dolomite sink, amplified after 40 Ma by weakening Mg-silicate removal (Fig. 2a–d).

Controls on Mg-silicate and dolomite fluxes
Marine Mg-silicate formation is considered to be dominated by seawater–seafloor reactions, with clay authigenesis contributing negligibly under Phanerozoic conditions of low dissolved silica and cooler bottom-water temperatures34. For instance, the modern ocean hosts an authigenic clay sink of only 0.02 Tmol/yr5, representing ~0.6% of the total Mg-silicate sink. Mg-silicate flux during seafloor alteration is mainly controlled by seafloor spreading, with faster spreading rates producing more fresh ocean crust available for hydrothermal alteration and Mg-silicate formation7. The tight correlation between FHT and FLT further supports a dominant driver of seafloor spreading on both fluxes (Supplementary Figs. 2–3), with the authigenic clay sink playing only a minor role.

Mg-silicate flux exhibits five prominent peaks at ca. 550, 420, 335, 240, and 110–40 Ma (Fig. 2c). Four of five peaks lag the corresponding maxima in reconstructed global seafloor production rate by 10–30 Myr, consistent with the observation that most hydrothermal alteration of oceanic crust occurs within 30 Myr of crustal accretion6,35. The exception at ca. 420 Ma, which lags the production maximum by ~50 Myr (Fig. 2c), may reflect uncertainties in fluid-inclusion [Mg] data and/or prolonged seafloor–seawater reactions, which can persist for up to 80 Myr35. Notably, the protracted Mg-silicate flux between ~140 and 40 Ma closely matches previously reconstructed hydrothermal fluxes during Pangea dispersal (140–35 Ma; Fig. 2f–g)36. Moreover, Fsil is negatively correlated with the continental fragmentation index during Pangea early assembly (500–320 Ma) and positively correlated during dispersal (200–0 Ma; Fig. 2c, g; Supplementary Fig. 12). These relationships indicate that Mg-silicate formation was significantly enhanced during active tectonic phases (assembly and breakup), whereas it was attenuated during periods of supercontinent stasis. These relationships highlight a strong coupling between the supercontinent cycle, seafloor production, and marine Mg-silicate formation.

Controls on dolomitization rate remain unclear. High sea level29, warm climates13,37, large-scale restriction38, the dissolution-reprecipitation process39, oceanic anoxia40, and microbial activity41 have all been proposed to promote dolomite formation. Microbial activity, elevated temperatures, and restricted hydrological conditions are suggested to accelerate dolomite nucleation and crystal growth37,39,41, while anoxic conditions can further stimulate microbial metabolisms that catalyze dolomitization. High sea level promotes expansive shallow platforms and accelerates reaction kinetics through enhanced evaporation and elevated temperatures, particularly in warm, semi-restricted settings29. Sea level and climate are fundamentally modulated by tectonics: supercontinent assembly leads to a lower sea level due to thermal contraction of aging seafloor and crustal thickening during continental collision, whereas breakup creates younger seafloor and thinner crust, raising sea level42,43. Warm climates typically accompany early supercontinent breakup, when diminished continental collision reduces silicate weathering, and increased arc volcanism, mantle plumes, and large igneous provinces drive rapid pCO₂ increases26,44,45. Large-scale restriction is controlled by global tectonic configuration46, while oceanic anoxia and microbial activity are influenced by atmospheric O₂ and broader environmental conditions, which have also been linked to the supercontinent cycle47,48.

Dolomite abundance has been reconstructed using various approaches, providing an independent benchmark for comparison with our dolomite flux evolution. Our modeled dolomite flux reproduces the major features of sedimentary records, including high dolomite abundance in the early Paleozoic and across the Permian–Triassic boundary, reduced values during the Mississippian, and a long-term decline since the early Mesozoic (Fig. 2d, 5)29,40. We further standardized all dolomite abundance curves through age alignment and long-term trend extraction using Generalized Additive Models and Locally Estimated Scatterplot Smoothing, to directly compare dolomite flux with those dolomite abundance curves. Spearman’s rank correlations (ρ) and bootstrap resampling of the aligned trends show strong agreement between our dolomite flux and three of the five sedimentary records (Fig. 5, Supplementary Fig. 13; ρ > 0.67). The remaining two records show weaker agreement due to intrinsic limitations, one with extremely low temporal resolution, and the other exhibiting a temporal pattern distinct from the broader dataset (Fig. 5).

Correlations between dolomite flux and environmental factors suggest that the dominant controls on dolomitization varied over the Phanerozoic. Dolomite flux is generally high in the early Paleozoic and declines through the Mesozoic, mirroring the secular evolution of atmospheric pO₂ (Fig. 5), suggesting a potential role of oceanic anoxia and microbial mediation on secular dolomite formation. Moreover, dolomite flux shows strong positive correlations with modeled pCO₂ (550–340 and 250–0 Ma) and reconstructed sea level (460–340 and 100–0 Ma) during most Phanerozoic intervals (Fig. 5, Supplementary Fig. 11), consistent with dominant control of climatic and eustatic variations on dolomite formation. However, the pronounced rise in dolomite flux between 340 and 250 Ma coincides with a low pCO₂ and low sea level interval (Fig. 5), suggesting other factors also influence dolomitization. This elevated dolomite flux may be attributed to the restriction of the Paleotethys Ocean during Pangea assembly38, which is also favored by arid continental climates and widespread evaporite deposition46 during this period.

The role of Phanerozoic supercontinent cycles on seawater Mg/Ca
Previous studies have emphasized the influence of tectonics and climate on Phanerozoic seawater Mg/Ca8,36,49, yet direct constraints on how these processes modulated marine Mg fluxes have been lacking. Our model quantitatively constrains variations in Mg-silicate and dolomite fluxes and reveals their tight links to specific tectonic–climatic processes (Fig. 2). Across the Phanerozoic, high Mg/Ca ratios broadly coincide with low sea level and coolhouse conditions associated with supercontinent assembly and major collisional phases, whereas low Mg/Ca ratios correspond to high sea level and greenhouse climates during supercontinent dispersal (Fig. 2). Our results show that this inverse relationship between seawater Mg/Ca and sea level reflects the combined impact of tectonics and climate on Mg-silicate formation and dolomitization. Collectively, these findings imply that the supercontinent cycle exerted a first-order control on Phanerozoic seawater Mg/Ca.

The supercontinent cycle regulates Phanerozoic seawater Mg/Ca as follows: during supercontinent early assembly (e.g. Pangea) or major collisional events (e.g. Cenozoic Alpine-Himalayan collisions), attenuated seafloor spreading and cooler climates suppress Mg-silicate formation and dolomitization, thereby increasing seawater Mg/Ca and shifting global oceans from calcite to aragonite seas (Fig. 6a). During periods of supercontinent early breakup, enhanced seafloor spreading and warmer climates stimulate both Mg-silicate formation and dolomitization, lowering Mg/Ca and returning the oceans to calcite seas (Fig. 6b). During intervals of stable supercontinent dispersal or sluggish plate reorganization, both tectonic and climatic conditions remain relatively steady, yielding limited variability in Mg fluxes and sustaining prolonged calcite seas (Fig. 6c). The supercontinent cycle also drives long-term variations in calcium carbonate deposition by modulating global shelf area50. Supercontinent breakup expands flooded margins and enhances CaCO₃ burial, whereas assembly contracts shelf space and suppresses CaCO₃ deposition. Because this process consumes more Ca during breakup and dispersal but significantly less during assembly, and contributes only marginally to the marine Mg cycle, its influence on Phanerozoic seawater Mg/Ca is likely important only during seawater [Mg] fluctuation periods (Fig. 2a).

Beyond Mg/Ca, the supercontinent cycle has been invoked to regulate other major and trace elements in seawater. Phanerozoic seawater sulfate shows a broadly similar evolution to Mg, which potentially reflects enhanced removal during high-sea-level intervals when microbial sulfate reduction8 and Ca-driven sulfate precipitation intensified49. Seawater Sr and Li concentrations also appear tied to the supercontinent cycle, decreasing during assembly and rising during breakup and dispersal, likely due to their sensitivity to changes in hydrothermal fluid fluxes51,52. These trace elements, strongly modulated by hydrothermal activity51,52, are therefore particularly responsive to tectonic reorganization associated with supercontinent cycling.

By integrating elemental and isotopic records into a marine Mg-cycle model, our study provides quantitative estimates of long-term fluxes through the major Mg sinks in seawater. This advance allows seawater Mg/Ca variability to be more directly attributed to specific geological processes. Our results offer a framework for interpreting Phanerozoic seawater chemistry variations and demonstrate that these variations reflect coupled interactions among ocean chemistry, sedimentation, climate, and the supercontinent cycle.

Methods
Model inputs
Required parameters in our model include Δdol-SW, ΔHT-SW, ΔLT-SW, and input ranges of Friv, FHT, FLT, and Fdol. Temperature-dependent Δdol-SW has been experimentally constrained20: Δdol-SW = -0.1554(±0.0096) × 10⁶/T². Calculated Δdol-SW is -1.8 ± 0.1‰ at 20 °C. Considering the variation of Earth’s surface temperature from ~9 to ~33 °C in the Phanerozoic21, we allow Δdol-SW to vary with temperature and propagate uncertainty in the equation (from -1.95 ± 0.1‰ to -1.66 ± 0.1‰). During high-temperature seafloor-seawater reactions, Mg is removed quantitatively from fluid to crust with nil isotope fractionation6,16, and we therefore set ΔHT-SW to 0‰. During low-temperature Mg-silicate formation, previous studies consistently suggest a positive ΔLT-SW, varying from 0.1‰ to 1.1‰ in authigenic clays5, with average values of 0.8‰ to 1‰ in altered oceanic crust17,18, and 0.25‰ to 7‰ in hydrothermal fluids19. Most studies adopt 0.6–0.8‰, and we therefore set 0.8 ± 0.2‰ to encompass a relatively large range in previous constraints. Reference ranges of these isotope fractionation in previous studies and this study are summarized in Supplementary Table 1.

Ranges of Friv, FHT, FLT, and Fdol for the implementation of Monte Carlo simulations in our model are grounded in previous estimates on modern and Phanerozoic models. Previous estimates on Friv, FHT, FLT, and Fdol range from 4 to 6.1 Tmol/yr, 0.24 to 2.6 Tmol/yr, 0.6 to 2.5 Tmol/yr, and 0.75 to 3 Tmol/yr, and we set a range of 4.5–6 Tmol/yr, 1–3 Tmol/yr, 0.6–2.5 Tmol/yr, and 1–4 Tmol/yr for input ranges of Friv, FHT, FLT, and Fdol in our model (Supplementary Table 2). Since some studies suggest lower FHT than that of FLT19, we also considered a scenario using ranges of 1–3 Tmol/yr for FLT and 0.6–2.5 Tmol/yr for FHT, which resulted in identical results to our primary outputs (Supplementary Fig. 10). To verify whether these inputs can effectively cover the actual ranges of these fluxes over the whole Phanerozoic and whether these prior input ranges would affect our model outputs, we further included variable upper and lower ranges for these fluxes at each time step to quantify their effect on our results in the following sensitivity tests.

Sensitivity test
We evaluate two types of controls on the modeled riverine, Mg-silicate, and dolomite Mg fluxes: (1) local sensitivity at each individual time step, and (2) secular drivers across the Phanerozoic. The local sensitivity quantifies which parameters constrain flux magnitudes within a given time step, whereas the secular analysis identifies which boundary conditions primarily explain the long-term evolution of these fluxes through Earth’s history. For both analyses, we use Partial Rank Correlation Coefficients (PRCC) to quantify the independent (partial) monotonic control of each variable on each flux, and we use Δrange metrics (ΔY or ΔF) to quantify the corresponding amplitude of that control in physical units (Tmol/yr).

Monte Carlo sampling was used to identify flux combinations satisfying mass–isotope balance, whereas Latin Hypercube Sampling was employed to systematically explore parameter sensitivities around these solutions. For each time step, we generated 3×10⁴ Latin Hypercube samples (LHS) covering the full range of model parameters. LHS provides uniform, stratified coverage of the parameter space, enabling efficient yet comprehensive exploration of model behavior. Local perturbations follow the data-driven uncertainties: seawater [Mg] is drawn from 68% confidence bounds, and δ²⁶MgSW is perturbed by ± 0.2‰, which is derived from current uncertainty in the reconstructed Phanerozoic curve13. Δdol–SW is drawn via the temperature-dependent coefficient (0.1554 ± 0.0096) × 10⁶/T², ΔLT–SW is sampled in 0.6–1.0‰, and δ²⁶Mgriv in ± 0.1‰. Physically plausible Mg fluxes are allowed by first drawing a lower and upper bound for each flux in Tmol/yr, then sampling uniformly within that interval: Friv is allowed to vary between ~4.5 and 6.0 Tmol/yr (Friv_min = 4.5–5.0; Friv_max = 5.5–6.0), Fdol between ~1 and 4 Tmol/yr (Fdol_min = 1.0–1.5; Fdol_max = 3.5–4.0), FHT between 1 and 3 Tmol/yr (FHT_min = 1.0–1.5; FHT_max = 2.5–3.0), and FLT between 0.6 and 2.5 Tmol/yr (FLT_min = 0.6–1.2; FLT_max = 2.0–2.5). These accepted ensembles form the foundation for both the local and secular analyses, while parameters that exhibit genuine secular trends are evaluated in the secular framework. Local sensitivity therefore reflects uncertainty within each geological interval, whereas the secular analysis evaluates the effect of boundary-condition evolution through time.

(a) Partial Rank Correlation Coefficients (PRCC)
PRCC quantifies the monotonic influence of each parameter (Xᵢ) on a modeled flux (Y) while removing covariation with all other parameters. For the local analysis, at each time step, all variables are rank-transformed, Xᵢ and Y are separately regressed against all remaining parameters, and the Spearman correlation (ρ) between their residuals defines the local PRCC:

PRCC(Xᵢ, Y) = ρ(rXi, rY) (3)

The sign indicates whether Xᵢ increases or decreases Y, and |PRCC| measures the strength of that relationship. |PRCC| values are interpreted following common practice, with higher absolute values indicating stronger monotonic influence. Local PRCCs reveal, for each time step, which parameters co-vary with the flux to satisfy Mg mass–isotope balance. For each time step, we additionally apply a nonparametric bootstrap by resampling the accepted Monte Carlo ensembles with replacement (1000 iterations) and recalculating PRCC, yielding a bootstrap distribution for that interval. The bootstrap median provides a stable PRCC estimate for that specific time step. We then integrate all local PRCC results across ~550 time steps to compute the median and interquartile range (IQR) for each parameter–flux pair.

(b) Effect size (Δrange).
To complement the correlation-based PRCC, we compute the physical amplitude of each parameter’s effect using Δrange metrics. In the local analysis, all accepted ensembles within a given time step are sorted by parameter X and divided into quantile bins; the difference between median fluxes in the highest and lowest bins define the local effect size (ΔY):

ΔY = median(Y|Xhigh) - median(Y|Xlow) (4)

ΔY (in Tmol/yr) represents the characteristic amplitude by which that parameter can shift the flux within that specific interval. Within each time step, ΔY is also bootstrapped by resampling accepted ensembles and recalculating ΔY 1000 times, producing a bootstrap ΔY distribution whose median represents the interval-specific effect size. Local ΔY values are likewise summarized across ~550 time steps (median ± IQR), where the reported 550 Ma statistics are derived from the set of per-interval bootstrap medians, reflecting long-term variability in effect magnitude.

(c) Secular drivers (20 Myr binned Deming regression with bootstrap)
To evaluate the leverage of boundary conditions on long-term flux evolution over the Phanerozoic, time-series data were aggregated into 20 Myr bins, and binned mean fluxes (F: Friv, Fsil, Fdol) were related to binned boundary conditions (X: seawater [Mg], δ²⁶MgSW, δ²⁶Mgriv, Δdol–SW). We first applied a Deming (total least-squares) regression, which accounts for uncertainties in both X and F, to estimate the secular slope:

F = a + bX, b = dF/dX (5)

ΔF = b × (Xmax - Xmin) (6)

Where b is the long-term Deming slope (dF/dX), ΔF is the flux change implied by the secular range of X, and the coefficient of determination (R²) is computed from fitted versus observed fluxes to evaluate explanatory power. The secular analysis is limited to parameters that exhibit genuine secular trends and can therefore drive multi-million-year flux evolution, whereas time-invariant parameters such as ΔLT-SW or prior ranges only affect local solutions and are not evaluated in the secular framework. We then apply a nonparametric bootstrap that resamples the 20 Myr bins with replacement (N≈27 per iteration), recalculating the Deming slope and ΔF for each realization. The resulting bootstrap medians and IQRs quantify the robustness of long-term relationships to temporal sampling.

The combined framework thus links instantaneous and long-term controls: local PRCC and ΔY capture short-term parameter sensitivity within each interval, their time-integrated medians reveal persistent controls across the Phanerozoic, and the secular regression (ΔF and R²) quantifies the actual magnitude and explanatory power of slow boundary-condition variations.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "I am honored to have the opportunity to review this manuscript. Building on previous work on the role of supercontinent cycles in regulating seawater chemistry, Zhang et al. combine newly published data on Phanerozoic seawater magnesium concentrations and Mg isotopes, and apply an inverse model to further constrain the sources and sinks of seawater Mg fluxes. Overall, the manuscript contributes to advancing our understanding of the elemental evolution of ancient seawater and its driving mechanisms. However, my main concerns are that (1) the study interprets changes in the Phanerozoic seawater Mg/Ca ratio solely in terms of Mg fluxes, and (2) the correlation between supercontinent cycles and Mg flux does not appear to be particularly strong. I therefore recommend that the manuscript undergo major revision. My specific comments are as follows:\nLine 14: Ensure consistent usage of either “Mg-isotope” or “Mg isotope.”\nLine 14: Clarify the meaning of “dynamics.”\nLine 28: Specify “modern seawater.”\nLine 30: The minimum y-axis value in Fig. 1a is 1, so Mg/Ca ratios <1 during the Jurassic–Cretaceous and Cambrian–Carboniferous intervals cannot be distinguished from the plot.\nLine 37: Consider listing the major views on the dominant processes controlling Phanerozoic seawater Mg/Ca.\nLine 65: Mg isotope data are a key parameter in the model, but both the number of evaporite data points and the fitted seawater Mg isotope curve differ somewhat from those recommended by Xia et al. (2024). Please verify.\nLine 66: Add the missing reference.\nLine 70: Revise to “The modern global riverine weighted average δ²⁶Mg,” which is more accurate.\nLines 83–84: Provide a more detailed explanation of the sensitivity test in Fig. 3 so that readers can better understand it.\nLines 96–99: Clarify the rationale for these interval selections. Why 550–500 Ma rather than 550–490 Ma? Similarly, why 400–290 Ma rather than 400–280 Ma?\nLines 99–101: Specify that these are 68% confidence intervals.\nLines 119–120: As mentioned above, explaining Mg/Ca changes solely from Mg fluxes is problematic. For example, around 400 Ma, seawater Mg/Ca ratios were close to 1 (~1.5). It is difficult to reconcile how, during 400–340 Ma, reduced dolomite removal (Mg:Ca removal molar ratio = 1:1) combined with enhanced silicate fluxes (removing Mg but introducing Ca) would result in a rise in Mg/Ca. This suggests that additional processes, not mentioned in the manuscript, also influenced seawater Mg/Ca ratio.\nLines 137–138: The 50 Ma peak lags behind maxima in global seafloor production by more than 10–30 Myr. Although you explain this later, it is still inaccurate here to state “each.”\nLine 142: Correct to “Fig. 1b.”\nLines 142–144: Should seafloor spreading rates correlate with the magnitude of Mg-silicate fluxes? If so, why does the 120 Ma peak not correspond to a significantly elevated Mg-silicate flux?\nLines 162–164: Although the study compares dolomite flux with Phanerozoic dolomite abundance, two recent studies (Li et al., 2021; Husson & Coogan, 2023) suggest that dolomite abundance was highest in the Cambrian, which is inconsistent with the present results.\nLines 164–166: The correlations appear weak. For example, Fig. 1 indicates that 250–120 Ma corresponds to supercontinent breakup and dispersal, which coincided with rising sea level, yet dolomite flux clearly declined—contrary to the manuscript’s interpretation.\nLines 166–168: Explain why only selected intervals were used to evaluate correlations rather than the entire Phanerozoic dataset.\nLines 181–183: The statement is overly absolute. During the Cambrian, high Mg/Ca ratios coincided with high CO₂ levels.\nLine 256–257: Consider whether “f×FLT” would be more appropriate.", "reviewer_2": "Zhang et al. developed a Mg isotope based inverse model to constrain changes in Mg fluxes of rivers, marine Mg-silicate formation, and dolomitization in the Phanerozoic oceans. Based on their modeling results, they argue that the Phanerozoic marine Mg cycle, including seawater Mg contents and Mg/Ca ratios, has been affected by supercontinent cycles as the first-order control.\nThis is an interesting study, and the novelty of the research lies in the fact that these authors used the existing Mg records (Mg/Ca, [Mg], and δ26Mg) to deconvolve possible changes in major Mg source (rivers) and sink (Mg-silicates and dolomites) fluxes over the entire Phanerozoic. Although their main conclusion that major ion chemistry of past seawater has been controlled by the tectonic forcing (i.e., supercontinent cycles) is not entirely new, the model-reconstructed Mg fluxes of rivers, Mg-silicates, and dolomites are valuable and significant, if true. This study is of broad interest and suitable for Nature Communications.\nThe main concern I have for this study is that this manuscript does not provide sufficient validation for their model-generated flux estimates. The principle of the kind of isotope mass balance model developed in this study is well known in literature. At face value, this kind of model can refine elemental source and sink fluxes based on combined elemental and isotopic mass balance equations. However, these refined elemental flux estimates depend on how well end-member Mg isotope compositions (or the relevant Mg isotope fractionation factors) are constrained, as well as how confident we think the previous estimates cover the “true” values of these fluxes. Although these authors presented sensitivity test results, it is challenging to assess if the tested parameter space indeed has covered what is possible in nature. This is particularly true for Mg isotope fractionation factors and how they may have varied over space and time. When it comes to their results, it is surprising to see that riverine Mg flux has remained almost constant over the entire Phanerozoic. If this is true, does it mean the riverine flux has been almost constant for other major ions? I am not suggesting the result to be false, but it is certainly counterintuitive, given that weathering must have varied considerably during the Phanerozoic. I think independent evidence to back up the modelled fluxes is essential but inadequate in the current manuscript.", "reviewer_3": "In this paper, Zhang et al. investigate the drivers behind changes in Mg concentrations throughout the Phanerozoic. To do this, they constructed an isotope mass balance model and used the known ranges in size of the key fluxes and their isotopic values in a Monte Carlo ensemble, to try and reproduce an estimated seawater curve of Phanerozoic Mg concentrations ([Mg]) and isotopic values (δ26Mg ), discarding any model runs which did not match the curve +/- a 1% error. From their results they conclude that supercontinent cycles have driven changes in [Mg].\nI think the paper shows just how much information about the Earth system can be gleaned from the use of simple isotope mass balance modelling, so in that regard it is a pretty neat piece of work. But, I feel like the work is lacking in detail in some places, meaning that it is difficult to exactly work out (at the moment) how the model has been constructed. The work could also be more comprehensive in its exploration of some of the uncertainties regarding e.g. isotopic values, but also limitations of the modelling etc. I will detail my concerns, grouping them into major comments and then some minor, easily fixable suggestions.\nMajor comments\nWhen it comes to the results/discussion sections regarding changes to Mg/Ca and their drivers the sole focus is on variations in Mg. However, you need to include what is going on with [Ca] as inferred by either fluid inclusion data or modelling using Ca isotopes because there are some disparities between the timing and relative magnitude of changes between [Mg] and Mg/Ca, but if the input from riverine is largely invariant as suggested by your modelling, then Mg : Ca exchange is going to be 1:1 so then Fig 1ab / Fig 2a needs a bit more discussion, as this then also relates to the timings of the shift between aragonite and calcite seas you mention. This also then relates to the title you have given the paper – really it only focuses on [Mg] at the moment, barely scratching the surface of [Ca] never mind any other seawater chemistry.\nI can understand why you combine low temperature (LT) and high temperature (HT) Mg removal via silicate formation for ease, but I’m guessing this means that your Fsil then has the same isotopic fractionation as the LT flux? I’m having to guess here because, unless I’ve completely missed it, this isn’t explicitly stated. But this means that your results might not be correct because if you keep these fluxes separate, a large HT flux with no isotopic fractionation and a smaller LT flux with some isotopic fractionation (and vice versa - e.g. what happens if, contrary to the rest of the studies, the Shalev paper referenced in the supplement is correct and the HT removal flux is lower than LT) is going to be different than if these are combined with a single fractionation value. Plus splitting these out would provide useful insights into whether LT or HT is key in driving Mg cycle changes. In either case, on L80 you say that you use equations 1 and 2 you compute the predicted [Mg] and δ26Mg but obviously this isn’t strictly true because you’re combining LT and HT, so you should write out new equations to illustrate this. It would also be incredibly helpful to have the isotopic values and their ranges used in a table as you have done with the flux sizes in the supplement.\nI think that it is important to include on L53 but also in the discussion section that Mg has a long residence time at present. think this is important to note, because the Xia et al. (2024) paper ignores some of the more rapid changes in δ26Mg , saying that this might be due to restricted basin effects because of the long residence time of Mg. However, other than for the Cenozoic and late Devonian - Permian, [Mg] is quite different and fluctuations are going to be driven by changes to residence time and steady state. Just because the residence time is long now doesn't mean it’s always been that long!\nI think there also needs to be more discussion regarding the riverine flux of Mg and its possible effect on [Mg] through the Phanerozoic – just because the mean from your model is largely invariant, it does not mean that there can’t be substantial variations per timestep between the bounds of the confidence intervals, which roughly matches (in terms of Tmol/yr) the size of the variations between the confidence interval bounds for the dolomite flux.\nThere could be more in the way of additional sensitivity tests, here are two as an example:\n1. What are the results like if the δ26Mg is different through time? You don’t really say what value is used, is it -1.09 +/- 0.19‰? If it is, you’re then not fully including the uncertainties of either reference, e.g. the Tipper paper says ± 0.05‰ for modern δ26Mg and I couldn’t really see where the 0.19‰ in Xia et al. (2024) comes from, albeit I’ve had to just skim through the paper, but there is surely some uncertainty attached to that too. Added to this, the modern day value is based on a flux weighted average of river water - how is that going to change when you have more carbonate weathering, e.g. when sea level drops? This doesn’t necessarily need modelling, but it needs to be considered when deciding on what range to use for δ26Mg and for your discussion too – how might it affect your results? Bluth and Kump (1991) have paleogeographic maps which can be used to aid some of the discussion.\n2. What happens to the results if the lines based on [Mg] and δ26Mg are different, considering some of the spread in that data in Fig. 1?\nThere needs to be a bit more in the discussion as to how changing lithology will also affect the proportions of Mg vs Ca released from the continents during weathering (e.g. L106-107) – e.g. anoxic dolomite can be a source/sink of Mg but not of Ca (Berner, 2004 in AJS), as well as any factors governing some variations in mobility.\nPlease include in the Methods section in which programming language the model been constructed. Or is it just in Excel? The mean and CI results of the model are in Excel, but the code should be publicly available too, so that what you’ve done can be properly examined.\nThere could be some greater justification behind some of the ranges used for the isotope values. For example, on L246-249, what is the reason for using +/- 0.2‰ specifically? Is there some defined temperature effect on the isotopic fractionation as there is for elements (e.g. lithium isotope fractionation during its uptake when clays form is sensitive to temperature)? Also, on L253-256, why do you not use up to the full isotope fractionation of 1.6‰, other than that other studies have not used this?\nMinor comments\nL43 – Change DSM3 in the equation to ‘standard’ and then outside of the equation state that the standard is DSM3. Also note that this is an equation and renumber the others.\nThroughout – there are some instances where ‘the’ is missing, e.g. L147-148 ‘Moreover, the silicate flux is negatively correlated with the continental fragmentation index during Pangea assembly but positively…’.\nFigure 1 – please include the geological timescale on the bottom of the figure, as you have done for Figure 2.\nFigure 2b – please include a value for the y-axis where it meets the x-axis. I cannot tell exactly how low the riverine Mg flux goes.\nL256 – I’m guessing that this should be Fsil = FHT + FLT = f x FLT ? And not f x Fsil ?\nL257 and 258 – please add formation after ‘silicate’\nL267-280 – Some of this text is more of a discussion and should sit there, rather than in the Methods section. But also, specifically on L276, its variation is very limited in the context of the modelling work you have conducted."}, "author_response": "-Reply to the Editor:\nComment: While all reviewers find promise in your approach, further analysis of your data will\nbe required to demonstrate the robustness of your conclusions to the level necessary for eventual\npublication in Nature Communications. We will need to see further analysis and discussion of the\nCa data and how this also relates to supercontinent cycling and seawater chemistry more broadly.\nIndependent validation for the modelled flux estimates, alongside sensitivity analysis and a full\ntreatment of uncertainties and limitations (including those related to currently simplified systems)\nwill also be required. Finally, we will need to see a significantly expanded description of your\nmethods and model, fully justifying your approach and assumptions.\nI note that you initially submitted this in short format, please take full advantage of Nature\nCommunications' more generous limits (5000 words, 10 figures/tables main text, unlimited\nmethods and supplementary information), which should allow plenty of scope to expand on all\nnecessary points.\nResponse: Thank you very much for forwarding the reviewers’ reports and providing us the\nopportunity for a thorough revision. We have used the reviewers’ constructive feedback to make\nsubstantial improvements in the manuscript’s scope, structure, and analytical transparency. We\nbelieve the revised manuscript satisfactorily addresses all concerns raised and provides a more\nrigorous and well-justified work. Reviewers' comments appear in black, and our responses are\nprovided in blue font. All revisions in the manuscript are marked to indicate which reviewer\ncomment each change addresses. Below is a concise overview of the main revisions:\n(1) We revised our figures, conducted analyses, and added discussions to address how Ca\ndata is related to supercontinent cycles and seawater chemistry (L139-157, 275-280).\nFirst, fluid inclusion seawater [Ca] data have been added to Fig. 1 to show that seawater [Mg]\nco-varies with seawater Mg/Ca while oppositely varying with [Ca]. Second, since seawater\n[Mg] and [Ca] were independently reconstructed, strong correlations between seawater [Mg]\nand Mg/Ca during seawater [Mg] increase and decrease periods demonstrate that long-term\nseawater [Mg] changes can account for the primary secular patterns in seawater Mg/Ca that\ndefine transitions between aragonite and calcite seas. Third, we also admit that an independent\nmarine Ca cycle via calcium carbonate burial likely also contributed to seawater Mg/Ca\nvariation during seawater [Mg] and Mg/Ca relatively stable periods, which requires future\nconstraints from the marine Ca cycle.\n(2) Independent validations for the modelled flux estimate and sensitivity analysis are also\nadded and conducted (L74-85, 102-125, 362-480, new Figs 3-4). Full parameter ranges and\nassociated uncertainties are now discussed in both the main text and the Methods section, and\nare summarized in the Supplementary Tables 1-2. Moreover, we conducted systematic\nsensitivity tests, evaluating the temporal and secular variations of all outputs relative to each\nparameter. Our sensitivity test results show that the direction of flux changes is primarily\ngoverned by δ²⁶Mg , whereas the magnitude is mainly controlled by seawater [Mg], with the\nSW\nprescribed range of F additionally modulating its own variability. Thus, our conclusions are\nriv\nrobust and not sensitive to the exact parameter ranges. Step-by-step sensitivity test results are\nshown in Supplementary Figs. 2-7, and integrated temporal and secular test results are\npresented in new Figs. 3-4 of the main text.\n(3) We expanded the description of our methods and modeling approach. We re-conducted\nour model with separated high- and low-temperature Mg-silicate fluxes (supplementary Fig.\nS2). We also conducted additional situations that assigned higher low-temperature fluxes,\nwhich produced results nearly identical to our main findings (supplementary Fig. S10).\nMoreover, we expanded our methods part (L74-125), and discussed the model limitations,\nincluding the Ca-dominated periods and the relatively invariant riverine Mg flux (L139-176).\n(4) All other moderate and minor issues raised by the reviewers have also been addressed.\nWe added discussions comparing our reconstructed dolomite flux with dolomite abundance\nrecords (Figs. 5, S13), updated our fitted curve to incorporate all available fluid inclusion\nseawater [Mg] data (revised Fig. 1), and performed a thorough check and revision throughout\nour manuscript.\n-Reply to Reviewer #1:\nGeneral Comment: I am honored to have the opportunity to review this manuscript. Building on\nprevious work on the role of supercontinent cycles in regulating seawater chemistry, Zhang et al.\ncombine newly published data on Phanerozoic seawater magnesium concentrations and Mg\nisotopes, and apply an inverse model to further constrain the sources and sinks of seawater Mg\nfluxes. Overall, the manuscript contributes to advancing our understanding of the elemental\nevolution of ancient seawater and its driving mechanisms. However, my main concerns are that (1)\nthe study interprets changes in the Phanerozoic seawater Mg/Ca ratio solely in terms of Mg fluxes,\nand (2) the correlation between supercontinent cycles and Mg flux does not appear to be\nparticularly strong. I therefore recommend that the manuscript undergo major revision. My specific\ncomments are as follows\nResponse: We sincerely thank Reviewer #1 for the thoughtful and constructive evaluation of our\nmanuscript. Main concerns are carefully considered and addressed as follows.\n(1) We demonstrate that first-order drivers of Phanerozoic Mg flux variations are also\nresponsible for seawater Mg/Ca (L139-157). Seawater [Mg] and [Ca] were independently\nreconstructed by measuring Mg and Ca concentrations in fluid inclusions of primary marine\nhalite (Horita et al., 2002; Lowenstein et al., 2001), hence, Mg/Ca ratios arise from two\nindependently determined quantities. To demonstrate that seawater [Mg] faithfully tracks\nseawater Mg/Ca, we plotted seawater [Mg] and Mg/Ca. Their high correlations (R² = 0.93–\n0.99; Fig. S11) during seawater [Mg] increase and decrease periods (550-390, 310-200, 90-0\nMa) demonstrate that they are controlled by the same drivers, namely seafloor alteration and\ndolomitization, which both remove Mg while adding Ca to the oceans. This coupling allows\nMg cycling alone to reproduce the broad structure of seawater Mg/Ca evolution. Similar\napproaches have been widely adopted in reconstructions of the Cenozoic Mg/Ca rise\n(Gothmann et al., 2017; Higgins and Schrag, 2015; Hu et al., 2022) and in explanations for the\nextensive dolomitization at the Permian–Triassic boundary (Hu et al., 2021). By contrast, low\ncorrelation between seawater [Mg] and Mg/Ca during seawater [Mg] fluctuation periods (500-\n490, 200-90 Ma; R² = 0.01–0.53; Fig. S11). Since seawater [Mg] varied marginally while\nseawater [Ca] varied significantly during those periods (Fig. 1b in the main text), Ca-related\nprocesses likely also influenced seawater Mg/Ca during those periods. For example, the rise\nof pelagic calcification since the Mesozoic likely contributed substantially to the observed\ndecline in seawater [Ca] (Ridgwell, 2005).\n(2) We provide evidence that Mg fluxes are highly related to the supercontinent cycle, and\nthe weak correlation suggests indirect control and the influence of multiple factors.\nFirst, Mg-silicate flux is positively correlated with the continental configuration index (Figs.\n2, S11), consistent with enhanced seafloor spreading and hydrothermal alteration during\nsupercontinent breakup and reduced hydrothermal activity during assembly. Moreover, the\nalignment of major peaks in seafloor production and Mg-silicate flux further supports a\nmechanistic link between Mg-silicate flux and (Fig. 2 in the main text). However, a direct\ncorrelation between Mg-silicate flux and seafloor production rate is absent due to temporal\nlags between spreading, crust formation, and subsequent alteration (L201-213).\nSecond, dolomite flux is positively correlated with modelled pCO₂ (550–340 and 250–0 Ma)\nand reconstructed sea level (460–340 and 100–0 Ma) during most Phanerozoic intervals (Figs.\n5, S11), consistent with the control of climate and sea level on dolomitization. Moreover,\ndolomite flux is generally high in the early Paleozoic and decreases since the Mesozoic,\nconsistent with the secular evolution pattern of pO ; elevated dolomite flux during the late\nPangaea assembly. These complicated model of dolomite flux is consistent with previous\nstudies that dolomitization is influenced by multiple factors, including sea level (Holland,\n2005), warm climates (Arvidson and Mackenzie, 1999; Xia et al., 2024), large-scale restriction\n(Hu et al., 2021; Warren, 2010), oceanic anoxia (Li et al., 2021), and microbial activity (Shang,\n2023), and all of which evolved significantly over the Phanerozoic. Even so, our dolomite\nflux curve highly aligns with previous dolomite abundance reconstructions (L214-253).\nComment 1: Line 14: Ensure consistent usage of either “Mg-isotope” or “Mg isotope.” Clarify\nthe meaning of “dynamics.”\nResponse: We now consistently use “Mg isotope” throughout the manuscript and have removed\nthe term “dynamics” to avoid ambiguity.\nComment 2: Line 28: Specify “modern seawater.”\nResponse: Thanks, and it is specified now.\nComment 3: Line 30: The minimum y-axis value in Fig. 1a is 1, so Mg/Ca ratios <1 during the\nJurassic–Cretaceous and Cambrian–Carboniferous intervals cannot be distinguished from the plot.\nResponse: We have adjusted the y-axis range in Fig. 1a so that Mg/Ca values <1 intervals can\nnow be clearly distinguished.\nComment 4: Line 37: Consider listing the major views on the dominant processes controlling\nPhanerozoic seawater Mg/Ca.\nResponse: We expanded the text and listed that fluctuations in the mid-ocean ridge hydrothermal\nflux, the rate of dolomitization, and clay authigenesis have been previously proposed to account\nfor Phanerozoic seawater Mg/Ca variation (L33-40).\nComment 5: Line 65: Mg isotope data are a key parameter in the model, but both the number of\nevaporite data points and the fitted seawater Mg isotope curve differ somewhat from those\nrecommended by Xia et al. (2024). Please verify.\nResponse: Thank you for pointing this out. We have re-examined the dataset and updated the\nfigure accordingly. The evaporite Mg isotope data now match those compiled by Xia et al. (2024)\nafter excluding non-marine samples (identified by 87/86Sr ratios). The fitted seawater Mg isotope\ncurve also follows Xia et al. (2024), except that we did not include the slight rise in the past ~9\nMyr because most previous studies reveal a constant late-Cenozoic seawater Mg isotope\ncomposition curve (Gothmann et al., 2017; Higgins and Schrag, 2015; Hu et al., 2022). This\nclarification has been added to the revised manuscript (L304-310).\nComment 6: Line 66: Add the missing reference.\nResponse: Reference is added.\nComment 7: Line 70: Revise to “The modern global riverine weighted average δ²⁶Mg,” which is\nmore accurate.\nResponse: Revised as suggested.\nComment 8: Lines 83–84: Provide a more detailed explanation of the sensitivity test in Fig. 3 so\nthat readers can better understand it.\nResponse: We have expanded the description of our sensitivity analyses in the main text (L103-\n125), and added new Figs. 3–4 in the main text to display sensitivity test results. Additional details\nare also shown in the Methods (L386-480) and supplementary Figs. S4–9.\nBriefly, we assess two complementary sensitivity tests: (1) local sensitivity at each time step and\n(2) secular sensitivity across the Phanerozoic. The local results show which parameters restrict\nflux magnitudes within a given interval, whereas the secular analysis identifies which boundary\nconditions explain their long-term evolution. Both local and secular sensitivity tests show that\nδ²⁶Mg and seawater [Mg] dominantly control the variation direction and magnitude of all fluxes,\nSW\nand F is also influenced by its prior range. These analyses demonstrate that the inferred controls\nriv\non all three Mg fluxes mainly arise from the intrinsic mass-balance framework and seawater [Mg]-\nδ²⁶Mg variations, rather than from prior ranges of input-flux or isotope-fractionation factors.\nSW\nComment 9: Lines 96–99: Clarify the rationale for these interval selections. Why 550–500 Ma\nrather than 550–490 Ma? Similarly, why 400–290 Ma rather than 400–280 Ma?\nResponse: Thank you for raising this point. The interval definitions are based on the Generalized\nAdditive Model (GAM)–fitted seawater [Mg] curve. The boundaries where the first derivative of\nthe fitted curve switches sign are selected to divide these intervals. For example, seawater [Mg]\nshows a clear decline from 550 to ~500 Ma, followed by the onset of a slight rise thereafter;\ntherefore, 550–500 Ma represents a distinct decreasing interval. These intervals are consistent with\nvariations of seawater Mg/Ca and oscillations between aragonite and calcite seas (L127-136).\nComment 10: Lines 99–101: Specify that these are 68% confidence intervals\nResponse: We have specified that these are ranges of average Mg fluxes instead of 68%\nconfidence intervals in L136-138.\nComment 11: Lines 119–120: As mentioned above, explaining Mg/Ca changes solely from Mg\nfluxes is problematic. For example, around 400 Ma, seawater Mg/Ca ratios were close to 1 (~1.5).\nIt is difficult to reconcile how, during 400–340 Ma, reduced dolomite removal (Mg:Ca removal\nmolar ratio = 1:1) combined with enhanced silicate fluxes (removing Mg but introducing Ca)\nwould result in a rise in Mg/Ca. This suggests that additional processes, not mentioned in the\nmanuscript, also influenced seawater Mg/Ca ratio.\nResponse: Thank you for raising this important point. We agree that interpreting Mg/Ca variations\nrequires considering both the Mg and Ca cycles. We have expanded the discussion, incorporated\nseawater [Ca] data, and added analyses to address this issue (L139-157).\n(1) Our inverse model solves the marine Mg budget strictly according to seawater Mg mass\nbalance equations; hence, during 400–340 Ma, the total source minus sink is positive\n(Flux = F –F – F > 0, Fig. 2e in the main text), which explains increasing seawater\nsource-sink riv dol sil\n[Mg] and is also consistent with increasing Mg/Ca (Fig. 2a in the main text).\n(2) We demonstrate that variations in seawater [Mg] and Mg/Ca were driven by the same\ngeological periods in most Phanerozoic periods. Seawater [Mg] and [Ca] were reconstructed\nindependently using fluid inclusions in early-stage primary marine halite (Horita et al., 2002;\nLowenstein et al., 2001); hence, seawater [Mg] and Mg/Ca variations are controlled by the\nsame processes if they are correlated. Seawater [Mg] and Mg/Ca are positively correlated\nduring the seawater [Mg] increase and decrease periods (550-500, 390-310, 310-200, 90-0 Ma;\nR² = 0.93–0.99; Fig. S11), demonstrating that they are controlled by the same drivers. This is\nreasonable since two major geological processes, seafloor alteration and dolomitization, would\nboth remove Mg while adding Ca into the oceans (L139-149).\n(3) Ca-related processes are also discussed. Seawater [Mg] and Mg/Ca are not well correlated\nduring the seawater [Mg] fluctuation periods (500-390, 200-90 Ma; R² = 0.01–0.53; Fig. S11),\nsuggesting decoupled [Mg] and Mg/Ca. Since seawater [Ca] varied significantly during these\nperiods (Fig. 1b), suggesting independent Ca-related processes may also be responsible for\nseawater Mg/Ca variations. Calcium carbonate burial is an important process controlling\nseawater Mg/Ca by influencing Ca burial but not affecting Mg. For example, the rise of pelagic\ncalcification since the Mesozoic likely contributed substantially to the observed decline in\nseawater [Ca] (Ridgwell, 2005) (L149-157).\nComment 12: Lines 137–138: The 50 Ma peak lags behind maxima in global seafloor production\nby more than 10–30 Myr. Although you explain this later, it is still inaccurate here to state “each.”\nResponse: We have revised our manuscript accordingly in L1201-204.\nComment 13: Line 142: Correct to “Fig. 1b.”\nResponse: Revised as suggested.\nComment 14: Lines 142–144: Should seafloor spreading rates correlate with the magnitude of\nMg-silicate fluxes? If so, why does the 120 Ma peak not correspond to a significantly elevated\nMg-silicate flux?\nResponse: Yes, seafloor spreading rates will influence Mg-silicate fluxes because faster spreading\ngenerates more fresh oceanic crust, expanding the area available for seafloor–seawater reactions\nand thereby enhancing Mg uptake by alteration (Coogan and Gillis, 2018).\nThe 120 Ma seafloor production peak corresponds to a persistent increase of Mg-silicate flux\nbetween 140 and 40 Ma, according to our updated model outcomes that used full fluid-inclusion\nseawater [Mg] data fitted the seawater [Mg] curve. This protracted interval is consistent with\nprevious hydrothermal reconstructions (Müller et al., 2013), and likely reflects the cumulative\neffect of sustained high seafloor production during Pangea breakup (L193-200, 207-209).\nComment 15: Lines 162–164: Although the study compares dolomite flux with Phanerozoic dolomite abundance, two recent studies (Li et\nal., 2021; Husson & Coogan, 2023) suggest that dolomite abundance was highest in the Cambrian, which is inconsistent with the present results.\nResponse: Thank you for raising this important point, and we added comparisons, analyses, and\nexplanations (L230-241, new Fig. 5, and new supplementary Fig. S13) to show that our dolomite\nflux is broadly consistent with previous dolomite abundance records.\n(1) First, dolomite flux is a proxy different from previous dolomite abundance records.\nPrevious studies have reconstructed the proportion of dolomite relative to total carbonate,\nnamely dolomite abundance, which is different from the absolute dolomite burial flux. If total\ncarbonate burial in the early Cambrian was relatively low, high dolomite abundance would not\nnecessarily translate into a proportionally high dolomite flux. In this context, our modelled\ndolomite flux for the early Cambrian is consistent with the dolomite abundance record.\n(2) Our dolomite flux does not exclude the highest value in the Cambrian, considering the\nconfidence intervals (Fig. 2d in the main text).\n(3) The secular trend of our dolomite flux is highly consistent with previous dolomite\nabundance reconstructions. To verify the robustness of our dolomite abundance records, we\ncompared our results with all dolomite abundance curves. We standardized all dolomite\nabundance curves and then conducted bootstrapped resampling Spearman’s rank correlation\n(ρ). Results show that our dolomite flux curve is highly consistent with most dolomite\nabundance records, even though these records have different temporal resolutions, assumptions,\nand uncertainties (Figs. 5, S12; ρ > 0.67 for three curves). The remaining two records exhibit\nweaker agreement due to intrinsic data limitations, with one dataset having extremely low\ntemporal resolution (only six spots), and the other displaying a distinct stratigraphic pattern\nrelative to the other curves due to the limited sample dataset (Fig. 5 in the main text).\nComment 16: Lines 164–166: The correlations appear weak. For example, Fig. 1 indicates that\n250–120 Ma corresponds to supercontinent breakup and dispersal, which coincided with rising sea\nlevel, yet dolomite flux clearly declined—contrary to the manuscript’s interpretation.\nResponse: Yes, and we deleted this sentence and just stated that dolomite flux is correlated to\npCO and sea level during specific periods in the revised manuscript (L2466-249).\ncorrelations rather than the entire Phanerozoic dataset.\nResponse: Thanks for pointing this out, and only selected intervals were used to evaluate\ncorrelations because likely the controlling factors on dolomite flux varied over the Phanerozoic.\n(1) Many factors have been proposed to account for dolomite formation, but no consensus\nhas been reached. High sea level (Holland and Zimmermann, 2000; Turchyn and DePaolo,\n2019), warm climate (Arvidson and Mackenzie, 1999; Xia et al., 2024), ocean anoxia (Li et\nal., 2021), microbial activity (Shang, 2023), the dissolution-reprecipitation process(Kim et al.,\n2023), and large-scale restriction (Hu et al., 2021; Warren, 2010) have been proposed to\npromote dolomite formation. Microbial activity, elevated temperatures, and restricted\nhydrological conditions are suggested to accelerate dolomite nucleation and crystal growth,\nwhile anoxic conditions may further stimulate microbial metabolism that catalyzes\ndolomitization reactions, and high sea level promotes the development of extensive shallow\nmarine platforms and accelerates reaction kinetics via enhanced evaporation and elevated\ntemperatures, particularly in warm, semi-restricted settings (L214-229).\n(2) Our results show that controls on dolomite flux varied over the Phanerozoic. Overall,\ndolomite flux shows strong positive correlations with modelled pCO₂ (550–340 and 250–0 Ma)\nand reconstructed sea level (460–340 and 100–0 Ma) during most Phanerozoic intervals (Figs.\n2, S11). Moreover, dolomite flux is generally high in the early Paleozoic, while it decreases\nsince the Mesozoic, which is consistent with the secular evolution pattern of pO (New Fig. 5),\n2\nindicating the potential influence of oceanic anoxia and microbial activity (Li et al., 2021;\nShang, 2023). Elevated dolomite flux during the late Paleozoic–early Mesozoic may be\nattributed to the Paleotethys ocean restriction during the Pangaea assembly (Hu et al., 2021)\n(L244-269). Thus, different controls dominantly regulate dolomite formation during different\nperiods. For instance, when pO is high and relatively constant, dolomite formation is mainly\n2\ncontrolled by sea level and climate. When the supercontinent assembles, global continental\nconfiguration changes significantly and likely contribute to dolomite formation (L230-253).\nComment 18: Lines 181–183: The statement is overly absolute. During the Cambrian, high Mg/Ca\nratios coincided with high CO₂ levels.\nResponse: The statement has been revised.\nComment 19: Lines 256–257: Consider whether “f×FLT” would be more appropriate.\nResponse: Thank you and revised.\n-Reply to Reviewer #2:\nComment 1: Zhang et al. developed a Mg isotope based inverse model to constrain changes in\nMg fluxes of rivers, marine Mg-silicate formation, and dolomitization in the Phanerozoic oceans.\nBased on their modeling results, they argue that the Phanerozoic marine Mg cycle, including\nseawater Mg contents and Mg/Ca ratios, has been affected by supercontinent cycles as the first-\norder control.\nThis is an interesting study, and the novelty of the research lies in the fact that these authors used\nthe existing Mg records (Mg/Ca, [Mg], and δ26Mg) to deconvolve possible changes in major Mg\nsource (rivers) and sink (Mg-silicates and dolomites) fluxes over the entire Phanerozoic. Although\ntheir main conclusion that major ion chemistry of past seawater has been controlled by the tectonic\nforcing (i.e., supercontinent cycles) is not entirely new, the model-reconstructed Mg fluxes of\nrivers, Mg-silicates, and dolomites are valuable and significant, if true. This study is of broad\ninterest and suitable for Nature Communications.\nResponse: We thank Reviewer 2 for the constructive evaluation of our manuscript and for\nrecognizing the suitability of our work for this journal.\nComment 2: The main concern I have for this study is that this manuscript does not provide\nsufficient validation for their model-generated flux estimates. The principle of the kind of isotope\nmass balance model developed in this study is well known in literature. At face value, this kind of\nmodel can refine elemental source and sink fluxes based on combined elemental and isotopic mass\nbalance equations. However, these refined elemental flux estimates depend on how well end-\nmember Mg isotope compositions (or the relevant Mg isotope fractionation factors) are\nconstrained, as well as how confident we think the previous estimates cover the “true” values of\nthese fluxes. Although these authors presented sensitivity test results, it is challenging to assess if\nthe tested parameter space indeed has covered what is possible in nature. This is particularly true\nfor Mg isotope fractionation factors and how they may have varied over space and time.\nResponse: We appreciate the reviewer’s careful assessment and constructive concerns regarding\nmodel validation and parameter uncertainty. To address your concern, we (1) provide a detailed\ndescription of input parameters that involve wider uncertainties and added supplementary tables\nto show that the parameters we used cover wide ranges that were discovered in nature and used in\nprevious studies (L362-379 in the Methods, and supplementary Tables S1-S2), and (2)\nconducted systematic sensitivity tests to show that prior ranges of isotope fractionation or flux\nranges played minor roles in regulating model outcomes (L103-125, details in L386-480 of the\nMethods, added new Figs. 3–4 and supplementary Figs. S4–S9).\n(1) We adopt ranges of parameters that have been frequently used in previous studies and\neven wider ranges (L74-92, 362-379). All input parameters and their uncertainties, including\nthe full ranges of fluxes and Mg isotope fractionation factors, are derived from previous work\nspanning the Cenozoic to the Proterozoic (supplementary Tables S1–S2). These values have\nbeen widely applied and successfully tracked the marine Mg cycle during critical geological\nperiods (Gothmann et al., 2017; Higgins and Schrag, 2015; Hu et al., 2021; Hu et al., 2022;\nShalev et al., 2019; Xia et al., 2024). Uncertainties for Δ (± 0.1‰), Δ (± 0.2‰) and\ndol–SW LT–SW\nδ²⁶Mg (± 0.1‰) are included according to uncertainties from calibrated equation, natural\nriv\nsample variations, and lithology variations, and Δ (± 0.1‰) is set to vary with\ndol–SW\ntemperature. As such, the explored parameter space reflects the full set of values currently\nsupported by empirical and experimental constraints (L74-85, 362-379).\n(2) We conducted two complementary sensitivity tests: local sensitivity at each time step and\nsecular sensitivity across the Phanerozoic to evaluate the effect of model parameters on\noutcomes (L103-125, 386-480, Figs. 3-4, supplementary Figs. 4-9). The local results show\nwhich parameters restrict flux magnitudes within a given interval, whereas the secular analysis\nidentifies which boundary conditions explain their long-term evolution. Both local and secular\nsensitivity tests show that δ²⁶Mg and seawater [Mg] dominantly control the variation\nSW\ndirection and magnitude of all fluxes, and F is also influenced by its prior range. These\nriv\nanalyses demonstrate that the inferred controls on all three Mg fluxes mainly arise from the\nintrinsic mass-balance framework and seawater [Mg]-δ²⁶Mg variations, rather than from\nSW\nprior ranges of input-flux or isotope-fractionation factors.\nComment 3: When it comes to their results, it is surprising to see that riverine Mg flux has\nremained almost constant over the entire Phanerozoic. If this is true, does it mean the riverine flux\nhas been almost constant for other major ions? I am not suggesting the result to be false, but it is\ncertainly counterintuitive, given that weathering must have varied considerably during the\nPhanerozoic. I think independent evidence to back up the modelled fluxes is essential but\ninadequate in the current manuscript.\nResponse: We agree that weathering flux likely varied significantly over the Phanerozoic. We\nadded explanations that our model is insensitive to the magnitude of riverine flux and we do not\nexclude that riverine flux may also play important roles in regulating seawater Mg/Ca during\ncritical geological periods (L158-176).\n(1) First, the apparent invariant weathering flux arises likely because our model is\nintrinsically less sensitive to riverine Mg input, whereas changes in dolomitization and Mg-\nsilicate formation more efficiently satisfy elemental and isotopic mass balance in equations #1-2\nin the main text for the observed seawater [Mg] and δ²⁶Mg trajectories.\nSW\n(2) The full output ranges of our modelled F actually allow substantial variation in riverine\nriv\nMg flux, even though we mainly discuss the mean values.\n(3) Riverine input is inherently less effective at shifting seawater Mg/Ca because, unless\ncontinental lithology changes dramatically, rivers deliver Mg and Ca in broadly similar proportions,\nwhile dolomite formation and seafloor alteration simultaneously remove Mg and add Ca.\nMoreover, Phanerozoic continental lithologic composition appears to have remained broadly\nstable (Peters et al., 2018; Walton et al., 2023). Thus, riverine input would be less important than\nMg-silicate and dolomite formation during most Phanerozoic periods (L165-168).\n(4) We agree that continental weathering likely varied substantially through the Phanerozoic\nand we do not completely exclude the potential effect of riverine input on seawater Mg/Ca.\nOur model indicates that riverine input could have been particularly influential during intervals\nsuch as 550–480 Ma, 350–320 Ma, and 50–0 Ma—times associated with enhanced silicate\nweathering in the Early Paleozoic (Wei et al., 2024), mafic crust exposure during Pangea assembly\n(Yang et al., 2023), and intensified uplift during the Alpine–Himalayan orogeny (Raymo and\nRuddiman, 1992), yet the exact magnitudes of these variations remain difficult to quantify\nprecisely in our model. We believe better constraints on weathering flux would further clarify\nthe role of riverine flux on seawater chemistry in future studies (L158-164).\n-Reply to Reviewer #3:\nComment: In this paper, Zhang et al. investigate the drivers behind changes in Mg concentrations\nthroughout the Phanerozoic. To do this, they constructed an isotope mass balance model and used\nthe known ranges in size of the key fluxes and their isotopic values in a Monte Carlo ensemble, to\ntry and reproduce an estimated seawater curve of Phanerozoic Mg concentrations ([Mg]) and\nisotopic values (δ26Mgsw), discarding any model runs which did not match the curve +/- a 1%\nerror. From their results they conclude that supercontinent cycles have driven changes in [Mg].\nI think the paper shows just how much information about the Earth system can be gleaned from\nthe use of simple isotope mass balance modelling, so in that regard it is a pretty neat piece of work.\nBut, I feel like the work is lacking in detail in some places, meaning that it is difficult to\nexactly work out (at the moment) how the model has been constructed. The work could also be\nmore comprehensive in its exploration of some of the uncertainties regarding e.g. isotopic values,\nbut also limitations of the modelling etc. I will detail my concerns, grouping them into major\ncomments and then some minor, easily fixable suggestions.\nResponse: We thank Reviewer #3 for recognizing the value of our approach and for the positive\nassessment of the scientific contribution. We have addressed the reviewer’s concerns regarding\nmodel transparency, parameter uncertainties, and the need to more clearly articulate the model\nstructure and its limitations by expanding the manuscript as detailed below.\nMajor Comment 1: When it comes to the results/discussion sections regarding changes to Mg/Ca\nand their drivers the sole focus is on variations in Mg. However, you need to include what is going\non with [Ca] as inferred by either fluid inclusion data or modelling using Ca isotopes because there\nare some disparities between the timing and relative magnitude of changes between [Mg] and\nMg/Ca, but if the input from riverine is largely invariant as suggested by your modelling, then Mg :\nCa exchange is going to be 1:1 so then Fig 1ab / Fig 2a needs a bit more discussion, as this then\nalso relates to the timings of the shift between aragonite and calcite seas you mention. This also\nthen relates to the title you have given the paper – really it only focuses on [Mg] at the moment,\nbarely scratching the surface of [Ca] never mind any other seawater chemistry.\nResponse: We agree that variations in Ca as well as Mg contribute to understanding seawater\nMg/Ca evolution. To address the reviewer’s concerns, we have broadened our discussion to\nincorporate available constraints on seawater Ca by including fluid-inclusion [Ca] data and\nevaluating the co-evolution of seawater [Mg] and Mg/Ca (L139-157; updated Fig. 1 in the main\ntext; new supplementary Fig. S11).\n(1) We included seawater [Ca] data that are derived from primary marine halite fluid\ninclusions. Seawater [Mg] and Mg/Ca varied in the same direction while seawater [Mg] and\n[Ca] varied oppositely over the Phanerozoic (Fig. 1 in the main text), which is consistent with\nthe assumption that the same processes controlled these variations.\n(3) We further quantified the co-variation of seawater [Mg] and Mg/Ca and demonstrated\nthat first-order drivers of Phanerozoic Mg flux variations are also responsible for\nseawater Mg/Ca (L139-149). Seawater [Mg] and [Ca] were independently reconstructed by\nmeasuring Mg and Ca concentrations in fluid inclusions of primary marine halite (Horita et al.,\n2002; Lowenstein et al., 2001), hence, Mg/Ca ratios arise from two independently determined\nquantities. To demonstrate that seawater [Mg] faithfully tracks seawater Mg/Ca, we plotted\nseawater [Mg] and Mg/Ca. Their high correlations (R² = 0.93–0.99; Fig. S11) during seawater\n[Mg] increase and decrease periods demonstrate that they are controlled by the same drivers,\nnamely seafloor alteration and dolomitization, which both remove Mg while adding Ca to the\noceans. This coupling allows Mg cycling alone to reproduce the broad pattern of seawater\nMg/Ca evolution. Similar approaches have been widely adopted in reconstructions of the\nCenozoic Mg/Ca rise (Gothmann et al., 2017; Higgins and Schrag, 2015; Hu et al., 2022) and\nin explanations for the extensive dolomitization at the Permian–Triassic boundary (Hu et al.,\n2021).\n(4) We also admit that independent Ca-dominated processes may be important in regulating\nseawater Mg/Ca during calcite sea periods (L146-154). Weak correlation between seawater\n[Mg] and Mg/Ca during seawater [Mg] fluctuation periods (500-490, 200-90 Ma; R² = 0.01–\n0.53; Fig. S11). Since seawater [Mg] and Mg/Ca varied marginally while seawater [Ca] varied\nsignificantly during those periods (Fig. 1 in the main text), Ca-related processes likely also\ninfluenced seawater Mg/Ca during those periods, such as calcium carbonate burial, which is a\nmajor Ca sink while it does not affect Mg. For example, the rise of pelagic calcification since\nthe Mesozoic likely contributed substantially to the observed decline in seawater [Ca]\n(Ridgwell, 2005). Currently, the marine δ⁴⁴/⁴⁰Ca record is not yet as complete and with large\nuncertainties (Jenkyns, 2012), and it would exponentially enhance the complication of our\nmass balance model if we introduced Ca and Ca isotopes. Moreover, seawater Mg/Ca remained\nrelatively stable during those calcite sea periods (~500-490, 200-90 Ma; Fig. 1 in the main\ntext). Thus, we choose to stick to the marine Mg cycle, while we believe further constraints\nfrom the marine Ca cycle would help better understand the seawater Mg/Ca fluctuation during\nthese periods with minor seawater Mg/Ca variations.\nMajor Comment 2: I can understand why you combine low temperature (LT) and high\ntemperature (HT) Mg removal via silicate formation for ease, but I’m guessing this means that\nyour Fsil then has the same isotopic fractionation as the LT flux? I’m having to guess here because,\nunless I’ve completely missed it, this isn’t explicitly stated. But this means that your results might\nnot be correct because if you keep these fluxes separate, a large HT flux with no isotopic\nfractionation and a smaller LT flux with some isotopic fractionation (and vice versa - e.g. what\nhappens if, contrary to the rest of the studies, the Shalev paper referenced in the supplement is\ncorrect and the HT removal flux is lower than LT) is going to be different than if these are combined\nwith a single fractionation value. Plus splitting these out would provide useful insights into whether\nLT or HT is key in driving Mg cycle changes. In either case, on L80 you say that you use equations\n1 and 2 you compute the predicted [Mg] and δ26Mgsw but obviously this isn’t strictly true because\nyou’re combining LT and HT, so you should write out new equations to illustrate this. It would\nalso be incredibly helpful to have the isotopic values and their ranges used in a table as you have\ndone with the flux sizes in the supplement.\nResponse: Yes, in our primary model, the LT and HT processes are combined, and the effective\nF therefore carries the isotopic fractionation of the LT component—only the LT fraction of F\nsil sil\ncontributes to Mg isotope fractionation. To address your concern, we revised our model and\nconducted new simulations.\n(1) We separated HT and LT processes in our updated model. We revised our model and set\nindependent input ranges for HT and LT processes now (L86-92; supplementary Table S2).\nThe results are identical to our previous results and totally address your concern. Moreover,\nour modelled HT and LT processes varied highly similarly over the Phanerozoic (L97-99,\nsupplementary Figs. 2-3). Thus, to reduce complications, we still combined the two processes\nin our discussions.\n(2) We conducted our model with the situation that the HT flux is lower than that of the LT\n(L119-122). Results are nearly identical to those of our primary results (supplementary Fig.\nS10), indicating that either higher or lower HT flux does not affect our results. This is\nreasonable since our sensitivity tests also show that all fluxes are insensitive to prior ranges of\nthese fluxes (Figs. 3-4 in the main text). Moreover, since prevalent studies claim a higher HT\nflux than that of the LT (supplementary Table S2), we follow this in our primary model\nsetting.\nMajor Comment 3: I think that it is important to include on L53 but also in the discussion section\nthat Mg has a long residence time at present. think this is important to note, because the Xia et al.\n(2024) paper ignores some of the more rapid changes in δ26Mgsw, saying that this might be due\nto restricted basin effects because of the long residence time of Mg. However, other than for the\nCenozoic and late Devonian - Permian, [Mg] is quite different and fluctuations are going to be\ndriven by changes to residence time and steady state. Just because the residence time is long now\ndoesn't mean it’s always been that long!\nResponse: We highly agree with you that a long residence time today does not imply that Mg has\nalways behaved conservatively, and have added the relevant discussion in the revised manuscript.\nNotably, unlike Xia et al. (2024), our framework explicitly incorporates changes in the residence\ntime of Mg by including seawater [Mg] evolution in equation #1 (L102-103). The residence time\nof Mg is calculated as the total mass of Mg in the global oceans divided by riverine input (or\ntotal sink), which is displayed here as the following figure for your reference, but does not need to\nbe shown in the manuscript.\nFig. 1 The variation of Mg residence time in the ocean over the Phanerozoic in our model.\nMajor Comment 4: I think there also needs to be more discussion regarding the riverine flux of\nMg and its possible effect on [Mg] through the Phanerozoic – just because the mean from your\nmodel is largely invariant, it does not mean that there can’t be substantial variations per timestep\nbetween the bounds of the confidence intervals, which roughly matches (in terms of Tmol/yr) the\nsize of the variations between the confidence interval bounds for the dolomite flux.\nResponse: We agree that invariant mean riverine flux does not mean that riverine Mg flux does\nnot vary through the Phanerozoic, and added discussions to clarify (L158-176).\n(1) As you suggested, the broad confidence ranges of riverine Mg flux explicitly allow large\ntimestep-to-timestep fluctuations comparable in magnitude to those of the dolomite flux.\n(2) The apparent invariant arises likely because our model is intrinsically less sensitive to\nriverine Mg input, whereas changes in dolomitization and Mg-silicate formation more\nefficiently satisfy elemental and isotopic mass balance in equations #1-2 in the main text for\nthe observed seawater [Mg] and δ²⁶Mg trajectories.\nSW\n(3) Riverine input is inherently less effective at shifting seawater Mg/Ca because, unless\ncontinental lithology changes dramatically, rivers deliver Mg and Ca in broadly similar\nproportions, while dolomite formation and seafloor alteration simultaneously remove Mg and\nadd Ca.\n(4) We also discussed the potential effect of riverine input on seawater Mg/Ca. Our model\nindicates that riverine input could have been particularly influential during intervals such as\n550–480 Ma, 350–320 Ma, and 50–0 Ma—times associated with enhanced silicate weathering\nin the Early Paleozoic (Wei et al., 2024), mafic crust exposure during Pangea assembly\n(Yang et al., 2023), and intensified uplift during the Alpine–Himalayan orogeny (Raymo and\nRuddiman, 1992), yet the exact magnitudes of these variations remain difficult to quantify\nprecisely in our model. We believe better constraints on weathering flux would further clarify\nthe role of riverine flux on seawater chemistry in future studies.\nMajor Comment 5: There could be more in the way of additional sensitivity tests, here are two\nas an example:\n1. What are the results like if the δ26Mg is different through time? You don’t really say what\nriv\nvalue is used, is it -1.09 +/- 0.19‰? If it is, you’re then not fully including the uncertainties of\neither reference, e.g. the Tipper paper says ± 0.05‰ for modern δ26Mg and I couldn’t really see\nriv\nwhere the 0.19‰ in Xia et al. (2024) comes from, albeit I’ve had to just skim through the paper,\nbut there is surely some uncertainty attached to that too. Added to this, the modern day value is\nbased on a flux weighted average of river water - how is that going to change when you have more\ncarbonate weathering, e.g. when sea level drops? This doesn’t necessarily need modelling, but it\nneeds to be considered when deciding on what range to use for δ26Mg and for your discussion\nriv\ntoo – how might it affect your results? Bluth and Kump (1991) have paleogeographic maps which\ncan be used to aid some of the discussion.\n2. What happens to the results if the lines based on [Mg] and δ26Mgsw are different, considering\nsome of the spread in that data in Fig. 1?\nResponse: These suggestions are very good and have been included in the revised manuscript.\nSuggestion 1 – variation of δ26Mg .\nriv\n(1) You are correct that the modern δ26Mg has a weighted global average value of −1.09 ± 0.05‰,\nriv\nwhich is revised and adopted in our model (L81-85).\n(2) The data used in Xia et al. (2024) is derived from their Fig. S11 as the following Figure.\nFig. 2. The evolution of the fraction of carbonate weathering (a) on the continent and its δ26Mg (b), yielding a\nδ²⁶Mg variation from −1.01 ± 0.15‰ to −0.82 ± 0.12‰ over the Phanerozoic, from Figure S11 in Xia et al.\nriv\n(2024).\n(3) We used a δ²⁶Mg that varied from −1.09 ± 0.1‰ to −0.90 ± 0.1‰ over the Phanerozoic\nriv\n(supplementary Fig. S1), shifted by ~ −0.08‰ from that in Xia et al. (2024) to match the\nmodern δ²⁶Mg value of −1.09‰. This range is larger than nearly all previous studies\nriv\n(supplementary Table S1), and corresponds to a large variation of the carbonate fraction from\n0.32 ± 0.04 to 0.46 ± 0.04. Although it remains unclear how δ²⁶Mg actually varied in the past,\nriv\nour sensitivity results show that model outputs are only weakly sensitive to δ²⁶Mg (new Figs\nriv\n3-4), so uncertainties in this parameter do not propagate strongly into the reconstructed flux\nhistories.\nSuggestion 2 – different lines of seawater [Mg] and δ26Mg from Fig 1.\nSW\n(1) The Phanerozoic δ26Mg curve has only reported one line with an uncertainty of 0.2‰\nSW\nin Xia et al. (2024). Thus, it is unclear if there are any other trajectories. While according to\nour sensitivity test (new Figs 3-4), both at each time step and the secular trend over the whole\nPhanerozoic, δ26Mg is the most (or second) factor that determines the direct and magnitude\nSW\nof all fluxes. Since our model is based on the current only curve, we stick to our results. If there\nare more secular δ26Mg curves reconstructed in the future, it is necessary to compare the\nSW\nresults derived from different curves and discuss their differences.\n(2) Major transitions between aragonite and calcite seas are obvious and would remain no\nmatter how different the seawater [Mg] curve is. In the former version, we fitted the\nseawater [Mg] curve but ignored its secondary fluctuation during calcite sea periods. Now we\nchave updated the curve by considering all fluid inclusion data (Fig. 1b). Our model results is\nconsistent with the former version, except that we identified Mg-silicate flux peaks during the\ncalcite sea periods (Fig. 2; which is reasonable since the curves during these periods are\nrevised). We could certainly construct an alternative curve that passes through every individual\ndata point, but curves with excessively high temporal resolution would be unrealistic and\nunnecessary because (1) the inherent analytical uncertainties in reconstructed seawater [Mg]\nand the poorly constrained timing of evaporite fluid-inclusion formation, and (2) this is not the\ngoal of our study. Our goal is to elucidate the long-term drivers of significant seawater Mg/Ca\nfluctuations between aragonite and calcite seas rather than to reproduce short-term fluctuations,\nand the fitted curve we use already captures the major structure of seawater Mg concentration\nthrough time. Although this curve may be refined and other trajectories are mathematically\npossible in the future, the first-order transitions between “aragonite seas” and “calcite seas”\nare robust.\nMajor Comment 6: There needs to be a bit more in the discussion as to how changing lithology\nwill also affect the proportions of Mg vs Ca released from the continents during weathering (e.g.\nL106-107) – e.g. anoxic dolomite can be a source/sink of Mg but not of Ca (Berner, 2004 in AJS),\nas well as any factors governing some variations in mobility.\nResponse: We appreciate this constructive suggestion, and we added corresponding discussions\nto discuss the potential variations in lithology and its effect on riverine input flux(L165-176).\n(1) The weathering of different lithologies provides different Mg-Ca ratios. Weathering of\nsilicate rocks supplies both Mg and Ca with typical Mg/Ca ratios of ~0.5–0.8, whereas\ncarbonate weathering releases much more Ca than Mg (Mg/Ca ≈ 0.1), and dolostone\nweathering contributes roughly equal proportions of the two elements. Changes in the global\ndistribution of these lithologies, including the extent of carbonate platforms, flood basalt\nprovinces, and silicate terrains, would therefore influence the Mg:Ca ratio of riverine inputs\nand, ultimately, seawater composition. We also note that anoxic dolomite may act as a Mg\nsource or sink without affecting Ca, as highlighted by Berner (2004), further contributing to\nspatial and temporal heterogeneity in weathering-derived fluxes.\n(2) Quantifying these lithologic effects, however, remains difficult within the current model\nframework because global-scale lithologic exposure cannot yet be constrained with sufficient\ntemporal resolution, and our model is comparatively insensitive to riverine Mg flux relative to\nthe Mg-silicate and dolomite sinks. Nevertheless, available reconstructions suggest that the\noverall proportions of exposed carbonate versus silicate terrains changed only modestly\nthrough the Phanerozoic (Bluth and Kump, 1991; Walton et al., 2023). Since currently the\nvariation of weathering lithology over the Phanerozoic remains poorly constrained, we try to\nmake as less assumptions as possible. Future work coupling this Mg isotope framework with\nlithologic reconstructions and weathering–climate models may allow a more complete\nevaluation of how continental lithology shaped coevolving Mg–Ca fluxes.\nMajor Comment 7: Please include in the Methods section in which programming language the\nmodel been constructed. Or is it just in Excel? The mean and CI results of the model are in Excel,\nbut the code should be publicly available too, so that what you’ve done can be properly examined.\nResponse: We thank the reviewer for pointing this out. The numerical model was developed and\nimplemented in R (version 4.3.1), which we used for all Monte Carlo simulations and uncertainty\npropagation. This information has been added to the Methods section.\nMajor Comment 8: There could be some greater justification behind some of the ranges used for\nthe isotope values. For example, on L246-249, what is the reason for using +/- 0.2‰ specifically?\nIs there some defined temperature effect on the isotopic fractionation as there is for elements (e.g.\nlithium isotope fractionation during its uptake when clays form is sensitive to temperature)? Also,\non L253-256, why do you not use up to the full isotope fractionation of 1.6‰, other than that other\nstudies have not used this?\nResponse: Thank you for this thoughtful comment. We have expanded the justification for the\nisotope-fractionation ranges (L74-85, 362-374, Table S1).\n(1) The fractionation factor of dolomite is based on the experimentally determined\ntemperature-dependent Δ , including the uncertainty in the temperature-dependent\ndol–SW\nequation (~ ±0.1‰): Δ = −0.1554(±0.0096)×106/T2.\ndol-SW\n(2) The fractionation factor for Mg-silicate formation is mainly based on geological\nobservations. During the low-temperature Mg-silicate formation, previous studies\nconsistently suggest a positive Δ , varying from 0.1‰ to 1.1‰ in authigenic clays\nLT-SW\n(Dunlea et al., 2017), averaging 0.8‰ to 1‰ in altered oceanic crust (Huang et al., 2018;\nSantiago Ramos et al., 2020), and 0.25‰ to 7‰ in hydrothermal fluids (Shalev et al., 2019).\nMost studies adopt 0.6–0.8‰, and we therefore set 0.8 ± 0.2‰ to encompass a large range\nin previous ranges.\n(3) ± 0.2‰ now is ± 0.1‰ for Δ but varies with temperature, corresponding to a variation\ndol–SW\nfrom ~ −1.66 ± 0.1‰ to −1.95 ± 0.1‰ over the Phanerozoic with a variation of temperature\nfrom ~9 to ~33 ℃ {Scotese, 2021 #2782}, and ± 0.1‰ is the uncertainty in the equation:\n−0.1554(±0.0096)×106/T2 (Li et al., 2015). Uncertainty for Δ dol–SW is ±0.2‰ to cover a large\nrange in previous estimates (0.5‰ to 0.83‰ in most studies). We cannot set this value to\nvary with temperature because Mg-silicate mineral is a mixture and would have no\ntemperature-dependent equation. No temperature-dependent isotope fractionation during\nMg-silicate formation was also widely set in previous studies.\n(4) We did not include the larger value of 1.6‰ because this value is unreliable. First,\nShalev et al. (2019) reported a range of 0.25–7‰ for hydrothermal fluids, but the few\nvalues >1.19‰ occur where fluid [Mg] (52.3–52.4 mmol/kg) and δ²⁶Mg (−0.95‰ ± 0.12‰\nto −0.86‰ ± 0.12‰) are indistinguishable from their assumed initial seawater\ncompositions ([Mg] = 53.3 mmol/kg; δ²⁶Mg = −0.83‰). Under such conditions, small\nuncertainties in sampling, assumed seawater chemistry, or analytical error (±0.12‰)\nproduce disproportionately large and unstable fractionation estimates, larger than the entire\nrelevant fractionation axis (as in the following Fig. 3). These values should therefore be\nconsidered very carefully. Second, the 1.6‰ value used in Shalev et al. (2019) to represent\naverage seafloor-alteration fractionation at 15 °C hinges on these outlier points and is not\nrobust once they are excluded (as in the following Fig. 3). It is possible to only remove\npartial Mg from the seawater under temperatures higher than 15 ℃ with a moderate\nfractionation factor (0.8‰ ± 0.2‰). Third, studies of authigenic clays and altered oceanic\ncrust consistently show Mg-silicate fractionations of ~0.1–1‰, which aligns with the lower\nvalues in (Shalev et al., 2019) (0.25‰ to 1.2‰) and with all other published constraints,\nbut not with these anomalously high values. Thus, it is more reasonable to adopt a\nconservative and empirically grounded fractionation range consistent with the broader\nliterature.\nFig. 3. Mg isotopes in the low-temperature hydrothermal (LTH) systems, from Figure 3 in Shalev et al. (2019).\nMinor Comment 1: L43 – Change DSM3 in the equation to ‘standard’ and then outside of the\nequation state that the standard is DSM3. Also note that this is an equation and renumber the others.\nResponse: Revised accordingly. This equation is included in the text instead of being displayed\nas an equation now. We can move this equation out if it is required from the text editor.\nMinor Comment 2: Throughout – there are some instances where ‘the’ is missing, e.g. L147-148\n‘Moreover, the silicate flux is negatively correlated with the continental fragmentation index\nduring Pangea assembly but positively…’.\nResponse: We appreciate the reviewer’s attention to language precision. In this particular sentence,\n“silicate flux” and “continental fragmentation index” refer to general variables rather than specific\ninstances; hence, the definite article the is not required according to scientific writing conventions.\nWe have nonetheless rechecked the entire manuscript to ensure consistent and grammatically\ncorrect use of articles.\nMinor Comment 3: Figure 1 – please include the geological timescale on the bottom of the figure,\nas you have done for Figure 2.\nResponse: Added as suggested.\nMinor Comment 4: Figure 2b – please include a value for the y-axis where it meets the x-axis. I\ncannot tell exactly how low the riverine Mg flux goes.\nResponse: Added accordingly.\nMinor Comment 5: L256 – I’m guessing that this should be Fsil = FHT + FLT = f x FLT ? And\nnot f x Fsil ?\nResponse: Yes, revised accordingly.\nMinor Comment 6: L257 and 258 – please add formation after ‘silicate’\nResponse: Added accordingly.\nMinor Comment 7: L267-280 – Some of this text is more of a discussion and should sit there,\nrather than in the Methods section. But also, specifically on L276, its variation is very limited in\nthe context of the modelling work you have conducted.\nResponse: We moved some of them into the discussion part.\nReferences:\nArvidson, R.S., Mackenzie, F.T., 1999. The dolomite problem; control of precipitation kinetics by\ntemperature and saturation state. Am. J. Sci. 299, 257-288.\nBluth, G.J., Kump, L., 1991. Phanerozoic paleogeology. Am. J. Sci. 291, 284-308.\nCoogan, L.A., Gillis, K.M., 2018. Low-temperature alteration of the seafloor: impacts on ocean\nchemistry. Annu. Rev. Earth Planet. Sci. 46, 21-45.\nDunlea, A.G., Murray, R.W., Santiago Ramos, D.P., Higgins, J.A., 2017. Cenozoic global cooling\nand increased seawater Mg/Ca via reduced reverse weathering. Nat. Commun. 8, 844.\nGothmann, A.M., Stolarski, J., Adkins, J.F., Higgins, J.A., 2017. A Cenozoic record of seawater\nMg isotopes in well-preserved fossil corals. Geology 45, 1039-1042.\nHiggins, J.A., Schrag, D.P., 2015. The Mg isotopic composition of Cenozoic seawater – evidence\nfor a link between Mg-clays, seawater Mg/Ca, and climate. Earth Planet. Sci. Lett. 416, 73-\n81.\nHolland, H., Zimmermann, H., 2000. The dolomite problem revisited. International Geology\nReview 42, 481-490.\nHolland, H.D., 2005. Sea level, sediments and the composition of seawater. Am. J. Sci. 305, 220-\n239.\nHorita, J., Zimmermann, H., Holland, H.D., 2002. Chemical evolution of seawater during the\nPhanerozoic: Implications from the record of marine evaporites. Geochim. Cosmochim. Acta\n66, 3733-3756.\nHu, Z., Li, W., Zhang, H., Krainer, K., Zheng, Q.-f., Xia, Z., Hu, W., Shen, S.-z., 2021. Mg isotope\nevidence for restriction events within the Paleotethys ocean around the Permian-Triassic\ntransition. Earth Planet. Sci. Lett. 556, 116704.\nHu, Z., Shi, Z., Li, G., Xia, Z., Yi, L., Liu, C., Li, W., 2022. The Cenozoic seawater conundrum:\nNew constraints from Mg isotopes in island dolostones. Earth Planet. Sci. Lett. 595, 117755.\nHuang, K.-J., Teng, F.-Z., Plank, T., Staudigel, H., Hu, Y., Bao, Z.-Y., 2018. Magnesium isotopic\ncomposition of altered oceanic crust and the global Mg cycle. Geochim. Cosmochim. Acta\n238, 357-373.\nJenkyns, C.L.B.G.M.H.H.C., 2012. Explaining the Phanerozoic Ca isotope history of seawater.\nGeology 40, 843-846.\nKim, J., Kimura, Y., Puchala, B., Yamazaki, T., Becker, U., Sun, W., 2023. Dissolution enables\ndolomite crystal growth near ambient conditions. Science 382, 915-920.\nLi, M., Wignall, P.B., Dai, X., Hu, M., Song, H., 2021. Phanerozoic variation in dolomite\nabundance linked to oceanic anoxia. Geology 49, 698-702.\nLi, W., Beard, B.L., Li, C., Xu, H., Johnson, C.M., 2015. Experimental calibration of Mg isotope\nfractionation between dolomite and aqueous solution and its geological implications.\nGeochim. Cosmochim. Acta 157, 164-181.\nLowenstein, T.K., Timofeeff, M.N., Brennan, S.T., Hardie, L.A., Demicco, R.V., 2001.\nOscillations in Phanerozoic seawater chemistry: evidence from fluid inclusions. Science 294,\n1086-1088.\nMüller, R.D., Dutkiewicz, A., Seton, M., Gaina, C., 2013. Seawater chemistry driven by\nsupercontinent assembly, breakup, and dispersal. Geology 41, 907-910.\nPeters, S.E., Husson, J.M., Czaplewski, J., 2018. Macrostrat: A Platform for Geological Data\nIntegration and Deep-Time Earth Crust Research. Geochem., Geophys., Geosyst. 19, 1393-\n1409.\nRaymo, M.E., Ruddiman, W.F., 1992. Tectonic forcing of late Cenozoic climate. Nature 359, 117-\n122.\nRidgwell, A., 2005. A Mid Mesozoic Revolution in the regulation of ocean chemistry. Mar. Geol.\n217, 339-357.\nSantiago Ramos, D.P., Coogan, L.A., Murphy, J.G., Higgins, J.A., 2020. Low-temperature oceanic\ncrust alteration and the isotopic budgets of potassium and magnesium in seawater. Earth\nPlanet. Sci. Lett. 541, 116290.\nShalev, N., Bontognali, T.R.R., Wheat, C.G., Vance, D., 2019. New isotope constraints on the Mg\noceanic budget point to cryptic modern dolomite formation. Nat. Commun. 10, 5646.\nShang, H., 2023. Dichotomous effects of oxidative metabolisms: A theoretical perspective on the\ndolomite problem. Global Planet. Change 222, 104041.\nTurchyn, A.V., DePaolo, D.J., 2019. Seawater Chemistry Through Phanerozoic Time. Annu. Rev.\nEarth Planet. Sci. 47, 197-224.\nWalton, C.R., Hao, J., Huang, F., Jenner, F.E., Williams, H., Zerkle, A.L., Lipp, A., Hazen, R.M.,\nPeters, S.E., Shorttle, O., 2023. Evolution of the crustal phosphorus reservoir. Sci. Adv. 9,\neade6923.\nWarren, J.K., 2010. Evaporites through time: Tectonic, climatic and eustatic controls in marine and\nnonmarine deposits. Earth-Sci. Rev. 98, 217-268.\nWei, G.-Y., Zhao, M., Sperling, E.A., Gaines, R.R., Kalderon-Asael, B., Shen, J., Li, C., Zhang, F.,\nLi, G., Zhou, C., Cai, C., Chen, D., Xiao, K.-Q., Jiang, L., Ling, H.-F., Planavsky, N.J., Tarhan,\nL.G., 2024. Lithium isotopic constraints on the evolution of continental clay mineral factory\nand marine oxygenation in the earliest Paleozoic Era. Sci. Adv. 10, eadk2152.\nXia, Z., Li, S., Hu, Z., Bialik, O., Chen, T., Weldeghebriel, M.F., Fan, Q., Fan, J., Wang, X., An,\nS., Zhang, F., Xu, H., Chen, J., Ji, Z., Shen, S., Lowenstein, T.K., Li, W., 2024. The evolution\nof Earth’s surficial Mg cycle over the past 2 billion years. Sci. Adv. 10, eadj5474.\nYang, J., Cawood, P.A., Yuan, X., Yuan, D., Zhou, Y., Liu, A., Liu, J., Du, Y., 2023. Enhanced\nDenudation of the Emeishan Large Igneous Province and Precipitation Forcing in the Late\nPermian. J. Geophys. Res.: Solid Earth 128, e2023JB027430."}, {"round": 2, "reviewer_comments": {"reviewer_2": "This is a manuscript that I reviewed previously. I carefully read all reviewers' comments from the previous round of review and these authors' responses and the changes they made. I feel these authors have done a thorough job in addressing all the concerns to the best they can. I am happy to recommend acceptance for publication.", "reviewer_3": "I appreciate the considerable amount of additional work that the authors have put into this paper, which I think have improved the strength of the results and the conclusions which arise from this.\nI have a couple of further, minor, suggestions for improvement:\nL58: This still needs to be made clear that it is the modern day residence time.\nL133-136: This sentence is not the easiest to follow, it may be that there are a few missing words. I suggest reading this again and doing some minor re-writing to make it clearer as to what is meant.\nL280: I would clarify that you mean the minor [Mg] cycle fluctuations which you've sort of given name to in this work, and reference your Fig. 2a panel.\nFigure 2, panels c and f: It looks like the arrows in c are offset in time by the requisite lag from the peaks in seafloor production rate in f. However, I think it would be better having the arrows (or some other symbol) in panel c being equivalent in time to those in f, so that readers looking at panel c could easily see the difference between peak seafloor production and the peaks in the Mg-silicate flux.\nRemarks on code availability:I have successfully downloaded, installed and run the code. As far as I can see, the model reproduces the main text results.\nIn the Methods section of the paper I would give an indication of roughly how long it takes for the model to run. As with most models, unless you write some code to force the command window to display information about the number of iterations the model has processed, it is not easy to know whether the model is operating as it should or has encountered some difficulties executing a particular line. This is especially the case when running any reasonably sized Monte Carlo ensemble, but is even more necessary here when the number of iterations is 1,000,000 and takes quite a long time to complete."}, "author_response": "-Reply to Reviewer #2:\nComment: This is a manuscript that I reviewed previously. I carefully read all reviewers'\ncomments from the previous round of review and these authors' responses and the changes they\nmade. I feel these authors have done a thorough job in addressing all the concerns to the best they\ncan. I am happy to recommend acceptance for publication.\nResponse: We thank Reviewer 2 for another round of the evaluation of our manuscript and the\npositive feedback.\n-Reply to Reviewer #3:\nComment:I appreciate the considerable amount of additional work that the authors have put into\nthis paper, which I think have improved the strength of the results and the conclusions which arise\nfrom this.\nI have a couple of further, minor, suggestions for improvement:\nResponse: We thank Reviewer #3 for constructive suggestions and careful assessment. Detailed\nresponses and modifications are listed below.\nMinor Comment 1:L58: This still needs to be made clear that it is the modern day residence time.\nResponse: Revised accordingly.\nMinor Comment 2: L133-136: This sentence is not the easiest to follow, it may be that there are\na few missing words. I suggest reading this again and doing some minor re-writing to make it\nclearer as to what is meant.\nResponse: We thank the reviewer for this suggestion and rewrote this sentence.\nMinor Comment 3:L280: I would clarify that you mean the minor [Mg] cycle fluctuations which\nyou've sort of given name to in this work, and reference your Fig. 2a panel.\nResponse: We agree and revised as suggested.\nMinor Comment 4:Figure 2, panels c and f: It looks like the arrows in c are offset in time by the\nrequisite lag from the peaks in seafloor production rate in f. However, I think it would be better\nhaving the arrows (or some other symbol) in panel c being equivalent in time to those in f, so that\nreaders looking at panel c could easily see the difference between peak seafloor production and\nthe peaks in the Mg-silicate flux.\nResponse: We have modified the arrows in panel c and use the length of the arrow to represent\nthe time lag between peaks of seafloor production and Mg-silicate.\nRemarks on code availability:I have successfully downloaded, installed and run the code. As far\nas I can see, the model reproduces the main text results.\nIn the Methods section of the paper I would give an indication of roughly how long it takes for the\nmodel to run. As with most models, unless you write some code to force the command window to\ndisplay information about the number of iterations the model has processed, it is not easy to know\nwhether the model is operating as it should or has encountered some difficulties executing a\nparticular line. This is especially the case when running any reasonably sized Monte Carlo\nensemble, but is even more necessary here when the number of iterations is 1,000,000 and takes\nquite a long time to complete.\nResponse: We have added a progress bar in the code."}]

Decision letter