Dowody z rdzeni lodowych świadczące o dużych erupcjach wulkanicznych na początku ocieplenia Dansgaard-Oeschger / Johannes Lohmann i Andersa Svensson

2.1  Grenlandzkie zapisy izotopów o wysokiej rozdzielczości i stos

Bierzemy pod uwagę zapisy δ ¹⁸ O o wysokiej rozdzielczości z czterech głębokich rdzeni lodowych Grenlandii, a także stos pochodzący z nich. Zapisy zostały zmierzone przy różnej rozdzielczości głębokości, przy czym każdy pomiar wykonano na materiale masowym o ciągłych odstępach głębokości. Pomiary nie są zatem próbkami punktowymi, ale średnimi w ciągłych odstępach. Wykorzystujemy zapisy δ ¹⁸ O z rdzenia lodowego NGRIP ( NGRIP Members ,  2004 a ; Gkinis i in. ,  2014 ) , a także rdzeni lodowych GRIP ( Johnsen i in. ,  1997 ) , GISP2 ( Stuiver i Grootes ,  2000 ) i NEEM ( Rasmussen i in. ,  2013 ; Gkinis i in. ,  2021 ) . Wszystkie rdzenie są zsynchronizowane z roczną warstwową chronologią rdzeni lodowych z Grenlandii 2005 (GICC05) ( Svensson i in. ,  2006 , 2008 ; Rasmussen i in. ,  2013 ; Seierstad i in. ,  2014 ) . Rekordy δ ¹⁸ O są przetwarzane w następujący sposób. Punkty środkowe przedziałów głębokości są interpolowane liniowo do skali czasowo-głębokiej GICC05, co daje nierównomiernie rozmieszczone szeregi czasowe. Następnie te szeregi są nadpróbkowane do 1-rocznej siatki równoodległej przy użyciu interpolacji najbliższego sąsiada. W ten sposób charakter pomiarów jako ciągłych średnich głębokości i oryginalnych wartości pomiarów są zachowywane, a wszystkie rekordy są umieszczane na tej samej równoodległej siatce czasowej. Na koniec uśredniamy cztery rdzenie lodowe, aby uzyskać stos ze znacznie zmniejszonym szumem o wysokiej częstotliwości. Szacujemy początki DO zarówno ze stosu, jak i z pojedynczych rdzeni lodowych, z wyjątkiem GISP2, ze względu na porównywalnie niską rozdzielczość surowych danych w stadiach: zaczynając od GS-2, rozdzielczość próbki spada z 8 lat do około 30 lat w GS-15.2. Dla porównania, rozdzielczości stadiów dla tych samych okresów czasu NGRIP, GRIP i NEEM wynoszą odpowiednio 3 do 6, 4 do 9 i 4 do 8 lat. Rozważany przez nas okres czasu rozciąga się od 11 700 lat b2k (lata przed 2000 r. n.e.) do 60 000 lat b2k (GS-18). Nie bierzemy pod uwagę zdarzeń DO 2.1 i 2.2, ponieważ nie uzyskano dopasowania bipolarnego dla tego przedziału.

2.2  Zapisy rdzeni lodowych z ostatniego wulkanizmu lodowcowego

Podstawowym zestawem danych wulkanicznych użytym w tym badaniu jest ostatni zapis lodowcowy bipolarnego wulkanizmu uzyskany przez SVE20. Zawiera on N = 82 erupcji wulkanicznych, dla których znaczące osadzanie się siarczanów można było zidentyfikować synchronicznie w rdzeniach lodowych Grenlandii i Antarktydy. Z niewieloma wyjątkami, wszystkie te erupcje mają nieznane źródło. Podczas gdy zestaw danych obejmuje przedział 11,7–60 tys. lat, istnieje luka w danych od 16,5 do 24,5 tys. lat z powodu trudności w synchronizacji rdzeni Grenlandii i Antarktydy w wyniku wysokiego poziomu szumu zanieczyszczeń w późnym zlodowaceniu. Pomijając tę ​​lukę, daje to czas trwania Δ T = 40  300 lat badanych w tym badaniu.

Aby przeanalizować statystyki tego zestawu danych i ograniczyć liczbę bipolarnych erupcji o podobnej skali, które mogą w nim brakować, bierzemy również pod uwagę zapis LIN22 obejmujący ten sam okres uzyskany za pomocą zautomatyzowanej metody wykrywania erupcji wulkanicznych z ciągłych pomiarów siarczanów w rdzeniach lodowych z Grenlandii i Antarktydy. Na podstawie tych pomiarów Lin i współpracownicy wyprowadzili dwie oddzielne listy erupcji wulkanicznych i ich całkowitego osadzania siarczanów, biorąc pod uwagę erupcje wykryte indywidualnie odpowiednio na Grenlandii lub Antarktydzie. Listy te zawierają erupcje z osadzeniem siarki przekraczającym 10 kg km ^-2 i 20 kg km ^-2 odpowiednio na Antarktydzie i Grenlandii, co odpowiada mniej więcej połowie osadzania z erupcji Tambora z 1815 r. n.e. Co ważne, listy zawierają również bipolarną erupcję z SVE20. W przypadku tych bipolarnych erupcji, na podstawie osadzania się siarki na obu biegunach Lin i współpracownicy obliczyli ładunek aerozolu stratosferycznego jako miarę globalnego wymuszania klimatycznego, którego używamy do oceny, jak duże są bipolarne erupcje w porównaniu do bardziej znanych niedawnych erupcji (sekcja  2.4 ). Ponadto, na podstawie względnego osadzania się na Grenlandii w porównaniu z Antarktydą, Lin i współpracownicy podają przybliżone szacunki szerokości geograficznej erupcji, co jest omówione bardziej szczegółowo w sekcji  3.7 .

Erupcje pozostałe z zestawów danych LIN22 po usunięciu erupcji bipolarnych z SVE20 będą dalej nazywane zestawami danych unipolarnych . W przedziale 11,7–60 tys. lat minus luka danych bipolarnych 16,5–24,5 tys. lat daje to odpowiednio 740 i 498 erupcji w zestawie danych grenlandzkich i antarktycznych. Porównując osady siarczanowe z tych erupcji unipolarnych z bipolarnymi, możemy wyprowadzić górny szacunek, ile niezidentyfikowanych erupcji bipolarnych o podobnej skali może brakować w zestawie danych SVE20, zakładając, że zestaw danych LIN22 nie zawierałby żadnych regionalnych (stosunkowo blisko miejsc występowania rdzeni lodowych) erupcji z dużym (unipolarnym) osadem, ale małym globalnym wpływem.

2.3  Określenie początku zdarzenia DO

Określamy czas rozpoczęcia ocieplenia DO z zapisów δ ¹⁸ O o wysokiej rozdzielczości z pojedynczych rdzeni lodowych, a także z zapisu skumulowanego. Dla każdego cyklu DO wykrywamy, kiedy dane zaczynają znacząco odchylać się w górę od średniej stadialnej przed nagłymi przejściami. Metoda jest szczegółowo wyjaśniona w Załączniku A. Zidentyfikowane początki dla stosu izotopowego są pokazane razem z zapisem izotopowym i najbliższą bipolarną erupcją wulkaniczną poprzedzającą początek na Rys. S1 w Suplemencie. Szacowane czasy rozpoczęcia w poszczególnych rdzeniach są porównywane ze stosem na Rys. S2, gdzie względne czasy rozpoczęcia są podane w odniesieniu do najwcześniejszego z czterech niezależnych ustaleń początku.

To porównanie pozwala nam ocenić wiarygodność szacunków początku ze stosu i uzyskać oszacowanie niepewności czasowej. Początek jest bardzo spójny dla zdarzeń 3, 4, 5.2, 6, 7, 8, 10, 16.1 i 16.2, dla których zakres czterech niezależnych szacunków początku wynosi 25 lat lub mniej. W przypadku zdarzeń 1, 5.1, 9, 11, 13, 14, 15.1, 17.1 i 17.2 początek stosu grupuje się z dwoma innymi rdzeniami lodowymi w okresie krótszym niż 20 lat, podczas gdy początek pozostałego rdzenia lodowego jest bardziej oddalony od tego klastra i może być uważany za wartość odstającą. Niektóre wartości odstające mogą być spowodowane niedoskonałą synchronizacją, jak w przypadku zdarzeń 5.1 i 17.1. Tutaj punkty dopasowania używane do synchronizacji NEEM z skalą czasową GICC05 są rozmieszczone kilkaset lat przed i po początku DO, a zatem prawdopodobnie występuje przesunięcie w czasie w momencie początku. W innych przypadkach kształt początku w jednym z rdzeni lodowych nie jest tak nagły jak w innych lub jest bardzo zaszumiony, co prowadzi do późnego wykrycia przez nasz algorytm. Może również wystąpić wczesne wykrycie, gdy przed początkiem występuje duży skok w szumie o wysokiej częstotliwości, jak w przypadku GI-14 w NEEM. W przypadku pozostałych zdarzeń 12 i 15.2 przejście DO wydaje się składać z dwóch kroków we wszystkich rdzeniach, a ze względu na wysoki poziom szumu nasz algorytm nie jest w stanie wykryć pierwszego kroku w NGRIP, który naszym zdaniem jest prawdziwym początkiem. Podsumowując, początki pochodzące ze stosu należy uznać za najbardziej reprezentatywne. Po pierwsze, dzieje się tak dlatego, że stos jest lepiej przystosowany do definiowania początków ze względu na ulepszony stosunek sygnału do szumu. Po drugie, początki stosu są zgodne z szacunkami czasowymi z poszczególnych rdzeni. Odrzucając wartości odstające, stwierdzamy, że rozrzut indywidualnych szacunków początków (w tym szacunków ze stosu), tj. różnica między najwcześniejszym a najpóźniejszym szacunkiem początku, wynosi średnio 12,9 roku. Może to służyć jako dobry szacunek niepewności początków określonych ze stosu izotopowego. Dlatego możliwe jest oszacowanie reperkusji klimatycznych erupcji wulkanicznych w skali dekadowej.

