Вулканология и сейсмология, 2021, № 3, стр. 18-28

1. ВВЕДЕНИЕ

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

Детальный мониторинг лавовых куполов велся на нескольких вулканах, например на Сент-Хеленс в США [Swanson et al., 1987], Пинатубо на Филиппинах [Daag et al., 1996], Унзен в Японии [Nakada et al., 1999], Сантьягуито (Санта-Мария) на Гватемале [Harris et al., 2003], Мерапи и Синабунг в Индонезии [Voight et al., 2000; Nakada et al., 2019], Суфриер-Хиллс на Монтсеррате [Watts et al., 2002; The Eruption …, 2014] и Колима в Мексике [Zobin et al., 2015]. Мониторинг позволяет картировать пространственное и временное развитие лавовых куполов и определять морфологические изменения в процессе их роста, а также определять изменение объема во времени (расход лавы).

На морфологию лавовых куполов влияют реология магмы и скорость ее поступления из канала вулкана (расход лавы, РЛ). Вязкость магмы зависит от температуры и объемной доли кристаллов, которая в свою очередь определяется кинетикой кристаллизации (характерным временем роста кристаллов, ВРК [Tsepelev et al., 2020]). При малых значениях ВРК, то есть при быстрой кристаллизации лавы, обелискообразные структуры развиваются при низких РЛ и блинообразные структуры при высоких РЛ; при больших ВРК купола имеют форму лавовых выжимок или образуют блинообразные структуры. Было показано, что охлаждение не играет существенной роли при развитии лавового купола. Это связано с тем, что если изменение количества кристаллов зависит только от температуры, то вязкость лавы увеличивается лишь в приповерхностном слое купола, а толщина температурного пограничного слоя остается небольшой по сравнению с высотой купола [Tsepelev et al., 2020]. В теле купола значительное увеличение вязкости происходит за счет кристаллизации, вызванной потерей летучих. Таким образом, по известной реологии, зависящей от ВРК, и для известного РЛ можно смоделировать эволюцию лавового купола. В то же время интерес представляет решение обратной задачи: по форме лавового купола и известному расходу восстановить реологические параметры лавы.

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

2. ПОСТАНОВКА ЗАДАЧИ И МЕТОД ЧИСЛЕННОГО МОДЕЛИРОВАНИЯ ЛАВОВЫХ КУПОЛОВ

Рассматривается двумерная осесимметричная модель двухфазной несмешиваемой несжимаемой жидкости, которая аппроксимирует лаву (одна фаза) и воздух (другая фаза), разделенные подвижной границей – поверхностью лавового купола. Влияние воздушной фазы на рост лавового купола незначительно из-за большого отношения плотностей и вязкостей фаз. В модельной области $\Omega $ (рис. 2) движение лавы описывается следующим набором уравнений вместе с начальными и граничными условиями [Ismail-Zadeh, Tackley, 2010; Tsepelev et al., 2019, 2020].

Мы используем уравнения Навье–Стокса с начальным условием ${\mathbf{u}}(t = 0,{\mathbf{x}}) = 0$ и уравнение неразрывности для описания динамики лавы

где ${\mathbf{x}} = ({{x}_{1}},{{x}_{2}}) \in \Omega $ – декартовы координаты; t $ \in [0,\vartheta ]$ – время; $\vartheta $ – конечный момент времени (длительность модельных экспериментов); ${\mathbf{u}} = ({{u}_{1}}(t,{\mathbf{x}}),\,\,{{u}_{2}}(t,{\mathbf{x}}))$ – скорость; $\rho $ – плотность; $\eta $ – вязкость; $p = p({\mathbf{x}})$ – давление; ${\mathbf{g}} = (0,g),$ g – ускорение силы тяжести; $\nabla ,$ ^T и $\left\langle {\centerdot ,\centerdot } \right\rangle $ обозначают вектор градиента, транспонированную матрицу и скалярное произведение векторов соответственно. Мы пренебрегаем температурной зависимостью физических параметров среды и силами поверхностного натяжения. Модельные плотность и вязкость представлены как $\rho = {{\rho }_{L}}\alpha (t,{\mathbf{x}})$ + + ${{\rho }_{A}}(1 - \alpha (t,{\mathbf{x}}))$ и $\eta = {{\eta }_{L}}\alpha (t,{\mathbf{x}})$ + ${{\eta }_{A}}(1 - \alpha (t,{\mathbf{x}}))$ соответственно. Здесь ${{\rho }_{A}}$ – плотность воздуха; ${{\rho }_{L}}$ – плотность лавы; ${{\eta }_{A}}$ – вязкость воздуха и ${{\eta }_{L}}$ – вязкость лавы.

