Cfd-подход к моделированию динамики лавовых потоков в вулканологии

Понимание поведения потоков лавы имеет основополагающее значение для оценки вулканических опасностей и разработки эффективных стратегий смягчения последствий. Во время извержений расплавленная порода может перемещаться со скоростью от метров в час до десятков километров в час, разрушая инфраструктуру, изменяя ландшафты и угрожая близлежащим популяциям. Вычислительная динамика жидкостей (CFD) обеспечивает надежную основу для моделирования сложной динамики движения лавы, позволяя ученым прогнозировать пути потока, скорости охлаждения и потенциальные зоны воздействия в широком диапазоне условий. В этой статье исследуются основные подходы CFD, используемые в вулканологии, физические процессы, которые они захватывают, и проблемы, которые остаются в достижении точных прогнозов в реальном времени.

Управляющие уравнения и реологические модели

По своей сути CFD решает уравнения Навье-Стокса для сохранения массы, импульса и энергии. Для потоков лавы эти уравнения должны быть адаптированы для учета неньютоновской реологии, сильной температурной зависимости вязкости и фазовых изменений (кристаллизации и газоразрешения). Обобщенное уравнение импульса для несжимаемой жидкости:

ρ (∂v/∂t + v· ⁇ v) = - ⁇ p + ⁇ ·τ + ρg

где ρ - плотность, v - скорость, p - давление, τ - девиаторный тензор напряжения, g - гравитационное ускорение. Ключевая сложность заключается в конститутивном отношении, связывающем напряжение со скоростью деформации. Лава ведет себя как вязкопластический материал: деформируется только тогда, когда приложенное напряжение превышает напряжение выхода τ0. Двумя распространенными реологическими моделями, используемыми в CFD, являются:

• Бингемская модель: τ = τ0 + μp γ? для τ > τ0, где μp — это вязкость пластика, а γ? — скорость сдвига. Эта модель фиксирует поведение потока пробки лавы — жесткую шапку, движущуюся по срезанному базальному слою.
• Модель Хершела-Булкли: τ = τ0 + K γ·^n, где K — индекс консистенции, а n — индекс потока (n < 1 for shear-thinning, n > 1 для сдвига). Эта более гибкая модель лучше представляет собой зависящее от температуры и скорости деформации поведение сложных магм.

Сама вязкость развивается по мере охлаждения лавы и зарождения кристаллов. Уравнение Роско-Эйнштейна часто используется для обновления эффективной вязкости как функции кристаллической фракции φ: μ eff = μ liquid · (1 - φ/φ max)^(-2,5 φ max). Над критической кристаллической фракцией (~60%) суспензионная переходная сила переходит в твердое поведение, создавая прочность выхода, которая может остановить поток. Температура отслеживается с помощью энергетического уравнения, которое включает в себя термины для адвекции тепла, проводимости, скрытого тепла кристаллизации и радиационного охлаждения на поверхности. Эти связанные уравнения делают проблему вычислительно интенсивной, но они необходимы для захвата остановки потоков лавы.

Численные методы моделирования потока лавы

Метод конечного объема (FVM)

Метод конечного объема делит вычислительную область на небольшие объемы управления и решает интегральную форму уравнений сохранения. FVM по своей сути консервативен - потоки массы, импульса и энергии точно сбалансированы по клеточным граням - что делает его подходящим для моделирования с высоким разрешением лавы, распространяющейся по сложной топографии. Популярные коды CFD с открытым исходным кодом, такие как OpenFOAM и FEniCS, были адаптированы для вулканологических применений. Исследователи использовали FVM для моделирования деталей лавовых каналов и долей в масштабе метра, захватывая образование дамб и переход от канализированного к листовому потоку. Метод хорошо обрабатывает сложные граничные условия; например, состояние без скольжения может применяться на поверхности земли, в то время как состояние свободного скольжения или сдвига-стресса применяется на вершине потока, если лава изолирована коркой. Однако FVM может бороться с резкими интерфейсами (например, между лавой и твердой землей) и часто требует адаптивной сетчатой очистки для решения фронта потока.

Метод конечных элементов (FEM)