2.4  Statystyki występowania erupcji bipolarnych

Aby sprawdzić, czy zdarzenia DO występują preferencyjnie krótko po erupcjach wulkanicznych, musimy znać częstotliwość występowania tych ostatnich w czasie i sprawdzić ważność naszej hipotezy zerowej. Zestaw danych SVE20 obejmuje przedział czasowy od 11,7 do 60 tys. lat minus luka danych od 16,5 do 24,5 tys. lat, co daje łącznie 20 641 lat warunków stadialnych z 42 erupcjami wulkanicznymi i 19 659 lat warunków międzystadialnych z 40 erupcjami wulkanicznymi. Daje to zadziwiająco podobne częstotliwości występowania erupcji bipolarnych λ = 2,0348 erupcji na tysiąc lat dla stadiów i λ = 2,0347 dla interstadiów, co odpowiada okresowi powrotu wynoszącemu blisko 500 lat. W okresach stadialnych rdzenie lodowe Grenlandii wykazują znacznie wyższy poziom różnych sygnałów zanieczyszczeń, co uważa się za wynik zmian w cyrkulacji atmosferycznej wzdłuż szlaku transportu, zmian wiatru u źródła i zmian w cyklu hydrologicznym w porównaniu z okresami międzystadialnymi ( Ruth i in. ,  2003 ) . Podczas gdy wyższy poziom tła w zapisach siarczanów stadialnych z pewnością obejmuje wiele mniejszych erupcji regionalnych, a także małe erupcje bipolarne, ma on również potencjał, aby utrudnić wykrywanie dużych erupcji bipolarnych, które są brane pod uwagę w tym badaniu. Jednak pomimo wyższych poziomów siarczanów tła stadialnego, podobne wskaźniki występowania wskazują, że nie ma systematycznego niedoliczania erupcji w stadiach w porównaniu z interstadiami.

Ponieważ zbiór danych SVE20 nie został uzyskany przy użyciu zautomatyzowanej metody, pozostaje obawa, że ​​może brakować znacznej liczby erupcji bipolarnych o podobnej wielkości i że może występować systematyczne odchylenie w kierunku preferowanego wykrywania erupcji w pobliżu nagłych przejść DO. Dopasowanie bipolarnych w SVE20 osiągnięto poprzez jednoczesne liczenie warstw w rdzeniach Grenlandii i Antarktydy. To liczenie warstw nie było w pełni ciągłe, ale obejmowało prawie cały okres, przy czym brakowało 16,5% rekordu (patrz tabela S2 w SVE20). Ponadto metoda dopasowania bipolarnego z założenia nie jest ciągła, ale pociąga za sobą identyfikację wzorców typowo trzech do pięciu erupcji zarówno na Grenlandii, jak i na Antarktydzie, które można dopasować na podstawie liczby rocznych warstw znalezionych między erupcjami. Jednakże biorąc pod uwagę średni odstęp między erupcjami wynoszący 500 lat, można łatwo zauważyć, że te wzorce obejmowałyby większość całego okresu czasu, nawet gdyby wybrano tylko wzorce skupione wokół 20 początków DO. Mimo to, nawet gdyby cały okres został zbadany tak dokładnie, jak to możliwe przez SVE20, niektóre duże erupcje bipolarne mogły zostać pominięte w objętych przedziałach czasu i musimy wykazać, że potencjalne niedoliczenie erupcji bipolarnych nie narusza przydatności zestawu danych i hipotezy zerowej, a także solidności naszych wyników. Jest to omówione poniżej poprzez porównanie wielkości i okresu powrotu zaobserwowanych w zestawie danych SVE20 z dobrze znanymi historycznymi erupcjami, a także poprzez wyprowadzenie górnego szacunku liczby potencjalnie niezidentyfikowanych erupcji bipolarnych o tej samej wielkości i przetestowanie solidności naszych wyników względem niej (sekcja  3.2 i 3.5 ). Obawy dotyczące systematycznego odchylenia katalogu bipolarnego SVE20 w stronę okresów nagłych przejść omówiono w rozdziałach  3.1 i  3.2 , gdzie porównuje się wielkość erupcji w pobliżu początków DO z erupcjami w innych miejscach.

W porównaniu do ponad 80 erupcji bipolarnych uzyskanych przez ostatnie 2500 lat przez Sigl i in. ( 2015 ) , jak również nowszej analizy rdzeni lodowych obejmujących cały holocen, które dają jedną erupcję z odciskiem bipolarnym co 35 lat ( Sigl i in. ,  2022 ) , erupcje rozważane tutaj są rzeczywiście dość rzadkie. Jest to jednak oczekiwane ze względu na przerzedzanie się warstw rdzeni lodowych wraz z głębokością, niższe wskaźniki akumulacji w okresie lodowcowym i znacznie wyższe poziomy tła sygnałów zanieczyszczeń w czasie ostatniego zlodowacenia ( Mayewski i in. ,  1997 ; Schüpbach i in. ,  2018 ) . Dlatego też zapisy siarczanowe z holocenu pozwalają na wykrycie znacznie mniejszych erupcji bipolarnych. Oczekujemy zatem, że w naszym zestawie danych będą obecne tylko erupcje o bardzo dużej skali. W zestawie danych Sigl et al. ( 2015 ) stwierdzono, że pięć erupcji w ciągu ostatnich 2500 lat było większych niż erupcja Tambora z 1815 r. n.e. pod względem ich bipolarnego osadzania siarczanów. Wielkość tych zdarzeń występujących raz na 500 lat można porównać do bipolarnego osadzania siarczanów podczas rozważanych tutaj erupcji lodowcowych, które zostało oszacowane przez LIN22 (patrz również sekcja  2.2 ). Stwierdzono, że 69 erupcji bipolarnych ma ładunek aerozolu stratosferycznego większy niż erupcje Tambora. Tak więc większość erupcji w zestawie danych SVE20 mieści się w kategorii zdarzeń występujących raz na 500 lat pod względem ich wielkości. Bardzo dobrze odpowiada to faktowi, że częstotliwość występowania erupcji wynosi rzeczywiście raz na 500 lat, jak pokazano powyżej. Innym znanym kandydatem do potwierdzenia częstotliwości erupcji w zestawie danych SVE20 jest erupcja Samalas z 1257 r. n.e. W zestawie danych SVE20 jest 29 erupcji, które według szacunków LIN22 mają większy ładunek aerozolu niż erupcja Samalas. Ten podzbiór erupcji daje okres powrotu wynoszący 1390 lat. Z danych w Sigl et al. ( 2015 ) erupcja Samalas wyłania się jako największa erupcja ostatnich 2000 lat lub druga co do wielkości erupcja ostatnich 2500 lat. Tak więc nasz okres powrotu jest również zgodny z ostatnimi szacunkami rdzeni lodowych dotyczącymi częstotliwości erupcji tych jeszcze większych erupcji.

Pożądane jest niezależne porównanie z zapisami wulkanizmu pochodzącymi ze źródeł innych niż rdzenie lodowe. Na przykład skompilowany zapis Rougier et al. ( 2018 ) daje okresy powrotu dla erupcji M6 i M7 wynoszące odpowiednio 110 i 1200 lat, przy czym M7 oznacza przedział wielkości erupcji obejmujący zarówno Tamborę, jak i Samalas. Jednak liczby te opierają się na szacunkach masy wybuchowej znanych erupcji, która nie jest bezpośrednio skalibrowana do polarnego osadzania siarki w rdzeniach lodowych. Ponadto, chociaż szacunki M6 i M7 autorstwa Rougier et al. ( 2018 ) opierają się wyłącznie na danych z ostatnich 2000 lat, dane te są prawdopodobnie nadal niekompletne w porównaniu z zapisem rdzeni lodowych, ponieważ zawierają one jedynie erupcje o znanym źródle. W badaniu Sigl et al. ( 2015 ) , 12 z 25 największych bipolarnych erupcji z ostatnich 2500 lat nie można było dopasować do erupcji o znanym źródle. Tak więc ilościowe porównanie danych z SVE20 i LIN22 jest obecnie możliwe tylko w odniesieniu do młodszego zapisu rdzenia lodowego. To porównanie, jak podano powyżej, wskazuje, że zestaw danych SVE20 jest stosunkowo kompletny, biorąc pod uwagę, że częstotliwość erupcji tej wielkości i towarzyszące jej polarne osadzanie siarki pozostały mniej więcej stałe. Mówiąc, że zestaw danych jest stosunkowo kompletny, mamy na myśli, że obejmuje on wystarczająco dużą część bipolarnych erupcji powyżej pewnego progu wielkości, co odpowiada erupcjom z okresami powrotu wynoszącymi jeden na 500 lat i dłuższymi. Natomiast nie mamy na myśli, że stanowi ona kompletny katalog wszystkich bipolarnych erupcji dowolnej wielkości, które można wykryć w bardziej rozdzielczych i lepiej zsynchronizowanych zapisach rdzeni lodowych, takich jak te z holocenu ( Sigl i in. ,  2022 ) . Szacunkowa górna granica liczby dużych bipolarnych erupcji potencjalnie niezidentyfikowanych przez SVE20 jest podana w sekcji  3.2 .

