Физика Земли, 2021, № 2, стр. 130-138

ВВЕДЕНИЕ

Высоковязкие риолитовые лавы образуют потоки, которые варьируются от десятков до сотен метров в толщину и обычно имеют незначительную длину (первые километры). Часто риолитовые потоки ассоциируются с континентальными горячими точками, поднимающимися струями горячих мантийных пород (или плюмами). Наиболее известным примером служит Йеллоустонская горячая точка в США, расположенная в настоящее время под штатом Вайоминг. Сейсмическая томография позволяет выделить канал низких скоростей под Йеллоустонским национальным парком, который и интерпретируется как мантийный плюм [Sigloch et al., 2008]. Эта горячая точка, вероятнее всего, отвечает за образование расплава для суперизвержений и крупных лавовых потоков и за создание, таким образом, вулканической провинции, покрывающей территорию нескольких штатов западнее Вайоминга в результате движения Североамериканской плиты над горячей точкой в течение последних 17 млн лет (см., напр., [Morgan, 1971; Smith et al., 2009; Camp et al., 2017]).

Крупные риолитовые потоки Йеллоустона в геологическом прошлом являются слабо изученными объектами. С использованием аналитического решения задачи течения вязкой несжимаемой жидкости на плоской поверхности [Huppert, 1982], в работе [Loewen et al., 2017] было показано, что риолитовый поток в районе Саммит-Лейк (Summit Lake; рис. 1а), образовавшийся в результате извержения лав 124000 лет назад, продолжался в течение приблизительно 2–5 лет при температуре порядка 800°C и высоких расходах магмы >100 м³ ⋅ с^–1. Такие высокие расходы, тем не менее, связаны с небольшими скоростями излияния (<1 см ⋅ с^–1), поскольку извержение происходило из трещины длиной около 6 км и большой площадью поперечного сечения, что позволило извержению оставаться эффузивным (невзрывным).

Хотя с помощью аналитического решения были оценены несколько важных физических параметров лавовых потоков в районе Саммит-Лейк Йеллоустона [Loewen et al., 2017], особый интерес представляет определение влияния вязкости лавовых потоков на их морфологию и продолжительность растекания лав. Несмотря на то, что потоки риолитовой лавы обычно не растекаются далеко от жерла из-за их высокой вязкости, в отдельных случаях лава может продвигаться на расстояние нескольких километров, если объем извергающейся магмы и толщина лавового потока достаточно велики или из-за более высокой температуры магмы, которая снижает ее вязкость. Лавовый поток при этом остывает медленно, его вязкость не увеличивается, и потому он может распространяться на большие расстояния.

Численные модели по оценке скорости теплопотерь при кондуктивном охлаждении крупных риолитовых лав показали, что охлаждение потоков толщиной 100–300 м происходит очень медленно и в дальнейшем тормозится образованием панциря (корки) и выделением скрытого тепла кристаллизации [Manley, 1992]. Модели показывают, что потоки лавы большой толщины могут оставаться активными в течение нескольких десятилетий. В то же время, по мере охлаждения лавового потока предел текучести его охлажденной корки будет оказывать влияние на его динамику [Balmforth, Craster, 2000] и приводить к более сложным его морфологиям [Fink, Griffiths, 1998]. Однако в работе [Loewen et al., 2017] было показано, что большой по объему риолитовый поток Саммита-лейк не имеет сложных морфологических структур, и поэтому он может быть смоделирован с помощью гравитационной изотермической модели вязкой жидкости.

Выбор лавового потока в районе Саммит-Лейк Йеллоустона связан с его гигантским объемом (порядка 50 км³) и большой толщиной (100–250 м), что при относительно коротком времени излияния потока приводит к его незначительному остыванию и образованию тонкой корки лавового потока. Хотя корка утолщается в случае нелинейного теплообмена (конвективного переноса и теплового излучения) по сравнению со случаем только конвективного переноса тепла на границе раздела с воздухом [Tsepelev et al., 2019], это увеличение толщины корки не оказывает существенного влияния на продвижение лавового потока большой толщины. Это позволяет в качестве первого приближения использовать предположение об изотермичности потока. В данной работе мы рассматриваем трехмерную модель изотермического течения лавы в районе Саммит-Лейк Йеллоустона с целью определения влияния топографии и вязкости лавовых потоков на процесс растекания лавы и продолжительность ее растекания.