Метод конечных элементов дискретизирует область на элементы (треугольники, четырехугольники, тетраэдры) и приближает решение с использованием базисных функций. FEM особенно привлекателен для моделирования потока лавы, поскольку он может обрабатывать нерегулярные, неструктурированные сетки, которые соответствуют реальным цифровым моделям возвышения (DEMs). Метод позволяет локальную уточнение сетки в областях с высоким градиентом, таких как фронт потока и близкие препятствия. Более поздние реализации используют вариационный многомасштабный (VMS) подход для стабилизации адвекционного доминирования импульса и энергетических уравнений, которые возникают, когда потоки лавы быстро спускаются по крутым склонам. FEM был использован в моделировании базальтовых потоков на горе Этна и Килауэа, где топографическая сложность (конусы скопий, трещины разлома, лавовые трубки) значительно влияет на направление потока. Ограничение стандартного FEM заключается в том, что он не является строго консервативным для массы и энергии; однако смешанные формулировки и использование прерывистых методов Галеркина могут

Методы без сетки: гидродинамика сглаженных частиц (SPH)

SPH - это лагранжианская безмешниковая техника, в которой жидкость представлена набором движущихся частиц, каждая из которых несет физические свойства. Частицы взаимодействуют через сглаживающее ядро, а управляющие уравнения решаются в их лагранжевой форме. SPH идеально подходит для моделирования больших деформаций, потоков свободной поверхности и взаимодействий со сложными границами - все характерно для лавы. Поскольку положения частиц адвентируются потоком, SPH может естественным образом захватывать расщепление доли лавы, образование нескольких каналов и размещение обломков. Метод успешно применяется для моделирования извержения Этны 2001 года и извержения нижней Восточно-Рифтовой зоны 2018 года, воспроизводя наблюдаемое продвижение потока с хорошей точностью. Основными недостатками SPH являются высокая вычислительная стоимость (число частиц может превышать миллионы) и трудности при наложении точных граничных условий (например, отсутствие скольжения на земле). Кроме того, модели SPH часто требуют тщательной калибровки вязкости и напряжения выхода для соответствия наблюдениям.

Сравнение и гибридные подходы

Каждый численный метод имеет сильные и слабые стороны. FVM консервативен и надежен для крупномасштабных операционных прогнозов; FEM превосходит по гибкости и адаптивности сетки; SPH является естественным для потоков свободной поверхности. Некоторые современные коды объединяют эти методы: конечный объем или конечный элементный решатель для основной массы потока в сочетании с подходом частица-в-ячейке (PIC) или адаптивная детализация сетки (AMR) на фронте потока. Кроме того, использование погруженных граничных методов позволяет FVM/FEM-решателям включать топографические данные без генерации сеток, оснащенных корпусом, сокращая время предварительной обработки.

Основные физические процессы в динамике потока лавы

Охлаждение, кристаллизация и эволюция вязкости

По мере продвижения лавы она теряет тепло в атмосферу и в землю. Верхняя поверхность быстро охлаждается, образуя тонкую кору, которая изолирует внутреннюю часть, позволяя потоку перемещаться на большие расстояния. Внутри потока температура остается выше жидкого в течение некоторого времени, но по мере ее падения кристаллы зародышевываются и растут, увеличивая объемную вязкость. Переход от ньютоновской к Bingham или Herschel-Bulkley жидкости контролируется кристаллической фракцией. Модели CFD обычно параметризируют это через связь между температурой и кристалличностью (фракция твердой фазы), используя данные экспериментальной петрологии или термодинамических моделей, таких как MELTS. Скрытое тепло, выделяемое во время кристаллизации, частично компенсирует радиационное охлаждение, стабилизируя длину потока. Неспособность объяснить эти эффекты приводит к предсказаниям нереалистично коротких или тонких потоков.

Газоразрешение и пузырьковый поток

