Results and Discussion
Tworzenie klatratów metanu badano wykonując sześć niezależnych i równoległych symulacji MD (Runs 1-6). Migawki z typowego procesu tworzenia się klatratu (przebieg 5) pokazane są na Rys. 1A. Czas nukleacji definiujemy jako moment, w którym liczba cząsteczek wody w największym krysztale hydratu zaczyna wzrastać. Wyraźnie widać, że czasy nukleacji dla sześciu symulacji są dość różne i wahają się od 183 do 550 ns (Bieg 1: ∼183 ns, Bieg 2: ∼208 ns, Bieg 3: ∼456 ns, Bieg 4: ∼273 ns, Bieg 5: ∼150 ns, i Bieg 6: ∼550 ns), jak pokazano na Rys. 1B. Można się zastanawiać, czy zauważalne różnice w czasach nukleacji wynikają z różnych szybkości rozpuszczania się klastra metanu w wodzie. Jak widać na Rys. 1C, klaster metanu ma prawie taką samą szybkość rozpuszczania i osiąga stałą liczbę ∼ 120 cząsteczek metanu (co odpowiada ułamkowi molowemu 3.3%) w ciągu 100 ns, co wskazuje, że to samo stężenie metanu prowadzi do zarodkowania w tych sześciu układach. Stąd, różnice w czasach zarodkowania są nadal zastanawiające, mimo że zachowanie to zostało przypisane głównie stochastycznemu zdarzeniu ewolucji struktury wody w trakcie zarodkowania klatratu.
Statystyki agregatowe multiringów wodnych. (A) Ewolucja struktur pierścieni wodnych z 4-7 cząsteczkami wody. (B) Ujęcia konkretnych TWRA: 2-1-0 i 3-0-0. TWRAs są zdefiniowane jako struktury hydratów łączące trzy typy pierścieni k (k = 5, 6, 7 były badane ze względu na znaczące zmiany ich liczby podczas tworzenia się klatratu), które są oznaczone jako i-j-l (np. 3-0-0, 2-1-0, 1-2-0, 0-2-1…), gdzie i, j, oraz l reprezentują liczby pierścieni odpowiednio 5, 6, oraz 7. Ewolucja: (C) dwóch specyficznych i (D) innych TWRA względem parametru porządkowego liczby 5 pierścieni w przebiegu 5. Szare słupki reprezentują etap nukleacji. N0 to liczba TWRA w konfiguracji początkowej. Ni to liczba TWRA przy określonej liczbie 5 pierścieni. Metoda obliczania Ni jest pokazana w Dodatku SI, Rys. S4.
Następnie przeanalizowano ewolucję struktur wodnych względem parametru porządkowego liczby 5 pierścieni, aby scharakteryzować ścieżkę nukleacji/wzrostu klatratów. Tutaj, struktury klatkowe, takie jak 512, 51262, 51263, 51264 i 4151064, nie zostały wybrane do prześledzenia ewolucji hydratów, ponieważ wiele nie dobrze uformowanych struktur klatkowych, takich jak struktury podobne do klatki, które nie są zamknięte, lub struktury amorficzne (21), są zaangażowane w proces zarodkowania/wzrostu katratu. TWRA (definicja w legendzie na Rys. 2), które mogą opisywać uporządkowane struktury wody w fazie ciekłej, a jednocześnie pełnią rolę ważnych fragmentów klatek wodnych, są wykorzystywane do śledzenia następujących ścieżek ewolucji w procesie tworzenia się klatratów: Woda ciekła → struktury amorficzne → klatki klatratów.
Ponieważ TWRAs 2-1-0 i 3-0-0 (pokazane na Rys. 2B) są głównymi agregatami multiringowymi wody w obrębie uformowanych klatek klatratów (takich jak 512, 51262, 51263 i 51264), są one traktowane jako struktury specyficzne, ponieważ ich liczba stale wzrasta podczas procesu nukleacji/wzrostu (Rys. 2C). W przypadku innych niespecyficznych TWRA obserwuje się różne procesy ewolucji, jak pokazano na Rys. 2D (przebieg 5) i w Dodatku SI, Rys. S5 (pozostałe przebiegi). Liczba TWRA 0-0-3 i 0-1-2 maleje wraz ze wzrostem parametru porządku, co wskazuje, że obie struktury TWRA od początku formowania się klatratu nieustannie przekształcają się w inne TWRA. Liczby TWRA 1-2-0 i 2-0-1 rosną najszybciej i osiągają swoje maksymalne wartości w okolicach czasu nukleacji. Ta cecha sugeruje, że tworzenie się struktur 2-1-0 i 3-0-0 jest prawdopodobnie kluczowym procesem dla zarodkowania klatratu. Porównując ewolucję różnych TWRA (zaznaczonych niebieską strzałką na Rys. 2D), zaproponowano następującą ścieżkę zarodkowania/wzrostu: 0-0-3 → 0-1-2 → (1-0-2, 0-2-1) → (0-3-0, 1-1-1) → (1-2-0, 2-0-1) → (2-1-0, 3-0-0), przy czym 1-2-0 i 2-0-1 są postrzegane jako struktury pośrednie podczas procesu nukleacji, a 2-1-0 i 3-0-0 jako struktury charakterystyczne, jak pokazano na Rys. 3. Zaproponowana ścieżka zarodkowania/wzrostu może być dalej weryfikowana w oparciu o inne różne parametry porządkowe, różne stężenia CH4 (ułamki molowe 2.8 i 7.9%, jak pokazano w Dodatku SI, Rys. S5F i S7, odpowiednio), oraz różne symulacyjne układy początkowe, takie jak zarodkowanie zainicjowane równomiernie rozłożonym w roztworze CH4 (Dodatek SI, Rys. S7) lub zarodkowanie zainicjowane mieszanymi cząsteczkami gościa CH4/H2S (Dodatek SI, Rys. S8). Te niezależne symulacje sugerują, że zidentyfikowana ścieżka zarodkowania/wzrostu jest dość ogólna w procesie tworzenia się klatratów. Co ciekawe, stwierdzamy, że ścieżka zarodkowania/wzrostu na Rys. 3 wydaje się być procesem zrzucania warstwy hydratacyjnej z cząsteczek metanu. Tutaj każda główna zmiana strukturalna jest związana z utratą jednej lub dwóch cząsteczek wody z warstwy hydratacyjnej metanu.
Ścieżka zarodkowania/wzrostu hydratów. Schematyczny schemat ścieżki ewolucji TWRA (zaznaczonej grubą niebieską strzałką) w zarodkowaniu hydratów. Liczby w nawiasach oznaczają liczbę cząsteczek wody w odpowiadających im strukturach TWRA. Tworzenie struktur TWRA 2-1-0 i 3-0-0 (zaznaczone czerwonym przerywanym prostokątem) jest prawdopodobnie kluczowym procesem dla zarodkowania klatratów.
Na podstawie powyższej analizy zidentyfikowaliśmy cztery kluczowe struktury TWRA, tj. dwie pośrednie TWRA (1-2-0, 2-0-1) i dwie specyficzne TWRA (2-1-0, 3-0-0), jako charakterystyczne struktury elementarne w procesie zarodkowania/wzrostu klatratów. Następnie zbadano czynniki odpowiedzialne za powstawanie tych czterech lokalnych struktur. Jak pokazano na Rys. 4A, CRk reprezentuje k liczbę ściśniętych regionów zajętych przez cząsteczki metanu. Duża liczba ściśniętych regionów odzwierciedla agregacje wielu cząsteczek metanu. Na Rys. 4 B-E (wyniki rundy 5), liczby różnych TWRA (1-2-0, 2-0-1, 2-1-0, 3-0-0) utworzonych w CR0, CR1, CR2, CR3 i CR>3 liczbie ściśniętych regionów są zliczane w zależności od czasu symulacji (wyniki innych rund są pokazane w Dodatku SI, Rys. S9, które wykazują ten sam trend). Można zauważyć, że pośrednie struktury TWRA (1-2-0, 2-0-1) powstają głównie przed czasem nukleacji (t < 150 ns). Co ważniejsze, obie struktury pośrednie mogą być wytworzone tylko przy dużej ilości CR3 i CR>3. Gdy skompresowanych regionów jest mniej niż 3, obie struktury pośrednie prawie w ogóle nie powstają. To samo zachowanie można zaobserwować dla dwóch specyficznych TWRA (2-1-0, 3-0-0), szczególnie w okresie nukleacji (100-200 ns w przebiegu 5). W rezultacie obserwuje się silną zależność pomiędzy tworzeniem kluczowych TWRA a agregacją cząsteczek metanu.
Dystrybucja TWRAs w warstwie hydratacyjnej metanu. (A) Rozpuszczony metan jest klasyfikowany według liczby utworzonych regionów ściśniętych (CR0, CR1, CR2, CR3, i CR>3). Cząsteczki metanu i wody są reprezentowane odpowiednio przez białe kulki i czerwone kropki. Warstwy hydratacyjne centralnej i otaczającej cząsteczki metanu są reprezentowane przez czerwoną i cyjanową ramkę drucianą. Obszar ściśnięty pomiędzy czerwonym i cyjanowym obramowaniem jest oznaczony żółtą powierzchnią. Liczby TWRAs CRk (k = 0-3 lub k > 3) utworzonych w różnych ściskanych regionach są pokazane odpowiednio w (B) 1-2-0, (C) 2-0-1, (D) 2-1-0, i (E) 3-0-0. Etapy nukleacji zaznaczono pionowymi zielonymi paskami.
Aby wyjaśnić dlaczego agregacja cząsteczek metanu może promować powstawanie kluczowych TWRA, przeanalizowano profile funkcji rozkładu radialnego (RDF) CH4-H2O na różnych etapach tworzenia klatratu. Na Rys. 5A porównano RDF CH4-H2O w fazie ciekłej (Rys. 5 C, ①) i w rejonie agregacji metanu (tuż przed nukleacją, Rys. 5 C, ②). Widać, że zdolność hydratacji CH4 w fazie ciekłej jest bardzo niska (krzywa niebieska). Dlatego też struktury wodne w fazie ciekłej nie są znacząco uporządkowane przez rozproszony CH4. W związku z tym, kluczowe TWRA nie mogą być wytworzone. W miarę jak cząsteczki CH4 agregują się w roztworze, RDF CH4-H2O wyraźnie wzrasta (czarna krzywa). Sugeruje to, że kompresja cząsteczek CH4 może sprawić, że cząsteczki wody staną się uporządkowane, tworząc obfite kluczowe TWRA, co następnie aktywuje nukleację klatratu. Wynik ten wskazuje, że agregacje cząsteczek CH4 odgrywają znaczącą rolę w promowaniu tworzenia się klatratów poprzez generowanie bogatych struktur wodnych wokół CH4, co również wyjaśnia szeroko akceptowany model blob (21).
Agregacja metanu i zrzucanie zbędnych cząsteczek H2O. (A i B) RDF pomiędzy CH4 i H2O w różnych okresach nukleacji/wzrostu klatratu. (C) Schematyczny schemat tworzenia się klatratów. Zielone kulki reprezentują metan, a czerwone słupki, krzywe wiggly i linie reprezentują cząsteczki wody.
Furthermore, Fig. 5B porównuje RDF dla CH4-H2O tuż przed nukleacją (czarna krzywa na Rys. 5B, czerwone krzywe wiggly na Rys. 5 C, ②) z RDF po utworzeniu klatratu (czerwona krzywa na Rys. 5B, czerwone linie na Rys. 5 C, ③). Wartości RDF sugerują, że metan uczestniczy w procesie dehydratacji podczas zarodkowania/wzrostu klatratu. Zachowanie to jest zgodne ze ścieżką dehydratacji TWRAs widoczną na Rys. 3. Na podstawie powyższej analizy można stwierdzić, że w procesie tworzenia się klatratów zachodzi prawdopodobnie następujący proces (Fig. 5C): Początkowo cząsteczki metanu agregują się w skompresowane regiony CR≥3. Następnie, niektóre cząsteczki metanu w CR≥3 zderzają się ze sobą, w wyniku czego cząsteczki wody oddalają się od warstwy hydratacyjnej metanu, co odzwierciedla pozorny proces dehydratacji w ścieżce zarodkowania/wzrostu, promując jednocześnie powstawanie kluczowych TWRA. Ten zidentyfikowany mechanizm tworzenia hydratów nazywamy mechanizmem ściskania/ zrzucania warstwy hydratacyjnej.
Aby dalej zwalidować mechanizm ściskania/ zrzucania warstwy hydratacyjnej, skonstruowano modelowy układ agregacji metanu z 46 cząsteczkami metanu i 257 cząsteczkami wody w pudełku ∼2 × 2 × 2 nm3. Temperatura wody jest utrzymywana na stałym poziomie 250 K, podczas gdy dla cząsteczek metanu rozważane są trzy różne temperatury, aby ocenić wpływ kompresji/wyrzucania warstwy hydratacyjnej na tworzenie się klatratów. Jak widać na Rys. 6A, gdy temperatura cząsteczek metanu wzrasta, szybkość tworzenia się klatratów znacznie wzrasta (v275 K > v250 K > v225 K). Zachowanie to można wytłumaczyć mechanizmem ściskania/wyrzucania warstwy hydratacyjnej: Ponieważ cząsteczki metanu w wyższych temperaturach mogą powodować silniejsze nakładanie się/kompresję warstwy hydratacyjnej, a następnie szybsze odrywanie (zrzucanie) warstwy hydratacyjnej podczas separacji cząsteczek metanu, może nastąpić szybsze tempo nukleacji/wzrostu. Wynik ten jest zgodny z poprzednim raportem, w którym stwierdzono, że wyższa temperatura może promować proces zarodkowania dla przesyconego roztworu metanu (36). Ponadto, zbudowaliśmy idealny model do odtworzenia nukleacji opartej na ściskaniu/wyrzucaniu warstwy hydratacyjnej (Rys. 6B). Wirtualne miejsce (bez masy, ładunku i oddziaływań) jest stałe i łączy się z otaczającymi je cząsteczkami metanu (kolor różowy) poprzez wirtualne wiązania (b = 0.61 nm). Centralny metan (kolor niebieski) wibruje harmonicznie w pobliżu miejsca wirtualnego, naśladując proces ściskania/usuwania warstwy hydratacyjnej. Gdy centralny metan oscyluje, wypycha niektóre cząsteczki wody z własnej powłoki hydratacyjnej, a następnie prowadzi do ścieżki nukleacji/wzrostu. Ogólnie rzecz biorąc, mechanizm kompresji/wyrzucania warstwy hydratacyjnej oferuje odpowiedź na pytanie, dlaczego blob jest wymagany do nukleacji klatratów i, co ważniejsze, na pytanie, dlaczego ewolucja struktur klatratowych jest związana z dynamicznym zachowaniem cząsteczek metanu.
Mechanizm ściskania/wyrzucania warstwy hydratacyjnej. (A) Procesy tworzenia klatratów dla cząsteczek metanu w różnych temperaturach. (B) Mechanizm kompresji/wyrzucania warstwy hydratacyjnej oparty na procesie nukleacji dla modelu idealnego.
Następnie, rozmiar krytycznego jądra i skala czasowa nukleacji zostały oszacowane w oparciu o mechanizm kompresji/wyrzucania warstwy hydratacyjnej. Biorąc pod uwagę, że agregacja i kompresja cząsteczek metanu w CR>3 działa jako kluczowy krok w zarodkowaniu klatratu, średnia odległość pomiędzy centralną cząsteczką metanu a otaczającymi ją cząsteczkami w CR>3 jest analizowana w celu określenia krytycznego jądra, jak pokazano na Rys. 7A. Przed zarodkowaniem, średnia odległość silnie oscyluje wokół 0.56 nm. Po rozpoczęciu nukleacji, oscylacja ta staje się znacznie słabsza. Ostatecznie, ta średnia odległość wzrasta do ∼0.6 nm, wskazując, że struktury amorficzne zakończyły konwersję do postaci jądra krytycznego. I w tym momencie liczba regionów CR>3 zbliża się do pięciu we wszystkich sześciu symulacjach, jak pokazano na Rys. 7B. W rezultacie, różnica w czasie nukleacji powinna być przypisana raczej prawdopodobieństwu formowania się regionu CR>3 niż losowej fluktuacji termicznej struktur wodnych. Co więcej, pozycje tych pięciu regionów CR>3 znajdują się blisko siebie (Rys. 7C), utrzymując kompresję/wyrzucanie warstwy hydratacyjnej pomiędzy cząsteczkami metanu, jednocześnie promując nukleację. Sugerujemy zatem, że jądro krytyczne składa się z około pięciu regionów CR>3, posiadających ∼13 cząsteczek metanu i ∼100 cząsteczek wody. Wartość pięć jest również zweryfikowana metodą forward flux sampling oraz analizą prawdopodobieństwa committer (43⇓⇓-46). Metody i szczegółowe wyniki są przedstawione w dodatku SI.
Jądro krytyczne w nukleacji klathratu dla przebiegów 1-6. (A) Ewolucja średniej odległości pomiędzy centralną cząsteczką metanu a otaczającymi ją cząsteczkami w CR>3 w sześciu niezależnych symulacjach MD. Czerwone pionowe paski reprezentują okres nukleacji. Pomarańczowy pionowy pasek reprezentuje przypadek, gdy jądro powstaje w tym okresie czasu, ale następnie ulega rozpadowi. (B) Ewolucja liczby regionów CR>3 w sześciu niezależnych symulacjach MD. (C) Migawki jądra krytycznego w sześciu niezależnych symulacjach MD.
Po określeniu liczby cząsteczek metanu w jądrze krytycznym, skala czasowa zarodkowania klatratu może być oszacowana przy użyciu następującego modelu analitycznego. Dla układu o objętości Vbox, N cząsteczek metanu jest rozpuszczonych w roztworze. Zakładając, że gdy n cząsteczek metanu połączy się w jądro krytyczne o objętości Vn, następuje inicjacja nukleacji. W związku z tym, N cząsteczek metanu można podzielić na dwa regiony, z których n cząsteczek jest rozmieszczonych w objętości Vn, a kolejne N – n cząsteczek jest rozmieszczone w pozostałej objętości Vbox – Vn. Dla każdej losowo zadanej konfiguracji całego układu, takie prawdopodobieństwo rozmieszczenia cząsteczek metanu (P0) spełnia rozkład dwumianowy dany równaniem (1):P0=∑nN(CNn(VnVbox)n(1-VnVbox)N-n),gdzie P0 można również traktować jako prawdopodobieństwo wystąpienia nukleacji. Zatem prawdopodobieństwo bez wystąpienia nukleacji (Pno) powinno byćPno=1-P0.Prawdopodobieństwo bez wystąpienia nukleacji po i-tej niezależnej ramce (Pnoi) jestPi,no=Pnoi, podczas gdy prawdopodobieństwo wystąpienia nukleacji po i-tej niezależnej ramce (Pyesi) jestPi,yes=1-Pi,no.Tutaj Pyesi jest oparte na kolejnych i niezależnych ramkach. W realistycznym procesie dynamicznym, każda niezależna ramka może trwać bardzo krótko. Jak widać na funkcji autokorelacji prędkości (VACF) dla układu klatratu, czas korelacji dla każdej ramki wynosi ∼0.1 ps (VACF maleje do 0; Dodatek SI, Rys. S10). W rezultacie, prawdopodobieństwo zarodkowania Pyesi można obliczyć jako ti/(0.1 ps), gdzie ti jest czasem MD dla i-tej ramki. Biorąc pod uwagę, że liczba cząsteczek metanu w jądrze krytycznym wynosi 13 (n ∼ 13), a objętości jądra krytycznego i pola symulacji wynoszą w przybliżeniu n × 0.15 nm3 (n = 13) i ∼(4.7)3 nm3 (dla Vn, Vbox), odpowiednio, prawdopodobieństwo nukleacji (NP = Pi,yes) jest obliczane w funkcji czasu ewolucji układu z różnymi liczbami rozpuszczonych cząsteczek metanu, jak pokazano na Rys. 8. Generalnie, gdy NP zbliża się do ∼0.5, prawdopodobieństwo zarodkowania jest większe. Dla układów rozpatrywanych w niniejszej pracy (krzywa czerwona), przewidywana skala czasu zarodkowania oparta na powyższym wzorze analitycznym dobrze zgadza się z wynikami symulacji. Ponadto, stwierdzamy, że czas zarodkowania jest silnie skorelowany z koncentracją metanu. Gdy stosunek molowy metanu do wody zmniejsza się z 3.9 do 2.8%, czas nukleacji wzrasta 10-krotnie (z ∼ 100 ns do ∼ 1000 ns). Względna niepewność czasu zarodkowania jest bardzo mała przy wysokich stężeniach metanu. Na przykład, maleje ona od ∼ 2500 ns przy stosunku molowym 2.8% (od ∼500 ns do ∼3,000 ns; NP: 0.5-1) do ∼50 ns przy stosunku molowym 3.9% (od ∼50 ns do ∼100 ns; NP: 0.5-1). Wynik ten jest zgodny z wynikami wcześniejszych badań (33, 47). W celu dalszej weryfikacji wzoru analitycznego użytego do oszacowania skali czasu nukleacji, przeprowadzono dodatkowe symulacje (patrz czarne kwadraty na Rys. 8; stosunek molowy 2.8%), wybrano również inną istotną pracę (33) (niebieskie kwadraty; stosunek molowy 3.9%). Wszystkie te wyniki były zgodne z naszymi szacunkami.
Skala czasowa nukleacji klatratu metanu. Oszacowane prawdopodobieństwo nukleacji (Pi,yes) versus czas ewolucji MD dla układów o różnych stosunkach molowych metan/woda. Niebieskie, czerwone i szare obszary reprezentują zakresy czasu dla zarodkowania o wysokim prawdopodobieństwie dla układów o trzech różnych stosunkach molowych metan/woda wynoszących odpowiednio 3.9%, 3.3% i 2.8%. Dla porównania, obserwowane czasy nukleacji w niniejszej pracy oraz z poprzedniej pracy (33) są również przedstawione jako kwadraty o różnych kolorach.
