Early‐Warning Signals for Marine Anoxic Events

Predicting which marine systems are close to abrupt transitions into oxygen‐deficient conditions (“anoxia”) is notoriously hard but important—as rising temperatures and coastal eutrophication drive many marine systems toward such tipping points. Rapid oxic‐to‐anoxic transitions occurred regularly within the eastern Mediterranean Sea on (multi)centennial time scales, and hence, its sedimentary archive allows exploring statistical methods that can indicate approaching tipping points. The here presented high‐resolution reconstructions of past oxygen dynamics in the Mediterranean Sea reveal that early‐warning signals in these deoxygenation time series occurred long before fast transitions to anoxia. These statistical indicators (i.e., rise in autocorrelation and variance) are hallmarks of so‐called critical slowing down, signaling a steady loss of resilience of the oxygenated state as the system approaches a tipping point. Hence, even without precise knowledge of the mechanisms involved, early‐warning signals for widespread anoxia in marine systems are recognizable using an appropriate statistical approach.


Introduction
Deoxygenation of our oceans and seas is one of the most profound consequences of global environmental and climatological change (Breitburg et al., 2018;Keeling et al., 2010;Levin, 2018). Increasing water temperatures result in lower gas solubility, enhanced respiration processes, and stratification of the water column, while eutrophication increases primary productivity and related aerobic respiration of organic matter that settles to the seafloor (Keeling et al., 2010;Kemp et al., 2009). As a result, parts of the ocean become anoxic (zero oxygen). Oxygen is fundamental for aerobic life and biodiversity and plays an essential role in the biogeochemical cycling of carbon, nitrogen, and other important elements. Deoxygenation of our oceans will, therefore, have large consequences for marine organisms and biogeochemical cycling. Yet, mechanisms involved in the commonly rapid transitions into anoxia are still poorly understood (Kemp et al., 2009;Levin, 2018), hampering prediction of these potential large-scale deoxygenation events in marine systems. Once shifted into an anoxic state, such systems may remain locked in an oxygen-depleted status quo-presumably through self-sustaining feedback mechanisms, related to processes such as sedimentary phosphorus regeneration maintaining cyanobacteria blooms (Vahtera et al., 2007)-making restoration virtually impossible. This is, for example, currently observed in the Baltic Sea (Carstensen et al., 2014). Hence, it would be useful to have diagnostic indicators that can probe whether a future transition into an anoxic state is likely for a specific marine system. Here we analyze past widespread anoxic events and associated organic-rich sediments ("sapropels") in the eastern Mediterranean Sea (Figure 1), in search for such early-warning signals.
The eastern Mediterranean sapropel deposits mark some of the largest-scale transitions from oxic to anoxic conditions and back again in Earth's recent past, that is, the late Cenozoic. These transitions may hold important information to better understand large-scale anoxic events (e.g., Cretaceous oceanic anoxic events), as they involve similar reorganizations of biogeochemical cycles, such as related to nitrogen fixation (Bale et al., 2019;Bauersachs et al., 2010;Higgins et al., 2010), nitrogen removal (Rush et al., 2019), and sedimentary phosphorus regeneration (Slomp et al., 2002(Slomp et al., , 2004. Eastern Mediterranean sapropels allow to study these processes with superior spatial and chronological control, making them compelling archives to study widespread anoxic transitions (see Rohling et al., 2015, for an extensive review). Excessive monsoonal flooding from North Africa inhibited intermediate-and deep-water formation and (indirectly) stimulated nutrient input, likely causing enhanced stratification and increased organic matter production that ultimately led to these oxygen-deficient intervals (Rohling et al., 2015;Rossignol-Strick, 1985). Changes in insolation, caused by variability in the Earth's axis and orbit relative to the Sun (i.e., orbital forcing), seem to be responsible for these regular changes in Mediterranean hydrology (Lourens et al., 2001;Rossignol-Strick, 1985). Importantly, the sudden transitions from oxic to anoxic conditions (and back) occurred within decades to centuries (Marino et al., 2007), while the monsoonal changes involved were rather gradual in comparison (Weldeab et al., 2014;Ziegler et al., 2010). This suggests that these abrupt transitions into anoxia corresponded to critical transitions that happened as the system reached a tipping point where change became self-propagating due to positive feedback mechanisms (Lenton, 2013).
In general, as a system gradually approaches a tipping point, recovery rates from small stochastic perturbations become increasingly sluggish. This phenomenon is called critical slowing down (Scheffer et al., 2009. In natural time series where systems are subject to stochastic perturbations such critical slowing down is often reflected in a rising variance and temporal autocorrelation, as systems that slow down resemble more and more their previous state (Carpenter & Brock, 2006;Ives, 1995;Scheffer et al., 2009Scheffer et al., , 2012. This implies that in systems that gradually approach a tipping point, such changes in statistical indicators may be interpreted as signals for changes in the resilience of the system (Scheffer et al., 2009).
To see if indicators of critical slowing down are detectable in the periods leading up to marine anoxic events, we reconstruct past oxygen conditions in the Mediterranean Sea using redox-sensitive trace elements. Specifically, we focus on molybdenum (Mo) and uranium (U) in sediments, as Mo has been found to be an excellent indicator of strong anoxic conditions ([O 2 ] = 0 mg L −1 ) with free sulfides, while U is a proxy for suboxic ([O 2 ] = 0.2-2 mg L −1 with Fe 2+ > 0 and H 2 S = 0) conditions (Algeo & Liu, 2020;Tribovillard et al., 2006). Recent technological advances in X-ray fluorescence (XRF) core scanning , combined with the multivariate log-ratio calibration approach (Weltje et al., 2015), allow accurate measurement of these trace elements at resolutions (i.e., 1 mm in depth, relating to~10-50 years in time) necessary to observe potential early-warning signals.

Sites, Analyses, and Chronology
We reconstructed past oxygen dynamics for two relatively deep cores, MS66 (33°1.9′N, 31°47.9′E; 1,630 m water depth) and 64PE406E1 (33°18.1′N, 33°23.7′E; 1,760 m water depth), and one shallower core, MS21 (32°20.7′N, 31°39.0′E; 1,022 m water depth), in the eastern Mediterranean Sea (Figure 1). All cores contain darker-colored, organic-rich sediment intervals (i.e., sapropels), with MS21 containing Sapropels S1 and S3, MS66 containing Sapropels S1 to S5, and 64PE406E1 containing Sapropels S1 to S9 (see supporting information Figure S10 for high-resolution pictures). Note that the existence of Sapropel S2 (~55 kyr BP) is still controversial, and in any case, this sapropel interval is mostly missing or very weakly expressed in eastern Mediterranean sediments (Bar-Matthews et al., 2000) and thus is not included in our analyses. Moreover, the top of core 64PE406E1 contains less than 1 cm of Sapropel S1, which seems disturbed due to piston coring, and therefore, 64PE406E1 Sapropel S1 was excluded from the early-warning-signal statistics. We minimized the effects of bioturbation on our statistical results by avoiding clearly bioturbated parts of the sediment for our geochemical analyses (see section 4.1 for a detailed discussion on the impact of bioturbation).
Bulk geochemical analyses were performed on all three cores with "trace-metal settings"  using an Avaatech XRF core scanner. Conversion of XRF core scanning data into concentrations was done using 440 calibration samples from conservative geochemical techniques implemented in the multivariate log-ratio calibration approach (Weltje et al., 2015). The obtained 1-mm sampling rate from XRF core scanning relates to an~10-50 year resolution, which allows evaluation whether increased variability and autocorrelation are reliable early-warning indicators for tipping points toward marine anoxia. Concurrent variability in barium, an export-productivity proxy (De Lange et al., 2008;Dymond et al., 1992), in these and neighboring cores (Grant et al., 2016;Ziegler et al., 2010) was used to intercalibrate age models. A more detailed description of the samples and geochemical analyses (Text S1) and chronological approach (Text S2) can be found in the supporting information.

Early-Warning Signal Statistics
All analyses were implemented in Rv3.4.3. (R project for Statistical Computing). The leading indicators were calculated using scripts that are based on the "earlywarnings" package in R . We selected the part of the time series that preceded the fast transition to an anoxic state (i.e., before the abrupt increase in the anoxia-proxy Mo to high values). The exact positions of transitions were sometimes placed slightly before the published ages for Sapropels S1 to S9 formation, to avoid including data points that were part of the transition itself, because inclusion of points that are part of the transition would bias the estimate of the autocorrelation at lag-1 (Dakos et al., 2008). These transitions toward anoxia with H 2 S occured typically within 500-1,000 years, in line with time intervals found for Sapropel S5 based on lower-resolution proxies (Marino et al., 2007). Sapropels and related anoxia are thought to be mainly driven by precessional forced monsoon variability (Lourens et al., 2001;Rohling et al., 2015;Rossignol-Strick, 1985), and therefore, we focus on 10 kyr (about half a precession cycle) before each transition.
We estimated the autocorrelation and variance (i.e., resilience indicators or potential "early-warning signals") within rolling windows of predetermined size up to the end of the records. The rolling window size was set at half the size of the records. Prior to these analyses, the time series were detrended by subtracting a Gaussian kernel smoothing function. The likelihood of obtaining our computed trends in these resilience indicators by chance were tested using 1,000 surrogate time series. A detailed description of detrending and calculating resilience indicators (Text S3) and testing of their significance (Text S4) can be found in the supporting information.

Geochemical Proxies and Rapid Anoxic Events
Repetitive episodes of eastern Mediterranean anoxia, sapropels numbered S1 to S9, and their link to marine organic matter production are clearly indicated by proxies for anoxic conditions (Mo), suboxic conditions (U), and export productivity (barium to aluminum [Ba/Al]), which all are increased during sapropel deposition (Figures 2b-2d). Ba/Al is used as export-productivity proxy, given that biogenic barite is enriched in decaying organic matter settling to the seafloor (Dymond et al., 1992). The elevated Ba/Al values indicate that marine productivity increased during these intervals, most likely due to additional nutrients by Nile River discharge, nutricline shoaling, and/or recycling of nutrients from the seafloor (Rohling et al., 2015). The association of sapropels with insolation-driven monsoon intensifications is indicated by the close correlation with titanium-to-aluminum ratios (Ti/Al; Figure 2a). In this area, particles relatively high in Ti are mainly of aeolian origin, while Nile-derived clay is relatively high in Al (Lourens et al., 2001). Hence, Ti/ Al values indicate dry/wet phases in North Africa, with the lower values associated with sapropels pointing to an increased monsoon intensity (Figure 2a; note that the Ti/Al y axis is descending). The intimate link of these monsoonal changes to orbital forcing is shown by the correspondence to insolation (Lourens et al., 2001;Rossignol-Strick, 1985) (Figure 2a).
The similarity of anoxia and productivity variability in the overlapping intervals of the three cores (Figures 2b-2d) suggests that these changes truly reflect a regional signal. The high-resolution anoxia records, based on Mo, indicate a rapid oxic-to-anoxic transition at the start of each sapropel (Figure 2c).

Geophysical Research Letters
Comparison of Mo with U shows that suboxic conditions generally occurred over a wider interval prior and after the sapropels, for example, clearly visible in the~250-160 kyr interval (Figure 2d). Importantly, tipping-point behavior seems especially apparent for the transitions toward anoxic conditions with free H 2 S, and therefore, we focus on Mo with the early-warning indicators (see section 2.2). Yet, for comparison, the early-warning indicators were also calculated for U and Ba/Al. Figure 3 shows the pre-anoxic-event autocorrelation and variance for anoxia proxy Mo in the three cores analyzed. Increasing variance is observed prior to all analyzed anoxic events in the relatively deep cores 64PE406E1 and MS66, while a rising temporal autocorrelation occurs prior to most of the analyzed anoxic events in these two deeper cores (Figure 3). Similarly, variance shows positive and high Kendall's τ values indicating strong rising trends before all events, while autocorrelation also shows mostly positive Kendall's τ values, except for Sapropel S6 (~0) and weak anoxic intervals S1 and S9 (Figure 3 and Table S1) in these two cores. Early-warning signals (i.e., increasing autocorrelation and variance) are generally absent or only weakly expressed in the two sapropels (S1 and S3) observed in the relatively shallow core MS21 (Figure 3). We note that statistical results are almost identical when calculated for the Mo/Al ratio (see Figure S3).

Early-Warning Signal Statistics
We explored the likelihood that these trends in variance and autocorrelation would be found just by chance, without an underlying critical slowing down of the system. We tested the occurrence of trends in a large number of randomized time series with the same power spectra and length as the original records, but with random phases in the Fourier modes (Ebisuzaki, 1997). These analyses (Table S2) indicated that the probability of finding the increase in variance detected prior to the abrupt transition in the deep cores (64PE406E1 and MS66) by chance is very low. For autocorrelation, the probability of finding these trends is higher. The probability of finding this whole set of p values by chance, however, appears to be very small for both variance (Fisher's combined probability test p < 10 −13 ) and autocorrelation (p < 0.003) (see Table S2).
We also analyzed the autocorrelation and variance for suboxia proxy U and the export-productivity proxy Ba/Al (see Figures S8 and S9, respectively). The results from the early-warning statistics for suboxia proxy U are in line with results from Mo. The data from U show a clear rising variance before most events (except for Sapropel S5 in MS66), while a rising autocorrelation occurs also for most events (except S4 in MS66 and S8 and S9 in 64PE406E1) in U, albeit with lower Kendall's τ values. On the other hand, Ba/Al shows mixed trends in the variance and temporal autocorrelation. For some events, the early-warning signals show a strong rising trend (S4, S6, S8, and S9), but the rise in Ba/Al occurs already thousands of years before the abrupt transitions observed in Mo for these intervals. These strong increases bias the estimation of the autocorrelation because of serial correlation (Dakos et al., 2008). Moreover, most of the other events do not show any rises in autocorrelation and variance in Ba/Al. As such, this suggests that the increasing variance and autocorrelation in anoxia/euxinia (Mo) are not directly reflecting changes in productivity only.
A systematic sensitivity analysis for bandwidth and window size was performed to test the robustness of our statistical results. This sensitivity analysis confirms that the observed increases in variance and autocorrelation in Mo are robust and, in fact, that more significant trends could have been obtained by choosing another bandwidth and window size for each specific record (see Text S5 for details). There is no need to interpolate the data as the spacing between the data points is almost regular. Interpolating in this case would have resulted in unwanted aliasing effects that strongly affect autocorrelation and variance by consistently raising the trends (see Text S6 for details).

Early-Warning Signals for Widespread Anoxic Events
The results suggest that the onsets of sapropel intervals with widespread, long-lasting anoxia with free H 2 S (euxinia), as indicated by Mo (Tribovillard et al., 2006), were indeed associated with rapid transitions at tipping points. Moreover, these results suggest a gradual loss of resilience when the system approaches such tipping points, resulting in a slowdown of the dynamics long (i.e., centuries to millennia) before the actual transitions (Figure 3). Slowing down in this case is expressed as a decelerating recovery back to oxygenated conditions after short-lived euxinic intervals. Our data are in line with studies using benthic faunal responses, which show that deoxygenation started several thousand years prior to the onset of sapropels (Schmiedl et al., 2003). Yet, our high-resolution trace metal profiles additionally show that the intensity of low-oxygen conditions prior to sapropels fluctuated on short time scales. High-frequency changes in redox conditions occurred thus not only during these sapropel events (Casford et al., 2003;Jilbert et al., 2010) but likely also prior to these events. In many settings, presapropel sediments are overprinted by a syndepositional diagenetic downward sulfidization process (e.g., Passier et al., 1996Passier et al., , 1997 resulting in a distinct gray interval underlying most sapropels. This gray band was initially interpreted as a distinct change in depositional conditions (i.e., the so-called "proto-sapropel" concept; Anastasakis & Stanley, 1984;Maldonado & Stanley, 1976;Murat & Got, 1987). Our findings point to mainly (sub)oxic conditions with short temporal changes to lower oxygen conditions and not a distinct change in depositional conditions. These short intervals with lower oxygen concentrations prior to the actual sapropels probably had an increasing intensity, as variance increases before all sapropels at depths >1,600 m.
Comparison of the autocorrelation and variance signals shows that increasing variance is the more consistent early-warning signal for rapid deoxygenation. The latter observation might indicate a phenomenon called flickering, in which the system starts jumping back and forth between alternative stable states (Scheffer et al., 2009;Wang et al., 2012), for example, when short periods of suboxia to euxinia occurred prior to the ultimate sudden shift to sapropel formation. The decrease in autocorrelation prior to the sudden shift might indicate the necessity of even higher-resolution records prior to some sapropels, as Wang et al. (2012) showed that lower-resolution subsampling of modeled time series of fast transitions in a similar setting (a eutrophic lake state) can conceal the autocorrelation increase while preserving the rising variance signal. Similarly, in a setting like the Mediterranean Sea, the ultimate resolution that can be achieved will be limited by bioturbation, which could be an additional explanation for the stronger rising variance signal relative to the autocorrelation signal in our records.
An impact of bioturbation on our statistical results cannot be fully excluded, as intervals underneath sapropels must have been oxygenated for significant amounts of time. As such, part of the observed signal that we record may be due to an increased preservation when bioturbation diminished toward a sapropel. The impact of bioturbation seems generally low closer toward the sapropels (see Figure S10 for core material pictures), while most of the observed increases in variance and autocorrelation persist within these intervals. Moreover, the results for the different events in the two deep cores (MS66 and 64PE406E1) are quite reproducible, which increases our confidence that we mainly capture the initial environmental conditions with our proxies.
The short anoxia intervals prior to sapropels that are recorded in these cores are likely an expression of a fast expanding and retracting area of (sulfidic) anoxia with its center in the deeper part of the basin-quite analogous to the spatiotemporal development of strong anoxia observed in the Baltic Sea from 1906 to 2012 (Carstensen et al., 2014). Reducing conditions ([O 2 ] = 0.2-2 mg L −1 ) generally prevailed for much longer time intervals and over much larger areas prior to sapropels, as indicated by suboxic conditions indicator U (Figure 2d), and in accordance with results from benthic foraminifer assemblages (Schmiedl et al., 2003). These reducing conditions may have been present as an oxygen-deficient "blanket" above the seafloor (Casford et al., 2003) that could rapidly turn anoxic under influence of (stochastic) environmental change. The early-warning signals for widespread strong anoxia are more apparent at depths >1,600 m-suggested by the absence of rising autocorrelation and variance in the record of shallow core MS21 (~1,000 m) for Sapropels S1 and S3 (Figure 3). This concurs with reconstructions showing the strongest anoxia at the greatest depths during sapropels (De Lange et al., 2008), and by inference, short-lived euxinic intervals seem to also develop earlier at larger depths. Similarly, the lack of increase in autocorrelation at~1,600-1,800 m during the relatively weak anoxic Sapropels S1 and S9 suggests that these signals for slowing down were either not well recorded (i.e., we miss resolution, see previous paragraph) or simply absent at these depths. This absence may be likely as especially the deep Mediterranean Sea below 1,800 m was devoid of oxygen during the Sapropel S1 interval (De Lange et al., 2008), thus explaining a lack of early-warning signals at core sites shallower than that depth, which may also apply to the relatively weak anoxic Sapropel S9.

Positive Feedback Mechanisms and Rapid Oxic-to-Anoxic Transitions
Bistability between oxic versus anoxic states with rapid transitions in marine systems, such as the eastern Mediterranean Sea, is likely driven by strong positive feedback mechanisms. A potential mechanism is the self-enforcing feedback between anoxic conditions and productivity. Increased productivity leads to anoxic conditions due to increased organic carbon burial and remineralization. In turn, anoxic conditions cause an increase of sedimentary phosphorus (P) efflux back to the water column (Kemp et al., 2009;Slomp et al., 2002Slomp et al., , 2004. Models that include these feedbacks are bistable and show critical slowing down in model simulations (Wang et al., 2012). Sedimentary phosphorous regeneration is well documented during sapropels in the Mediterranean Sea (Slomp et al., 2002(Slomp et al., , 2004. The mixing of the Mediterranean Sea remained much shorter (i.e., maximum of~1,000 years for the strong anoxic Sapropel S5; Andersen et al., 2018) than the duration of these events, allowing sufficient time for phosphorous to reach the surface waters and stimulate primary productivity, such as nitrogen-fixing cyanobacteria blooms (Bale et al., 2019;Bauersachs et al., 2010). At the same time, also other positive feedback mechanisms may have been at play. For instance, models of the thermohaline circulation show a bistablility due to the salt-advection feedback (Stommel, 1961), showing critical slowing down in model simulations (Dakos et al., 2008;Rahmstorf et al., 2005). This may be important for the Mediterranean thermohaline circulation as well, although needing independent confirmation.
The critical slowing down that we observe in our deoxygenation proxy record prior to the anoxic events does not necessarily allow unraveling which positive feedback mechanisms are most relevant. Still, positive feedback processes are likely important in shifting an ocean basin swiftly into widespread anoxia and sustaining such a state. Importantly, even though our statistical results apply directly to semi-enclosed seas, they may also be relevant for coastal environments and on global spatial scales (albeit on longer time scales) as the biogeochemical and physical oceanographic processes and feedbacks are quite similar (Breitburg et al., 2018;Keeling et al., 2010;Kemp et al., 2009;Levin, 2018).

Conclusions
We reconstructed variability in seafloor (de)oxygenation in the eastern Mediterranean Sea over the past 250 ka in high resolution. These records reveal that the repeated shifts into anoxia ("sapropels") in this basin had the character of critical transitions. Such transitions happen when a system reaches a tipping point where change becomes self-propelling due to positive feedback mechanisms. These geochemical results show, for the first time, empirical evidence that oxygen fluctuations in a marine system show characteristic slowing down well in advance of a shift to anoxia, in line with model data. As such, the statistical indicators associated to characteristic slowing down (i.e., rise in autocorrelation and variance) seem useful early-warning signals to predict the likelihood of a fast transition to an anoxic state.