Многие потоки лавы содержат растворенные летучие вещества (в основном H2O, CO2, SO2), которые растворяются при падении давления вблизи вентиляционного отверстия. Пузыри могут значительно снизить объемную плотность и вязкость пены, увеличивая скорость потока и расстояние выдоха. Для захвата этого поведения необходимы модели двухфазного потока (жидкость + газ) или модели смеси. Взаимодействие между пузырьками и жидкой фазой часто описывается с использованием уравнения Рэлея-Плессета для роста пузыря и модели скорости скольжения (например, модель дрейфа-потока Зубера-Финдлея). В то время как большинство операционных кодов CFD рассматривают поток как однофазный с эффективной плотностью и вязкостью, исследовательские коды все чаще включают двухфазные эффекты, чтобы лучше соответствовать наблюдаемой инфляции и дегазации активных потоков.

Топографическое взаимодействие и формирование Леви

Топография перед извержением — в том числе долины, хребты и ранее существовавшие лавовые поля — сильно контролирует направление и скорость потока. Модели CFD используют топографические данные высокого разрешения (например, 1 м воздушно-капельной лидарной или спутниковой стереоизображения) для создания цифровых моделей возвышения (DEM), с помощью которых строится вычислительная сетка. Взаимодействие между лавой и топографией приводит к образованию рукавов, которые являются самообразующимися поднятыми ободками, ограничивающими поток. Образование Леви — это процесс обратной связи: медленно движущиеся, высоковязкие края потока охлаждаются и затвердевают, создавая границу, которая предотвращает боковое распространение. Модель CFD должна точно захватывать теплообмен и реологию на полях для воспроизведения этого явления. Несколько исследований показали, что только путем включения зависящего от температуры напряжения выхода модели производят реалистичные ширины и высоты леви.

Ввод данных и граничные условия

Надежное моделирование CFD требует высококачественных входных данных. Наиболее важными параметрами являются:

• Скорость извержения (скорость выпота): Разряд лавы за единицу времени, обычно измеряемый в м3/с. Это можно оценить по спутниковым тепловым изображениям, полевым измерениям или историческим средним значениям. Скорость выпота определяет интенсивность потока и является основным контролем длины потока.
• Начальная температура лавы:] Для базальтовых магм температура жидкого вещества составляет ~1200°C; для андезитических или дацитических лав она ниже (~900-1100°C). Начальная температура влияет на скорость охлаждения и кинетику кристаллизации.
• Топография (DEM): Для точного предсказания траектории рекомендуется горизонтальное разрешение не менее 10 м. Степерные наклоны (≥30°) могут вызывать быстрое ускорение и переход от ламинарного к турбулентному потоку.
• Материальные свойства: Плотность (обычно 2500—2800 кг/м3), теплоемкость, теплопроводность, излучательность для радиационного охлаждения, и реологические параметры (удельный стресс, индекс консистенции, индекс потока). Лабораторные измерения на природных образцах являются золотым стандартом, но на практике эти значения настроены против прошлых извержений.

На вентиляционном отверстии накладываются условия фиксированной скорости (или постоянной скорости потока массы). На поверхности земли используется условие без скольжения, часто совмещенное с границей теплового потока (например, конвективный перенос тепла на подстилающую породу). Верхняя поверхность моделируется как свободный скольжение или сегментированная граница: там, где существует кора, условие без скольжения может применяться локально. Для крупномасштабных симуляций вычислительная область должна простираться достаточно далеко, чтобы содержать весь путь потока, который может составлять десятки километров. Это увеличивает размер сетки и вычислительное время.

Проверка и тематические исследования

Kīlauea 2018 Извержение в Восточно-восточной рифтовой зоне

Извержение вулкана Килауэа на острове Гавайи в 2018 году произвело один из самых разрушительных потоков лавы в новейшей истории, уничтожив более 700 структур. Обширные полевые наблюдения, изображения дронов и спутниковые данные (например, из Обсерватории вулканов Гавайев и USGS) обеспечивают отличный ориентир для моделей CFD. Исследователи применили подходы FVM и SPH для моделирования потока из трещины 8 вентиляционных отверстий, включая DEM с высоким разрешением и измеренные скорости выпота. Модель SPH Dietterich et al. (2022) воспроизвела общую геометрию доли потока и время прорыва леви с ошибкой менее 10% в охвате области. Модель подчеркнула критическую роль ранее существовавших потоков лавы Пуу О'О в управлении лавой 2018 года, деталь, которую упрощенные (например, самые крутые) модели не смогли захватить.

