Enhanced terrestrial nutrient release during the Devonian emergence and expansion of forests: Evidence from lacustrine phosphorus and geochemical records

The Devonian was a watershed moment in Earth history, with substantial changes in the biologic, ecologic, and atmospheric composition. The expansion and radiation of land plants occurred on an unprecedented scale during this period (Algeo et al., 1995; Stein et al., 2012, 2020; Giesen and Berry, 2013; Berry and Marshall, 2015; Morris et al., 2015; Lenton et al., 2016; Davies et al., 2021). Coincident with this emergent terrestrial biosphere, the Devonian hosted numerous marine extinction events (Fig. 1), including the Late Devonian mass extinction, one of the “big five” with the loss of ~40% of marine families and 60% of genera (McGhee, 1996). Additionally, the Devonian saw a pronounced decrease in atmospheric CO₂ to near contemporary levels (Lenton et al., 2016). Despite these key biological innovations and planetary transitions, much is still unknown about the specific feedbacks within the Devonian ecosphere, particularly with respect to the export of plant-mobilized nutrients, such as phosphorus (P), from the terrestrial to the marine realm. Whereas it is generally accepted that the colonization of land plants and the evolution of roots would have had a marked impact on nutrient weathering, the magnitude, timing, and duration have yet to be defined. Such insights are critical to putting together all of the pieces of the Devonian puzzle, particularly as terrestrial events such as the expansion and radiation of land plants have been implicated in the marine extinctions that occurred throughout the period (Algeo et al., 1995; Algeo and Scheckler, 1998).

Phosphorus is a critical biologically limiting nutrient on land (Walker and Syers, 1976; Vitousek and Farrington, 1997; Vitousek et al., 2010; Filippelli and Souch, 1999) and in aquatic systems (Krom et al., 1991; Correll, 1999; Blomqvist et al., 2004; Planavsky et al., 2010). Unlike other essential nutrients such as nitrogen and carbon, phosphorous has no natural stable gaseous form. Phosphorus can be transported long distances through the suspension of dust particles in the atmosphere; however, this delivery method constitutes a minute fraction of the total P flux (Chadwick et al., 1999; Filippelli, 2008). Submarine groundwater discharge can potentially deliver large amounts of P, but this varies greatly depending on the sorption capacity of aquifers as well as their redox conditions and remains an extremely difficult quantity to constrain, even in modern systems (Slomp and Van Cappellen, 2004; Santos et al., 2021). P is generally removed from groundwater through sorption to iron oxides or coprecipitation with Al, Ca, and Fe (Slomp and Van Cappellen, 2004). While this is an important source of P, the resultant nutrient delivery to aquatic systems via groundwater discharge is generally around 10% that of typical riverine flux (Santos et al., 2021). The most significant source of P to the oceans and other bodies of water is through the weathering of continental materials and ultimate delivery by way of rivers (Holland, 1978; Filippelli and Delaney, 1994; Filippelli, 2002).

Early in landscape development, P in the mineral phase is the primary source for biologic uptake (Filippelli, 2002, 2008). Because plants cannot directly access mineral-bound P, they liberate P through the acidification of root pore spaces via degradation of organic matter (Filippelli, 2002; Filippelli et al., 2010) and the release of organic exudates from roots (Jurinak et al., 1986; Filippelli, 2002; Filippelli et al., 2010). Phosphorus is lost in large amounts from the mineral phase during initial landscape development, particularly in young volcanic landscapes (see McPeek et al., 2007). This refractory form of P is used quickly, and over time sources of P shift to more labile and occluded forms (Filippelli, 2002, 2008; Filippelli et al., 2010). In well-developed soils, a significant amount of P is locked up in soil organic matter. In addition to shifting forms of P over time, total P within a system also decreases, with mineral-bound P being depleted and only organic and occluded fractions remaining (Fig. 2; Filippelli and Souch, 1999; Filippelli, 2002, 2008). It is thus expected that a landscape with little to no vegetation would transition through these phases stepwise, beginning with predominantly mineral sources of P and shifting to more occluded sources. Likewise, it would be expected that the overall release of P in a landscape would peak very quickly upon plant colonization and stabilize over time, as has been noted (Filippelli and Souch, 1999; McPeek et al., 2007; Filippelli et al., 2010).

Although first appearing in the Middle Ordovician (Wellman et al., 2003; but see also Strother and Foster, 2021), it was not until the Middle Devonian that land plants (embryophytes) were thought to be widespread and occupied continental interiors. Prior to the Middle Devonian, land plants were confined to areas immediately adjacent to bodies of water (or in wetland environments), had limited root systems, and had yet to develop the necessary biological innovations that would lead to arborescence—the development of the first trees (Algeo et al., 1995, 2001; Algeo and Scheckler, 1998, 2010; Stein et al., 2007). Continental interiors would have been bereft of significant vegetation and well-developed (and weathered) soil profiles. From a P-weathering perspective, most continental interiors would likely have been very “young,” meaning a significant amount of P was likely available in mineral form. By the Middle–Late Devonian, land plants began to diversify and spread into continental interiors and uplands, populating the once barren landscape with vegetation (Algeo and Scheckler, 1998; Algeo et al., 2001; Piombino, 2016). The diversification and expansion of plants led to the first appearance of trees and, by the Tournaisian, the widespread propensity of plants to produce seeds (Algeo et al., 1995, 2001; Algeo and Scheckler, 1998). The earliest of these trees appeared in the Eifelian and Givetian (Giesen and Berry, 2013), the most notable of which was Archaeopteris, a large progymnosperm that grew up to 30 m in height (Algeo et al., 1995). Archaeopteris first appeared in the earliest Givetian, achieving hegemony by the Frasnian (Meyer-Berthaud et al., 1999; Meyer-Berthaud et al., 2010; Stein et al., 2020). While not seed-bearing, Archaeopteris is the first known plant to display characteristics common in modern seed plants, possessing leaves and generating root systems that were simultaneously deep, laterally extensive, and complex (Algeo et al., 2001; Stein et al., 2020). Although possibly only occupying niche biomes on inception (Stein et al., 2020), such root systems likely would have had enormous impacts on geochemical cycling in those regions.

The development of significant and complex root systems in the Devonian is perhaps one of the single most critical events in the evolution of life in the Paleozoic and would have cascading impacts on Earth’s biosphere. The earliest evidence of roots associated with vascular land plants appeared in the Early Devonian (late Pragian) and were shallow and simplistic in comparison to modern analogues (Algeo et al., 2001). Even by the Middle Devonian, with the exception of archaeopteridaleans, most major plant groups (aneurophytaleans, cladoxylopsids, and lycopsids) had relatively simple root structures with minimal root penetration that would have had limited soil formation capabilities (Stein et al., 2012, 2020). More substantial root systems did not appear until the Middle to Late Devonian and were notably associated with Archaeopteris (Algeo et al., 2001; Kenrick and Strullu-Derrien, 2014; Stein et al., 2020). Significant root development led to the creation of modern soil weathering processes, nutrient cycling, and extensive soil development (Algeo and Scheckler, 1998; Kenrick and Strullu-Derrien, 2014; Morris et al., 2015). Unlike its contemporaries such as Wattieza (a cladoxylopsid whose life cycle terminated shortly after reproducing), Archaeopteris were likely long-lived (based on ring counts, Stein et al., 2007) and thus established significant soil profiles (Algeo and Scheckler, 1998) with a measurable impact on local biogeochemical cycling and nutrient loads into the oceans. Indeed, many studies have shown that roots dramatically alter soil formation and thus the soil biogeochemistry (e.g., Filippelli and Souch, 1999; Filippelli, 2008). The development of extensive root systems, stress-resistant seeds, and an overall increase in size subsequently contributed to the expansion of land plants into continental interiors and uplands (Algeo and Scheckler, 1998, 2010). The vast increase in land plants has been implicated in the drawdown of atmospheric CO₂ through photosynthesis and elevated rates of silicate rock weathering through increased soil residence times (Algeo et al., 1995; Algeo and Scheckler, 1998; Falkowski et al., 2004; Quirk et al., 2015; Kenrick and Strullu-Derrien, 2014; Berry and Marshall, 2015; Morris et al., 2015; Lenton et al., 2016). Furthermore, this expansion likely had a marked impact on the weathering and export of P-bearing minerals to the oceans. Some have suggested that this export happened so quickly and on such a global scale that it drove immense algal blooms in the shallow Devonian seas, which resulted in the creation of vast anoxic zones that ultimately drove the multiple biotic crises (Algeo et al., 1995, 2001; Algeo and Scheckler, 1998, 2010; Falkowski et al., 2011).