Функция $\alpha (t,{\mathbf{x}})$ принимает значение единицы для лавы и нуля для воздуха в каждой точке x и в каждый момент времени t, и эта функция переносится со скоростью u согласно уравнению адвекции

с начальным условием $\alpha (t = 0,{\mathbf{x}}) = 0,$ означающим, что в начальный момент времени вся расчетная область заполнена воздухом.

Мы предполагаем, что вязкость лавы (измеряемая в Па · с) зависит от объемной доли кристаллов [Costa et al., 2009] как

где $\varphi = {\phi \mathord{\left/ {\vphantom {\phi {{{\phi }_{*}}}}} \right. \kern-0em} {{{\phi }_{*}}}},$ $\phi $ – объемная доля кристаллов; ${{\phi }_{*}}$ – специфическое значение объемной доли кристаллов; B – теоретическое значение коэффициента Эйнштейна, который определяется из уравнения Эйнштейна как $B = {{({{\eta }_{L}}(\phi ) - 1)} \mathord{\left/ {\vphantom {{({{\eta }_{L}}(\phi ) - 1)} \phi }} \right. \kern-0em} \phi }$ [Mardles, 1940] (экспериментально установлено, что коэффициент Эйнштейна колеблется от 1.5 до 5 [Jeffrey, Acrivos, 1976]); $\delta = 7.24,$ $\gamma = 5.{\text{76}}$ и $\xi = 4.63 \times {{10}^{{ - 4}}}$ [Lejeune, Richet, 1995; Costa et al., 2009]; ${\text{erf}}( \cdot )$ – функция ошибок. Объемная доля кристаллов определяется из эволюционного уравнения, описывающего упрощенную кинетику роста кристаллов при кристаллизации вследствие дегазации магмы с начальным условием $\phi (t = 0,{\mathbf{x}}) = 0.$ Здесь ${{\phi }_{{eq}}}$ – объемная доля кристаллов в равновесии, которая зависит от количества воды, растворенной в магме, и от температуры; $\tau $ – ВРК. Чем меньше ВРК, тем быстрее процесс кристаллизации сходится к своему равновесному состоянию. ВРК называют временем релаксации, которое требуется для уменьшения разности в е (~2.72) раз между фактическим ($\phi $) и равновесным (${{\phi }_{{eq}}}$) значениями объемных долей кристаллов по отношению к начальной разности (${{\phi }_{{in}}} - {{\phi }_{{eq}}}$), где ${{\phi }_{{in}}}$ – начальная объемная доля кристаллов в магме. Для $\alpha = 1$ и u = 0, ВРК может быть найдено аналитически как $\tau = - t{{\left( {\ln \left[ {{{({{\phi }_{{eq}}} - \phi )} \mathord{\left/ {\vphantom {{({{\phi }_{{eq}}} - \phi )} {({{\phi }_{{eq}}} - {{\phi }_{{in}}})}}} \right. \kern-0em} {({{\phi }_{{eq}}} - {{\phi }_{{in}}})}}} \right]} \right)}^{{ - 1}}}.$ Значения $\tau $ могут варьироваться от нескольких часов до нескольких месяцев в зависимости от температуры и водонасыщенности магмы, количества ранее существовавших кристаллов в магме и ее состава [Tsepelev et al., 2020]. Заметим, что, хотя вязкость зависит также от петрологического (химического) состава лавы и содержания летучих в лаве (ее водонасыщенности), эти зависимости вязкости не рассматриваются в данной работе.