Гора Этна 2021 Пароксизмальный эпизод

В течение февраля-марта 2021 года на горе Этна на Сицилии наблюдался ряд впечатляющих лавовых фонтанов, каждый из которых генерировал кратковременный, но быстро движущийся поток лавы, который продвигался на несколько километров вниз по Валле-дель-Бове. Наблюдения итальянского Национального института геофизики и вулканологии (INGV) предоставили данные высокого временного разрешения о скоростях выпота и фронтальном продвижении потока. Моделирование конечных элементов с использованием кода LavaSIM (на основе неглубоководного приближения уравнений Навье-Стокса) смогло предсказать окончательную длину потока в пределах 5%, когда в качестве входа использовалась история скорости выпота. То же исследование показало, что поток развивался пульсирующим образом из-за взаимодействия между охлаждением и скоростью подачи, явление, которое не может быть воспроизведено более простыми статистическими моделями.

Проблемы и ограничения

Несмотря на значительный прогресс, CFD-моделирование потоков лавы по-прежнему сталкивается с рядом препятствий:

• Расчетная стоимость: 3D-моделирование полного извержения с высоким разрешением может занять дни или недели на суперкомпьютерах, что делает прогнозирование в реальном времени непрактичным. Даже модели со средним значением глубины 2D (например, уравнения с мелководьем) требуют тщательного проектирования сетки, чтобы избежать непомерного времени выполнения.
• Реологическая неопределенность:] Реология естественной лавы не до конца понята, особенно для богатых кристаллами и пузырчатых магм. Лабораторные эксперименты на переплавленных образцах могут не представлять поведение быстро охлаждающихся, газозаряженных потоков. Неопределенность в напряжении выхода и вязкости может привести к большим изменениям прогнозируемой длины потока.
• Многомасштабная физика:] Динамика пузырьков и рост кристаллов происходят в масштабах от микрометров до миллиметров, в то время как поток продвигается на километры.Объединяя эти масштабы в одном моделировании (многомасштабное моделирование) является активной темой исследования.
• Отсутствие данных в реальном времени: Во время продолжающегося извержения скорость выпота и температура лавы могут быстро меняться. Модели CFD, которые предполагают постоянный вход, не могут захватывать фазы воска и убывания. Включение тепловых и позиционных данных в реальном времени со спутников (например, VIIRS, Sentinel-2) в структуру ассимиляции данных является многообещающим, но сложным путем.
• Изменение топографии: Сам поток лавы изменяет топографию, создавая новые земли (например, образование дельты, если лава попадает в море). Большинство моделей CFD предполагают фиксированную топографию, которая может стать неточной для длительных извержений.

Будущие направления

Машинное обучение - улучшенные симуляции

Суррогатные модели, основанные на глубоком обучении (например, сверточные нейронные сети или графовые нейронные сети), могут быть обучены на библиотеке предыдущих симуляций CFD для получения почти мгновенных прогнозов. Эти «эмуляторы» могут предсказать путь потока и толщину с учетом скорости выпота и топографии без решения полных уравнений Навье-Стокса. Первые такие модели для потоков лавы были разработаны Mossoux et al. (2021) [FLT: 1] и показаны с точностью, сопоставимой с физическими моделями, но с ускорением от трех до четырех порядков величины. Гибридные подходы, которые используют физико-информированные нейронные сети (PINN) для обеспечения законов сохранения также изучаются.

Усвоение данных и оперативное прогнозирование

Ассимиляция данных в реальном времени - объединение спутниковых наблюдений с прогнозами модели CFD - может уменьшить неопределенность и улучшить навыки прогнозирования. Фильтр Ансамбля Калмана и фильтры частиц были применены для корректировки реологических параметров на лету по мере продвижения потока. Программа USGS по оценке рисков вулканов инвестирует в операционные инструменты CFD, которые могут работать на скромных вычислительных кластерах и предоставлять обновленные карты опасности каждый час. Проблемы включают задержку спутниковых наблюдений (обычно 10-30 минут) и необходимость надежной автоматизации сетчатой генерации и конфигурации решателя.