Метод постановки начальных условий в задаче численного моделирования цунами Текст научной статьи по специальности «Физика»

96

ВМУ. Серия 3. ФИЗИКА. АСТРОНОМИЯ. 2009. № 2

Метод постановки начальных условий в задаче численного

моделирования цунами

М. А. Носов", C.B. Колесов

Московский государственный университет имени М. В. Ломоносова, физический факультет, кафедра физики моря и вод суши. Россия, 119991, Москва, Ленинские горы, д. 1, стр. 2. E-mail: " tiosov@phys.msu.ru

Предложен метод постановки начальных условий в задаче о волнах цунами, согласно которому начальное возвышение водной поверхности в очаге определяется из решения трехмерной задачи в рамках потенциальной теории течения жидкости с учетом вертикальной и горизонтальной деформации дна, а также батиметрии в области источника. Работоспособность метода демонстрируется на примере очагов цунами, произошедших на Центральных Курильских островах в 2006 и 2007 гг.

PACS: 91.30.Nw.

Ключевые слова: волны цунами, численное моделирование, начальные условия, подводное землетрясение, остаточные деформации дна, генерация цунами, потенциал скорости течения, уравнение Лапласа.

Статья поступила 11.06.2008, подписана в печать 19.11.2008.

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

Традиционный способ постановки начальных условий нашел широкое применение в численном моделировании реальных событий, так как он достаточно адекватно воспроизводит основной механизм генерации цунами — вытеснение воды деформацией дна. В настоящей работе предлагается и обосновывается модифицированный способ постановки начальных условий, который физически более корректно описывает процесс генерации цунами. Модифицированный способ сопоставляется с традиционным на примере цунами, вызванных землетрясениями 15.11.2006 и 13.01.2007 в районе Центральных Курильских островов.

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

Ь — горизонтальная протяженность источника, § — ускорение силы тяжести). Во-вторых, при наклонном дне горизонтальные компоненты деформации также приводят к вытеснению воды и, следовательно, вносят вклад в начальное возвышение.

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

быстро затухающей функцией \fchkH, где к — волновое число [3]. Следовательно, движения дна не могут создать на поверхности возмущения с длиной волны А <Н. Для неглубоких землетрясений, и особенно в тех случаях, когда разлом выходит на поверхность дна, пространственное распределение деформаций дна может характеризоваться горизонтальными неоднородностями, масштаб которых меньше глубины океана. Поэтому прямой перенос деформаций дна на водную поверхность приводит к искусственному насыщению спектра цунами коротковолновыми компонентами, которые не существуют в реальности. В численных моделях адекватное описание коротковолновых компонент требует уменьшения шагов по пространству и времени, что резко увеличивает время счета. Резонансный отклик шельфа и мелкомасштабных бухт на коротковолновые компоненты, затрудняет интерпретацию расчетных данных, а в ряде случаев способствует развитию неустойчивости численных схем. Еще раз подчеркнем, что все эти проблемы связаны со спектральными компонентами, которые не свойственны реальным цунами, а искусственно привносятся в модель чрезмерно упрощенным способом постановки начальных условий.

На необходимость учета «сглаживающего эффекта» водного слоя в тех случаях, когда глубина океана сопоставима с горизонтальной протяженностью источника, указывал еще К. Каджаура [6]. На настоящее время нам известны лишь две работы, в которых сделаны попытки учесть этот эффект в расчетах реальных цунами. В работе [7] начальное возвышение рассчитывается с применением аналитической формулы, полученной для случая плоского горизонтального дна. Применение этой формулы для бассейна переменной глубины, очевидно, является не вполне корректным. В работе [8] утверждается, что начальное возвышение рассчитано из решения трехмерного уравнения Лапласа, но описание математической модели, по которой выполнен расчет, отсутствует.

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

АР = 0, (1)