W dalszej części sprawdzamy, czy czasy erupcji bipolarnych są zgodne ze stacjonarnym procesem Poissona. Najpierw musimy sprawdzić, czy rozkład czasów oczekiwania t między erupcjami jest zgodny z rozkładem wykładniczym o prawdopodobieństwie skumulowanym . Tak jest w istocie, co widać w teście Andersona–Darlinga (Kołmogorowa–Smirnowa) z p = 0,60 ( p = 0,66 ). Odpowiednia funkcja rozkładu empirycznego wraz z P _λ ( T ≥ t ) przy użyciu λ = 2,0348 jest pokazana na rys. S3a. Po drugie, potwierdzamy, że własność bezpamięciowości czasów oczekiwania t jest zachowana. Osiąga się to za pomocą dwustronnego testu hipotezy bootstrapowej na korelacji Spearmana kolejnych czasów oczekiwania, uzyskując p = 0,871 dla korelacji danych r _S = 0,019 . Na koniec testujemy założenie stałej szybkości λ w czasie, dzieląc zapis wulkaniczny na krótkie, ciągłe segmenty Δ T lat i testując, czy liczba n erupcji w każdym z nich jest zgodna z procesem Poissona przy ustalonym λ = 2,0348 . Robimy to, obliczając skumulowaną funkcję rozkładu Poissona:

Wybierając Δ T =2  kyr, z wyjątkiem nieco krótszych wartości na marginesach badanych przedziałów czasu, odkrywamy, że 2 z 21 segmentów leży poza 90% obszarem ufności oznaczonym przez P =0,05 i P =0,95 (pełne wyniki można znaleźć na rys. S3b). Te dwa segmenty występują w najmłodszej części zapisu i mogą być powiązane z wcześniej zgłoszonym wzrostem aktywności wulkanicznej podczas deglacjacji ( Zielinski i in. ,  1997 ) . Ostatnia analiza danych z rdzeni lodowych przeprowadzona przez LIN22 (patrz rys. 4 w nim) sugeruje, że wzrost aktywności podczas deglacjacji ogranicza się do erupcji z dużym osadem siarki na Grenlandii. Nie ustalono jednak, czy jest to statystycznie istotne, i nie możemy odrzucić możliwości, że po prostu odzwierciedla to większe szanse na wykrycie sygnałów wulkanicznych ze względu na lepsze zachowanie sygnału w młodszych częściach zapisu.

W przypadku rozważanych tu erupcji bipolarnych spodziewamy się 2,1 fałszywych wyników pozytywnych podczas testowania 21 niezależnych hipotez przy 90% ufności. Tak więc, nawet przy dwóch segmentach w kwestii, dane nie dostarczają dowodów na to, że tempo erupcji zmieniało się znacząco w czasie. Zamiast tego jest zgodne ze stacjonarnym procesem Poissona. Wykonaliśmy tę samą analizę, wybierając poszczególne stadiały i interstadiały jako segmenty danych zamiast segmentów przy ustalonym Δ T . Wyniki przedstawiono na rys. S4 i ponownie nie widzimy żadnych dowodów sugerujących, że tempo erupcji nie jest stałe w czasie lub nie jest stronnicze dla stadiów lub interstadiów. Tak więc, podczas gdy może występować wyższa częstotliwość erupcji w kierunku deglacjacji, w kontekście całego przedziału czasu 12–60 tys. lat temu mogło to nastąpić przypadkowo i nie jest uzasadnione uwzględnianie wyższej częstotliwości w naszej hipotezie zerowej. Niemniej jednak testujemy solidność naszych wyników w zestawieniu z wyższą częstotliwością erupcji (sekcja  3.5 ) i ponieważ stanowi to tylko krótki fragment całego okresu, potencjalnie zwiększona częstotliwość aktywności wulkanicznej podczas deglacjacji prawdopodobnie nie zmieni naszych wniosków.

Na koniec zauważamy, że brakujące luki w danych nie wpływają na przedstawione tutaj wyniki. Biorąc pod uwagę, że częstotliwość erupcji jest w przybliżeniu stała, nie ma potrzeby ciągłego okresu czasu, ponieważ do obliczeń wchodzi tylko tempo erupcji, a to ostatnie nie zmienia się znacząco przy wyborze różnych punktów, w których rekord jest odcinany. Bardziej konkretnie, ilość „pustej” przestrzeni bez zdarzeń DO, która pojawia się na granicach segmentów danych (lub luk), nie ma znaczenia w naszej analizie. Dzieje się tak, ponieważ pytamy, jak często zdarzeniu DO poprzedza erupcja, a nie jak prawdopodobne jest, że erupcja wywoła zdarzenie DO. Tylko w tym drugim przypadku ilość „pustej przestrzeni DO” na początku i na końcu przedziału danych miałaby wpływ na wyniki. Tak więc, zakładając, że statystyki erupcji i statystyczna relacja między zdarzeniami DO a erupcjami nie różniły się znacząco w luce danych, nasze wyniki nie są obciążone tym ostatnim.

2.5  Symulacje globalnego modelu oceanu

Wykorzystując globalny model oceanu, badamy, czy statystyczne powiązanie dużych erupcji wulkanicznych i cykli DO ujawnione tutaj można wyjaśnić przez wyzwolenie przejść atlantyckiej cyrkulacji południkowej (AMOC) z powodu krótkotrwałego ochłodzenia wulkanicznego. Te symulacje modelu będą służyć jako dowód zasady, a nie kompleksowe badanie wpływu wymuszenia wulkanicznego na zmienność atmosfery, lodu i oceanu DO w realistycznych warunkach granicznych lodowców. To ostatnie jest bardzo trudne, ponieważ ostatecznie główne czynniki napędzające zdarzenia DO pozostają tylko częściowo zrozumiane, a niewiele modeli jest w stanie odtworzyć zmienność klimatu podobną do cykli DO. Tak więc nasze symulacje jedynie ilustrują hipotezę, która zostanie przetestowana w dalszych realistycznych badaniach modelowania.

Tutaj używamy modelu oceanicznego Veros ( Häfner i in. ,  2018 ) , który jest modelem różnic skończonych równań pierwotnych skonfigurowanym w globalnej obecnej konfiguracji z określonym obecnym wymuszaniem atmosferycznym. Te obecne warunki różnią się znacząco od warunków brzegowych ostatniego okresu lodowcowego, a w konsekwencji siła i stabilność AMOC mogą być inne. Jednak prawdziwa siła i stabilność lodowcowego AMOC nie są dobrze ograniczone i znacznie się zmieniają w czasie. Zamiast próbować wyprowadzić spójne glacjalne warunki brzegowe, dla których wyłania się pożądany dynamiczny reżim AMOC, dostrajamy warunki brzegowe tak, aby system wykazywał rodzaj niestabilności AMOC, który jest powszechnie uważany za odpowiedzialny za przejścia między warunkami stadialnymi i międzystadialnymi podczas cykli DO. Konkretnie, używamy odpowiedniej zmiany warunku brzegowego zasolenia poprzez zwiększony dopływ słodkiej wody północnego Atlantyku (NA), w ramach którego model obsługuje wiele stabilnych stanów AMOC, tj. stany silnego, jak również zapadniętego przewracania. Stany te służą jako analogie do stanów stadialnych i międzystadialnych ostatniego zlodowacenia. Przechył między tymi stanami modelu został wcześniej zbadany przez Lohmanna i Ditlevsena ( 2021 ) . W porównaniu z ostatnim badaniem wprowadziliśmy drobne zmiany w konfiguracji modelu w celu poprawy realizmu i właściwości numerycznych. Po pierwsze, podczas gdy używamy tej samej rozdzielczości z 90 komórkami siatki podłużnej i 40 równoleżnikowej (gdzie rozdzielczość równoleżnikowa wzrasta z 5,3 ^∘ na biegunach do 2,1 ^∘ na równiku), a także 40 warstwami pionowymi (zwiększającymi grubość z 23 m na powierzchni do 274 m na dole), interpolacja batymetrii jest wykonywana przy użyciu filtra Gaussa o różnej szerokości pasma, co prowadzi do bardziej realistycznych mas lądowych. Po drugie, stosuje się współczesne wymuszanie naprężenia wiatru ERA-40. Wymiana ciepła i słodkiej wody z atmosferą jest modelowana przez warunki brzegowe, w których temperatura powierzchni morza i zasolenie są relaksowane w kierunku współczesnych pól klimatologicznych ( Uppala i in. ,  2005 ) w danej skali czasowej relaksacji. Wybór dłuższej skali czasowej dla zasolenia w porównaniu do relaksacji temperatury umożliwia dodatnie sprzężenie zwrotne adwekcji soli, które uważa się za odpowiedzialne za istnienie wielu stabilnych reżimów AMOC ( Weijer i in. ,  2019 ).. Tutaj używamy skali czasowej relaksacji wynoszącej 2 lata dla zasolenia i 30 dni dla temperatury. Należy zauważyć, że model nie obejmuje aktywnego lodu morskiego ani składnika atmosferycznego. Obecność tych składników może znacząco zmienić stabilność AMOC i jego reakcję na erupcje wulkaniczne. Nasz proponowany mechanizm jest zatem zależny od reakcji lodu morskiego i atmosfery, aby nie unieważnić bistabilności AMOC i gęstej formacji wodnej w północnym Atlantyku jako wulkanicznego wyzwalacza przejścia AMOC. Kolejną konsekwencją ram wyłącznie oceanicznych jest to, że nie są uwzględniane żadne szybkie procesy klimatyczne. Tak więc przejścia między różnymi stanami cyrkulacji nie są tak nagłe, jak można by się spodziewać ze względu na obecność szybkich, wzmacniających sprzężeń zwrotnych, takich jak nagłe zmiany zasięgu lodu morskiego ( Li i in. ,  2005 ) . Aby uzyskać więcej szczegółów na temat modelu oceanicznego, odsyłamy czytelnika do Lohmann i Ditlevsen ( 2021 ) oraz Häfner i in. ( 2018 ) .