Perturbations in the P cycle have been implicated as a driver of some of the Devonian marine crises, but there is no direct evidence. Nearly all of the extinction events in the Devonian were characterized by widespread anoxia, with the exception of the Taghanic Crisis, which had differing regional impacts and varying levels of hypoxia/anoxia (successive, linked fluctuations in climate and sea level are noted as potential causal factors, though a collapse in terrestrial vegetation occurred toward the end of the event; Marshall et al., 2011; House, 1985, 2002; McGhee, 1996; Racki, 2005; Aboussalam and Becker, 2011; Marshall et al., 2011; Becker et al., 2016, and references therein). For these vast ocean anoxic zones to form, it is reasonable to assume an external eutrophication pulse (Meyer and Kump, 2008).

Although P is readily mobilized in young landscapes, its short residence time in the modern ocean (20–50 k.y.; Ruttenberg, 2003) makes it unclear whether a terrestrial nutrient pulse (or pulses) on a global or basin-wide scale is sufficient to drive significant changes in the marine system. Attempts have been made to provide such evidence utilizing nutrient signals from marine records; however, this approach is problematic due to profound changes in the marine P cycle caused by ocean stratification and expansive anoxic zones in shallow seas characteristic of this period (Racki, 2005; Becker et al., 2016). Consequently, many questions remain about the mechanisms of enhanced P input to ancient oceans.

Here, we present a novel approach to constraining the timing and source of nutrient pulses in the Devonian associated with terrestrialization. Through the analysis of Devonian lacustrine sequences, versus marine sediments or paleosols as has been attempted previously, we can gain proximal insights into a higher temporal resolution of landscape changes and nutrient dynamics within the Devonian during the interval when plant evolution and expansion was well underway. Lacustrine sequences have long been used to study nutrient dynamics and can record vital information about weathering, nutrient flux, and a host of other important parameters. Although depositional environments vary greatly, lake environments are generally low-energy in nature and thus can provide a high-resolution sediment record with fine-scale laminations that can be narrowed down to specific decadal, and even sometimes annual and seasonal, variations. Constraining land plant-associated nutrient export will provide crucial evidence to link the timing and pace of land plant expansion with the marine biotic crises in the Devonian.

The target for this study was lacustrine and near-lacustrine sequences proximal to critical events and time transects within the Middle to Late Devonian (Eifelian, Givetian, and Frasnian). Where possible, outcrops spanning a biotic crisis with associated marine anoxia were chosen. These include the Choteč and Kačák events in the Eifelian, the Taghanic Crisis in the Givetian, and the Kellwasser Event in the Frasnian. All samples are archived at the University of Southampton, UK. Five distinct sites from the paleocontinent of Euramerica were chosen, including three from the Orcadian Basin (Scotland and Orkney) and two from the Devonian Basin in East Greenland. Study sites include Geanies in Easter Ross, Scotland; Hoxa Head on South Ronaldsay, Orkney; Hegglie Ber on Sanday, Orkney; Ella Ø (island), and Heintzbjerg in East Greenland (Fig. 1). All sites contain continuous lacustrine sequences except for the Heintzbjerg samples, which are primarily fluvial. In most cases, samples were taken at various points throughout a given lake cycle. For the Geanies set, samples were selected from deep lake facies (the rest of the lake cycle was not sampled; however, the sequence is an uninterrupted succession of lake cycles).

Throughout much of the Devonian, the Orcadian Basin was the site of a series of lacustrine systems formed in extensional half graben basins with distinct depocenters (McClay et al., 1986; Marshall and Hewett, 2003). The depths and areal extents of these lakes varied based on climatic conditions but they were generally shallow playa lakes with periods of deep permanent lakes (Marshall et al., 2007). During intervals of enhanced precipitation, these lakes occasionally merged to form extensive lacustrine systems, such as the one that formed the distinctive Achanarras fish bed referred to as a mega-lake (Marshall et al., 2007). The Orcadian Lake (also referred to as the Achanarras Lake) is a particularly long-lived example, which at one time occupied much of the eastern Highlands of Scotland and the Moray Firth and extended as far north as the Shetland Islands and east to the margins of Norway (Trewin and Thirlwall, 2002; Marshall et al., 2007). At its largest extent during the Eifelian and early Givetian, the Orcadian Lake would have been ~27,000 km² in size, or roughly the size of Belgium (Andrews and Hartley, 2015). The northern Scotland and Orkney samples (Geanies, Hoxa Head, and Hegglie Ber) were all collected from Orcadian Basin lacustrine laminates and contain mostly mudstone-limestones, siltstones, and sandstones (Donovan et al., 1974) with fluvial sandstones and conglomerates near the margins of the Orcadian Basin along with significant features such as the Achanarras fish bed (Andrews and Hartley, 2015; Johnstone and Mykura, 1989; Trewin and Thirlwall, 2002).

The paleolatitude was between ~20°S and 30°S, in an arid to semi-arid zone (Torsvik and Cocks, 2004). Although the area experienced many phases of transgression and regression (most likely due to climatic forcing), in general the depositional environment is characterized as low-energy, with fine-grained facies providing the potential for high-resolution data collection and analyses (Andrews and Hartley, 2015). The Geanies samples represent a series of cyclic lacustrine deposits near Tain, Scotland, in Easter Ross, of which 12 distinct lake cycles were sampled. The oldest samples in the Geanies sequence lie in the Balintore Formation (the first nine lake cycles), continue into Jessie Port (a fish bed equivalent to the distinctive Achanarras fish beds in Caithness correlative with the onset of the Kačák Event and spanning two lake cycles), and terminate in the latest Eifelian in the Geanies Formation (the final lake cycle of the sequence) (Marshall et al., 2007, 2011; Fig. 3). The Hoxa Head samples represent a similar cyclic lacustrine sequence with 21 lake cycles collected from an outcrop on northwestern South Ronaldsay, Orkney. They are sourced from the Rousay Flagstone Formation (the uppermost section of the Upper Stromness Flagstone Formation; Fig. 3) and are early Givetian (Speed, 1999; Trewin and Thirlwall, 2002), occurring above the Kačák Event (Fig. 1). The Hegglie Ber sequence originates from southwestern Sanday, Orkney (Fig. 1). This continuous but thin sequence also contains cyclic lacustrine deposits with at least three distinct lake cycles. It was collected from the Eday Flagstone Formation (Fig. 3) and is mid-Givetian in age (Marshall, 1996; Marshall et al., 2011).

The lacustrine sequences in Ella Ø formed under similar conditions as lakes that were comparable in size in the Orcadian Basin (McClay et al., 1986; Marshall and Stephenson, 1997; Trewin and Thirlwall, 2002). Devonian deposits comprise the entire southeastern corner of the ~145 km² island and represent the fringes of lakes preserved amongst the topography of an unconformity surface (Marshall and Stephenson, 1997, and references contained therein). Samples are from the Solstrand Formation and comprised primarily of eastward-draining braid-plain systems with a lacustrine bed present at three separate localities (Marshall and Stephenson, 1997; Larsen et al., 2008). The sample set presented here represents a single but complete lake cycle (Marshall and Stephenson, 1997) of Givetian age and occurs between the end of the Kačák Event and the beginning of the Taghanic Crisis (Fig. 1). During monsoon periods, it is possible that the outflow from this lake linked to that from the Orcadian Lake and may be correlative with the third lake cycle within the Eday Group (Marshall and Hemsley, 2003).

The Heintzbjerg samples represent a continuous 1100 m sequence on nearby Ymer Ø from the Sofia Sund, Midnatspas, Zoologdalen, and Andersson land formations consisting of fluvial deposits and abandoned channel lakes. They range in age from the Frasnian to early Famennian and encompass a portion of the Lower Kellwasser Event (LKW) as well as the entirety of the Upper Kellwasser Event (UKW) (Larsen et al., 2008; Fig. 1). This Devonian basin stage is marked by a shift in drainage direction to the south and is comprised of mostly fluvial systems with rare aeolian intercalations (Larsen et al., 2008).

Age constraints for these samples were obtained via precise palynological calibrations and in some cases in combination with fossil fish assemblages (see Marshall and Stephenson, 1997; Marshall and Hemsley, 2003). Precision dating of the Kellwasser events by Percival et al. (2018) enabled the employment of age control points for the Heintzbjerg samples, specifically bounding the end of the LKW, the beginning of the UKW, and the Frasnian/Famennian boundary and thus established a 525,000-year span between the oldest and youngest of these events. This allowed the determination of both sedimentation rates and P accumulation rates for the Heintzbjerg sequence using bulk density data (with an average value of 2.33 g/cm³) estimated from published U.S. Geological Survey data (Manger, 1963) for similar Devonian sequences. For the remaining sample sets, relative age control points were estimated by counting lake cycles present in the sequence, with one lake cycle defined as beginning at the base of one organic matter-rich layer and ending before the appearance of the next organic matter-rich layer. Lake cycles in the Orcadian Basin and the east Greenland Devonian Basin recur with Milankovitch cyclicity, specifically with the eccentricity parameter approximately every 100,000 years (Astin, 1990; Kelly, 1992; Olsen and Larson, 1993; Olsen, 1994; Marshall, 1996; Andrews and Trewin, 2010; Andrews et al., 2016). This allowed calculation of individual sedimentation rates and P accumulation rates for each lake cycle in each sample set.