На границе модели $\Gamma = {{\Gamma }_{1}} \cup {{\Gamma }_{2}} \cup {{\Gamma }_{3}} \cup {{\Gamma }_{4}} \cup $ ${{\Gamma }_{5}} \cup {{\Gamma }_{6}}$ задаются следующие условия (см. рис. 2). На границе ${{\Gamma }_{1}}$ задаются условия симметрии, то есть условие непроницаемости $\left\langle {{\mathbf{u}},{\mathbf{n}}} \right\rangle = 0$ и условие свободного скольжения $\left( {\nabla {\mathbf{u}} + \nabla {{{\mathbf{u}}}^{T}}} \right){\mathbf{n}}$ – ‒ $\left\langle {\left( {\nabla {\mathbf{u}} + \nabla {{{\mathbf{u}}}^{T}}} \right){\mathbf{n}},{\mathbf{n}}} \right\rangle {\mathbf{n}} = 0.$ Предполагается, что лава плотности ${{\rho }_{L}}$ и вязкости ${{\eta }_{L}}$ поступает в модельную область через границу ${{\Gamma }_{2}}$ при заданном РЛ Q. На границах ${{\Gamma }_{3}},$ ${{\Gamma }_{4}}$ и ${{\Gamma }_{5}}$ задается условие прилипания u = 0. Воздух удаляется из модельной области через границу ${{\Gamma }_{6}}$ согласно заданному РЛ. Предполагается, что объемная доля кристаллов равна $\phi = {{\phi }_{{in}}}$ на ${{\Gamma }_{2}}$ и $\phi = 0$ на ${{\Gamma }_{6}}.$ Значения параметров модели, используемые при численном моделировании, представлены в табл. 1.

Для численного решения задачи (1)–(5) используется метод конечных объемов, реализованный в программном обеспечении Ansys Fluent. Для определения положения границы лавового купола используется метод объема жидкости (VOF) [Hirt, Nichols, 1981]. Используется неявная схема интегрирования для решения уравнений (1)–(5) с соответствующими граничными и начальными условиями. Давление дискретизируется по схеме второго порядка PRESTO! [Peyret, 1996]. Для аппроксимации оператора Лапласа мы используем численную схему второго порядка точности, а монотонные схемы используются для дискретизации конвективных членов в уравнениях [Ismail-Zadeh, Tackley, 2010]. Используется SIMPLE метод [Patankar, Spalding, 1972] для решения уравнений (1)–(2), где параметры релаксации выбираются равными 0.01 и 0.3 для скорости и давления соответственно. Для функции $\alpha $ и объемной доли кристаллов параметр релаксации выбирается равным 0.5. Учитывая большой скачок между вязкостями лавы и воздуха, выбор параметров релаксации является критическим. Временной шаг выбирается в диапазоне от 0.1 до 1 с для оптимизации скорости и обеспечения сходимости решения системы линейных алгебраических уравнений (СЛАУ), полученной после дискретизации уравнения (1). Точность в решении СЛАУ для функции $\alpha $ и объемной доли кристаллов составляет 10^–6.

3. СОЗДАНИЕ БАЗЫ МОРФОЛОГИЧЕСКИХ ФОРМ И СИНТЕЗАЦИЯ НАБЛЮДАЕМОГО КУПОЛА

Создание базы морфологических форм модельных лавовых куполов рассмотрим на примере построения одного купола. Модельная область дискретизируется шестигранными ячейками (примерно 70000 ячеек). Зададим следующие значения ВРК $\tau = 5 \times {{10}^{4}}$ с, РЛ Q = 0.7 м³ с^–1 и радиуса жерла вулканического извержения r = 15 м. Используя параметры, указанные в табл. 1, средствами Ansys Fluent численно решаем задачу (1)–(5) на промежутках времени, указанных в табл. 2.