(9Т7 (дг)

дп

дС

2 = 0,

, г = -Н (х,у),

(2) (3)

где И — потенциал скорости течения, п — нормаль к поверхности дна. Искомое начальное возвышение рассчитывается по формуле £(х,у,т) = —g~lFt(x,y,0,т). Если процесс деформации дна оказывается продолжительным (т~1(£#)-1/2), то в начальных условиях задачи распространения наряду с начальным возвышением следует использовать начальное распределение скоростей течения

У(х,у,т) =

1

Н(х, у)

йг ЧР(х,у, г, т).

-Н(х,у)

В том случае когда деформацию дна можно считать мгновенной, эволюционная задача (1)-(3) легко сводится к более простой — статической:

аР = о,

Р = 0, 2 = 0, дР

— ^(т],п), г = -Н(х,у),

(4)

(5)

(6)

где Р = Начальное возвышение теперь опре-

о

деляется из решения задачи (4)-(6) по формуле €(х,у) = рг(х,у, 0).

II. В ноябре 2006 и январе 2007 г. восточнее о. Симу-шир (Центральные Курильские острова) произошли два сильных землетрясения с эпицентрами, локализованными в районе Курило-Камчатского желоба (15.11.2006, ма. = 8.3 и 13.01.2007, = 8.1). Оба сейсмических события вызвали волны цунами, которые наблюдались

по всему Тихому океану и были зарегистрированы множеством береговых и глубоководных станций [8, 9]. По данным удаленных мареографов, ноябрьское цунами характеризовалось большей амплитудой волн по сравнению с январским событием. Ближайшее к очагам цунами побережье Центральных Курильских о-вов не оборудовано мареографами и, кроме того, является необитаемым и труднодоступным регионом. Экспедиционное обследование этого побережья, проведенное летом 2007 г., выявило высоты заплесков до 20 м [10]. Какому из двух рассматриваемых событий соответствовали большие высоты волн, достоверно установить не удалось.

Для расчета остаточных деформаций дна океана, вызванных землетрясением 15.11.2006 (13.01.2007), мы использовали данные о структуре подвижки, которые доступны в цифровом формате на сервере Геологической службы США (ШвЗ). Область разрыва представлена прямоугольником, вытянутым по горизонтали, с размерами 400 х 137.5 км (200 х 35 км). Ориентация плоскости разрыва определяется углом простирания 220.23° (42.36°) и углом падения 14.89° (57.88°). Верхняя кромка области разрыва располагалась под поверхностью дна на глубине 5.9 км (4.6 км), нижняя — на глубине 40.6 км (33.4 км). Область разрыва разбита на 220 (175) элементарных прямоугольных подобластей с размерами 20х 12.5 км (8x5 км), для каждой из которых определен свой вектор подвижки вдоль разрыва. Величина подвижки вдоль плоскости разрыва достигала 8.9 м (20.3 м).

Векторное поле деформации дна г], входящее в формулу (6), мы рассчитывали как суперпозицию вкладов от элементарных подобластей с использованием формул Окада. Результаты расчета представлены на рис. 1. Амплитуда вертикальной деформации дна составила: поднятие 2.7 м (1.9 м), опускание 0.6 м (7.7 м). Горизонтальная деформация дна достигала 3.8 м (3.4 м). В обоих случаях наиболее значительные деформации приурочены к боль-

Рис. 1. Остаточные деформации дна в очагах цунами 15.11.2006 (а) и 13.01.2007 (б). Белыми изолиниями показаны вертикальные деформации (сплошная линия — поднятие, пунктирная — опускание). Цифры у изолиний — величина деформации в метрах. Черные стрелки — вектор горизонтальной деформации. Распределение глубин

показано серым тоном. Шаг изобат 1 км

ВМУ. Серия 3. ФИЗИКА. АСТРОНОМИЯ. 2009. № 2

3

X, км

Рис. 2. Начальное возвышение в источнике цунами 15.11.2006. (а) Сплошные белые линии — поднятие, пунктирные — опускание. Цифры у изолиний — величина отклонения водной поверхности в метрах, (б) Профиль начального возвышения вдоль черной пунктирной линии АВ: 1 — расчет с учетом всех трех компонент вектора деформации дна, 2 — расчет с учетом только вертикальной компоненты деформации, 3 — вертикальная

деформация дна

А Л А б В

/

— 1 1 3 /

1 /

/ 1/

50

100

150 х, км

х, км

Рис. 3. Начальное возвышение в источнике цунами 13.01.2007. Обозначения аналогичны рис. 2

шим глубинам и крутым склонам Курило-Камчатского желоба.

Начальное возвышение определялось из задачи (4)-(6), которая решалась численно методом установления [11]. В расчетах использовались шаги по пространству Ах « Ау « 1000 м, Аг « 200 м. Батиметрия расчетной области была заимствована из цифрового атласа ОЕВСО.

Результаты расчета начального возвышения показаны на рис. 2 и 3. Для события 2006 г. максимальные зна-

чения отклонения водной поверхности (поднятие 2.5 м, опускание 0.5 м) незначительно отличаются от соответствующих характеристик вертикальной деформации. Но для события 2007 г. амплитуда отклонения водной поверхности (поднятие 0.8 м, опускание 4.4 м) оказывается приблизительно в два раза меньше, чем амплитуда вертикальной деформации дна. Такое существенное различие возникает из-за того, что во втором случае ширина области деформации дна сопоставима с глубиной океана.

Особенности «тонкой структуры» начального возвы-

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

Из рис. 2, б видно, что двугорбая «тонкая структура», присущая вертикальной деформации дна в области максимума, не проявляется в смещении водной поверхности. Кроме того, из различия кривых 1 и 2 следует, что горизонтальная деформация, т. е. смещение северо-западного склона желоба в юго-восточном направлении, действительно вносит заметный вклад в начальное возвышение.

Существенное отличие вертикальной деформации дна и смещения водной поверхности наблюдается на рис. 3, б. Это означает, что для цунами января 2007 г. традиционный подход, т. е. перенос остаточных деформаций дна на водную поверхность, приводит к двукратной переоценке амплитуды. Близость кривых 1 и 2 говорит о том, что в данном случае роль горизонтальной деформации оказалась незначительной.

Показательными характеристиками очага цунами являются объем воды, вытесненный деформацией дна, и потенциальная энергия начального возвышения. Эти характеристики мы рассчитывали по следующим формулам: Vq = JJ £(х, у, т) dx dy, Wq = Щ- JJ £2(х, у, т) dx dy, где р — плотность воды. Расчеты проводились с учетом всех трех компонент вектора деформации дна, с учетом только вертикальной компоненты и по традиционной методике, т.е. £, = г]г. Результаты сведены в таблице. Видно, что вклад горизонтальной компоненты деформации дна в вытесненный объем превышает 30% для