Zaczynamy od przebiegu kontrolnego przy użyciu obecnych warunków początkowych i brzegowych. Biorąc pod uwagę grubą rozdzielczość, daje to w miarę realistyczną cyrkulację oceaniczną na dużą skalę, z maksymalną siłą AMOC około 18 Sv i maksymalną siłą północnoatlantyckiego i subpolarnego wiru odpowiednio 45 i 25 Sv (patrz rys. S9). Następnie stopniowo wprowadzamy wymuszanie słodkiej wody w NA, aby znaleźć punkt krytyczny odpowiadający załamaniu AMOC i ustanowić reżim bistabilności, w którym przejścia AMOC między alternatywnymi stanami stabilnymi mogłyby zostać wywołane przez wymuszanie wulkaniczne. Ze względu na warunki brzegowe relaksacji zasolenia, dodany stały strumień słodkiej wody (ujemny strumień zasolenia) jest równoważny zmianie pola wymuszającego zasolenie ϕ _i w komórce siatki i . W ten sposób stopniowo zmieniamy pole wymuszające w komórkach siatki między 296 a 0 ^∘  W i między 50 a 75 ^∘  N (patrz rys. S8). Odpowiada to obszarowi około A = 1,5  miliona km2 ^. Przyjmując stałe zasolenie odniesienia S _ref = 35  g kg ⁻¹ , zmianę w polu wymuszającym zasolenie można przekształcić w równoważne całkowite wymuszanie słodkowodne: (w Sv; 1 Sv  ≡  10 ⁶  m ³  s ⁻¹ ). Zwiększamy anomalię słodkowodną F w małych przyrostach, przy czym jeden przyrost składa się z 200-letniego liniowego wzrostu F , po którym następuje 100-letni okres relaksacji przy stałej wartości F . Przy wartości F = 0,49  Sv anomalia została zmniejszona z powrotem do 0 przy użyciu tych samych przyrostów. Całkowity czas symulacji tego eksperymentu przejściowego wyniósł 49 200 lat. Załamanie AMOC nastąpiło przy wyższej wartości F w porównaniu z odrodzeniem. Ujawnia to histerezę i bistabilność, które stanowią podstawę do przetestowania naszej hipotezy bezpośredniego wyzwalania przejść między współistniejącymi stanami AMOC jako odpowiedzi na ochładzanie wulkaniczne. Aby to zrobić, zaczynając od długich symulacji równowagi przy ustalonych wartościach F w obszarze bistabilności, przeprowadzamy symulacje, w których krótkotrwałe ochładzanie wulkaniczne jest wprowadzane przez tymczasową zmianę wartości warunku granicznego temperatury atmosferycznej. W przypadku symulacji przedstawionych tutaj, wszystkie erupcje zostały zainicjowane jesienią, ale wyniki nie są wrażliwe na porę roku erupcji (patrz rys. S11). Wybieramy strefowo jednorodne ochładzanie wulkaniczne, ale rozważamy różne profile południkowe, które emulują erupcje wulkaniczne na różnych szerokościach geograficznych. Więcej szczegółów podano w sekcji  3.7 .

3.1  Opóźnienia czasowe między początkiem zdarzenia DO a poprzedzającymi je erupcjami wulkanicznymi

Dokładne czasy początku ocieplenia DO oszacowaliśmy na podstawie indywidualnych, wysokorozdzielczych zapisów rdzeni lodowych Grenlandii δ ¹⁸ O, a także na podstawie stosu uzyskanego z nich. Ze względu na lepszy stosunek sygnału do szumu, ułożony zapis pozwala na najdokładniejsze określenie początku, ze średnią niepewnością czasu początku wynoszącą około 13 lat (sekcja  2.3 ). Dlatego skupimy się na omówieniu wyników uzyskanych ze stosu δ ¹⁸ O. Rysunek  1 przedstawia nasze szacunki czasów początku ocieplenia DO i wieku N = 82 bipolarnych erupcji z SVE20 wraz ze ułożonym zapisem Grenlandii δ ¹⁸ O. Istnieje znaczna liczba zdarzeń DO zainicjowanych w krótkim czasie po wybuchu wulkanu (patrz rys.  1 b–h). Opóźnienia czasowe między początkiem ocieplenia M = 20 DO a najbliższymi poprzedzającymi bipolarnymi erupcjami wulkanicznymi pokazano na rys.  2 a. W przypadku następujących siedmiu interstadiałów początek następuje w ciągu 50 lat po dwubiegunowej erupcji wulkanicznej: GI-1, GI-3, GI-4, GI-8, GI-9, GI-13 i GI-15.1. Jeśli zamiast tego przyjmiemy tolerancję 20 lat, znajdziemy pięć zdarzeń: GI-1, GI-3, GI-8, GI-13 i GI-15.1.

Na rys.  2 b pokazano opóźnienie najbliższej erupcji występującej po początku ocieplenia DO. Widać, że opóźnienia te są znacznie większe w porównaniu z erupcjami poprzedzającymi początki DO. Wskazuje to, że identyfikacja erupcji bipolarnych na krótko przed początkiem DO nie odzwierciedla jedynie faktu, że rdzenie lodowe są już dobrze zsynchronizowane blisko ociepleń DO w porównaniu z innymi miejscami w zapisie, co potencjalnie ułatwia znalezienie dopasowania bipolarnego z wymaganym poziomem pewności. Gdyby tak było i nie było związku przyczynowo-skutkowego między erupcjami bipolarnymi a ociepleniami DO, można by oczekiwać mniej więcej równego rozkładu opóźnień czasowych dla erupcji występujących przed i po początku DO. Można się nawet spodziewać, że opóźnienia czasowe dla erupcji po początku DO będą krótsze, ponieważ erupcje bipolarne powinny być łatwiejsze do zidentyfikowania ze względu na wyższą szybkość akumulacji międzystadialnej, a tym samym rozdzielczość proxy.

3.2  Wielkość erupcji i górna granica częstości występowania

Aby sprawdzić, czy jest to istotne, musimy oszacować liczbę erupcji, które miałyby miejsce na krótko przed początkiem DO przez przypadek, gdyby zdarzenia DO i erupcje bipolarne nie były skorelowane. Częstotliwość występowania dużych erupcji z SVE20 wynosi w przybliżeniu raz na 500 lat (patrz sekcja  2.4 ). Aby upewnić się, że nie jest to znaczące niedoszacowanie z powodu dużej liczby erupcji potencjalnie brakujących w zestawie danych, pokazujemy, że erupcje bipolarne miały w rzeczywistości wielkość, która odpowiada ich częstotliwości występowania raz na 500 lat, wykorzystując zapisy rdzeni lodowych obejmujące historyczne erupcje, a także ciągły zestaw danych dotyczących erupcji wulkanów lodowcowych z LIN22. Ponadto wyprowadzamy górny szacunek liczby równie dużych erupcji, które potencjalnie mogły mieć miejsce, a nie zostały zidentyfikowane w SVE20. Ten górny szacunek jest następnie później używany do testowania odporności naszych wyników na niedoszacowanie współczynnika występowania.

Podczas gdy ze względu na osadzanie się znacznych ilości siarczanów na obu biegunach jest oczywiste, że wszystkie zdarzenia w zestawie danych SVE20 były dużymi erupcjami, ich skala może się znacznie różnić. Niedawno Lin i współpracownicy zebrali zapis erupcji wulkanicznych z ciągłych pomiarów siarczanów w rdzeniach lodowych z Grenlandii i Antarktydy w przedziale 60–9 tys. lat (patrz sekcja  2.2 ) i obliczyli ich całkowitą depozycję siarczanów (LIN22). W swoim badaniu Lin i współpracownicy dokonali również ilościowej oceny erupcji bipolarnych z SVE20 jako podzbioru. W ten sposób możemy użyć zestawów danych LIN22 do porównania skali erupcji bipolarnych ze wszystkimi innymi erupcjami w tym samym przedziale czasowym, jak również w ciągu ostatnich 2500 lat ( Sigl i in. ,  2015 ) , które można wykryć w rdzeniach lodowych. Jak szczegółowo opisano w sekcji.  2.4 , wielkość pięciu największych erupcji w ciągu ostatnich 2500 lat odpowiada wielkości erupcji bipolarnych z SVE20, co potwierdza zaobserwowaną częstotliwość raz na 500 lat. Aby sprawdzić, czy nadal mogą występować erupcje o podobnej wielkości, które nie zostały zidentyfikowane jako bipolarne, w dalszej części porównano osady siarczanowe z erupcji bipolarnych z pozostałymi 740 (498) erupcjami wykrytymi w rdzeniach lodowych Grenlandii (Antarktydy) (zwanych dalej unipolarnymi), których nie można było zidentyfikować jako erupcji bipolarnych w SVE20.