Рассмотрим прямоугольную область ${{\Omega }_{{ABCD}}} \subset \Omega $ с вершинами A = (0, 30), B = (100, 30), C = (100, 100), D = (0, 100). Обозначим L = $ = L(\tau ,Q,r,t) \subset {{\Omega }_{{ABCD}}}$ – область, содержащую лавовый купол, и $L{\kern 1pt} ' = {{\Omega }_{{ABCD}}}\backslash L$ – область, содержащую воздух. Назовем границу между лавой и воздухом $F = F(\tau ,Q,r,t)$ = $\left\{ {({{x}_{1}},{{x}_{2}}) \in \partial L\backslash \partial {{\Omega }_{{ABCD}}}} \right\}$ границей (или морфологической формой) лавового купола.

По результатам вычислительного эксперимента для момента времени ${{t}_{k}}$ находим форму купола F_k и помещаем ее в базу морфологических форм. Для организации тестового исследования морфологические формы были рассчитаны также для значений $\tau = \{ 1.8 \times {{10}^{4}},$ $6 \times {{10}^{4}},$ $7 \times {{10}^{4}},$ $8 \times {{10}^{4}},$ $5 \times {{10}^{5}}\} $ с и r = 5 м (см. табл. 2). В дальнейшем база данных может пополняться модельными формами куполов для различных параметров ВРК, РЛ и радиуса жерла вулканического извержения. На рис. 3 приведены несколько морфологических форм куполов, полученных при моделировании для различных параметров задачи. Были выбраны моменты времени, в которые площадь области L (лавового купола) примерно одинакова.

Так как основная цель работы – найти в базе данных такую форму купола, которая наилучшим способом аппроксимирует природный лавовый купол, в данном исследовании был создан также синтетический купол, в качестве возможного природного лавового купола. Для этого из базы морфологических форм выбирается любой купол, например купол с морфологической формой ${{F}_{{135}}} = F(\tau = 6 \times {{10}^{4}}\,\,{\text{с}},$ $Q = 0.7\,\,{{{\text{м}}}^{3}}\,{{{\text{с}}}^{{ - 1}}},$ $r = 15\,\,{\text{м}},$ $t = 5 \times {{10}^{4}}\,\,{\text{с}})$ (рис. 4а). Вдоль границы этого купола случайным образом вносится “зашумление”, в результате которого и получается синтетическая форма купола F* (см. рис. 4б).

Данное зашумление производится следующим образом. В некоторой точке $a = ({{a}_{1}},{{a}_{2}}) \in F$ границы лавового купола задается случайное смещение этой точки как $a + \delta = ({{a}_{1}} + {{\delta }_{1}},{{a}_{2}} + {{\delta }_{2}}),$ где ${{\delta }_{1}},$ ${{\delta }_{2}}$ – нормально распределенные случайные величины c нулевым математическим ожиданием и стандартным отклонением, равным 1 [Венцтель, 1969]. Для каждой смещенной точки рассматривается круг ${{C}_{{a + {{\delta }}}}}\left( {\left| \varepsilon \right|} \right)$ с центром в точке $a + \delta $ и радиусом |ε|, где $\varepsilon $ – равномерно распределенная случайная величина, принимающая значения на отрезке $\left[ { - \omega ,\omega } \right],$ $\omega $ – постоянная величина; в данном случае ω = 1. Далее строится следующая область:

Таким образом, синтетический купол с границей F* получается с помощью объединения купола с границей F с множеством случайных кругов в случае неотрицательных ε и с помощью усечения синтетического купола множеством случайных кругов при отрицательных ε. Такой шум при искусственном искажении границы купола позволяет имитировать реальные искажения, которые могут быть следствием как неравномерности роста и частичным обрушением купола, так и погрешностей измерения его формы.