Samples were collected from whole rocks at the archives of the University of Southampton. Between 5 g and 10 g of sample was removed with care taken to avoid collecting weathered surfaces, washed with Milli-Q water, dried, and then powdered by hand using a porcelain mortar and pestle.

Approximately 100 mg of each homogenized sample was acidified with 6N HCl and then digested in a microwave system for 30 min using a variation of Environmental Protection Agency method 3052. Samples were then spun in a centrifuge at 10,000 rotations per minute for 10 min, and the supernatant was extracted. In preparation for inductively coupled plasma–optical emission spectrometry (ICP-OES) analysis, samples were diluted in accordance with expected concentrations of the analyte. In general, a dilution ratio of 1:10 was used for P and Ti, and 1:50 for Al, Ca, and Fe. ICP-OES analysis was performed on a PerkinElmer Optima 7000 DV using three replicates per sample, with analytical precision within 5%.

Total digestion was chosen over a sequential extraction for three reasons. The first was the uncertainty about the diagenetic impacts on the various P fractions in nearly 400-m.y.-old rocks. Second, while some of the smaller lakes analyzed may have had intervals in which they were likely oligotrophic and oxic (based on the presence of fossilized fish scales; Marshall and Stephenson, 1997), lakes within the Orcadian Basin had a varied redox history with evidence of intervals of anoxia (Trewin and Thirlwall, 2002; Andrews and Hartley, 2015). This is important, as redox cycling in anoxic environments could provide misleading P geochemical results (as addressed in detail below). Third, total element ratios enable the comparison of total P to detrital elements such as Al, Ca, and Ti as well as redox cycling via total Fe. Comparing total P to detrital inputs provides an estimate of P enrichment (either from authigenic accumulation of P or from elevated levels of P within the detrital signal as would be expected during soil development in a young landscape) while eliminating the bias introduced by sediment focusing, turbidity currents, in situ carbonate and aluminum hydroxide production, etc. Titanium serves as a reliable proxy for terrigenous input, with elevated P/Ti values suggesting an increase in P content of the detrital fraction. Elevated P/Ti on its own, however, is not definitive due to complex geochemical cycling within the water column as well as potential organic acid-facilitated Ti mobility in immature sediment (see Pe-Piper et al., 2011). The P signal must be decoupled from authigenic P; thus, comparing P to Al, Ca, and Fe is necessary. Ca is used as an indicator of in situ carbonate production and Al as an indicator of both terrigenous input and potential photochemical production of Al hydroxides. Using a suite of elements provides a more accurate portrayal of detrital inputs as Al is prone to some alteration during settling primarily due to biogenic scavenging, but this is most common in deep marine settings (Murray and Leinen, 1996; Pattan and Shane, 1999). Comparing total P to total Fe is useful in addressing concerns stemming from the redox condition during deposition. Under oxic conditions, orthophosphates within the water column are adsorbed onto Fe-oxyhydroxides that become stored within the sediment and may go on to form authigenic phosphate minerals (Filippelli, 2002). Under anoxic conditions, P bound to Fe-oxyhydroxides readily dissolves, releasing P into the water column to actively participate in the continued biologic cycling of P (Filippelli, 2002). The fate of both P and Fe are heavily intertwined; thus, the P/Fe ratio is useful in separating the P signal from reductive dissolution of these Fe-oxyhydroxides. An elevated P/Fe ratio indicates an increase in terrestrial delivery of P that is independent of redox cycling. While useful, as with any geochemical proxy, P/Fe on its own cannot be used to accurately determine redox conditions (this is addressed in section 2.4.3). To complement the analyses and assist with interpretations, correlations between P and various other analytes were performed and statistical significance was calculated using GraphPad Prism (only significant results are presented in section 4, but all statistical analyses are included in Supplemental Table S1¹). Additionally, as both Al and Ti are prone to some alteration, correlations between Al and Ti were calculated for each site. Although no geochemical proxy can be universally reliable, significant correlation between Al and Ti serves as an indicator of their utility as terrestrial input proxies in this study (Supplemental Table S1).

Total C and total inorganic C (TIC) was measured using an ELTRA CS-580 Carbon/Sulfur Analyzer with an attached TIC module. The amount of total organic carbon (TOC, also referred to as C_org) was calculated as the difference between total C and TIC. The C_org values were used in combination with P_tot determinations to determine the C_org:P_tot ratio, one of multiple paleoredox proxies used. Analytical precision was within 5% based on replicate analyses of standard reference materials.

Powdered samples were analyzed for Al, Cu, Cr, Mo, Ni, Pb, Rb, Si, Sr, Ti, U, and Zr using an Olympus X-ray fluorescence (XRF) instrument. These elemental distributions were used to calculate ratios of Rb/Sr, Sr/Cu, Si/Al, Ti/Al, and Zr/Al and metal enrichment factors (EFs) for Cu_EF, Cr_EF, Mo_EF, Ni_EF, Pb_EF, and U_EF.

In lake sediments, the Rb/Sr ratio has been used as both a chemical weathering indicator as well as a proxy for mass accumulation rate from terrestrial input (see Xu et al., 2010). Elevated values of Rb relative to Sr are indicative of material that has been weathered with respect to its source rocks and are interpreted here as elevated terrestrial input. Si/Al ratios are similarly used as a weathering proxy, but unlike Rb/Sr they are generally indicative of physical weathering intensity. As quartz is generally resistant to chemical weathering and is coarser than clay minerals, higher values of Si/Al will reflect a greater degree of physical weathering (Calvert and Pedersen, 2007). Similarly, Ti/Al and Zr/Al ratios have been used as proxies for grain size, specifically silt-to-clay ratio as both are transported in the silt/fine sand portion (see Calvert and Pederson, 2007; Wei and Algeo, 2020). Sr/Cu ratios have been used in both paleolacustrine and paleomarine settings as paleoclimate indicators (e.g., Vosoughi Moradi et al., 2016; Tang et al., 2020; Li et al., 2020; Pan et al., 2020). Strontium is readily released from source rocks via chemical weathering, and under hot and arid conditions it tends to concentrate in lake waters (Liang et al., 2014). Copper will precipitate as Cu₂S in organic-rich environments and is generally found in larger quantities in biologically productive aquatic environments associated with warm and humid environments (Liang et al., 2014). Thus, low ratios of Sr/Cu are suggestive of warm and humid conditions (values between 1 and 10 are generally indicative of these conditions, although some authors have suggested values of 1–5), and higher Sr/Cu values indicate hot and arid conditions (generally for Sr/Cu values >10; again, some authors argue that values >5 should be used) (Liang et al., 2014; Cao et al., 2015; Jinhua et al., 2018; Li et al., 2020). Sr/Cu tends to vary inversely with Rb/Sr such that a high Sr/Cu ratio accompanied by a low Rb/Sr ratio indicates hot and arid conditions with low weathering and vice versa (Cao et al., 2015; Shen et al., 2020).

Given that geochemical proxies can be highly variable between and within formations (Algeo and Liu, 2020), a suite of seven proxies was selected in combination with P/Fe to provide a robust approach to characterizing paleoredox. Multiple proxies allow for comparison among lakes in the various sequences from distinct geographic regions and potentially different elemental inventories. Many of these proxies were designed for marine environments, and their utility for paleolacustrine environments is unclear. Also, because TOC for the different sample sites varied (and sometimes significantly within a given sequence), it was necessary to have proxies effective in both anoxic (organic-rich) and oxic-suboxic (organic-poor) environments as both could be encountered in the Devonian sequences.

Most samples (but not all) were expected to be relatively organic-poor; thus, Pb_EF was chosen as these values generally display a negative covariation with TOC (see Algeo and Liu, 2020). C_org:P_tot, Cu_EF, Cr_EF, Mo_EF, Ni_EF, and U_EF were chosen due to their varying effectiveness under oxic and anoxic conditions (i.e., Cr and U are mobilized by oxic weathering and sequestered into sediments in anoxic environments, and Mo enrichments are highest under euxinic conditions, etc.). We chose C_org:P_tot (total P versus reactive P) for ease of comparison with other paleomarine and paleolacustrine studies. This serves as a proxy for C_organic:P_reactive and is generally effective in areas where the detrital P fraction is minor (Algeo and Ingall, 2007), though we acknowledge that the detrital fraction may be significantly larger in some of our sequences. The remaining proxies used (Cu_EF, Cr_EF, Mo_EF, Ni_EF, Pb_EF, and U_EF) represent enrichment factors for the various elements. As in the study of Liu and Algeo (2020), the detrital fraction of each element was estimated. The enrichment factor (relative to the detrital fraction) was calculated using the following formula:

