**P/Al and Phosphorus Accumulation Data**

P/Al data was taken directly from26 (Fig. 2) and scaled as described below. Due to limitations in the sample repository at the University of Southampton from which the Heintzbjerg samples were sourced, the initiation of the LKW was not covered in the Heintzbjerg sample set and thus, 26 were unable to establish an age control point at the base of the LKW (nor were they able to establish an age control point for the significant number of samples in the Famennian). Age control points were available for the end of the LKW, base of the UKW and the Frasnian-Famennian boundary and were used to calculate average sedimentation rates and phosphorus (P) accumulation rates between each of those sections26. Although the P/Al data is informative of terrestrial phosphorus export, P accumulation rate is a more robust parameter from which to gauge phosphorus weathering by land plants and subsequent landscape stabilization (see26,49,50). In order to use P accumulation data as an input into our model however, it was critical that sedimentation rates, and thus P accumulation rates, be calculated for the entire sequence to capture not only the LKW, but also the post-UKW geochemical response.

Recent refinements in the age of the Frasnian-Famennian boundary from similar sequences16,27 make it possible to use relative dating techniques to established estimated ages for the entire sequence based on this anchor point, thus allowing calculation of sedimentation rates and P accumulation rates for the entire Heintzbjerg sequence. The field of astrochronology has long used Milankovitch cycles to gauge the passage of time in sedimentary sequences. Climatic oscillations recorded in the sedimentary record can be tied to changes in the Earth’s orbital parameters. Specifically, the ellipticity of Earth’s orbit (referred to as eccentricity, with periods of 405 kyr, 123 kyr and 95 kyr), tilt (referred to as obliquity, with a modern period of 41 kyr and varying periods in deep time51) and precession (with a modern period of 20 kyr and also varying in deep time51). Meyers52 (further refined in53) utilized a statistical optimization method called TimeOpt in the *Astrochron* package in the R platform54 in order to correlate climatic changes in the stratigraphic record with Milankovitch cycles to establish relative age of a sedimentary sequence. We employ a similar approach here, but utilize the timeOptTemplate function which was designed for more complex sedimentation models53 which we expected to encounter given the wet/arid climatic shifts noted in the Heintzbjerg sedimentary sequence. The timeOptTemplate function allows the use of proxy or lithology-specific templates which greatly enhance the spectral alignment and overall fit for complex sedimentary sequences. We evaluated the utility of multiple climate-related proxies using this method along with geochemical data from26 ultimately selecting Rb/Sr and Log Ti as proxies which produced the best spectral alignment and overall fit (e.g., 16). Eccentricity periods of 94.9, 98.9, 123.8, 130.7 and 405 kyr55 along with precession periodicities for the Devonian of 16.85 and 19.95 kyr51 were utilized. The sedimentation rates estimated using timeOptTemplate were compared with sections calculated by26 and found to be of similar magnitude (Midnatspas average sedimentation rate = 43 cm kyr-1 (see 26) compared with 71.9 cm kyr-1 (timeOptTemplate); UKW average sedimentation rate = 272 cm kyr-1 (see 26) compared with 228.6 cm kyr-1 (timeOptTemplate). Age estimates were then determined based on the Frasnian-Famennian boundary (371.870 Ma16) as an anchor point and compared with age estimates for similar published sequences16,27. Using timeOptTemplate estimations, we calculated the base of the UKW of our sequence at 372.06 Ma (compared with 372.00 – 372.02 Ma based on16,27) and the end of the LKW at 372.49 Ma (compared with 372.45 Ma based on16), a mere 40 kyr difference for each and lending confidence in the validity of the timeOptTemplate results.

Based on the relative confirmation of the timeOptTemplate results, new P accumulation rates were calculated for the entire Heintzbjerg sequence. Once again, these results were compared with those sections calculated in26 and found to be strikingly similar in overall trends, but notably with somewhat differing magnitudes (Extended Data Fig. 2). The differences in magnitude however, do not change the overall interpretation by26. Thus, these new P accumulation rates were used as input into our model as described below.

**Earth system box model**

The basic model design in this study was based on the Carbon Oxygen Phosphorus Sulfur Evolution (COPSE) model41,42, which has been extensively tested and validated against geologic records during the Phanerozoic. We incorporate refinements to this model from Ozaki and Reinhard40 and extend the biogeochemical framework for the global methane (CH4) cycle without altering the basic behavior of the model. While Ozaki and Reinhard40 include the mass exchange between the surface system (atmosphere-ocean-crust) and the mantle, we ignore this in this study. We initialized the model at 600 Ma and ran the model forward in time until the end Famennian. Specific model parameters for each experiment are outlined in the following two sections.

**Enhanced terrestrial nutrient export model**

To test the enhanced P export hypothesis, we introduced a forcing factor, *f*pP, which represents the amplification factor for the terrestrial phosphorus weathering flux:

where *J*Pw, *J*silw, *J*carbw, and *J*orgw denote the rates of phosphorus weathering, silicate weathering, carbonate weathering, and the oxidative weathering of organic matter on land, respectively, and * represents the reference value (pre-industrial flux; see40,41). To assess the uncertainty in the timing and amplitude of pulses in phosphorus weathering, we explored the following different scenarios: (i) In an attempt to reproduce the 2-3‰ positive excursions in d13C, the following forcing factors were assumed; (LKW) *f*pP = 1→2→1 over 372.72-372.54-372.36 Ma. (UKW) *f*pP = 1→2→1 over 372.1-371.91-371.73 Ma. (ii) As *i* but decreasing the forcing factor to 1.5 for both the LKW and UKW. (iii) As *i* but decreasing the forcing factor to 1.5 for the LKW only. (iv) *f*pP based on the P/Al data from Heintzbjerg, Greenland (Fig. 3, left panel), with a scaling factor of 0.5 (see below). (v) *f*pP based on P accumulation data from Heintzbjerg (Fig. 3, right panel) and applying a scaling factor of 0.15 to account for the hypothesized episodic expansion of land plants as predicted by26,46, assuming a maximum value of *f*pP = 3.35.

To convert the sedimentary data of P (P/Al or P accumulation rate) into the forcing factor of *f*pP in (iv) and (v), we first calculate the typical ‘non-event’ value by averaging data between LKW and UKW (Midnatspas Fm.), *X*avg. Then, *f*pP is calculated, as follows:

where *X* denotes the sedimentary data of P (P/Al or P accumulation rate) and *a* represents the scaling factor. When *a* = 1, *f*pP (and global P loading flux) simply reflects the variation of P flux observed at Heintzbjerg. But, in this case, our model predicts an extremely large positive excursion of d13C (10‰), inconsistent with geologic records. The scaling factor is therefore introduced to represent the relationship between the Heintzbjerg P records to global P export flux. The value of a was tuned so that the model produces a reasonable positive excursion of d13C of 2-3‰.

**Enhanced volcanic activity model**

To understand the biogeochemical consequences of the enhanced volcanic activity, we introduced a forcing factor, *f*LIP, which is multiplied by the degassing rate of CO2 via carbonate and organic carbon metamorphisms:

where *J*carbm and *J*orgm denote carbon dioxide (CO2) degassing flux via metamorphism of sedimentary carbonate (C) and organic carbon (G), respectively (* represent reference values found in41). *f*G and *f*C represent the forcing factors of degassing flux and pelagic carbonate deposition which are used in the original model. By varying the value of *f*LIP, we explore the possible impacts of CO2 degassing rate on global biogeochemistry (Fig. 4): *f*LIP = 1→2→1 over 372.72-372.54-372.36 Ma. (UKW) *f*LIP = 1→5→1 over 372.1-371.91-371.73 Ma.