Na obu biegunach rozkład depozycji siarczanów odpowiadających erupcjom bipolarnym, jak pokazano na histogramach na rys.  3 a–b, jest wyraźnie przechylony w stronę większych wartości w porównaniu do pozostałych (unipolarnych) erupcji. W populacji erupcji bipolarnych, siedem erupcji mających miejsce w ciągu 50 lat przed wystąpieniem DO ma średnio nieco mniejsze depozycje, ale nie jest to statystycznie istotne. W przypadku Grenlandii średnie depozycje zdarzeń przed wystąpieniem DO wynoszą 107,9 kg km ^-2 w porównaniu do 157,3 kg km ^-2 dla całej populacji bipolarnej. Test hipotezy bootstrap dla średniej pokazuje, że ta różnica w średniej nie jest istotna ( p =0,16 ). W przypadku Antarktydy średnia wynosi 51,1 kg km ^-2 w porównaniu do 63,0 kg km ^-2 dla populacji bipolarnej ( p =0,28 ). Na podstawie względnych depozycji na Grenlandii i Antarktydzie, LIN22 wyprowadził również całkowity ładunek aerozolu stratosferycznego z erupcji bipolarnych, który determinuje ich globalne wymuszanie klimatyczne. Rozkład tego ładunku aerozolu pokazano na rys.  3 c. Średnia dla erupcji przed wystąpieniem DO (126,2 Tg) jest ponownie nieznacznie niższa niż średnia dla populacji (159,9 Tg), ale nieznacznie ( p =0,21 ). Tak więc pod względem wielkości erupcje przed wystąpieniem DO można uznać za reprezentatywną próbkę populacji bipolarnej i nie ma dowodów na stronniczość systematycznego dopasowywania mniejszych erupcji bipolarnych w pobliżu wystąpień DO.

Podczas gdy erupcje w zestawach danych unipolarnych wykazują wyraźnie niższe średnie poziomy depozycji zarówno na Grenlandii, jak i na Antarktydzie w porównaniu z erupcjami bipolarnymi, podzbiory erupcji pozostają w zestawach danych unipolarnych z poziomami depozycji, które pasują do populacji bipolarnej. Jest prawdopodobne, że niektóre z nich są w rzeczywistości bipolarne, chociaż nie było to możliwe do ustalenia na podstawie obecnych zapisów rdzeni lodowych i chronologii. W ten sposób wyprowadzamy konserwatywny szacunek górnej granicy liczby erupcji bipolarnych o tej samej charakterystycznej wielkości, które mogły być niezidentyfikowane. Mówiąc ogólnie, należy zidentyfikować maksymalną liczbę erupcji w zestawie danych unipolarnych z depozycją siarczanów porównywalną z populacją bipolarną, co można zrobić na kilka sposobów. Proste oszacowanie uzyskuje się, znajdując podzbiór największych erupcji w dowolnym z zestawów danych unipolarnych, który ma średnią depozycję równą siedmiu erupcjom, które miały miejsce w ciągu 50 lat przed zdarzeniami DO (dalej populacja DO). Z zestawu danych z Antarktydy obejmującego łącznie 498 erupcji wynika, że ​​podzbiór 84 największych erupcji ma średnią depozycję odpowiadającą średniej populacji DO. W przypadku większego zestawu danych z Grenlandii (łącznie 740 erupcji) podzbiór ten jest znacznie większy, a 212 największych erupcji ma średnią depozycję odpowiadającą średniej populacji DO.

Aby uzyskać pojedynczą uczciwą ocenę liczby potencjalnie niezidentyfikowanych dużych erupcji bipolarnych, należy wybrać mniejszą liczbę z dwóch. Wynika to z faktu, że w przypadku erupcji bipolarnej o rozmiarze rozważanym w naszym badaniu, wymagane jest wystarczająco duże osadzanie na obu biegunach, a zatem podzbiór musi spełniać wymaganą średnią depozycję na obu biegunach, aby być porównywalnym z erupcjami bipolarnymi z SVE20. Jeśli weźmiemy pod uwagę 212 największych erupcji w zestawie danych Antarktydy, otrzymamy podzbiór ze średnią depozycją, która jest znacznie mniejsza niż średnia populacji DO. Zatem nie może być tak wielu dodatkowych erupcji w zapisie rdzenia lodowego, które mają charakterystyczną wielkość porównywalną z danymi z SVE20. Zamiast tego rozbieżność między zestawami danych z Grenlandii i Antarktydy oznacza, że ​​zestaw danych z Grenlandii zawiera dużą liczbę regionalnych (np. islandzkich) erupcji z dużym grenlandzkim osadem ze względu na bliskość źródła, ale bardzo małym lub nieobecnym osadem antarktycznym, a zatem globalnym wpływem, który nie jest porównywalny z bipolarnymi erupcjami z SVE20. Potwierdza to analiza wulkanizmu w rdzeniach lodowych z ostatnich 2000 lat, opracowana przez Sigl i in. ( 2013 ) , w której rdzenie lodowe Grenlandii i Antarktydy są bardzo dobrze zsynchronizowane, a sygnały wulkaniczne są bardzo dobrze zachowane ze względu na dużą grubość warstwy. Prowadzi to do bardzo małej liczby potencjalnie niezidentyfikowanych bipolarnych erupcji. Biorąc pod uwagę cztery erupcje z największym grenlandzkim osadem siarczanowym w tym okresie, trzy erupcje nie mają żadnego powiązanego osadu antarktycznego i w związku z tym nie są klasyfikowane jako bipolarne. Dwa z nich mają znane pochodzenie islandzkie, przy czym erupcja Laki z 1783 r. n.e. jest największą ze wszystkich. Pozostała z czterech erupcji ma depozycję antarktyczną mniejszą niż 10 kg km ^-2 , co jest progiem użytym w zestawie danych LIN22. Te cztery erupcje mają miejsce w ciągu 1200 lat, co wskazuje na częstotliwość dużych grenlandzkich erupcji unipolarnych, która jest jeszcze wyższa w porównaniu z erupcjami bipolarnymi w zestawie danych SVE20. Należy zauważyć, że nie oznacza to, że takie erupcje nie miały wpływu na klimat. Na przykład wpływ klimatyczny erupcji Laki z 1783 r. n.e. jest nadal aktywnie badany ( Zambri i in. ,  2019 ; Edwards i in. ,  2021 ) , co jest trudnym zadaniem w przypadku pojedynczych erupcji ze względu na naturalną zmienność klimatu. Jednakże wyraźne globalne ochłodzenie wulkaniczne wydaje się mniej prawdopodobne w przypadku braku osadzania się siarki na Antarktydzie, a w przypadku wybuchu Laki argumentowano, że duże osadzanie się siarczanów na Grenlandii jest w dużej mierze wynikiem transportu w troposferze lub niższej stratosferze (Lanciki i in. ,  2012 ) . W przeciwieństwie do dużych unipolarnych erupcji na Grenlandii, cztery erupcje z największymi depozycjami siarki na Antarktydzie, do których należą erupcje Tambora z 1815 r. n.e. i Samalas z 1257 r. n.e., mają duże towarzyszące im depozycje grenlandzkie.

To dodatkowo potwierdza wykorzystanie zestawu danych antarktycznych do ograniczenia liczby potencjalnie brakujących erupcji bipolarnych, oprócz wyżej wymienionego wymogu metodologicznego dużej depozycji na obu biegunach. W rezultacie, postępując zgodnie z tą metodologią, znajdujemy górną granicę szacunkową 84 potencjalnie niezidentyfikowanych erupcji bipolarnych. Na podstawie tego szacunku, zestaw danych bipolarnych SVE20 używany w naszym badaniu może nie obejmować 50% erupcji o podobnej wielkości w maksymalnym zakresie. Nazywamy to górną granicą, ponieważ, podczas gdy pozornie mniej wyraźne w porównaniu do zestawu danych z Grenlandii, zestaw danych z Antarktydy zawiera również regionalne erupcje z silnym depozycją antarktyczną, ale słabym lub nieobecnym depozycją grenlandzką i globalnym wpływem na klimat. Jest to częściowo spowodowane lokalnym wulkanizmem antarktycznym, który można wykryć w rdzeniach lodowych Antarktydy w całym ostatnim zlodowaceniu ( Narcisi i in. ,  2017 ) .