4. МЕТОДЫ РАСПОЗНАВАНИЯ ОБРАЗОВ

Сравнивать форму синтетического купола F* с формами лавовых куполов F_k из базы данных будем с помощью методов, используемых в теории распознавания образов [Salomon, 2007]. Представим лавовый купол в виде бинарного изображения. Для этого введем равномерное прямоугольное разбиение области ${{\Omega }_{{ABCD}}}$ на $I \times J$ ячеек ${{\Omega }_{{ij}}}\left( {{{\Omega }_{{ABCD}}} = \bigcup\nolimits_{i = 0,j = 0}^{I - 1,J - 1} {{{\Omega }_{{ij}}}} } \right).$ Каждой форме купола F поставим в соответствие прямоугольную матрицу $P(F) = \{ {{p}_{{ij}}}\} _{{i = 0,j = 0}}^{{I - 1,J - 1}}$ размера $I \times J,$ ячейки которой содержат 0, если в соответствующей ячейке ${{\Omega }_{{ij}}}$ находится более 50% воздуха, и 1 во всех остальных случаях. Оценивать близость формы синтетического купола F* и произвольной формы купола F_k из базы морфологических форм будем с помощью следующих функционалов качества.

Функционал на основе симметрической разности:

где $S( \cdot )$ (м²) – площадь области и k₁ (м^–2) – масштабирующий множитель.

Функционал на основе метрики пикового отношения сигнала к шуму PSNR (peak signal-to-noise ratio) [Salomon, 2007]:

где k₂ – масштабирующий множитель и k₃ – положительная константа. В случае полного совпадения P(F) и P(F*) пользователю выдается сообщение, содержащее номер купола, на котором достигнуто это условие.

Функционал на основе индекса структурного сходства SSIM (structure similarity) [Wang et al., 2004]:

где $\mu (P),$ $\mu (P*)$ – математическое ожидание, $\sigma (P,P*)$ – ковариация, ${{\sigma }^{2}}(P),$ ${{\sigma }^{2}}(P*)$ – дисперсия, k₄ – масштабирующий множитель, c₁ = 0.01, c₃ = = 0.03. Здесь рассматривается вероятностная модель представления изображений, а именно, изображение исследуется как поле случайных величин, и значение в каждой точке этого поля представляет собой реализацию случайной величины.

Значения функционалов вычисляются на каждом элементе базы данных, и полученный набор значений упорядочивается по убыванию. При этом то, что тот или иной функционал дает более точную оценку, основывается на информации о синтетическом куполе. Например, на рис. 5 выведены значения функционалов на элементах, представленных на рис. 3 и показано, что все три функционала достигают минимума на куполе 135. Заметим, что функционалы J₁, J₂ и J₃ оценивают количественное отклонение модельного и синтетического куполов, в то время как функционал J₃ – также структурные особенности форм куполов, хотя приводит к сложным и ресурсоемким вычислениям.

5. РЕЗУЛЬТАТЫ: ОПРЕДЕЛЕНИЕ ВЯЗКОСТИ ЛАВЫ ПО ИЗВЕСТНОЙ ФОРМЕ КУПОЛА

Используя описанные выше функционалы, сравним синтетический купол F* со всеми куполами F_k из базы морфологических форм. Масштабирующие множители подберем так, чтобы значения всех функционалов на элементах базы лежали на отрезке [0, 1]: k₁ = 1/30, k₂ = 1/25, k₃ = 33, k₄ = 1/14. На рис. 6 и 7 на горизонтальной оси указаны номера k куполов F из базы, на вертикальной оси указаны значения функционалов J_m(F*, F), m = 1, 2, 3. Из рис. 7 видно, что набор элементов, на которых достигаются наименьшие значения рассматриваемых функционалов, практически идентичен, наиболее близкие к синтетическому куполу F* купола – F₄₁, F₁₃₃, F₁₃₄ и F₁₃₅. На рис. 8 изображены синтетический купол и четыре наиболее близких модельных купола из базы данных. Время образования купола, расход лавы и размер жерла помогут с выбором близкого модельного купола.