where X is the element and (X/Al)_Detrital are the average concentrations of the element and Al of upper continental crust as reported by McLennan (2001). Average upper continental crust values used are 8.04 wt% for Al, 25 ppm for Cu, 83 ppm for Cr, 1.5 ppm for Mo, 44 ppm for Ni, 17 ppm for Pb, and 2.8 ppm for U (McLennan, 2001). All seven paleoredox proxies were employed for each separate sequence and reported here. To determine the effectiveness of each proxy at each locality, correlations among each of the seven proxies to each other and TOC were calculated, and only those whose correlations were statistically significant were considered viable for redox determinations.

Results are presented here by section from oldest to youngest. Each section contains a figure comparing the P input relative to metal ions (P/Al, P/Ca, P/Fe, and P/Ti; all of which are P to metal ion ratios and henceforth abbreviated as P/Me), P accumulation rate, the weathering proxy Rb/Sr, the paleoclimate proxy Sr/Cu, elemental enrichment factors, C_org and C_org:P_tot.

Results for the Eifelian samples at the Geanies locality are shown in Figure 4. At the base of the section, P/Ti is elevated in lakes 9 and 10 and gradually decreases to less than 1 within lakes 12–15 (Fig. 4A). A gradual increase in P/Ti begins again in lake 15, reaching values of ~2 in the upper half of lake 16. P/Fe generally tracks P/Ti in lakes 9 and 10 but remains at background values like the other P/Me ratios. Significant perturbations are observed in P/Al in lakes 14, 17, and 20, aligning with maxima in Rb/Sr (Fig. 4C). Phosphorus accumulation rate is relatively low at the base of the sequence and increases to greater than 300 µmol P cm^–2 kyr^–1 for lakes 11–14, peaks during lake 15, and peaks again in the upper half of lake 16 into lake 17. Lakes 10, 12, 16, and 19 have the highest concentrations of C_org (>3 wt%; Fig. 4E). The Rb/Sr values remain relatively low for most of the sequence but notably increase in the upper half of lake 14, decrease to minimum values in lake 16, and return to relatively higher values in lakes 17, 18, and 20. The Sr/Cu values are generally antithetical to the Rb/Sr values, with significant maxima in lakes 9, 14, 16, and 19 (Fig. 4C). Paleoredox proxies exhibit high enrichment factors in lakes 9, 10, 12, 16, and 19 (Fig. 4D). Notably, Cu_EF does not trend like the other EFs in lakes 9 and 10 but is in relative agreement with the other enrichment factors in the remainder of the sequence. Additionally, Cr_EF does not appear to show elevated values in lakes 12 and 19. Peaks in C_org:P_tot are observed during lakes 10, 16, and 19.