Oszacowanie jest ponadto ostrożne, ponieważ rozkład depozycji siarczanów tego potencjalnego podzbioru 84 erupcji jest silnie skośny. Dzieje się tak, ponieważ ogon rozkładu unipolarnego zanika znacznie szybciej niż rozkład bipolarny (patrz rys.  3 b). Mniejsze zdarzenia depozycji są zatem reprezentowane częściej niż w populacji bipolarnej. Ponieważ nasza analiza opiera się wyłącznie na populacji zdarzeń w ramach danego rozkładu wielkości, bardziej sprawiedliwe może być zadanie następującego pytania: jaka jest maksymalna liczba zdarzeń bipolarnych, które mogłaby zawierać próbka unipolarna i które podążają za rozkładem wielkości znanej populacji bipolarnej? Ponieważ taki podzbiór zdarzeń będzie pochodził głównie z ogona rozkładu unipolarnego, a zdarzenia w ogonie są stosunkowo rzadkie, przyjęto podejście probabilistyczne w celu uzyskania niepewności w oszacowaniu górnej granicy. Metoda jest szczegółowo opisana w Załączniku B i wykorzystuje losowe próbkowanie ze zwracaniem w celu wielokrotnego generowania maksymalnych podzbiorów z populacji unipolarnej, które podążają za rozkładem depozycji bipolarnej. Rys.  4a pokazuje, że te wybrane podzbiory rzeczywiście podążają za tym rozkładem z dobrym przybliżeniem. Na rys.  4b pokazujemy rozkład liczby zdarzeń w podzbiorach. Średnia liczba zdarzeń w zestawie danych antarktycznych wynosi 37,5, z przedziałem ufności [20,0, 55,0] opartym na 10. i 90. percentylu. Co ciekawe, używając zamiast tego zestawu danych z Grenlandii, znajdujemy nieco mniejszą liczbę 32,8 erupcji średnio (patrz rys. S7) w przeciwieństwie do znacznie większej liczby erupcji ze średnią depozycją populacji DO podaną powyżej. Dzieje się tak, ponieważ rozkład depozyów z Grenlandii w zestawie danych bipolarnych jest znacznie bardziej skośny w porównaniu ze zbiorem danych z Antarktydy (porównaj rys.  4a i S7a). Zatem średnia jest zdominowana przez największe erupcje, a podzbiór skonstruowany na podstawie średniej będzie zawierał znacznie więcej małych erupcji, aby zrównoważyć duże. W ten sposób stwierdzamy, że zgodnie z oczekiwaniami uwzględnienie obserwowanego rozkładu depozycji odpowiadającej wydarzeniom bipolarnym zamiast dopasowania tylko do średniej depozycji prowadzi do mniejszego szacunku. Jednakże wykorzystamy większą górną granicę 84 brakujących erupcji, aby przetestować solidność naszych wyników.

Nieco naiwne podejście polegające na wykorzystaniu jednobiegunowych osadów do wnioskowania o wielkości potencjalnych erupcji bipolarnych ostatecznie nie może ujawnić prawdziwej częstotliwości erupcji z bipolarnym odciskiem i dużym globalnym wymuszaniem klimatycznym. Dlatego chcemy podkreślić, że liczby podane powyżej nie są naszymi najlepszymi szacunkami rzeczywistej liczby niezidentyfikowanych erupcji bipolarnych, ale raczej ich górną granicą, zgodnie z którą zakładamy, że większość dużych zdarzeń depozycji na poszczególnych biegunach to rzeczywiście erupcje z bipolarnym odciskiem. Ponadto liczby podane tutaj, tj. okres powrotu 500 lat, jak również szacunek erupcji potencjalnie brakujących w zestawie danych SVE20, nie odnoszą się do erupcji bipolarnych dowolnej wielkości, ale do erupcji bipolarnych o charakterystycznym (dużym) rozmiarze erupcji bipolarnych w danych SVE20. Wiadomo, że te pierwsze występują znacznie częściej ( Sigl i in. ,  2022 ) . Rzeczywista częstotliwość dużych erupcji bipolarnych pozostaje kwestią otwartą, dopóki nie pojawią się nowe dane. Dla naszego celu ograniczenia maksymalnej częstotliwości zdarzeń w ramach obecnie dostępnych zestawów danych szacunki podane tutaj są wystarczające.

3.3  Porównanie danych z hipotezą zerową o losowo występujących erupcjach wulkanicznych

Biorąc pod uwagę liczbę erupcji bliskich początkom DO znalezioną w sekcji  3.1 , zadajemy teraz pytanie, czy jest prawdopodobne, że zaobserwowano to przypadkowo w przypadku, gdy bipolarne erupcje wulkaniczne wystąpiły całkowicie losowo i nieskorelowane z występowaniem zdarzeń DO. Aby odpowiedzieć na to pytanie, testujemy hipotezę zerową, że erupcje wulkaniczne występują jako stacjonarny proces Poissona niezależnie od zdarzeń DO; tj. nie występują one preferencyjnie krótko przed początkami DO. Opis erupcji wulkanicznych jako stacjonarnego procesu Poissona oznacza, że ​​występują one ze stałą szybkością w czasie i że kolejne erupcje nie są skorelowane ( De la Cruz-Reyna ,  1991 ) . Nasze dane dotyczące erupcji z N = 82 i Δ T = 40,3  tys. lat są zgodne z tym założeniem, dając λ = 2,03 erupcji na tysiąc lat (sekcja  2.4 ). Można również założyć, że czasy oczekiwania między kolejnymi wybuchami DO nie są skorelowane. Dwustronny test hipotezy bootstrapowej na korelacji Pearsona (Spearmana) kolejnych czasów oczekiwania daje p = 0,60 ( p = 0,32 ) dla korelacji r = 0,13 ( r = 0,26 ). W tym przedziale glacjalnym, który obejmuje Marine Isotope Stage 3 z dość regularnymi cyklami DO, empiryczny rozkład czasów oczekiwania nie pasuje idealnie do funkcji wykładniczej. Test Andersona–Darlinga (Kolmogorowa–Smirnowa) prawie odrzuca rozkład wykładniczy przy 95% ufności z p = 0,060 ( p = 0,051 ). Jednak wynikające z tego potencjalne odchylenie wybuchów DO od procesu Poissona nie jest istotne, ponieważ średni odstęp między wybuchami jest duży w porównaniu z oknem tolerancji i odstępem między erupcjami. Zatem zbiór danych jest odpowiedni dla naszej hipotezy zerowej.

Zgodnie z hipotezą zerową oczekiwana wartość opóźnienia czasowego niezależnie występującego zdarzenia DO do najbliższej poprzedzającej erupcji bipolarnej wynosi  lat. Kwantyle opóźnienia czasowego wynoszą . Tak więc istnieje 5% (10%) szansy na zaobserwowanie opóźnienia czasowego mniejszego lub równego 25,2 (51,8) lat. Linie poziome na rys.  2a pokazują, że w obu przypadkach istnieje znaczna liczba zdarzeń z mniejszym opóźnieniem czasowym. Aby pokazać, że jest to rzeczywiście istotne, powiemy w dalszej części, że istnieje dopasowanie zdarzenia, jeśli erupcja wulkaniczna następuje w oknie opóźnienia czasowego τ lat przed wystąpieniem DO, które nazywamy tolerancją. Liczbę dopasowań zdarzeń w danych dla tolerancji τ należy porównać z liczbą razy, w których jedno lub więcej zdarzeń wulkanicznych zostałoby przypadkowo znalezione w oknie czasowym τ lat podczas losowego próbkowania M okien procesu Poissona o szybkości λ . W tym przypadku liczba zdarzeń n występujących w ciągu τ lat jest podana przez rozkład Poissona. Prawdopodobieństwo zaobserwowania jednego lub więcej zdarzeń w ciągu τ lat wynosi

Losowe wybieranie M niezależnych okien i obserwowanie, czy zawierają one jakiekolwiek zdarzenia, jest równoważne sekwencji prób Bernoulliego z prawdopodobieństwem sukcesu P ( n ( τ )≥1) . Tak więc prawdopodobieństwo znalezienia jednego lub więcej zdarzeń w k spośród M okien jest podane przez rozkład dwumianowy:

Wreszcie prawdopodobieństwo znalezienia co najmniej k z M okien zawierających jedno lub więcej zdarzeń reprezentuje wartość p dla naszej hipotezy zerowej. Jest to skumulowane prawdopodobieństwo poprzedniego wyrażenia:

Dla k =7 dopasowań znalezionych przy τ =50 lat daje to prawdopodobieństwo p =0,002 . Tak więc możemy odrzucić hipotezę zerową przy poziomie ufności 99%. Dla k =5 dopasowań znalezionych przy τ =20  lat otrzymujemy p =0,0009 . Wartość oczekiwana dopasowań wynosi , co daje tylko 1,9 (0,8) zdarzeń poprzedzonych jedną lub większą liczbą erupcji w ciągu 50 (20) lat zgodnie z hipotezą zerową. Natomiast opóźnienia czasowe erupcji następujących po początkach DO (rys.  2 b) dobrze odpowiadają tej hipotezie zerowej.

Pobierać

3.4  Odporność wyników na zmienne wybory maksymalnego opóźnienia czasowego