события 2006 г. и почти достигает 20% для события 2007 г. Такой заметный вклад связан с обширностью площади, охваченной горизонтальной деформацией. По абсолютной величине вытесненного объема воды первое событие является более мощным. Интересно, что оценка энергии дает совсем иной результат: более мощным оказывается цунами января 2007 г. Значение энергии довольно чувствительно к методу расчета, что особенно заметно для события 2007 г. Обратим внимание на тот факт, что традиционная методика завышает значение энергии почти в 2 раза.

Заключение. Отметим, что начальное возвышение, рассчитанное предложенным методом, конечно, сглаживает мелкомасштабные пространственные неоднородности деформации дна, но оно содержит в себе компоненты с длиной волны А > ~Я. Такие волны, в отличие от длинных (Х^>Н), подвержены фазовой дисперсии, а следовательно, их распространение, особенно на трансокеанские расстояния, не может быть адекватно описано теорией длинных волн. Тот факт, что второе событие (13.01.2007), имея большую амплитуду и энергию в источнике, проявилось как относительно слабое на удаленных мареографах, можно объяснить «размазыванием» энергии в пространстве из-за фазовой дисперсии.

Работа выполнена при финансовой поддержке РФФИ (грант 07-05-00414).

Списож литературы

1. Titov V.V., Gonzalez F.I. // NOAA Technical Mémorandum ERL PMEL-112. 1997.

2. Kowalik Z., Knight W., Logan T., Whitmore P. // Science of Tsunami Hazard. 2005. 23, N 1, P. 40.

3. Левин Б.В., Носов M.А. Физика цунами и родственных явлений в океане. М., 2005.

4. Gisler G.R. // Annu. Rev. Fluid Mech. 2008. 40. P. 71.

5. Okada Y. // Bull. Seis. Soc. Am. 1985. 75, N 4. P. 1135.

6. Kajitira К. Il Bull. Earthquake Res. Inst. Tokyo Univ. 1963. 41, N 3. P. 535.

7. Tanioka Y, Seno T. // Geophys. Res. Lett. 2001. 28, N 17. P. 3389.

8. Rabitiovich A.В., Lobkovsky L.I., Fine I.V. et al. Il Adv. Geo-sciences. 2008. 14, P. 105.

9. Tanioka Y, Hasegawa Y, Kuwayama T. Il Adv. Geo-sciences. 2008. 14, P. 129.

10. Левин Б.В., Кайсмренко В.M., Рыбин А.В. и др. // Докл. РАН. 2008. 419, № 1. С. 118.

11. Годунов С.К., Рябенький B.C. Разностные схемы (введение в теорию). М., 1973.

Объем воды Ц>, вытесненный деформацией дна, и потенциальная энергия начального возвышения

а* b** с***

Дата Vo. км3 Щ), 1014 Дж Vo. км3 Го, 1014 Дж Vo, км3 Го, 1014 Дж

15.11.2006 8.97 1.00 6.09 0.85 6.08 0.97

13.01.2007 ^6.37 1.23 ^5.20 1.36 ^5.20 2.36

* Начальное возвышение, рассчитанное с учетом всех трех компонент вектора деформации дна. ** Начальное возвышение, рассчитанное с учетом только вертикальной компоненты вектора деформации дна. *** Начальное возвышение, приравненное к вертикальной деформации дна.

Method of specification of the initial conditions for tsunami numerical simulation M. A. Nosov ' , S. V. Kolesov

Department of Marine and Inland Water Physics, Faculty of Physics, M. V. Lomonosov Moscow State University, Moscow 119991, Russia. E-mail: " nosov@phys.msu.ru.

A new method of initial conditions specification in the tsunami problem is suggested. The initial elevation of water surface is determined from the solution of the 3D problem in the framework of potential theory taking into account both horizontal end vertical components of the bottom deformations and the distribution of depths in the vicinity of the source. The tsunamigenic Central Kuril Islands earthquakes of 2006 and 2007 are taken as examples to demonstrate efficiency of the new method. PACS: 91.30.Nw.

Keywords: tsunami waves, numerical simulation, initial conditions, bottom earthquake, permanent bottom deformations, tsunami generation, fluid velocity potential, the Laplace equation. Received 11 June 2008.

English version: Moscow University Physics Bulletin 2(2009).

Сведения об авторах

1. Носов Михаил Александрович — д. ф.-м. п., профессор; тел.: 939-36-98.

2. Колесов Сергей Владимирович — мл. научи, сотр.; тел.: 939-36-98, e-mail: kolesov@phys.msu.ru.