ПОСТАНОВКА ЗАДАЧИ

Рассмотрим трехмерную численную модель течения двухфазной несжимаемой несмешивающейся вязкой жидкости, которая аппроксимирует экструзию лавы (одна фаза) из жерла вулкана в воздух (другая фаза) и последующий лавовый поток. Модельная область $\Omega = [0,{{l}_{1}}]\,\, \times \,\,[0,{{l}_{2}}]\,\, \times \,\,[f({{x}_{1}},{{x}_{2}}),f({{x}_{1}},{{x}_{2}}) + {{l}_{3}}]$ ограничена снизу топографией исследуемой местности $f({{x}_{1}},{{x}_{2}})$, где l₁ и l₂ – горизонтальные размерности и l₃ – высота модельной области. В этой области изучается течение вязкой ньютоновской неоднородной несжимаемой жидкости в поле силы тяжести. Это течение в декартовых координатах описывается нестационарным уравнением Навье–Стокса [Ismail-Zadeh, Tackley, 2010; Tsepelev et al., 2019] с начальным условием ${\mathbf{u}}(t = 0,{\mathbf{x}}) = 0$:

и уравнением неразрывности где: ${\mathbf{u}} = ({{u}_{1}},{{u}_{2}},{{u}_{3}})$ – вектор скорости; ${\mathbf{G}} = (0,0, - g\rho )$ – вектор внешних массовых сил; g – ускорение силы тяжести; p – давление; $\rho $ – плотность; $\eta $ – вязкость; ${\mathbf{x}} = ({{x}_{1}},{{x}_{2}},{{x}_{3}})$ – пространственная переменная; t – время. Перенос двухфазной жидкости с начальным условием $\alpha (t = 0,{\mathbf{x}}) = {{\alpha }_{0}}({\mathbf{x}})$ описывается уравнением: где $\alpha (t,{\mathbf{x}}) \in [0,1]$ определяет объемную долю жидкости в точке ${\mathbf{x}} \in \Omega $ в момент времени t. В начальный момент модельная область заполнена воздухом и поэтому ${{\alpha }_{0}}({\mathbf{x}}) = 0.$ Плотность и вязкость тогда описываются равенствами: где: ${{\rho }_{a}}$ и ${{\eta }_{a}}$ – плотность и вязкость воздуха; ${{\rho }_{l}}{\text{ и }}{{\eta }_{l}}$ – плотность и вязкость лавы.

Хотя тепловые эффекты влияют на образование лавовых корок и течение лавы, было показано, что при высоких расходах магмы толщина корки достаточно малая (~3–5%) по сравнению с толщиной лавового потока [Tsepelev et al., 2019]. Корка лавы растрескивается, несется потоком лавы под действием гравитационной силы до тех пор, пока ее толщина малая и не влияет значительно на продвижение потока [Tsepelev et al., 2016]. В данной модели рассматривается гравитационное течение лавы без учета тепловых эффектов.

На нижней модельной поверхности и боковых границах задается условие прилипания: ${\mathbf{u}} = 0$. На той части нижней поверхности, где расположено жерло, задаются условия: ${\mathbf{u}} = {{{\mathbf{u}}}_{0}},{\text{ }}\rho = {{\rho }_{l}},{\text{ }}\eta = {{\eta }_{l}}$, где ${{{\mathbf{u}}}_{0}}$ – скорость экструзии магмы. На верхней границе модельной области задаются условия: $\nabla {\mathbf{u}} = 0,{\text{ }}p = 0,{\text{ }}\rho = {{\rho }_{a}},{\text{ }}\eta = {{\eta }_{a}}$.