Nasza pewność odrzucenia hipotezy zerowej może zależeć od tolerancji τ . Aby sprawdzić, czy wyniki są solidne, bierzemy pod uwagę prawdopodobieństwo dwumianowe dla prawdopodobnego zakresu tolerancji. Tolerancje niższe niż 10–15 lat nie powinny być brane pod uwagę, ponieważ są tego samego rzędu, co szacowana niepewność w naszych czasach wystąpienia DO. Ponadto, przy tych tolerancjach dopasowania zdarzeń w danych stają się zbyt rzadkie, aby mogły być wiarygodnym szacunkiem rzeczywistego prawdopodobieństwa współwystępowania zdarzeń DO i erupcji bipolarnych. W przypadku tolerancji wynoszących sto lat lub więcej, bezpośredni wpływ klimatyczny nawet dużych erupcji wulkanicznych staje się mniej prawdopodobny. Rysunek  5a przedstawia liczbę dopasowań jako funkcję τ . Ciągła szara linia wskazuje oczekiwaną wartość dopasowań przy hipotezie zerowej. Szare (żółte) cieniowanie wskazuje odpowiadające im pasma ufności 90% (95%). W przypadku tolerancji dłuższych niż 13 lat dane stale znajdują się powyżej pasma ufności 95% i dlatego jest mało prawdopodobne, aby wystąpiły przy hipotezie zerowej. Dokładne prawdopodobieństwa obliczone z równania ( 4 ) podano na rys.  5 b. Odpowiadają one wartościom p hipotezy zerowej jako funkcji τ i widzimy, że dane leżą poniżej p = 0,01 dla wszystkich tolerancji większych niż 13 lat. Wyniki są również istotne, gdy weźmiemy pod uwagę czasy początku poszczególnych rdzeni lodowych (kolorowe krzywe). Ponadto wyniki analityczne dla prawdopodobieństwa zdarzeń występujących zgodnie z hipotezą zerową podaną w tej sekcji i poprzedniej są dokładne i odporne na (niewielkie) odchylenia naszych zestawów danych od założeń leżących u podstaw hipotezy zerowej. Pokazano to na rys. S6, gdzie wyniki analityczne porównano z bezpośrednimi symulacjami Monte Carlo, w których faktycznie obserwowane dane dotyczące początku DO i erupcji wulkanicznej są używane jawnie.

3.5  Odporność wyników na potencjalne niedoreprezentowanie erupcji dwubiegunowych

Co najważniejsze, pewność, z jaką możemy odrzucić hipotezę zerową, zależy od oszacowania częstości występowania erupcji. Jeśli zbiór danych SVE20 jest niekompletny, częstość występowania erupcji λ jest niedoszacowana. Zostało to szczegółowo omówione w sekcji  3.2 , a my znaleźliśmy 84 erupcje jako konserwatywny szacunek górnej granicy maksymalnej liczby potencjalnie niezidentyfikowanych erupcji. Po zwiększeniu N z 82 do 166, wartość oczekiwana prawie się podwaja, jak widać linią przerywaną na rys.  5a . Mimo to liczba dopasowań w danych pozostaje znacząca: Cieniowanie na rys.  5b pokazuje zakres prawdopodobieństw uzyskanych przy rozważaniu częstości poprzez zwiększenie N z 82 do 166. Nawet wtedy konsekwentnie znajdujemy p ≤0,1 , co pozwala nam odrzucić hipotezę zerową przy 90% pewności. Zatem nasze wyniki są odporne na potencjalne niedoszacowanie liczby dużych erupcji bipolarnych, biorąc pod uwagę, że połowa erupcji, które nie miały miejsca w pobliżu zdarzenia DO, nie została zidentyfikowana.

3.6  Kontrastujący wpływ wulkanów na fazy ocieplenia i ochłodzenia cykli DO

Może się wydawać sprzeczne z intuicją, że radiacyjne chłodzenie aerozolem erupcji wulkanicznych może wywołać nagłe przejścia ocieplenia cykli DO. Wpływ na nagłe fazy chłodzenia cykli DO może wydawać się na początku bardziej prawdopodobny. Wykonaliśmy tę samą analizę przejść chłodzenia, które oznaczają zakończenie interstadiów i nie znajdujemy znaczącego skupiska erupcji prowadzących do ochłodzeń (patrz rys. S5). Zauważamy jednak, że przejścia chłodzenia w zapisach izotopowych są na ogół znacznie mniej dobrze zdefiniowane. Podstawa statystyczna do oceny istotności tego nie jest zatem tak solidna w porównaniu z przypadkiem ociepleń DO. W rezultacie kwestia potencjalnego wulkanicznego wyzwalacza przejść chłodzenia DO pozostaje niejednoznaczna, jak omówiono poniżej.

Po pierwsze, około połowa interstadiów nie ma dobrze zdefiniowanego nagłego ochłodzenia, i można by zatem argumentować, że nie ma sensu szukać wyzwalacza zakończeń tych interstadiów. W Lohmann i Ditlevsen ( 2019 ) podano i przetestowano kryteria dotyczące kształtu piły cykli DO, tj. istnienia nagłego ochłodzenia, które jest odrębne od stopniowego ochłodzenia podczas interstadiów. Spośród 19 interstadiów badanych tutaj, 10 nie wykazuje wyraźnego nagłego ochłodzenia, zgodnie z Lohmann i Ditlevsen ( 2019 ) . Przy ograniczeniu próby do dziewięciu zdarzeń z wyraźnym nagłym ochłodzeniem, najbliższa odległość poprzedniej erupcji wulkanicznej wynosi 170 lat. Biorąc pod uwagę, że 10. percentyl pozostaje 51,8 roku (sekcja 3.3), jest to wyraźnie nieistotne statystycznie. Spośród wszystkich 19 przejść chłodzenia znajdujemy tylko jedno chłodzenie DO, dla którego poprzedzająca erupcja następuje w ciągu 100 lat. Na podstawie tych szacunków czasowych wydaje się niemożliwe znalezienie podzbioru zdarzeń, dla których można ustalić statystycznie istotny związek między zakończeniami interstadiów i erupcjami wulkanicznymi.

Jednak po drugie, w porównaniu do ociepleń DO istnieją większe niepewności w naszych szacunkowych czasach ochłodzeń DO, które zostały tutaj określone inną metodą wprowadzoną przez Lohmann i Ditlevsen ( 2019 ) w porównaniu do nagłych ociepleń. Jest to widoczne po stosunkowo dużym rozrzucie międzyrdzeniowym na rys. S5a i dotyczy to również interstadiów z wyraźnym nagłym ochłodzeniem. W rezultacie mogliśmy przegapić pojedyncze dopasowania zdarzeń, takie jak wybuch wulkanu występujący blisko zakończenia GI-16.2. Krótkie interstadiały GI-16.2 i GI-17.2 (w których znajdujemy wybuch wulkanu tuż przed zakończeniem) są wartościami odstającymi ( Lohmann i Ditlevsen ,  2019 ) . Tak więc wpływ wulkanu, który zakończy je przedwcześnie, może być prawdopodobny. Po trzecie, niektóre duże zdarzenia bipolarne występujące w pobliżu ochłodzeń mogły zostać przeoczone przez dane SVE20, ponieważ wcześniejsza synchronizacja rdzeni lodowych Grenlandii i Antarktydy nie jest tak dobra w pobliżu ochłodzeń DO, jak w przypadku ociepleń DO.

Mając na uwadze te niepewności co do istnienia i czasu nagłych ochłodzeń, obecne dane bipolarne sugerują, że takie wpływy wulkaniczne przed ochłodzeniami DO nie występowały częściej, niż można by oczekiwać przypadkowo. Jednak w naszym badaniu nie bierzemy pod uwagę potencjalnego wpływu dużych erupcji pozatropikalnych bez śladu bipolarnego, które mogą również wywierać znaczący długotrwały wpływ na klimat NH ( Kobashi i in. ,  2017 ; Toohey i in. ,  2019 ; van Dijk i in. ,  2022 ) . Aby wyjaśnić wpływ chłodzenia wulkanicznego na ocieplenie lub przejścia w ochłodzenie DO, należy rozważyć jego wpływ na główne części systemu klimatycznego, o których uważa się, że biorą udział w cyklach DO. W dalszej części badamy bardzo prosty mechanizm, który uwzględnia regionalny wpływ chłodzenia na AMOC przy użyciu globalnego modelu oceanu.

3.7  Eksperymenty z globalnym modelem oceanu w celu zbadania wpływu chłodzenia wulkanicznego na cykle DO

Przejścia między energicznym a zapadniętym lub bardzo słabym stanem cyrkulacji AMOC są powszechnie sugerowane jako wyjaśnienie znanych wzorców czasoprzestrzennych cykli DO. Czy przejścia nastąpiły jako autonomiczne samopodtrzymujące się cykle, w wyniku stopniowych zmian warunków tła klimatycznego, czy przynajmniej częściowo przez wymuszenie stochastyczne, nie jest obecnie wiadome. Niezależnie od długoterminowej dynamiki, mechanizmy te mają to wspólnego, że ze względu na nieliniowość alternatywne stabilne stany AMOC współistnieją dla szeregu warunków tła klimatycznego. Przejścia mogą wystąpić albo wtedy, gdy warunki tła zmieniają się tak bardzo, że jeden ze stanów traci stabilność (w punkcie krytycznym), albo gdy system jest szybko zaburzany przez wystarczająco duży szok. To drugie staje się bardziej prawdopodobne, im bliżej punktu krytycznego.