Results from the Givetian sequence at Hoxa Head, South Ronaldsay, Orkney, are shown in Figure 5. P/Me ratios show some variation, but magnitudes are low, except for P/Ti, which varies between 0.4 and 1.8 (Fig. 5A) throughout the record. P accumulation rates range from 13 µmol P cm^–2 kyr^–1 to near 100 µmol P cm^–2 kyr^–1 (Fig. 5B). Maxima in P/Ca align with Rb/Sr maxima at 16 m, 53 m, 65 m, 96 m, 102 m and between 120 m and 140 m. Concurrent maxima are observed in all P/Me ratios at 14 m, 63 m, 132 m, and 145 m. C_org ranges between 0 wt% and 7 wt%, with values above 5 wt% recorded in five samples (Fig. 5E). The maxima and minima in C_org coincide with those in Sr/Cu (which are opposite in magnitude to variations in Rb/Sr; (Fig. 5C)—trends that continue throughout the sequence. Trace element enrichments are observed at 14 m, 60 m, 68 m, 116 m, 130 m, 146 m, and 176 m (Fig. 5D). C_org:P_tot is in general agreement with maxima observed in EFs, with the highest values at 14 m, 22 m, 36 m, 60 m, 68 m, 85 m, 112–134 m, and 146 m (Fig. 5E).

Results from the Givetian sequence at Ella Ø, Greenland, are shown in Figure 6. P/Ti values are relatively high at the bottom of the sequence and at 0.85 m and 1.1 m (Fig. 6A). Phosphorus accumulation rates (Fig. 6B) and the remaining P/Me ratios have patterns similar to those of P/Ti. The exception is P/Al at 0.7 m and 1 m, which appears to trend out of phase with the other P/Me ratios in those samples. Notably, fossilized plant remains were discovered at various points in this sequence, including in proximity to the two latter P perturbations at 0.85 m and 1.05 m. C_org varies from below 1 wt% to ~2 wt%, with a pattern that is roughly similar to that of C_org:P_tot, with the exception being between 0.85 m and 1.0 m (Fig. 6E). Sr/Cu values are low at the bottom of the sequence and trend toward higher values at the top of the sequence, which is opposite to the patterns observed in Rb/Sr (Fig. 6C). The paleoredox proxies exhibit maxima at 0.85 m and 1.35 m for most EFs (Fig. 6D). Mo_EF records a peak at 0.85 m that is not substantially expressed in the other EFs (particularly not Cu_EF). The C_org:P_tot ratios are generally opposite the patterns exhibited in the EFs (Fig. 6E). The peak in the C_org:P_tot values occurs at 1.6 m, notably prior to the P/Ti peak, and is not in agreement with other paleoredox proxies.

Results from the Givetian sequence at Hegglie Ber, Sanday, Orkney, are shown in Figure 7. One sample at 10 m has relatively high P/Me ratios (Fig. 7A), P accumulation rates (Fig. 7B), and C_org concentrations (Fig. 7E). Notably, the P/Ti peak value is the largest measured in all of the Devonian sequences presented here. The Rb/Sr ratio and Sr/Cu are again broadly opposite one another (Fig. 7C). With the exception of Cr_EF, paleoredox EFs exhibit maxima at 10 m (Fig. 7D). The general trend for C_org:P_tot aligns with the EFs with a maximum value of 700 at 10 m (Fig. 7E).

Results from the Frasnian/Famennian sequence at Heintzbjerg, Greenland, are shown in Figure 8. At the base of the record during the LKW (lower blue shaded region), pronounced increases in P/Fe, P/Ca, P/Al (Fig. 8A), and Rb/Sr (Fig. 8C) ratios are observed at 15 m and 51 m. P/Ti values are elevated when compared to values between the end of the LKW and beginning of the UKW (Fig. 8A). Sr/Cu drops rapidly from maxima near 15 at the bottom of the sequence to less than 5 by the end of the LKW and then fluctuates between 2 and 8 for the remainder of the LKW and the intermission between the LKW and UKW (Fig. 8C). Although generally low (<0.65 wt%), C_org (Fig. 8E) maxima within the LKW coincide with those in the P/Me ratios. All P/Me ratios decrease to background levels between the end of the LKW and the beginning of the UKW. The Rb/Sr ratio (Fig. 8C) is elevated prior to the start of the UKW, with increases in both P accumulation (Fig. 8B) and C_org observed at UKW boundary. The most significant P/Ti perturbation occurs during the latter half of the UKW beginning around 600 m and is reflected across all P/Me ratios and P accumulation, after which all values (with the exception of a pronounced increase in P/Ca at the Frasnian/Famennian boundary) return to background levels in the transition from the latest Frasnian into the Famennian. Of note, P accumulation markedly increases at the onset of the UKW and peaks near 8000 µmol P cm^–2 kyr^–1 at 680 m, which is the highest rate observed across all sequences presented here. Sustained values above 1500 µmol P cm^–2 kyr^–1 are observed for the entire portion of the UKW in which P accumulation rates were calculated (P accumulation rates were not calculated for the LKW, the latest Frasnian, and the Famennian due to the lack of age control points for those portions of the sequence). C_org remains low and variable for the rest of the sequence after peaking near 1 wt% at the beginning of the UKW. Rb/Sr values fluctuate throughout most of the UKW, and the largest values within the UKW align with the beginning of the positive trend in P/Me and P accumulation just below 600 m before dropping in the latest Frasnian and rising sharply in the Famennian around 1000 m. Sr/Cu similarly fluctuates during the UKW, with minima also coincident with maxima in both P/Me and P accumulation and maxima of 12.6 at the bottom of the Zoologdalen Formation. The paleoredox proxies show significant maxima in EFs (Fig. 8D) during the LKW (10 m) and again during the UKW (720 m). Prominent enrichments are observed in Pb_EF at 240 m and 620 m. The peak at 240 m corresponds with modest increases in the other EFs, whereas the peak at 620 m does not correlate with other EFs. The C_org:P_tot data (Fig. 8E) are highly variable, but peaks are present during the LKW, the beginning of the UKW at 350 m, and at the end of the UKW beginning at 700 m and extending past 800 m into the Famennian.

Overall, the data from these sequences provide insight into the geochemical cycling of Devonian lacustrine environments. In general, increases in P concentrations in a lacustrine stratigraphic sequence, particularly if there is a consistent P/Me increase, indicate an increase in net P input to the system (e.g., Filippelli, 2002). Sediment concentrations of a given element can also reflect dilution or enrichment in a mixed input system, such as one where terrigenous allochthonous input and autochthonous in situ production both impact net sedimentation. This is typically constrained by normalizing to different elements—in the case of this study, normalizing to Ti (and Al) to constrain for terrigenous input and Ca to constrain for in situ carbonate production. Uniformly, our records show a mostly stable background P signal with brief intervals of higher P concentration and higher P/Me ratios. Combining P/Me ratios with redox data adds additional detail regarding the genesis of P enrichment within the lakes. Paleoredox determinations remain one of the largest challenges in studies of ancient marine and particularly ancient lacustrine systems. Local lithological variations, combined with changes in sedimentation rates as well as diagenetic impacts, may render certain proxies ineffective (e.g., Algeo and Liu, 2020); however, the general agreement of the majority of proxies employed here provides a reasonable argument that their interpretation is the most likely scenario. Overall, the striking agreement between the six EFs utilized and, to a lesser extent, agreement with C_org:P_tot at the different sample sites, is persuasive evidence that they are sentinels of actual redox conditions present in these Devonian lacustrine systems. Combining redox determinations with paleoclimate and weathering proxies was particularly effective for most sequences and provided additional context for the interpretation of both redox conditions and nutrient input into the lakes.

The weathering proxy Rb/Sr and paleoclimate proxy Sr/Cu generally display an antithetic relationship in each of the study sites presented here. This is important in establishing a relational link between weathering and climate in both the Orcadian Basin and the Devonian Basin in East Greenland. During wetter periods (low Sr/Cu values), higher amounts of precipitation and runoff would result in enhanced weathering (elevated Rb/Sr values). The reverse would be expected during more arid periods (elevated Sr/Cu values) with reduced precipitation and runoff resulting in a reduction in weathering rates (low Rb/Sr values). While there is sample-to-sample variation within the data, the inverse relationship between Rb/Sr and Sr/Cu is maintained almost universally at each of the five sites.

In the higher resolution records such as those of Geanies (Fig. 4), Hoxa Head (Fig. 5), and Heintzbjerg (Fig. 8), cyclic trends can be observed where wetter periods lead into arid periods before returning to wetter periods. These trends are consistent with what is known about regional climate in both the Orcadian Basin and the Devonian Basin in East Greenland. As discussed earlier, lake cycles within East Greenland and the Orcadian Basin are linked to Milankovitch cyclicity, with Milankovitch orbital parameters the genesis for climatic variation driving the cyclicity (Astin, 1990; Kelly, 1992; Olsen, 1994; Marshall 1996; Andrews and Trewin, 2010; Andrews et al., 2016). Deep and semi-permanent lakes formed during wetter periods, followed by progressive shallowing into playa lake environments as aridity returned (Marshall and Stephenson, 1997). This ~100 k.y. cycle is most easily identifiable in the Geanies sequence (with 12 distinct lake cycles) and particularly so in the upper half of the sequence. While the lake cycles at the Hoxa Head sequence are more difficult to separate, the samples presented here represent 21 separate lake cycles, and most can be identified by the alternation of Sr/Cu from wetter to arid. Even among the smaller data sets such as Ella Ø (Fig. 6) and Hegglie Ber (Fig. 7), this cyclicity can be observed. Ella Ø represents one lake cycle (thus ~100 k.y.), and a clear trend of increasing aridity can be observed with Sr/Cu values near 5 at the base of the record, which rise to 20 in the upper portion of the record. Hegglie Ber represents three distinct lake cycles (~300 k.y.), and while it is difficult to identify three lake cycles solely based on Sr/Cu ratios, there is a clear change from arid at the base of the record, to wet by 10 m, and back to arid by 25 m. Thus, the general agreement between the Sr/Cu data reported here and known regional climatic cyclicity, combined with the expected antithetic behavior of Sr/Cu and Rb/Sr, suggest that these data are generally reflective of the climate and weathering conditions present at each site.

The general agreement among the six distinct enrichment factors with both C_org and C_org:P_tot lend confidence that these proxies are effectively capturing redox conditions within the lakes at each of these sites. In the five records presented here, all exhibit distinct periods of anoxia, with the one exception being the Heintzbjerg record. In the Geanies sequence, significant periods of anoxia occur in lakes 9, 10, 16, and 19 (Figs. 4D and 4E). At Hoxa Head, five distinct periods of anoxia can be identified that are nearly in line with the largest increases in C_org at 12 m, 62 m, 68 m, 82 m, and 146 m (Figs. 5D and 5E). At Ella Ø, the presence of amorphous organic matter, and thus anoxia, throughout the sequence noted by Marshall and Stephenson (1997) is supported at various intervals by the EFs, specifically at 0.55 m, 0.8 m, 1.05 m, and 1.4 m (Figs. 6D and 6E). In this case, C_org:P_tot values seem out of phase with the EFs, which complicates redox determinations. However, elevated values of C_org near or above 2 wt% are present in much of the sequence, which suggests that this lake was likely anoxic/suboxic from ~0.55 m to 1.5 m. At Hegglie Ber, there is a clearly defined anoxic period at ~10 m, with the rest of the sequence reflecting conditions that are likely oxic (Figs. 7D and 7E). At Heintzbjerg, the organic matter content of this sequence was low, with C_org mostly below 1 wt% (Fig. 8E). This is not surprising for a fluvial sequence predominantly composed of sandstone, but higher values of C_org are observed during both the LKW and UKW. Interestingly, elevated values of C_org:P_tot are observed during the LKW, at the onset of the UKW, and at the end of the UKW. Aside from C_org:P_tot, redox proxies generally do not support the existence of anoxia in the channel lakes of this fluvial system, save for the bottommost portion of the sequence as well as at 720 m during the UKW. What is most significant across all the sites presented here is that there are distinct periods when each site hosted both oxic and anoxic systems. In some cases, such as the anoxia observed in lakes 16 and 19 at Geanies, climatic fluctuations likely played a significant role in the development of anoxia as increased aridity may have led to water column stagnation, stratification, and anoxia. Other cases, such as lake 10 at Geanies, are more complex and will be discussed in detail in the following sections.

In each sequence presented, there is evidence for varying levels of P enrichment throughout a given record. While the magnitude of P enrichment varies among the different sites, the relative agreement between P/Ti and other P/Me ratios in most cases suggests that P perturbations were indeed occurring. The Geanies sequence is our earliest evidence of elevated terrestrial P input. Enrichments in P/Ti and P/Fe in lake 9, and to a lesser extent lake 10, lead the sequence, both of which gradually decrease until reaching minima at the end of lake 14 (Fig. 4A). Concurrently, the P accumulation rate is lowest at the bottom of the sequence, building rapidly through lake 10, steadily through lake 12, and then remaining stable through lake 14 (Fig. 4B). The antithetic trends observed between P/Ti and the P accumulation rate align with expectations if a large and relatively rapid nutrient flux perturbed this lacustrine system. Similar trends are observed in the upper portion of lake 16, lower portion of lake 18, and in lake 19 (lakes 18 and 19 were formed concurrently with the Kačâk Event) accompanying climatic shifts to wetter conditions with increases in P/Me ratios and accompanying increases in P accumulation. Beyond the Geanies sequence, these patterns are again observed at Ella Ø at 0.15 m, 0.85 m, and 1.10 m (Fig. 6); Hegglie Ber at 10 m and, to a lesser extent, at 25 m (Fig. 7); and at Heintzbjerg in the first 60 m of the LKW and between 600 m and 700 m in the UKW (Fig. 8). The magnitudes of the P/Ti maxima at Geanies (~3.1), Ella Ø (~3.7), Hegglie Ber (~8.2), and Heintzbjerg (~1.1) are all significantly above background for each sequence and, in most cases, are sustained across multiple samples (the magnitude at Heintzbjerg is relatively small in comparison, but this is likely due to the fluvial nature of this sequence with the lithology being primarily sandstone with fine-grained material interspersed). In many of the cases above, anoxia or suboxia is concurrent with P/Ti maxima, which introduces the possibility of reductive dissolution of Fe-bound P minerals driving increases in P/Ti values. However, in each instance there is either an accompanying increase in P/Fe or an increase in P accumulation rate, either of which would be indicative of an external P source. Additionally, strong correlations of P vs. Ti, P vs. Al, and Ti vs. Al during nearly all of these elevated P/Ti intervals indicate that terrestrial origin is likely (Supplemental Table S1).

Unlike the other four sites, Hoxa Head is significantly different and shows no evidence of significant nor sustained P/Ti maxima. The magnitudes of the P/Ti and P accumulation rate maxima are significantly lower than those of most other sites, and there are no observable trends in the data. The lack of an observable trend in P/Me ratios (correlations run between P and Me returned no significant values for almost the entire sequence; Supplemental Table S1), combined with the low rates of P accumulation, suggest there is no significant P flux occurring throughout the entirety of this 2.1 m.y. record. While this record lacks a definable P flux, it is important to consider the temporal location of this sequence (Fig. 1). Located in the early Givetian, Hoxa Head is well removed from the Kačák Event (latest Eifelian) and the Taghanic Crisis (late Givetian). This critical point is considered in detail in section 4.6.

When designing this study, it was suspected that analyses of Devonian lacustrine and near lacustrine sequences would reveal much about both terrestrial nutrient flux and geochemical cycling within these systems. What was not necessarily expected was how similar these Devonian lakes were to extant systems regarding their internal geochemical response to external nutrient flux and climatic perturbations. Our record of P input into lakes and fluvial systems ranges from the Eifelian to the early Famennian. In each case presented here, P input and climatic perturbations (or both) are catalysts for a geochemical response within the lake. This is perhaps best demonstrated in the Geanies sequence, given that it hosts pronounced climatic fluctuations as suggested by the Sr/Cu data. At the beginning of lake 16, low Sr/Cu values are concurrent with low values of C_org, and there are no indications of anoxia. A shift from low values of Sr/Cu to high values signifies a shift to more arid conditions. Coincident with Sr/Cu maxima are indications of anoxia and a rapid rise in C_org to around 4 wt%. With the development of anoxia, P is mobilized from Fe-bound minerals, and the P accumulation rate drops (only slightly in this case) and P/Ti increases, which reflects an overall increase as P is liberated from the anoxic portion of the lake. This scenario also is observed in the middle of lake 19, at 56–60 m and 174–178 m at Hoxa Head, and between 1.0 m and 1.5 m at Ella Ø. This geochemical response to arid conditions is similar to what would be expected in a modern lake.

As conditions shift from arid to wet, a similar model for lake response can be constructed. Continuing with the example of lake 16, as aridity transitions to wetter conditions, more oxic conditions would prevail within the lake as runoff increases and water column stagnation decreases. However, as conditions become wetter, redox conditions alternate between reducing and oxidized (evidenced by the brief drop in EFs as well as C_org:P_tot) until P input and P accumulation both rise significantly, leading to the continuation of anoxic conditions. This likely occurred for one of two reasons. The first is that wetter conditions contributed to a larger and deeper lake, briefly increasing mixing, but ultimately promoting the development of stratification at depth and anoxia. The second possibility is that elevated terrestrial nutrient influx, as suggested by the P/Ti values, promoted eutrophication in the water column, sustaining anoxic conditions and increasing C_org preservation. The latter is the more likely for several reasons. First, elevated values of P relative to detrital input (P/Ti and P/Al) and redox cycling (P/Fe) in the upper portion of lake 16 are coincident with anoxia, which suggests an external P source. Additionally, P accumulation increases significantly in the upper portion of lake 16. Finally, sustained high values of C_org (around 3–5 wt%) persist through the end of lake 16, which supports an increase in net primary production (NPP). This model of terrestrial nutrient-driven eutrophication during wetter periods is observed in lakes 9, 10, and at the beginning of lake 19 (prior to the onset of arid conditions). It is also observed at Hoxa Head at 12 m, 62 m, and 145 m; at Ella Ø between 0.8 m and 1.1 m; and at Hegglie Ber between 4 m and 11 m. At the fluvial sequence at Heintzbjerg, a nearly identical response is observed between 20 m and 80 m during the LKW and again between 600 m and 700 m during the UKW. Although Heintzbjerg lacks a similar redox response, the turbulence of fluvial systems tends to promote exchange with the atmosphere and more oxic conditions. Nevertheless, like the lake systems presented here, a climatically driven change to wetter conditions was accompanied by an increase in terrestrial P flux and a significant increase in P accumulation rates.

The geochemical responses to changing climatic conditions and terrestrial nutrient flux allow the construction of a basic geochemical model for Devonian freshwater systems in the Orcadian Basin and the Devonian Basin in East Greenland (Fig. 9). Intensely arid periods in this region are characterized generally by low rates of weathering, stable/background values of terrestrial nutrient input, and frequently, stratification-driven anoxia. During wetter periods, elevated terrestrial nutrient flux appears to drive eutrophication, sustaining high rates of C_org deposition and anoxia. The robustness of this model extends across each sequence presented here, beginning in the late Eifelian and extending into the early Famennian, a period spanning roughly 15 m.y. This begs the question, what is causing the increase in terrestrial nutrient flux during some wet periods and not others?

A common thread connecting the most significant P/Ti perturbations in each sequence is their occurrence during the wettest conditions. This is most evident at lakes 9 and 10 at Geanies (Fig. 4), between 0.15 m and 0.85 m at Ella Ø (Fig. 6), at 10 m at Hegglie Ber (Fig. 7), and between 600 m and 700 m at Heintzbjerg (Fig. 8). At face value, this might suggest that increased runoff resulted in increased net export. Paradoxically, some P/Ti maxima occur during periods of reduced weathering such as in lakes 9, 10, 16, and 19 at Geanies; 62 m at Hoxa Head; 1.10 m at Ella Ø; and during the first 20 m at Heintzbjerg. Thus, the explanation that increased terrestrial nutrient delivery is simply a result of increased runoff cannot fully explain the data. However, wetter conditions would produce a more favorable environment for plant growth, which could drive elevated nutrient export from the land. The significant P/Ti and P/Fe maxima observed at Geanies in lake 9, followed by the overall decrease in P input (inferred by the gradual decrease in P/Ti) from lakes 9–15, is consistent with what would be expected as plants colonize a young landscape, using up the readily available mineral phosphate (Fig. 2; Filippelli, 2008). This possibility is further evidenced by the gradual increase in P accumulation through lake 15. Similar trends between P/Me and P accumulation rates can be observed in the Heintzbjerg sequence during the UKW, with the notable differences at Heintzbjerg being lower P/Ti maxima and much greater P accumulation rates. The former was addressed in section 4.4, but the latter could be the result of several different factors that will be explored below.

The most significant P/Ti perturbations at Heintzbjerg occur during the UKW, where P accumulation rates increase from a background value of between 100 µmol P cm^–2 kyr^–1 and 200 µmol P cm^–2 kyr^–1 to peak values between 4000 µmol P cm^–2 kyr^–1 and 8000 µmol P cm^–2 kyr^–1 (Fig. 8B). These values compare strikingly well, both in magnitude and behavior, to those obtained during studies of post-glacial Holocene lakes (e.g., Filippelli and Souch, 1999) and relatively well to P accumulation rates on continental margins with ranges of 80–1200 µmol P cm^–2 kyr^–1 for the margin of Baja California and 1100–8000 µmol P cm^–2 kyr^–1 for the St. Lawrence Seaway (see Filippelli, 1997, and references therein). By this stage in the Devonian, the progymnosperm Archaeopteris had achieved hegemony in Euramerica. During both the LKW and UKW, a warm and wet climate driven by high atmospheric CO₂ concentrations (a likely result of the eruption of large igneous provinces [LIPs]; Racki, 2020) would have resulted in a favorable environment for plant growth. The favorable environmental conditions would have encouraged further expansion of Archaeopteris in Euramerica, with its substantial root systems serving to accelerate soil development in areas that had previously only been colonized by much more primitive plants. The increased plant growth and soil formation would have served to mobilize large amounts of P from relatively under-weathered or un-weathered areas. The products of this mobilization would have been effectively transported to lakes and oceans via the enhanced runoff experienced due to the wetter climate. Thus, the prevalence of deeply rooting plants—combined with climatic and tectonic drivers—could explain the marked difference in P accumulation rates between Geanies and Heintzbjerg while simultaneously addressing the similarity between Late Devonian and Holocene values.

Despite the favorable comparison with Holocene and modern P accumulation rate values, the comparison breaks down somewhat when considering the landscape stabilization timeframe. Where the Pleistocene–Holocene transition analog stabilized over a period of 5–10 k.y. (Filippelli and Souch, 1999), the Late Devonian landscape at Heintzbjerg reflects a stabilization time frame closer to 34 k.y. The Middle Devonian landscape at Geanies reflects an even greater stabilization time frame on the order of 600 k.y. There are several reasons why this could be the case. First, the lakes studied by Filippelli and Souch (1999) were all small and from upland areas in headwater catchments. Thus, given their limited catchments, it is possible that landscape stabilization occurred relatively rapidly. Second, landscape colonization by plants may have been staggered throughout the catchments, and not all areas were colonized at the same time. Given Milankovitch-driven climatic variations in the region, it is possible that the arid/wet sequences, which produced the well-defined lake cycles in the Orcadian Basin and the Devonian Basin in East Greenland, also contributed to the episodic expansion of land plants throughout the respective drainage basins and resulted in extended periods of P accumulation. Wetter periods would serve to promote growth and expansion, slowing during the subsequent arid periods, and repeating during the next climate cycle. Within each cycle, the available mineral phosphate would weather quickly, and the landscape would stabilize, with less and less total P in the system during each cycle. This would result in diminished P export over a period of ~600 k.y. as seen at Geanies. Additionally, the ubiquity and biodiversity of vascular land plants in the Holocene compared to during the Devonian may have facilitated a much more rapid stabilization timeframe in the Holocene. A variation of this point can also be utilized to partially explain the difference in stabilization rates between Geanies and Heintzbjerg. The shallow and primitive root systems of Middle Devonian land plants, coupled with their limited reach into continental interiors, likely would have resulted in significantly slower weathering rates.

The above discussion can also serve as a partial explanation for the more short-lived P perturbations that occur in the upper half of the Geanies record (including the Kačák Event) as well as at Ella Ø and Hegglie Ber. As land plants continued to evolve and expand into continental interiors and weather deeper soil profiles, P would be mobilized from previously un-weathered or under-weathered areas. It is possible that Milankovitch-driven arid/wet cycles may have been one of the key factors driving this expansion. For example, during the Kačák Event (Geanies lakes 18 and 19), Milankovitch-driven climatic cycles resulted in an insolation maximum, driving enhanced warming in the southern low latitudes (up to 6 °C warmer than elsewhere on Earth; Suttner et al., 2021). This seasonal insolation maximum would have served to reinforce the continental high-pressure cell, subsequently drawing in more moisture-laden air that would then cool and precipitate out through orographic lifting over the nearby Caledonian Mountains (Marshall et al., 2007). The increased rainfall, or monsoon, would primarily manifest in the mountains, but runoff would eventually make its way through the drainage basin, feeding the lakes in the Orcadian Basin. The increased runoff resulted in the growth and interconnection of many lakes within the Orcadian Basin, which formed a mega-lake, with some portions of the basin sustaining lakes that persisted through the accompanying arid periods and likely for up to three lake cycles or roughly 300 k.y. (Marshall et al., 2007). Our data may support this interpretation, as elevated P/Ti values in lake 18 coincide with minima in Sr/Cu and elevated Rb/Sr (Fig. 4). Enhanced terrestrial nutrient export likely would result from seasonal monsoons. It is possible that the increased rainfall in the distant mountains fueled plant growth in that region, with the nutrient pulse eventually making its way to the lakes at Geanies. While the magnitudes of the P/Ti pulses seen near the Kačák Event are modest compared with others at Geanies or our other study sites, they likely would have been influenced by the geographic separation between the Caledonian Mountains and the Geanies region. A similar scenario can be applied to other lakes in the upper portion of the Geanies sequence (i.e., lake 16) and indeed to both Ella Ø and Hegglie Ber. The lakes at Ella Ø and Hegglie Ber, however, mark a notable step forward in time to a period when both macrofossil and microfossil evidence of Archaeopteris existed within the region (i.e., Marshall, 1996; Marshall and Stephenson, 1997). At Ella Ø, fossilized stems of Archaeopteris were found at several locations throughout the sequence (Fig. 6); they likely floated away from the shore and sank in deeper water (Marshall and Stephenson, 1997). While Hegglie Ber lacks similar macrofossil evidence of Archaeopteris, there is abundant spore evidence for its widespread presence within the Eday Flags (Marshall, 1996) along with macrofossil evidence of cladoxylopsida, another major Middle Devonian vascular land plant group (Berry and Hilton, 2006). Both Ella Ø and Hegglie Ber contain significant P/Ti perturbations (the two largest maxima observed in all five sequences) that occur during the wettest period of each record (Figs. 6 and 7, respectively). Like the upper lakes at Geanies, climate appears to be driving land plant expansion. At Ella Ø and Hegglie Ber, however, more advanced vascular land plants elicit a more dramatic nutrient release. Together, this evidence would seem to support the hypothesis that P weathering was episodic in nature and linked to both climate and land plant evolution.

The final piece of crucial evidence supporting the hypothesis of episodic plant expansion is the record from Hoxa Head. Unlike the other four sequences, Hoxa Head contains no substantial evidence of a terrestrial P pulse (Fig. 5). Indeed, the variability and lack of significant trends or maxima in P/Ti suggest a relatively stable terrestrial P flux. Given its location under an arid climate between two biologic crises, the lack of significant P perturbations within this record is perhaps its most compelling aspect. While large amounts of fossilized fish were present throughout the sequence, no fossilized plant remains were noted. If the evolution and radiation of land plants played a role in the extinction events of the Devonian, then we would not necessarily expect to find significant nutrient pulses so far outside of any known biotic crisis. Additionally, the lack of macro/microfossil evidence for the presence of Archaeopteris during this period in Orkney (Marshall, 1996) suggests that deeply rooting plants had not yet begun to colonize this region. There is abundant spore evidence that suggests colonization by aneurophytalean progymnosperms, among other vascular plants (Marshall, 1996). However, as discussed above, the less complex root systems possessed by Early Devonian vascular land plants may not have been sufficient to mobilize large amounts of P from the landscape, or the initial nutrient export may have occurred earlier in this region and left no evidence later in the record. Keep in mind, however, that this sequence occurs well after the Kačák Event. In effect, Hoxa Head may be reflective of a landscape that had long since stabilized with respect to terrestrial P export and is likely to have remained in equilibrium until a significant event, such as the appearance of more deeply rooting plants, perturbed the system. If nothing else, Hoxa Head serves as evidence that despite the occurrence of numerous events in the Devonian that significantly impacted biogeochemical cycling, there were periods of regional equilibrium that likely lasted millions of years.

The observed peaks in P at multiple locations throughout the Middle to Late Devonian indicate substantial global changes in the terrestrial P cycle. While these peaks may not necessarily coincide in time or in magnitude, land plant colonization was certainly not a single punctuated event, but likely staggered geographically, peaking at different times in different parts of Euramerica and other parts of the Devonian Earth. Additionally, the evolution of root systems throughout the Devonian occurred gradually, which leaves open the possibility that more deeply rooting plants could impact terrestrial nutrient export even in areas previously colonized by more primitive vascular land plants. Thus, temporal alignment of nutrient export events is not necessarily crucial. Additionally, these results point to the conclusion that episodes of large nutrient export into the ocean were transitory events that lasted briefly at each site as colonization took hold and the initial pool of mineral-bound P was depleted. This is analogous to what is seen in post-glacial records and consistent with soil weathering trends (see Filippelli and Souch, 1999).

As plants evolved and proliferated, large-scale changes to terrestrial P inventories, weathering regimes, and soil development were all but inevitable. As P is mobilized and used within a system (or catchment), the total amount of reactive P will decrease over time (Fig. 2; Filippelli, 2008). Thus, in a large but also geographically constrained region, such as the Orcadian Basin and the Devonian Basin in East Greenland, the amount of P in the region would be expected to decline over the roughly 15 m.y. represented by these study sites. A trendline fit across all five data sets displays a shallow but significant negative trend in P/Ti ratios from the Eifelian to the early Famennian (Fig. 10A), which indicates the predicted decrease in P over time. Similarly, as plants proliferate and colonize progressively more of the landscape, a transition would occur from a predominantly physical weathering regime to a primarily chemical weathering regime. Rb/Sr is most frequently used as a chemical weathering proxy, while Si/Al is a proxy for physical weathering (see section 2.4.3). The striking antithetic relationship between Rb/Sr and Si/Al seen in the composite records (Figs. 10C and 10D) confirms that this transition took place in the region encompassing the study sites over the course of 15 m.y. Finally, as land plants colonize a landscape, significant changes in the ratio of silt to clay would be expected. More mature landscapes would be expected to have a higher clay content. Ti/Al and Zr/Al have been utilized as a proxy for clay content, as lower values of both ratios are indicative of higher clay content. While linear fits to both Ti/Al and Zr/Al plots reveal no global trends in the data (Figs. 10E and 10F), locally weighted scatterplot smoothing (LOWESS) reveals significant smaller scale trends. There is a marked decrease in the silt/clay ratio concurrent with landscape stabilization noted in the lower half of the Geanies record, with fluctuations in the latter half of the record perhaps reflecting the noted climatic instability. Through the Hoxa Head sequence, the silt/clay ratio remains low for most of the sequence, which is indicative of a more mature and stable landscape as suggested by interpretations of the P/Me data. While interpretations of the remainder of the record are mixed, Ti/Al ratios reflect increasing clay content exiting the LKW while Zr/Al ratios reflect increasing clay content exiting the UKW, which potentially indicates an increase in soil maturity during each event. Based on this evidence, it is clear that the colonization of land plants had a significant long-term impact on P inventory, weathering, and soil development in this region.

Most studies attempting to address the impact of land plants on the Devonian biosphere utilized marine records and have been crucial to our understanding of events. The redox proxy C_org:P_tot utilized in this study has been extensively reported in Devonian marine records and is often used as a proxy for benthic nutrient fluxes. In a study performed by Algeo and Ingall (2007), median C_org:P_tot ratios were compiled from 91 Phanerozoic studies, including many covering the Devonian. Algeo and Ingall (2007) identified a significant negative trend in C_org:P_tot ratios from the Middle Devonian to the Late Devonian, where median values decreased from as high as ~1625:1 to as low as 200:1. A similar negative trend is observed from Geanies to Heintzbjerg with median values dropping from 314:1–63:1 (Fig. 10B). While lacustrine and marine records would not necessarily be expected to be in agreement, it is significant that a substantial decrease in C_org:P ratios observed in marine records is reflected across our terrestrial records as well. The decrease observed in marine records is likely directly related to the increase in atmospheric O₂ concentrations and reflects a transition of increasing ocean ventilation from the Middle to Late Devonian (Algeo and Ingall, 2007). Using this model, if Devonian seas were poorly ventilated, then it is likely that lakes shared a similar fate based on the C_org:P_tot ratios reported here. A poorly ventilated sea would more easily be pushed into anoxia, as would a poorly ventilated lake.

Earth’s oceans have been described as “on the edge of anoxia” (Watson, 2016) such that even small changes in marine NPP driven by any significant increase in P delivery to oceans have the potential to drive the formation of anoxic zones (Watson et al., 2017). This is particularly relevant in the Devonian given the prevalence of shallow, restricted, and poorly ventilated seas. The particulate load of rivers would presumably have had a higher proportion of P held in primary apatite minerals, which are essentially unreactive in the marine system (Ruttenberg and Berner, 1993; Filippelli and Delaney, 1996). Some of the particulates may be in the form of early soil-forming materials that could have contained P trapped in iron oxides and oxyhydroxides. Under anoxic conditions, these minerals readily dissolve, releasing P back into the water column where it can be utilized once again, driving an increase in NPP (Filippelli, 2002; Murphy et al., 2000). If plant colonization took place in the staggered manner suggested, the presence of significant P pulses staggered temporally and geographically (but in proximity to one another with respect to geologic time), whether transient or sustained, have the potential to cause a positive feedback loop in an anoxic environment and create self-sustaining eutrophication long after termination of the P pulse.

This presents an interesting possibility regarding sustained eutrophication in Devonian oceans. Given the residence times for P in modern oceans of 15–30 k.y. (Ruttenberg, 1993; Filippelli and Delaney, 1996) and assuming a similar value for Devonian oceans, it is possible that that the geographically distinct weathering events may be a driving factor in sustained marine eutrophication. As elevated nutrient input declines in one region, eutrophication would be self-sustaining for multiple residence times, which may be long enough for colonization to occur in a separate geographic region and thus drive another distinct eutrophication pulse. Additionally, with widespread eutrophication-driven anoxia, this residence time could be longer given an active P-Fe shuttle remobilizing deep P. This would, in effect, provide the catalyst to drive the sustained eutrophication events present in Devonian oceans. Ultimately, new P supply to marine environments is controlled by the weathering of continental rocks, which is significantly increased by biota (Lovelock and Watson, 1982; Watson et al., 2017). As demonstrated by Lenton et al. (2012), the presence of a biological agent such as moss increases P weathering by as much as a factor of 60. Thus, even if plant colonization increased by just a few percent of the land surface, this would effectively double the flux of P into the Devonian oceans (Lenton et al., 2012).

Our record of P input into lakes and fluvial systems ranges from the Eifelian to the early Famennian. In each case presented here, P input is a catalyst for a geochemical response within the lake (Fig. 9). Devonian lakes appear to respond similarly to modern lakes, which brings up another series of interesting possibilities. If Devonian lakes responded no differently than modern lakes, were land plants and forest development truly that important with respect to geochemical cycling in lakes? Or (perhaps an entirely different scenario, given the similarity with modern lakes), could it be that significant plant-enhanced P weathering began before the Middle Devonian and before ecosystems had stabilized to a pseudo modern state? Did vascular land plants impact global biogeochemical cycling far earlier than we suspect, causing Devonian lakes to look so much like modern lakes? After all, it was during the Devonian that atmospheric CO₂ and O₂ reached levels comparable to those of the Holocene. Given that land plants have been around since at least the Silurian, this scenario seems more than plausible. Were it not for the record at Heintzbjerg, this argument would be supported. Given the similarity in P response during the UKW to those of Holocene analogues, there is little doubt that we are seeing initial landscape stabilization in this region. The fact that this nutrient export is concurrent with a major Paleozoic extinction event provides a compelling argument that elevated nutrient export driven by land plant expansion certainly played a role in exacerbating at least one biotic crisis in the Devonian. If this were still occurring in the Frasnian, it could have occurred earlier in the Devonian in other regions and suggests that landscape stabilization had not occurred everywhere even by the Late Devonian. The timing of regional colonization by deeply rooting plants such as Archaeopteris likely would have played a major factor in determining how and when such pulses impacted the Devonian biosphere. With the recent discovery of mid-Givetian forests containing Archaeopteris (as well as significant root systems; Stein et al., 2020), this is certainly plausible. Given the success in deciphering geochemical cycling in Devonian lacustrine environments in this study, these methods can and should be applied to other Devonian lacustrine sequences in geographically distinct locations to support or refute the implication of land plants in the biotic crises of the period. Additionally, it would be prudent to investigate Early Devonian and even Late Silurian lacustrine sequences.

Large-scale changes in the terrestrial P cycle occurred at critical evolutionary intervals in the Devonian in at least one portion of Euramerica. Significant mobilization of P occurred in close temporal proximity to global biotic crises. These transient pulses sustained a long-term eutrophication feedback. Behavior similar to that of modern lacustrine systems was observed in all five study locations with respect to internal geochemical responses to external nutrient input and climatic variability. Sites in proximity to a biotic crisis displayed evidence of enhanced nutrient input, except for the Hoxa Head sequence, which is not temporally proximal to a biotic crisis. The Hoxa Head sequence provides a >2 m.y. record of steady-state lacustrine cycling free from significant external nutrient pulses and thus provides a basis for comparison with the other four sequences.

Land plant expansion on different continents was likely staggered, but if evidence on multiple continents occurring within a few million years of each other is shown for major transient P pulses such as those reported here for East Greenland, Easter Ross, and Orkney, we can more confidently develop a temporally relevant model of nutrient inputs. Terrestrial evolution of the soil pool and marine extinction events typically last many millions of years, which makes timing differences of a few millions of years among different continents largely irrelevant. Much work remains to be done to constrain the global picture, particularly with the addition of data for lacustrine and near-lacustrine sequences from Gondwana and other parts of Euramerica and paleo-China, but putting together these pieces of the paleo puzzle, section by section, should lead to answers to the long-standing questions of the role of plant evolution in the extraordinary biotic events of the Devonian.

For support of this research, we acknowledge the Geological Society of America (Graduate Student Research Grant to M. Smart), the National Science Foundation (EAR-1850878 to G. Filippelli and W.P. Gilhooly III), and the donors of the American Chemical Society Petroleum Research Fund (59949-ND2 to W.P. Gilhooly III and G. Filippelli). Also, we thank Alice Bosco-Santos, Brooke Vander Pas, Emeline Frix, and Jonathan Fislar for invaluable assistance in the laboratory and the multiple undergraduate assistants who contributed to this project.