Лавовые потоки в районе Саммит-Лейк Йеллоустона вносят вклад в современную топографию исследуемого региона. Для моделирования лавовых потоков в геологическом прошлом необходимо преобразовать современную топографию таким образом, чтобы скорректировать ее за счет снятия толщины излившихся лав. Далее описывается трансформация рельефа исследуемой области. Пусть $F({{x}_{1}},{{x}_{2}})$ – современная топография изучаемой области, максимальная толщина лавы $Z(x_{1}^{v},x_{2}^{v})$, равная 230 м [Loewen et al., 2017], достигается в жерле вулкана и уменьшается линейно с расстоянием от жерла, и $\left( {x_{1}^{v},x_{2}^{v}} \right)$ – координаты пересечения прямых A–A', B–B' и C–С', где располагается жерло вулкана. Тогда модельная топография вычисляется по формуле $f({{x}_{1}},{{x}_{2}}) = F({{x}_{1}},{{x}_{2}})$ – ${{Z(x_{1}^{v},x_{2}^{v})} \mathord{\left/ {\vphantom {{Z(x_{1}^{v},x_{2}^{v})} {(1 + {\text{ }}r({{x}_{1}},{{x}_{2}}))}}} \right. \kern-0em} {(1 + {\text{ }}r({{x}_{1}},{{x}_{2}}))}}$, где $r({{x}_{1}},{{x}_{2}}) = {{\sqrt {{{{\left( {{{x}_{1}} - x_{1}^{v}} \right)}}^{2}} + {{{\left( {{{x}_{2}} - x_{2}^{v}} \right)}}^{2}}} } \mathord{\left/ {\vphantom {{\sqrt {{{{\left( {{{x}_{1}} - x_{1}^{v}} \right)}}^{2}} + {{{\left( {{{x}_{2}} - x_{2}^{v}} \right)}}^{2}}} } D}} \right. \kern-0em} D}$, D = 20 км. Рисунок 1б иллюстрирует положение вычисленной таким образом модельной топографии.

Таким образом, задача моделирования лавовых потоков сводится к решению уравнений (1)–(4) с описанными выше начальными и граничными условиями в модельной области $\Omega $. При численных расчетах все переменные обезразмериваются с масштабом времени $t* = {{{{\eta }_{l}}} \mathord{\left/ {\vphantom {{{{\eta }_{l}}} {({{\rho }_{l}}gh)}}} \right. \kern-0em} {({{\rho }_{l}}gh)}}$, длины h и скорости ${h \mathord{\left/ {\vphantom {h {t{\text{*}}}}} \right. \kern-0em} {t{\text{*}}}}$.

ЧИСЛЕННЫЙ МЕТОД

