CC BY
217
УДК 519+531
МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ МЕГАЦУНАМИ
В.А. Симоненко, Н.А. Скоркин, В.П. Елсуков, А.С. Углов
С помощью конечно-разностных кодов: двухмерного вычислительного кода МЕЧ, предназначенного для расчета газодинамических задач; вычислительного кода ТМВ, предназначенного для решения задач, описываемых теорией мелкой воды, - проведено численное моделирование задачи о формировании и распространении мегацунами, образующегося в результате падения крупного небесного тела в воды мирового океана.
Ключевые слова: удар, астероид, комета, цунами, численное моделирование, теория мелкой воды.
Введение. В результате исследований последних десятилетий общепризнанным стало понимание того, что одну из главных опасностей для нашей цивилизации представляют столкновения с Землей малых космических тел, принадлежащих Солнечной системе, - астероидов, комет и их фрагментов. Такая опасность стимулирует дальнейшее расширенное изучение последствий катастрофических космических столкновений. Отдельного внимания заслуживают падения опасных космических объектов (ОКО) в океан. Прежде всего, вероятность таких падений примерно вдвое выше, чем на континенты. При этом значительная часть кинетической энергии передается толще воды и в дальнейшем переходит в энергию мегацунами. Еще одним существенным разрушительным фактором падений в океан может быть испарение большого количества воды. Образованные в результате удара пары воды устремляются в атмосферу и распространяются в верхних ее слоях. Появление большого количества водяного пара приводит к изменению оптических свойств атмосферы: уменьшается доля отраженного света, вступает в действие парниковый механизм. Он, в свою очередь, приводит к дополнительным испарениям больших водных масс. В глобальном масштабе это приводит к повышенным водным осадкам и крупномасштабным наводнениям. По-видимому, описания «всемирного потопа», присутствующие в ряде древних сказаний, являются дошедшим до нас свидетельством одного или нескольких таких явлений. В связи с этим представляет интерес более детальное исследование развития климатических последствий падений каменных и кометных тел в океан.
В представленной работе в основном рассматриваются процессы формирования и распространения в океане столкновительного цунами (мегацунами). При распространении цунами по океану фактически не происходит диссипации энергии. Уменьшение амплитуды волны происходит в основном за счет геометрического фактора - удаления фронта волны от источника. Благодаря большой длине волны (в направлении, перпендикулярном фронту), в открытом океане цунами не представляет опасности. Но при выходе на шельф из-за уменьшения массы вещества, вовлекаемого в движение, происходит увеличение амплитуды волны, повышение плотности энергии на ед. массы вещества. Скорость набегающего ротационного течения становится столь большой, что происходит интенсивная эрозия дна. В набегающий поток оказывается вовлеченным большое количество твердых частиц и каменных фрагментов. Вся эта масса накатывается на берег, приводя к громадным разрушениям и изменяя формы ландшафта в зоне затопления волной цунами. Хорошо наблюдаемым следствием воздействия грандиозных цунами являются характерные наносные образования, которые наблюдаются во многих местах побережья - прибрежные дюны, сложенные из пород, в основном поднятых со дна шельфовой зоны. Они имеют очертания шевронов, осевая линия которых ориентирована по направлению прихода и не согласуется с линиями преимущественных направлений ветра и линией максимального градиента шельфа. Такие прибрежные образования получили название шевронных дюн (рис. 1).
Прибрежные шевронные дюны являются вторичным признаком состоявшегося в прошлом достаточно мощного удара в океан, первичным - служит кратер на дне океана. В настоящее время выявлено несколько десятков подводных структур, которые могут быть импактными (ударными) кратерами. Одной из таких структур является кратер Беркла в Индийском океане [1] (30,87° ю.ш., 61,36° в.д.; диаметр 29,1 км). В ряде случаев положение таких структур достаточно хорошо можно определить по вторичным проявлениям на побережье. Так, в заливе Карпентария
с помощью спутниковой альтиметрии удается выделить две локализованные эллипсоидальные депрессии - Канмар (Капшаге, со средним диаметром 18 км) и Таббан (ТаЬЬап, со средним диаметром 12 км) [2]. Имеющиеся данные указывают, что соответствующее двойное падение произошло около 1500 лет тому назад.
Рис. 1. Шевронные дюны Фенамбози на южной оконечности Мадагаскара, высотой до 180 м,
в пяти километрах от океана
Более определенная информация имеется по кратеру Мауйка (МаИшка, 48,3° ю.ш., 166,4°
в.д.) [3]. Этот ударный кратер диаметром около 20 км располагается в шельфе на глубине 300 м на расстоянии 250 км от Новой Зеландии и 600 км от юго-восточного побережья Австралии. Совокупность данных указывает, что соответствующее столкновение произошло около 500 лет назад. На противоположной стороне Австралийского континента (северо-запад, область Кимберли) имеется еще один комплекс свидетельств воздействия в недалеком прошлом (по оценкам в 17 веке) гигантского цунами. Воздействием было охвачено около 1 500 км побережья Австралии. Однако до настоящего времени подводный кратер так и не обнаружен. Возможно, что столкновение произошло при косом падении кометного тела в области больших глубин.
Перечисленные примеры, если они не будут опровергнуты дальнейшими более детальными исследованиями, свидетельствуют, что прежние оценки частоты мощных столкновений являются заниженными. Они показывают, что в каждом тысячелетии происходят столкновения, оказывающие глобальное воздействие. Поэтому особое значение приобретает исследование сценариев развития таких событий, разработка средств своевременного обнаружения опасных объектов и способов предотвращения столкновений. В связи с этим большую ценность представляют более продвинутые технологии комплексного моделирования процессов, сопровождающих столкновения, несущих глобальную опасность. Данная работа посвящена именно таким методическим целям.
Постановка проблемы и способы решения. В силу существенных различий временных и пространственных масштабов комплексного явления, в настоящее время невозможно осуществить сквозное описание всех процессов с помощью единой физической модели и одного математического кода. Поэтому в представленной здесь технологии расчетов используются две модели и два кода; подробно - чуть ниже. В силу трудоемкости представленной технологии детальное описание с привлечением достаточно полной геологической и геофизической информации, учета трехмерного характера процессов следует осуществлять применительно к моделированию конкретных событий. В частности, особого внимания заслуживает изучение возможного падения ОКО в Индийский океан, которому соответствует кратер Беркла.
Исходя из методической направленности работы, при постановке задачи использован ряд упрощающих предположений. Чтобы не проводить трудоемкие трехмерные расчеты, была рассмотрена задача о вертикальном падении ОКО в океан, хотя вероятность такого события мала. В представленных расчетах для простоты было принято, что дно океана сложено однородными силикатными породами с начальной плотностью р = 2,67 г/см3; использовалось относительно простое уравнение состояния гранита, константы которого подбирались ранее для моделирования сильных взрывов на выброс. В работе рассматривалось падение каменного астероида. Предполагается, что тело имеет сферическую форму, массу М »1,4 Гт и скорость и = 22 км/с. При этом диаметр астроида (каменного тела) составляет 1 км, а свойства породы близки к свойствам гранита.
Качественная картина явления. Столкновение малых космических тел существенно зависит от многих факторов: размера тел, их состава, скорости сближения с Землей и угла входа в атмосферу, от свойств вещества в области падения. Движение космического тела в атмосфере
сопровождается образованием воздушной ударной волны большой интенсивности. В силу высокой скорости движения тела фронт головной ударной волны в воздухе будет отстоять на незначительном расстоянии от объекта, а часть воздуха, обтекающего тело, будет увлечена отошедшей ударной волной, фронт которой вытянут вдоль следа. Поэтому при рассмотрении удара тела о воду присутствием воздуха можно пренебречь. При столкновении такого тела с преградой в ней будут протекать следующие физические процессы: испарение воды, ударное нагружение, плавление и фазовые превращения пород океанического дна, и механическое разрушение их. На поверхности воды поднимется гигантская волна, которая способна распространяться на сотни и тысячи километров. Характерное время динамических процессов в окрестности удара составляет ~ 1 сек, а характерное время движения волны цунами в океане - несколько тысяч секунд.
Расчет столкновения космического тела с преградой был разбит на два этапа. На первом этапе с помощью конечно-разностной методики МЕЧ [4] с учетом поля силы тяжести был рассчитан процесс вертикального падения астероида в мировой океан, глубина которого в районе падения была принята равной 3 км. Расчеты проведены в газодинамическом приближении, т.е. не учитывалась прочность пород. При этом оценивалась зона испарения воды и астероида, а также возможность плавления материала скальной породы у дна океана. При описании материалов использовались уравнения состояния типа Ми-Грюнайзена с соответствующими параметрами.
К концу первого этапа продолжительностью ~12 с в системе «океан-скала» все нелинейные процессы, а именно: полиморфные переходы в воде и скале закончились, затухли ударные волны, сформировалась зона кратера, а сам астероид превратился в пар и множество мелких осколков. На поверхности воды сформировалась волна высотой ~7 км и диаметром воронки ~14 км. Плотность воды стала равной —1,0...1,01 г/см3, вектор скорости движения частиц среды в слое воды стал почти горизонтальным, а модуль вектора скорости - приблизительно постоянным по глубине воды. Разброс величины модуля скорости по глубине не превосходил 10 %. Остаточное давление паров воды в зоне кратера составляло примерно 18 атм. Так как длина волны, распространяющейся по поверхности воды, значительно больше глубины океана, то движение ее можно рассчитывать по теории мелкой воды (ТМВ). Иначе говоря, на втором этапе для описания распространения поверхностной волны в океане правомерно использовать классическую теорию мелкой воды.
Итак, для рассматриваемой проблемы необходимо было решить две задачи: 1) рассчитать начальную фазу удара и формирование поверхностной волны в океане; 2) решить задачу о движении поверхностной волны в океане и выходе ее на шельф. Профиль дна мирового океана и шельфа приведены далее на рис. 4, в связи с решением задачи о распространении цунами.
Приведем результаты расчетов первой задачи. Процесс соударения астероида с Землей при падении в океан показан на рис. 2, где на фоне изменяющейся геометрии взаимодействующих тел демонстрируется векторное поле скоростей (стрелки), а также изолинии удельной внутренней энергии, соответствующие испарению воды и гранита. Размерность длины - в км, времени - в сек. Ось Х - ось симметрии задачи, ось У направлена вдоль поверхности Земли. В воде выделена другим цветом дополнительная область - зона испарения.
В результате падения тела в океан образуется осесимметричная волна, высота которой в районе падения достигает почти 7 км и резко убывает при удалении от места падения. На дне океана при ударе астероидом образовался кратер, первичная глубина которого составляла —4 км, его диаметр - 8,8 км, а диаметр образовавшейся в воде воронки из смеси водяного пара и испаренного астероида - около 14 км. По результатам расчетов было установлено, что на большую высоту поднимается различная масса воды и водяного пара. При этом заметную массу водяного пара, поднявшуюся на высоту 7 км, можно наблюдать уже через 1,5 с после удара, а через 3 с этой высоты достигла вода в виде капель. Через 12 с после падения ОКО процесс поднятия больших масс водяного пара практически завершился. Дальнейшую динамику подъема воды на такую высоту в расчетах отследить не удалось, так как изначально не предполагалось считать задачу до больших времен. Для принятой постановки расчетов невозможно определить и размеры площади выпадения осадков от поднявшейся на большую высоту океанической воды. Эта отдельная задача должна решаться уже с учетом атмосферы Земли.
Рис. 2. Поле течения на моменты времени 1 с и 12 с. Стрелками указаны направления движения частиц среды
Перейдем ко второму этапу расчетов. В качестве исходных данных для расчетов распространения цунами использовался профиль поверхности воды, дна (левая часть рис. 4) и распределение скорости частиц воды, полученных из расчетов по вычислительному коду МЕЧ на момент времени 12 с.
75 700 750 800 850 900 950 1000 R, km
Рис. 4. Профили поверхности воды и дна океана на момент 12 с, полученные из расчетов по коду МЕЧ, и по коду ТМВ на различные моменты времени: 3500 с, 4000 с, 4500 с, 7000 с
Для описания цилиндрически расходящейся волны и ее воздействия на побережье с учетом профиля шельфа океана была разработана специальная программа ТМВ, в которой реализовано приближение мелкой воды.
Уравнения теории мелкой воды, описывающие распространение одномерной круговой поверхностной волны, имеют вид:
дг
д R
Здесь R - расстояние от начала координат (точка входа тела в воду), расположенного на невозмущенной поверхности океана; t - время; g - ускорение свободного падения; и = и(R,t)- скорость частиц на поверхности волны; г = h(R,t) - профиль поверхностной волны; £ = h(R) + h(R,t);
д и д u
----+ и---
д t д R
(1)
к(Я)- профиль дна мирового океана. Уравнения (1) являются уравнениями гиперболического типа и для них можно применять известные конечно-разностные аппроксимации.
Своеобразием рассматриваемой задачи является то, что распространение поверхностной волны осуществляется на очень большие расстояния —1000 км. Время распространения составляет величину порядка десятка тысяч секунд. Это означает, что при численном решении уравнений приходится делать большое число циклов по времени. Поэтому конечно-разностные уравнения, аппроксимирующие уравнения (1), должны обладать минимальным сглаживающим эффектом.
Были опробованы схемы Мак-Кормака, Лакса-Вендроффа, Неймана. Все они оказались неудовлетворительными. Заслуживающей внимания оказалась конечно-разностная 1 -схема Мо-ретти, которая при аппроксимации дифференциальных уравнений использует информацию о характеристиках уравнений (1). Схема Моретти не нуждается в использовании искусственной вязкости, дает монотонные профили, в ней практически нет численной диссипации. Однако прямое использование 1 -схемы Моретти в уравнениях (1) давало решения, физически неправдоподобные.
Поэтому была выполнена модификация 1 -схемы, которая сохранила все положительные свойства исходной схемы Моретти. С помощью этой модифицированной схемы была решена задача о распространении поверхностной волны на большие расстояния. Результаты расчетов представлены на рис. 4 (правая часть рисунка).
По данным расчетов получается, что в результате взаимодействия волны с донным уступом высота волны увеличилась с 60,6 м до 92,9 м, т. е. в 1,5 раза. Затем высота начала плавно уменьшаться. Данное обстоятельство требует внимательного изучения. С этой целью, прежде всего, рассмотрим вопрос о достоверности расчетов по коду для мелкой воды.
В цитированной работе К. Мейдера [5] приведены результаты численного решения задачи о распространении уединенной волны в океане глубиной 4,550 км. На расстоянии 460 км от скалистого берега, т.е. на левой границе счетной области задавалось граничное условие для скорости частиц воды в виде синусоидального источника волн:
Такой источник формирует уединенную волну высотой 1 м и шириной 140 км, распространяющуюся вправо со скоростью 210 м/с. На правой границе счетной области задавалось условие равенства нулю скорости частиц воды, т.е. берег в виде отвесной скалы.
Также была рассмотрен случай, когда граничное условие имело вид
Это означает, что данный источник формирует цуг волн с амплитудой от минус 1 м до 1 м и шириной 280 км, распространяющихся вправо со скоростью 210 м/с.
В цитированной работе [5] расчеты были проведены по двум вычислительным кодам: SWAN - код для теории мелкой воды и ZUNI - код для двухмерных уравнений Навье-Стокса. Геометрия счетной области имела вид, представленный рис. 5 (левая часть рисунка).
Глубина океана Y = 4,550 км. На расстоянии X = 283,500 км от источника волны начинается скальный уступ с уклоном 1:15. На расстоянии 344,250 км начинается шельф постоянной глубины 500 м.
На рис. 5 (правая часть рисунка), взятом из работы [5], представлена рассчитанная по кодам SWAN и ZUNI поверхность океана на два момента времени, когда в качестве источника волн на левой границе бралось соотношение (3).
Аналогичная задача была рассчитана и по коду ТМВ, разработанному авторами данной работы. Результаты расчетов вполне удовлетворительно согласуются с данными работы [5]. Чтобы в этом убедиться, достаточно сравнить рис. 5 (правая половина рисунка) и рис. 6.
Рассмотрим снова результаты расчетов, приведенных в работе [5]. По этим результатам получается, что уединенная волна (граничное условие (2)) в результате взаимодействия с подводным уступом с уклоном 1:15 увеличила свою высоту с 0,96 м до 1,50 м.
Дальнейшее ее распространение по шельфу до момента взаимодействия с правой границей счетной области (скальной стеной) характеризовалось постоянством высоты волны. Данное обстоятельство отмечено на рис. 5 прерывистой штриховой огибающей линией в верхнем правом углу рисунка - SWAN (shallow water, long wave).
(2)
и = 0,04666sin(0,0047131) м/с.
(3)
ч»
Я ton
0.0035
2as,500m 344)£S0a 419,000л
Рис. 5. Схема счетной области и рассчитанные профили волн на моменты времени 3000 и 3840 с,
согласно работе [5]
7 /<т
0.0030
0.0025
0.0020
0.0015
0.0010
0.0005
0.0000
-0.0005
-0,0010
-0.0015
-0.0020
1 1
time 3000s 3020s
time i
i
i
-- — * i \ / \ i A
/ \ \ jr 1 * / N i у V
/ / / \ \ 1 t V
/ \ V /' T 1
✓ - \ \ I A
v -/ \ I
V
100 200 300 400 К7 km
Рис. 6. Профили волн на различные моменты времени
time = 1440S
time= 1800s time= 2200s
/1
time = 2400s r,
time = 3000s / \
/ ; if 1
^ \ 4 I
/ \ / i 1
/ V f 1 I
l>
_ _ - ‘*4 /'■
r
R, кт
Рис. 7. Профили волн для шельфа с уклоном
Аналогичная задача с граничным условием (2) была рассчитана по вычислительному коду ТМВ для теории мелкой воды. Получено хорошее согласование результатов расчетов с данными из работы [5]: получено увеличение высоты волны с 1 м до 1,5 м и дальнейшее движение волны без изменения высоты. На рис. 7 показано увеличение высоты волны в случае шельфа переменной глубины: глубина шельфа при X = 459 км равнялась нулю.
Таким образом, показано, что результаты расчетов аналогичных задач по вычислительному коду ТМВ и апробированному (эталонному) коду SWAN хорошо согласуются друг с другом. Имея это в виду, обратимся снова к рис. 4, на котором по результатам расчетов по коду ТМВ имеет место быть уменьшение высоты волны при ее движении по шельфу, что якобы противоречит рис. 7. Из геометрических данных, приведенных на рис. 5 (левая часть рисунка), следует, что отношение высоты волны к глубине шельфа в рассмотренной задаче из работы [5] равно 1:500. А в задаче о мегацунами, рассматриваемой в данной работе, это отношение равно 1: 2,5.
Явно видна зависимость высоты волны, распространяющейся по шельфу, от его глубины. Чтобы это проверить, была рассчитана задача в постановке работы [5], но глубина шельфа была принята равной 3 м на расстоянии X = 283,500 км и равной нулю на правой границе X = 459,00 км. Результаты расчетов приведены на рис. 8. Усматривается уменьшение высоты волны при ее движении по шельфу.
Сопоставление результатов, полученных нами и в работе [5], показывает, что предложенная технология может быть применена. В то же время, можно ожидать, что реальные физические процессы будут иметь более сложный характер: возможна существенная эрозия шельфа, и предложенная модель должна быть существенно изменена, чтобы учитывать взаимодействие течения с грунтом, что можно сделать, например, в рамках модели многокомпонентной среды.
Рис. 8. Движение волны над неглубоким шельфом
Выводы
Использование предложенного в данной работе подхода к оценке параметров мегацунами, образующегося от падения астероида в океан, вполне приемлемо как для качественного, так и для количественного описания данного физического явления. Как показали полученные приближенные оценки, последствия от падения каменного астероида диаметром ~1 км могут быть разрушительными для побережья океана.
В работе представлены результаты двухмерных осесимметричных расчетов при вертикальном падении астероида, что является весьма редким случаем. В природе наиболее вероятным является косое падение, что приводит к эффективному увеличению глубины океана, а, следовательно, к дополнительному увеличению доли переданной ему энергии. Очевидно, что при таком ударе интенсивность волны цунами будет зависеть от угловой координаты с максимальной амплитудой по направлению движения ОКО. Однако выполнение 3-х мерных расчетов является весьма трудоемкой задачей и оправдано только тогда, когда можно опираться на рассмотрение реальных случаев падения ОКО в океан.
Работа выполнена при частичной поддержке грантом РФФИ, проект № 07-01-96011, и междисциплинарным интеграционным проектом СО РАН № 2006-113.
Литература
1. Abbott, D. Burckle Abyssal Impact Crater: Did this Impact Produce a Global Deluge? / D. Abbott, L. Burckle, W. Masse et al. // In the Atlantis Hypothesis: Searching for a Lost Land, Helio-topos Publications, St. P. Papmarinopoulous, Ed. - 2007. - P. 179-190.
2. Abbott, D. Impact Ejecta and Megatsunami Deposits from a Historical Impact into the Gulf of Carpentaria, Australia / D. Abbott, E. Tester, C. Meyers et al. // Geological Society of America Annual Meeting, Denver, CO. - Abstacts with Programs. - 2007. - V. 39. - P. 312.
3. Brayant, E. Cosmogenic Mega-Tsunami in the Australia Region: Are They Supported by Aboriginal and Maori Legends / E. Brayant, G. Walsh, D. Abbott // In Myth and Geology, Piccardi L., Masse W.B. (eds). - Geological Society of London Special Publication 273. - London. - 2007. -P.203-214.
4. Аврорин, Е.Н. Численное моделирование взаимодействия частиц кометы Галлея с космическим аппаратом / Е.Н. Аврорин, Н.Н. Анучина, В.В. Гаджиева и др. // Препринт ИПМ имени Келдыша. - 1985. - № 177.
5. Mader Charles L. Numerical Simulation of Tsunamis / L. Mader Charles // Journal of Physical Oceanography. - 1974. - V. 4. - P. 74-82.
Поступила в редакцию 25 января 2008 г.
MATHEMATICAL MODELING OF MAGATSUNAMIS
Using the finite-difference keys (the bidimensional calculating key MPC designed for the gasdy-namic problems solution and the calculating key LWT designed for solution the problems described by the shallow water theory) the authors carried out the computational modelling of a problem about forming and propagation of megatsunamis which are caused by a celestial object coming into the World's water.
Keywords: stroke, asteroid, comet, tsunami, computational modelling, shallow water theory.
Simonenko Vadim Aleksandrovich - Dr.Sc. (Physics and Mathematics), Professor, Acting Scientific, Advisor of the Russian Federal Nuclear Center of All Russian Research Institute of Technical Physics of the name E.I. Zababakhin, city of Snezhinsk.
Симоненко Вадим Александрович - доктор физико-математических наук, профессор, заместитель научного руководителя Российского федерального ядерного центра ВНИИ ТФ,
г. Снежинск.
e-mail: sva@snezhinsk.ru
Skorkin Nikolaj Andreevich - Dr.Sc. (Engineering), Professor, Branch of the South Ural State University in the city of Snezhinsk.
Скоркин Николай Андреевич - доктор технических наук, профессор, филиал ЮжноУральского государственного университета, г. Снежинск.
e-mail: n.a.scorkin@rambler.ru
Elsukov Vasilij Pavlovich - Head of the of the Russian Federal Nuclear Center of All Russian Research Institute of Technical Physics named after E.I. Zababakhin, city of Snezhinsk.
Елсуков Василий Павлович - начальник лаборатории Российского федерального ядерного центра - Всероссийского НИИ технической физики им. акад. Е.И. Забабахина.
e-mail: v.p.elsukov@vniitf.ru
Uglov Aleksandr Sergeevich - Senior Scientific Researcher of the of the Russian Federal Nuclear Center of All Russian Research Institute of Technical Physics named after E.I. Zababakhin, city of Snez-hinsk.
Углов Александр Сергеевич - старший научный сотрудник Российского федерального ядерного центра - Всероссийского НИИ технической физики им. акад. Е.И. Забабахина.
e-mail: a.s.uglov@vniitf.ru