Siła AMOC zależy od konwekcji w NA, gdzie gęsta woda tworzy się na powierzchni, opada w kierunku dna, a następnie płynie na południe. Gęstość powierzchni NA jest zależna od temperatury atmosfery, a także od budżetu zasolenia poprzez wymianę słodkiej wody. Osłabienie i ostateczne załamanie AMOC może nastąpić, jeśli do NA dostanie się zbyt dużo słodkiej wody, np. z topnienia Grenlandii. Na rys.  6a pokazujemy reakcję siły AMOC w globalnym modelu oceanicznym Veros, w którym wymuszanie słodkiej wody jest wprowadzane do NA, powoli zwiększane i ostatecznie ponownie zmniejszane (więcej szczegółów w sekcji  2.5 ). Ujawnia to reżim parametru wymuszania słodkiej wody F ze współistniejącymi stanami silnego i załamanego AMOC, które są analogiczne do warunków AMOC podczas ostatniego okresu lodowcowego: okresy stadialne były na ogół związane z załamanym AMOC, podczas gdy w interstadiałach stwierdzono silne AMOC ( Henry i in. ,  2016 ) . Stabilność stanów modelu potwierdzono długimi symulacjami równowagi (duże okręgi na rys.  6a ). Patrz rys. S9 dla ilustracji odpowiadających funkcji strumienia AMOC dla stanu zapadniętej i energicznej cyrkulacji przy tym samym wymuszeniu F i rys. S10 dla szeregów czasowych maksymalnej siły AMOC symulacji równowagi. Zakres bistabilności (około 0,02 Sv) jest dość wąski w porównaniu z innymi modelami wyłącznie oceanicznymi i modelami klimatycznymi o średniej złożoności, a w przeciwieństwie do wielu innych modeli obecny klimat nie znajduje się w reżimie bistabilnym ( Rahmstorf i in. ,  2005 ) . Jednak zakres bistabilności w tych ramach wyłącznie oceanicznych zależy od parametrów, które są w dużej mierze nieograniczone, takich jak skale czasowe relaksacji temperatury i zasolenia ( Lohmann i Ditlevsen ,  2021 ) . Szerokość i lokalizacja reżimu bistabilności nie wpływają istotnie na możliwość przetestowania naszej hipotezy.

Biorąc pod uwagę, że istnieją dwa alternatywne stany stabilne, bezpośredni wpływ termiczny chłodzenia wulkanicznego na AMOC może wywołać przejście z zapadniętego do energicznego AMOC (początek ocieplenia DO), ale nie odwrotnie. Konkretnie, chłodzenie atmosferyczne NA prowadzi do zwiększonej utraty ciepła przez ocean, a tym samym do tworzenia się gęstej wody, co może zainicjować lokalną konwekcję. W ten sposób może zostać wywołane spontaniczne odrodzenie zapadniętego AMOC, podczas gdy odwrotnie, spontaniczne zapadnięcie się energicznego AMOC staje się mniej prawdopodobne. Ponieważ zaburzenie jest stosunkowo krótkotrwałe, spontaniczne przejście byłoby prawdopodobne tylko wtedy, gdyby system był już blisko punktu krytycznego.

Testujemy to w modelu za pomocą dwóch zestawów symulacji zainicjowanych odpowiednio w stanie energicznym i zapadniętym AMOC. Warunki tła w odniesieniu do wymuszania słodkiej wody NA są dobierane tak, aby system był blisko punktu krytycznego w obu scenariuszach. Odpowiednie stany są oznaczone zieloną gwiazdką na rys.  6 a. Symulacje rozpoczęte w stanie częściowo zapadniętym przy F = 0,360 zostały również przeprowadzone (rys. S13b). Zaczynając od długich symulacji spin-up ze stałymi warunkami brzegowymi, wprowadzamy globalne ochłodzenie wulkaniczne, które osiąga szczyt 1 rok po erupcji i stopniowo zanika w ciągu 10 lat (patrz rys.  6 b). Rozważamy dwa różne wzorce przestrzenne z wyraźniejszym chłodzeniem NH i SH, odpowiednio (patrz rys.  6 c–d). Pierwszy scenariusz mógłby reprezentować erupcje północno-zachodnie i równikowe, które wykazały wyraźniejsze ochłodzenie północno-zachodnie w symulacjach modelowych ( Schneider i in. ,  2009 ; Pausata i in. ,  2020 ; Black i in. ,  2021 ; Zhuo i in. ,  2021 ) . Scenariusz z silniejszym ochłodzeniem północno-wschodnich półkul odpowiadałby raczej pozatropikalnej erupcji północno-zachodnich półkul.

Rysunek  7 pokazuje, że przejście AMOC można zaobserwować tylko w przypadku erupcji NH–równikowej rozpoczynającej się od stanu zapadniętego (zielona trajektoria na rys.  7 a). Jest to odporne na zmiany pory roku erupcji, wielkości i warunków początkowych (rys. S11 i S12). Erupcja SH nie prowadzi do ponownego wzrostu AMOC, prawdopodobnie z powodu niewystarczającego chłodzenia NA. Obie symulacje zainicjowane w stanie energicznego AMOC pokazują nagłe wzmocnienie AMOC, po którym następują tłumione oscylacje powracające do stanu pierwotnego. Zauważamy, że przejście na rys.  7 a nie jest bardzo nagłe, w przeciwieństwie do obserwowanych początków DO. Podczas gdy w wyniku początku formowania się NA głębokowodnej maksymalna siła AMOC początkowo gwałtownie wzrasta po zaburzeniu wulkanicznym, początkowa silna konwekcja NA zostaje zastąpiona bardziej stopniowym poszerzaniem się przewracającej się komórki po kilku dekadach. W terminologii układów dynamicznych początkowe zaburzenie AMOC wypycha układ z zapadniętego stanu stabilnego nie bezpośrednio do energicznego stanu stabilnego, ale do jego basenu przyciągania. Przejście w obrębie tego basenu przyciągania jest następnie bardzo stopniowe i może wykazywać przejściowe spowolnienie dynamiki, gdy odwiedzane jest sąsiedztwo jednego lub więcej stanów siodłowych, co naszym zdaniem jest specyficzne dla naszego modelu i jego konfiguracji (np. konfiguracja w Lohmann i Ditlevsen ,  2021 , zachowywała się znacznie bardziej gwałtownie). W terminologii fizycznej powolne przejście może również wskazywać na fakt, że model nie obejmuje aktywnych komponentów o szybkich skalach czasowych, tj. dynamicznego lodu morskiego i dynamiki atmosferycznej. Uważa się, że w szczególności dynamika lodu morskiego jest zaangażowana w nagłość zdarzeń DO ( Li i in. ,  2005 ; Dokken i in. ,  2013 ) .

Nadal uważamy, że symulacje modelu pomagają zilustrować potencjalne wyjaśnienie naszej obserwacji, że przejścia ocieplenia DO wydają się być związane z poprzedzającymi je erupcjami, podczas gdy w zestawie danych użytym tutaj nie mogliśmy znaleźć podobnych dowodów w odniesieniu do przejść chłodzenia DO. Obserwacja, że ​​tylko niektóre początki DO są poprzedzone erupcjami, jest zgodna z ideą, że erupcja musi nastąpić w momencie, gdy system klimatyczny jest już blisko punktu krytycznego, co jest również wspierane przez nasze symulacje modelu (patrz rys. S13a dla symulacji zainicjowanych dalej od punktu krytycznego).

Wiadomo, że szerokość geograficzna i półkula erupcji wpływają na jej wpływ klimatyczny ( Sun i in. ,  2019 ; Zhuo i in. ,  2021 ) . Jeśli rzeczywiście do wywołania ocieplenia DO wymagane jest duże ochłodzenie NA, można by argumentować, że nie wszystkie rozważane tutaj erupcje bipolarne powinny zostać uwzględnione w naszej analizie. LIN22 sklasyfikował erupcje bipolarne do dwóch pasm szerokości geograficznej, biorąc pod uwagę względne osadzanie na Grenlandii w porównaniu z Antarktydą, z następującymi wynikami. Spośród 82 erupcji bipolarnych 34 jest klasyfikowanych jako niskie szerokości geograficzne lub półkula południowa (niskie szerokości geograficzne poniżej 40 ^∘  N lub półkula południowa, LLSH), podczas gdy 48 jest klasyfikowanych jako wysokie szerokości geograficzne półkuli północnej (powyżej 40 ^∘  N, NHHL). Jeśli chodzi o erupcje zidentyfikowane tutaj jako występujące krótko przed początkiem DO, klasyfikacja daje pięć (LLSH) w porównaniu z dwoma (NHHL) erupcjami występującymi w ciągu 50 lat i cztery erupcje (LLSH) w porównaniu z jedną erupcją (NHHL) w ciągu 20 lat. Może się to wydawać sprzeczne z naszym proponowanym mechanizmem wyzwalającym poprzez chłodzenie NA. Należy jednak zauważyć, że klasyfikacja LLSH obejmuje wszystkie erupcje tropikalne, które nadal mają duże chłodzenie NH, jak omówiono pokrótce powyżej. Niewiele dużych pozatropikalnych erupcji SH jest znanych podczas ostatniego zlodowacenia (LIN22). Ogólnie rzecz biorąc, chociaż klasyfikacja według LIN22 jest pomocna, nie uważamy, że uzasadnia ona a priori faworyzowanie niektórych erupcji w naszej analizie statystycznej. Istnieje duża niepewność w klasyfikacji poszczególnych erupcji ze względu na dużą zmienność międzyrdzeniową depozycji siarczanów. Ponadto względne osadzanie siarczanów nie ogranicza w pełni faktycznego wpływu erupcji na klimat, który zależy również od tego, gdzie i jak długo aerozole były transportowane. Niemniej jednak klasyfikacja według LIN22 w odniesieniu do naszych wyników powinna być brana pod uwagę w przyszłych badaniach.