Для численного анализа модельной задачи мы используем программное обеспечение Ansys Fluent (https://www.ansys.com/products/fluids/ansys-fluent), основанное на методе конечных объемов. Модельная область разбивается на $174 \times 200 \times 70$ гексаэдров, составляющих конечные объемы. Численные коды используют многофазную нестационарную VOF модель, а решатель использует неявную по времени схему интегрирования уравнений (1)–(3) для совместного определения поля скоростей, давления и объемной доли жидкости. Для аппроксимации давления и лапласиана используются численные схемы второго порядка; для дискретизации конвективных членов – монотонные схемы (см., напр., [Ismail-Zadeh, Tackley, 2010]). Для совместного учета давления и скорости в уравнении (1) применяется аппроксимация PRESTO! и численный метод SIMPLE [Patankar, Spalding, 1972], где параметры релаксации выбраны равными 0.01 и 0.3 для скорости и давления соответственно. Учитывая нелинейность задачи, временной шаг выбирается в диапазоне от 1  до 10 с в зависимости от устойчивости и для оптимизации скорости сходимости решения линейных алгебраических уравнений (СЛАУ), полученных после дискретизации задачи. Для решения СЛАУ используется multi-grid метод (см., напр., [Ismail-Zadeh, Tackley, 2010]).

Основная трудность при численном моделировании данной задачи состоит в том, что отношение вязкости лавы к вязкости воздуха составляет ~14 порядков. Для осуществления устойчивого численного эксперимента была построена неадаптивная расчетная сетка, такая что в области границы лава–воздух сетка сгущалась в вертикальном направлении. Это позволило уменьшить численную диффузию в области этой границы. Также была выбрана неявная схема совместного интегрирования систем дифференциальных уравнений, которая позволяет выполнять стабильные вычисления с относительно большим временным шагом. Хотя неявные схемы интегрирования не требуют жесткого ограничения числа Куранта ($\left| {\mathbf{u}} \right|{{\Delta t} \mathord{\left/ {\vphantom {{\Delta t} {\Delta x}}} \right. \kern-0em} {\Delta x}}$, где $\left| {\mathbf{u}} \right|$ – магнитуда скорости, $\Delta t$ – шаг по времени и $\Delta x$ – пространственный шаг [Courant et al., 1928]), большие шаги по времени приводят к увеличению численной диффузии на границе лава–воздух и к нефизичности получаемых результатов, и поэтому шаг по времени выбирался достаточно малым. Явные схемы интегрирования, которые требуют значительно меньших затрат для численной реализации, не позволили осуществить стабильные расчеты даже при малых числах Куранта. Шаг по времени и параметр релаксации для скорости уменьшались эмпирически для обеспечения устойчивости вычислительного процесса.

Течения высоковязкой лавы являются ламинарными, так как число Рейнольдса достаточно малое (<10^–6) для параметров вязкости и плотности лавы, характерной длины и скорости течений в исследуемых моделях (см. табл. 1). Воздух не оказывает влияния на динамику течения лавовых потоков, так как его вязкость и плотность значительно ниже значений этих физических параметров для лавы, и граница между лавой и воздухом рассматривается как свободная. И хотя число Рейнольдса при этом достаточно большое в воздухе, программный продукт Fluent позволяет избегать турбулентных течений в воздушном слое.

РЕЗУЛЬТАТЫ

Предполагается, что магма в районе Саммит-Лейк Йеллоустона изливалась из жерла вулкана вдоль эллиптической трещины длиной 4 км (рис. 1а), расположенной на модельной топографии вдоль профиля С–С'. Здесь рассмотрены три численных эксперимента: 1) течение лавы по восстановленной выше топографии поверхности, 2) течение лавы по наклоненной подстилающей поверхности и 3) течение лавы с увеличенной на порядок (по сравнению с экспериментом 1) вязкостью. Значения модельных параметров представлены в табл. 1. Плотность лавы (${{\rho }_{l}} = 2350$ кг ⋅ м^–3) соответствует плотности риолитовых пород с высоким содержанием кремнезема [Loewen et al., 2017]. Вязкость лавы рассчитывается по формуле Хапперта [Huppert, 1982] ${{\eta }_{l}}\sim 0.023{{\rho }_{l}}g{{V}^{4}}{{Q}^{{ - 1}}}{{r}^{{ - 8}}}$, где V – объем лавового потока (~50–55 км³), Q – расход лавы (~2000–6000 м³ ⋅ с^–1) и r – радиус фронта лавового потока (~10–12 км). В данной работе принимаются два значения вязкости (${{\eta }_{l}} = 5.6 \times {{10}^{9}}$ Па ⋅ с и ${{\eta }_{l}} = 5.6 \times {{10}^{{10}}}$ Па ⋅ с), находящиеся в пределах вязкостей, определенных для лав с температурами между 750 и 850°С и концентрацией воды порядка 0.1 весовых процента [Farquharson et al., 2015; Romine, Whittington, 2015; Loewen et al., 2017].

Эксперимент 1. Эволюция лавового потока в первые 8 мес. с начала извержения представлена на рис. 2. Лава распространяется почти осесимметрично вплоть до 2 мес., а затем в юго-западном и северо-восточном направлениях, так как рельеф местности препятствует течению на север и юго-восток. Через 8 мес. современная область, покрытая риолитовой лавой, заполняется модельной лавой за исключением ее южной и западной оконечностей, однако распространение модельной лавы в северо-восточном направлении, заметном на рис. 2 с 4-го мес., не подтверждается наблюдениями. Течение лавы вдоль разреза B–B' является асимметричным: оно простирается на 2 км к северо-востоку и более чем на 22 км к юго-западу от предполагаемого жерла. На северо-востоке нет ни одной современной топографической высоты, которая могла бы ограничить поток лавы в этом направлении.

Если считать, что оконтуривание лавового потока произведено с учетом возможной эрозии и размываний [Christiansen, 2001], то можно предположить, что либо в прошлом существовала возвышенность, ограничивающая течение на северо-восток, либо средний угол наклона подстилающей поверхности в прошлом отличался от современного на незначительную величину (приблизительно на 1.3°). Это может объясняться общей деформацией Йеллоустонской кальдеры за счет изменения давления в магматическом очаге во время излияния лавовых потоков, хотя доказательств данного предположения в литературе не приводится.

Эксперимент 2. Для того чтобы проверить гипотезу об измененном угле наклона подстилающей поверхности, был проведен численный эксперимент с измененным ускорением силы тяжести так, чтобы создать условие течения на юго-запад без изменения рельефа местности (см. табл. 2 ). На рис. 3 сравниваются результаты расчетов численных экспериментов 1 и 2. В экс. 2 модельный лавовый поток направлен на юго-запад, и современная область, покрытая лавой на юге, практически полностью заполняется модельной лавой за 6 мес. (рис. 3б). Однако модельная лава выходит за пределы современной границы лавового потока на западе и востоке, что может быть связано с неточностью в построении модельной топографии поверхности или выбора расхода, или выбора длительности течения. При этом лава заполняет достаточно плотно современную лавовую область на северо-востоке, хотя лавовый поток распространяется медленно в северо-западном направлении.

Эксперимент 3. Вязкость лавы влияет на морфологию и процесс распространения лавового потока. Проводился численный эксперимент по моделированию течения лавы при ее вязкости, увеличенной на один порядок по сравнению с предыдущими экспериментами (см. табл. 1). На рис. 4 показано трехмерное изображение течения модельной лавы в случае экс. 1 (левая колонка; ${{\eta }_{l}} = 5.6 \times {{10}^{9}}$ Па ⋅ с) и экс. 3 (средняя колонка; ${{\eta }_{l}} = 5.6 \times {{10}^{{10}}}$ Па ⋅ с). В правой колонке сравниваются формы растекания этих лавовых потоков. На рис. 5 приведено сравнение результатов расчетов вдоль трех профилей А–А', B–B' и C–C' для экспериментов 1 и 3. Из-за высокой вязкости в экс. 3 течение в горизонтальном направлении замедляется и потому высота лавового потока увеличивается, изменяя его морфологию (рис. 5).

ОБСУЖДЕНИЕ И ВЫВОДЫ

Расчеты показали хорошее совпадение формы поверхности лавового потока и его толщины с измеренными величинами при условии, что средний угол наклона подстилающей поверхности отличается от современного на незначительную величину (~1.3°). Величина вязкости в численных экспериментах находится в пределах, оцененных в работе [Loewen et al., 2017] для лавовых потоков в районе Саммит-Лейк Йеллоустона, хотя оценки потенциальных вязкостей лавы имеют неопределенность, связанную с неопределенностью в оценках объема излившейся лавы.

Заметим, что при объеме в 50 км³ излившихся лав и модельном расходе лавы $3500{\text{ }}{{{\text{м}}}^{3}} \cdot {{{\text{с}}}^{{ - {\text{1}}}}}$ (скорости экструзии $\sim 2.7 \times {{10}^{{ - 3}}}{\text{ м}} \cdot {{{\text{с}}}^{{ - 1}}}$) понадобится приблизительно 5.5 мес. непрерывного извержения из модельной трещины. Однако максимальная ширина модельной трещины достаточно велика (400 м) и нереалистична. Поэтому результаты расчетов должны быть скорректированы по времени за счет ширины трещины. При уменьшении максимальной толщины трещины до 40 м и той же скорости экструзии потребуется ~4.6 лет для извержения такого же количества лавы. Это согласуется с оценками длительности лавового потока в районе Саммит-Лейк. В целом трехмерное моделирование лавовых потоков подтвердило оценки, сделанные на основе аналитической модели в приближении тонкого слоя [Loewen et al., 2017], и позволило более детально исследовать динамику лавовых потоков в районе Саммит-Лейк и их морфологию в зависимости от направленности течения и вязкости лавы.