Имея информацию о том, что синтетический купол F* был получен путем “зашумления” модельного купола F₁₃₅, можно сказать, что функционал на основе SSIM дает более точную оценку близости. Отметим, что для бинарного изображения функционалы J₁ и J₂ дают качественно схожие результаты и достаточно рассматривать один из них. Чтобы одновременно учесть качественную и количественную близость синтетического и модельного куполов можно применить линейную комбинацию описанных функционалов, тогда при подходящем выборе весовых коэффициентов можно достичь оптимальной качественной и количественной оценки. Например, рассмотрим функционал $J(F*,F)$ = $0.5\left[ {{{J}_{1}}(F*,F) + {{J}_{3}}(F*,F)} \right].$ При этом минимум такого функционала достигается для купола 135 (рис. 9). Хотя с помощью линейной комбинации функционалов можно точнее выбрать тот модельный купол, который будет минимизировать оптимальным образом природный купол, надо заметить, что на практике выбор весовых коэффициентов может представлять сложность, так как оцениваться будет природный, а не синтетический лавовый купол.

6. ОБСУЖДЕНИЕ И ЗАКЛЮЧЕНИЕ

Морфологические изменения лавовых куполов объясняются изменением вязкости магмы, вызванным дегазацией и кристаллизацией при подъеме магмы по каналу вулкана из очага извержения [Melnik, Sparks, 1999, 2005]. Быстрая скорость подъема магмы сокращает время, доступное для кристаллизации магмы, поэтому магма ведет себя как относительно маловязкая жидкость и после экструзии на поверхность растекается, образуя блинообразные морфологические структуры. С уменьшением РЛ или при малых ВРК кристаллическая магма переходит из текучего состояния к квазитвердому, образуя на поверхности лавовые выжимки и обелиски [Tsepelev et al., 2020].

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

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

Надо заметить, что рассматриваемая обратная задача является достаточно сложным объектом исследования как с теоретической, так и с вычислительной точек зрения [Самарский, Вабищевич, 2004; Кабанихин, 2009]. Обратные задачи являются, как правило, некорректными, т.е. неустойчивыми и могут не обладать единственностью решения. Рассматриваемая обратная задача является некорректной, она является неустойчивой (малые ошибки в исходных данных или ошибки округления могут привести к большим ошибкам в решении задачи), а также она может иметь неединственное решение. Например, при различных расходах лавы и различных характерных временах роста кристаллов может получаться похожая форма купола. Однако расход лавы можно практически определить по объему купола, и потому РЛ можно рассматривать как известную характеристику процесса. Классические методы не пригодны для решения некорректных задач. Альтернативными в этой ситуации являются различные подходы, связанные с математическими преобразованиями обратных задач [Тихонов, Арсенин, 1979]. Многие из таких подходов могут оказаться удовлетворительными при теоретическом исследовании задачи, но неудовлетворительными при ее численном решении. Другие подходы связаны с различными упрощениями обратных задач.

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

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

При применении данной методики к реальному извержению можно учесть отклонения расчетной формы купола от измеренной не только в некоторый момент времени, но и за весь период наблюдений, задать ограничения на форму кратера, подводящего канала и петрологию магмы. Реальные лавовые купола являются трехмерными объектами. И хотя представленный в работе подход является двумерным, он может быть распространен на трехмерный случай и использоваться для восстановления условий роста природных куполов.

ИСТОЧНИКИ ФИНАНСИРОВАНИЯ

Данное исследование поддержано грантом Российского научного фонда (РНФ-19-17-00027). Численные эксперименты проводились с использованием суперкомпьютера URAN в Институте математики и механики им. Н.Н. Красовского УрО РАН, г. Екатеринбург.