CC BY
15
ОБЗОР
ФИЗИКА ЗЕМЛИ, АТМОСФЕРЫ И ГИДРОСФЕРЫ
Численное моделирование цунами и рельеф дна
Е.А. Куликов1,0, В. К. Гусяков2,3, А. А. Иванова1, Б. В. Баранов1
1 Институт океанологии имени П. П. Ширшова РАН.
Россия, 117997, Москва, Нахимовский пр., д. 36.
2 Институт вычислительных технологий Сибирского отделения РАН.
3 Институт вычислительной математики и математической геофизики Сиб. отд. РАН. Россия, 630090, Новосибирск, пр. Академика Лаврентьева, д. 6.
E-mail: а[email protected] Статья поступила 31.05.2016, подписана в печать 24.06.2016.
В работе рассматривается эффект влияния качества данных батиметрии на точность расчетов волнового поля цунами. Представлен обзор, посвященный истории развития численных моделей распространения цунами. Особое внимание уделено моделям рельефа дна Мирового океана. Показано, что современные цифровые карты батиметрии, например GEBCO, не обеспечивают адекватное воспроизведение рельефа дна в численных моделях распространения волн, что может приводить к существенным ошибкам при оценке максимальных заплесков цунами на побережье.
Ключевые слова: цунами, численная модель, рельеф дна, данные батиметрии.
УДК: 551.466.62. PACS: 91.30.Nw.
Введение
Развитие методов оперативного прогноза цунами и получение долгосрочных оценок цунами-риска в значительной степени основывается на использовании технологии компьютерного моделирования. В оперативном прогнозе численные модели необходимы для заблаговременной оценки высот цунами на побережье по данным сейсмических наблюдений и расчетов параметров очага землетрясения. В задаче долгосрочного прогноза (цунамирайонирования побережья) компьютерное моделирование используется в сценарных расчетах для получения ожидаемых высот волн при землетрясениях различной магнитуды.
Современные системы инструментального мониторинга цунами (например, система глубоководных буев DART [1]) способны заблаговременно обнаруживать сигнал цунами на значительном удалении от берега. Однако чтобы определить время прихода и оценить высоту волны на конкретном участке побережья, требуется использование эффективных методов расчета распространения цунами в глубоком океане и на шельфе.
Не менее важной задачей является цунамирайо-нирование побережья, т. е. предварительная оценка долгосрочной цунами-опасности для конкретных участков берега. Современный подход к решению этой задачи базируется на методике PTHA (Probabilistic Tsunami Hazard Assessment) [2]. Он использует результаты численных расчетов распространения цунами от модельных очагов, приуроченных
к основным сейсмоактивным (цунамигенным) зонам акватории океана.
Эффективность применения численных моделей для решения этих главных научно-практических задач проблемы цунами зависит не только от корректности используемых физико-математических моделей и реализующих их численных алгоритмов, но и от ряда других факторов, одним из которых является точность и детальность батиметрических данных, используемых для построения расчетных сеток.
В настоящем исследовании представлены результаты анализа влияния качества данных батиметрии на точность расчетов волнового поля цунами.
1. Численное моделирование генерации и распространения цунами
Наиболее распространенной и широко применяемой физической моделью цунами является модель так называемых длинных волн, распространяющихся в слое однородной несжимаемой жидкости, лежащем на жестком дне. Поскольку для типичных сей-смогенных цунами, возбуждаемых в морских и океанических бассейнах очагами подводных землетрясений, длина волны много больше глубины бассейна, для математического описания распространения можно использовать приближение мелкой воды, в котором пренебрегается вертикальными ускорениями и соответственно вертикальными скоростями частиц жидкости. Как следствие, использование приближения гидростатики позволяет в трехмерных уравнениях движения исключить зависимость горизонтальных скоростей течения от глубины и перейти к интегрированным по вертикали двумерным урав-
2 ВМУ. Физика. Астрономия. № 6
нениям, например относительно уровня и полных потоков [3, 4]. Эта процедура значительно упрощает задачу моделирования распространения волны цунами в океане по сравнению с использованием полных трехмерных моделей (типа Навье-Стокса).
Численные методы для расчета распространения цунами на отдельных участках акватории океана начали использоваться с середины 1960-х гг., с момента появления в исследовательских центрах достаточно крупных ЭВМ (типа IBM-360, CDC 6400 и БЭСМ-6), пригодных для такого рода расчетов. Подробный анализ начального этапа развития численного моделирования цунами можно найти в работах [5, 6]. Хронологически первые работы в этой области были выполнены в Японии, в них были предприняты попытки промоделировать проявление Чилийского цунами 1960 г. в Токийском заливе [7], а также воспроизвести основные особенности трансокеанского распространения этого цунами [8]. В этих расчетах использовался программный код, созданный для моделирования приливных волн и реализующий численное решение нелинейной системы мелкой воды, записанной для вращающейся сферы. Расчеты трансокеанского распространения Чилийского цунами 1960 г. выполнялись на ламповой ЭВМ первого поколения IBM 704 на батиметрической сетке 25 х 30 узлов с пространственным шагом в 5 угловых градусов (порядка 550 км). Историю развития компьютерных моделей, специально разработанных для воспроизведения исторических событий возникновения и распространения волн цунами, по-видимому, следует начинать с работы И. Аиды [9]. Целью исследования была попытка воспроизвести основные особенности проявления Ниигатского цунами 16.06.1964 г. в Японском море и цунами Токачи-оки 16.06.1968 г. вблизи Хоккайдо. В модели применялась численная конечно-разностная схема, аппроксимирующая уравнения мелкой воды, без учета сил Кориолиса и донного трения, записанная в прямоугольной системе координат. При моделировании Ниигатского цунами 1964 г. в качестве начального возвышения свободной поверхности моря в источнике использовались данные повторной батиметрической съемки очаговой области, выполненной сразу после землетрясения [10, 11] и выявившей значительные постсейсмические деформации дна. Таким образом, в модели предполагалось, что подвижки дна происходили мгновенно и форма начального возвышения морской поверхности повторяет форму донных смещений. Дальнейшая эволюция этого возвышения определялась численно с учетом реальной конфигурации береговой линии и рельефа морского дна. Расчет выполнялся в сеточной области размерностью 20 х 30 узлов с шагом по пространству, равным 10 км. Глубины в узлах сетки задавались путем ручной оцифровки морских навигационных карт. Расчетное распределение высот вдоль берега сравнивалось с измеренными
заплесками, и было найдено их удовлетворительное совпадение в отношении максимальных значений и характера убывания высот при удалении вдоль берега от очаговой области.
Аналогичный подход был применен для моделирования цунами Токачи-оки 1968 г. В расчетах использовалась сетка размерностью 30 х 25 узлов с шагом по пространству 20 км. Для этого цунами характер подвижек дна был известен лишь в самых общих чертах, однако в распоряжении оказалось довольно большое количество мареографных записей на побережье. Было рассмотрено несколько вариантов распределения подвижек дна, задаваемых достаточно произвольно внутри очаговой области в форме эллипса, положение и размеры которого выбирались на основе данных об афтершоках и результатов расчетов обратных рефракционных диаграмм [12]. Далее путем сравнения рассчитанных мареограмм с наблюдаемыми была выбрана модель деформации дна с наилучшим соответствием наблюдаемым данным.
В США первой работой по численному моделированию цунами была работа [13], в которой была построена численная модель разрушительного Аляскинского цунами 28.03.1964 г. и промоделировано распространение этого цунами в северо-восточной части Тихого океана. В ней авторы применили нелинейную (первого порядка) систему уравнений движения в приближении мелкой воды. Эти уравнения преобразовывались в конечно-разностную форму и решались численно на сетке с шагом порядка 0.5°, покрывающей область очага и прилегающие к ней участки океана. Были вычислены волновые профили в нескольких точках вблизи очага, их совпадение с реальными данными не проверялось, поскольку отсутствовали записи этого цунами на глубокой воде. В одном пункте побережья вычисленные колебания сравнивались с наблюдениями и было обнаружено их удовлетворительное соответствие. В последующей работе [14] эта модель была улучшена путем введения поправок на кривизну Земли и записи основных уравнений в сферических координатах. Позднее она применялась для численного моделирования цунами от Чилийского 1960 г. и некоторых других сильнейших землетрясений, характер подвижек дна при которых более или менее достоверно известен.
В России первые работы по численному моделированию цунами были выполнены в Сибирском отделении РАН, в них исследовалось поведение сначала гипотетических, а потом и реальных исторических цунами в районе Южных Курил [15, 16, 17] в рамках нелинейной системы мелкой воды, записанной в декартовых координатах и поршневой модели возбуждения цунами.
Уже первые попытки применения численных методов для моделирования реальных цунами показали, насколько плодотворен такой подход. В дальней-
шем появилось большое число программных пакетов и комплексов для моделирования цунами, основанных на существенно более сложных моделях, включающих эффекты вращения Земли, амплитудной и частотной дисперсии, донного трения. При этом достаточно рано (к началу 1990-х гг.) в международном цунами-сообществе была осознана необходимость верификации и валидации численных моделей и пакетов программ, используемых для расчетов цунами. В США были организованы и проведены две специализированные школы-семинара, где была разработана и принята система тестов, которой должен удовлетворять каждый программный комплекс, используемый для расчетов цунами [18, 19].
В течение первых двух десятилетий (60-е и 70-е гг. прошлого века) широкое использование компьютерных моделей для изучения цунами сдерживалось не только недостатком оперативной памяти и быстродействием имеющихся в распоряжении исследователей ЭВМ, но также недостатком информации о структуре очага цунами. Данные о косейсмических смещениях дна в очаговой области, полученные путем повторных батиметрических промеров, были доступны лишь для очень немногих цунамигенных землетрясений. Для оценки положения и размеров очага использовался метод обратных рефракционных диаграмм, а начальная величина смещения уровня океана оценивалась по корреляционным формулам, устанавливающим связь между магнитудами землетрясения и цунами, которые имели очень приблизительный характер. С середины 1970-х гг., однако, стали появляться систематические данные о механизмах очагов сильных землетрясений, получаемые по методике CMT (Centroid Moment Tensor) [20]. C 1976 г. такие определения стали рутинными для всех землетрясений, начиная с магнитуды Mw = 5.5 [21]. Это дало возможность использовать полученные решения по механизмам подводных цунамигенных землетрясений для вычисления остаточных смещений дна в очаговой области.
Первой и достаточно знаковой работой в этом отношении была работа И. Аиды [22], в которой оценивалась надежность получения информации об очаге цунами на основе данных о механизме очага подводного землетрясения. В ней были построены компьютерные модели для пяти цунами, происшедших в 1952-69 гг. вблизи восточного побережья о-вов Хоккайдо и Хонсю (4.03.1952 г. Токачи-оки,
16.05.1968 г. Токачи-оки, 12.06.1968 г. Ивате-оки,
12.08.1969 г. Шикотан и 17.06.1973 г. Немуро-оки), для которых имелись данные об ориентации и размерах плоскости разрыва и величине подвижки. Остаточные смещения дна в очаговой области вычислялись на основе формул, полученных в работе [23]. Смещения использовались в качестве начальных условий для задачи возбуждения цунами. Дальнейшая эволюция начального возвышения воды рассчитывалась на основе линейной модели мелкой
воды, которая решалась численно на последовательности вложенных сеток, с шагом по пространству, меняющимся от 10 км на глубокой воде до 312.5 м в окрестностях береговых мареографных точек. Сравнение расчетных амплитуд с измеренными проводилось только для головной (первой) волны, поскольку предполагалось, что последующие колебания могли содержать значительные расчетные погрешности из-за недостаточной точности аппроксимации прибрежной батиметрии и конфигурации береговой линии. Точность воспроизведения в расчетах амплитуд первого колебания оказалась в пределах 50-100%, что в те времена считалось достаточно хорошим результатом, соответствие расчетных времен вступлений наблюдаемым оказалась в пределах 5-10 мин. Для каждой очаговой модели вводился так называемый корректирующий фактор К, показывающий среднее (по 4-6 мареографным точкам) отношение наблюдаемых и расчетных амплитуд. Для большинства рассмотренных очаговых моделей он оказался больше единицы, со средним значением 1.4 (с вариациями в пределах 0.82-4.45), т. е. расчетные амплитуды были в среднем меньше наблюдаемых. Общий вывод этой работы — при наличии данных по механизму очага и измеренной инструментально величины сейсмического момента землетрясения расчетные смещения дна в очаговой области вполне могут быть использованы для построения достаточно реалистичных численных моделей исторических цунами.
В дальнейшем использование расчетных смещений дна в очаговой области в качестве аппроксимации очага цунами (так называемая «поршневая модель» возбуждения цунами) стало стандартной схемой для построения компьютерных моделей исторических и гипотетических цунами. Такие смещения обычно рассчитываются по формулам, выведенным И. Окадой в работе [24]. Они представляют собой статические (остаточные) смещения поверхности однородного упругого полупространства, возникающие под действием внутреннего пространственного источника дислокационного типа, который задается шестью параметрами: длиной и шириной плоскости разрыва, ее углом падения и направлением подвижки, азимутом простирания, глубиной верхнего (или нижнего) края разрыва и величиной подвижки (дислокации). Справедливости ради стоит отметить, что работа [24] явилась лишь завершающей в целой серии работ, посвященных решению этой задачи. В различных подходах к постановке и методам получения решения эта задача рассматривалась в работах [25, 26, 27, 28, 29]. В России формулы, аналогичные выведенным в [24], были получены в работе [30]. При некотором различии в форме записи итоговых выражений для компонент смещений многочисленные проверки показали полную идентичность результатов расчетов по формулам,
полученным в этих работах, при любых наборах параметров модельных очагов.
До начала 2000-х гг. в качестве модели источника цунами использовалась однородная подвижка вдоль плоскости разрыва, моделирующей очаг землетрясения, смещения от которой рассчитывались по формулам, полученным И. Окада [24]. Развитие систем сейсмических наблюдений в направлении использования широкополосных цифровых приборов и применение методов теории обратных задач для восстановления тонкой структуры очага позволило во многих случаях использовать данные о реальном распределении подвижки по плоскости разрыва [31], при этом формулы И. Окады применяются для каждой из элементарных подвижек, на которые разбивается одна или несколько плоскостей разрыва в очаге реального землетрясения. По такой технологии в настоящее время выполняется большинство работ по численному моделированию реальных исторических цунами (обзор таких работ можно найти, например, в [4, 32, 33]) и проводится анализ свойств косейсми-ческих деформаций в очагах цунами [34, 35, 36].
В настоящее время качество численного моделирования процессов генерации и распространения цунами достигло весьма хорошего уровня. Глубоководные записи цунами, получаемые датчиками системы DART [37], как правило, воспроизводятся в расчетах с очень хорошей точностью с различиями по амплитуде первой волны в пределах 5-10% [38, 39]. Также обычно не составляет особых проблем воспроизвести в расчетах (вплоть до выполняемых в режиме реального времени) приборные записи головной волны цунами, полученные в отдельных бухтах, в особенности расположенных в дальней от очага зоне [40].
Однако адекватное воспроизведение в расчетах мареографных записей и заплесков на берегу в ближней от очага зоне по-прежнему является большой проблемой, возникающей даже в тех случаях, когда расчетные мареограммы хорошо воспроизводят записи глубоководных датчиков. Естественно, что моделирование процесса распространения цунами на мелководье и выхода волны на берег является существенно более сложной задачей по сравнению с моделированием распространения на глубокой воде. Здесь необходим учет амплитудной дисперсии, что возможно только в рамках нелинейной системы, донного трения, а также применение специальных численных методов для моделирования выхода волны на сухой берег. В современных моделях расчета наката все эти факторы обычно учитываются, однако проблема точного воспроизведения колебаний уровня на мелководье и максимальных заплесков на берегу остается.
Существенное (на десятки, а иногда на сотни процентов) различие расчетных высот от измеренных может быть связано не только с упомянутыми выше проблемами задания очага цунами, эффек-
тивностью численного алгоритма и детальностью используемой сетки, но также и с точностью аппроксимации реальной батиметрии морского дна. Данная работа посвящена оценке качества современных цифровых батиметрических карт Мирового океана и оценке возможных ошибок, возникающих при численном моделировании цунами из-за неадекватности используемой батиметрии.
2. Данные по батиметрии Мирового океана
В ранних работах по моделированию цунами, выполненных в 60-е и 70-е гг. прошлого века, сеточные массивы батиметрии создавались, как правило, самими исследователями путем ручной оцифровки морских навигационных карт [10]. Первым широкодоступным массивом глобальной батиметрии стал появившийся в 1988 г. массив ETOPO5, содержащий округленные до целых метров высоты суши и глубины Мирового океана на сетке с шагом 5 угловых минут. Массив был создан в США путем интеграции топографических и батиметрических данных из различных источников. Фактическая точность значений глубин варьировалась от нескольких метров до 100-150 м в областях со слабым покрытием судовыми треками [41].
Международная программа картирования морского дна (GEBCO — General Bathymétrie Chart of Oeeans), поддерживаемая Международной гидрографической организацией (IHO — International Hydrologieal Organization) с 1903 г., выпустила первую цифровую версию атласа GEBCO в 1994 г. [42]. Цифровая версия была подготовлена на основе опубликованного в 1978 г. пятого издания обычного атласа GEBCO, показывающего рельеф дна Мирового океана в виде изолиний глубин и отдельных точечных промеров и, по сути, представляла собой отсканированные копии контурных батиметрических карт GEBCO.
Следующим важным шагом стало появление в 1997 г. глобального сеточного массива высот и глубин с шагом в 2 угловые минуты [43]. Источником данных были результаты высокоточных спутниковых измерений высот рельефа суши и уровня свободной поверхности океана, затем зарегистрированные гравитационные аномалии по весьма сложному алгоритму пересчитывались в функцию, описывающую форму рельефа дна. Этот алгоритм хорошо работает для больших (более 1 км) абиссальных глубин, но на мелководье (в окраинных морях, на материковом склоне и шельфе) дает большие относительные погрешности в оценках глубин, достигавшие в ряде случаев 50-100%. Основным источником этих погрешностей являются неопределенности в толщине слоя рыхлых донных осадков.
Сеточный массив глубин GEBCO с разрешением в одну угловую минуту появился в 2008 г. (GEBCO_08). Глубины в узлах регулярной сетки рассчитывались путем интерполяции оцифрованных
изолиний глубин и данных точечных промеров, содержавшихся в юбилейном издании цифрового атласа GEBCO, посвященном столетию проекта [44].
В марте 2015 г. был выпущен глобальный сеточный массив GEBC0_2014, содержащий высоты суши и глубины дна в узлах 30-секундной сетки. Эта база GEBCO включает в себя данные из большого числа источников: это и национальные и региональные организации (Национальное управление океанических и атмосферных исследований — N0AA, Европейский научный фонд — ESF и др.), и научно-исследовательские институты, и отдельно взятые экспедиции, и измерения. Традиционно данные и карты GEBCO содержали в себе лишь батиметрию наиболее глубоких акваторий Мирового океана, начиная с глубин 200 м и более, недостаточно детально представленных на морских навигационных картах. Сейчас основная работа ведется над повышением качества данных в мелководных районах океана. Для этого используются электронные морские карты (ENCs), коллективно подготовленные государствами-членами Международной гидрографической организации (IH0). Многие гидрографические службы и гидрографические организации уже добавили значительное количество батиметрических данных своих прибрежных зон. Это привело к существенному улучшению качества рельефа некоторых мелководных участков. Благодаря этим данным GEBCO обеспечивает более полную батиметрическую модель всего Мирового океана [45].
Примером современной глобальной модели рельефа поверхности Земли с разрешением в 1 мин является массив ET0P01. Для его создания использованы данные многочисленных глобальных и региональных источников. Отличительной особенностью этой модели является то, что она была создана в двух вариантах: версия «Ice Surface», которая отображает поверхность ледяного покрова Антарктики и Гренландии, и версия «Bedrock», которая отображает коренную породу под двумя ледниковыми щитами [46].
Модель рельефа, основанная на региональных источниках батиметрических данных, разработана в рамках проекта IBCA0 (The International Bathy-metric Chart of the Arctic 0cean) и CRM (Coastal Relief Model). Целью проекта IBCA0 является создание массива цифровых данных, содержащего все доступные батиметрические данные к северу от 64° с.ш., для использования в работах, требующих детального и точного знания глубины и формы дна Северного Ледовитого океана. IBCA0 версии 3.0 включает в себя улучшенные данные, собранные арктическими странами, данные попутных измерений рыболовных судов, данные, полученные с помощью подводных лодок ВМС США и научно-исследовательских судов разных стран [47]. Сетка, построенная с помощью усовершенствованного алгоритма гридинга, имеет разрешение 500 м, что позволя-
ет различить гораздо больше деталей арктического морского дна, нежели при использовании IBCAO версии 1.0 (2.5 км) и версии 2.0 (2.0 км).
Национальный центр информации об окружающей среде (NCEI — National Centers for Environmental Information NOAA) создает цифровые модели рельефа побережья (DEM — digital elevation model) и батиметрии прилегающих участков морского дна (DBM — digital bathymetric model) с разрешением 3 угловые секунды для береговых районов всех штатов США, включая Аляску. Батиметрические, топографические и береговые данные, используемые в этом массиве, заимствованы из различных источников, таких как федеральные и местные органы власти, научные учреждения и частные компании. Эта объединенная цифровая модель рельефа прибрежной области США, подготовленная NCEI, с разрешением 3 угловые секунды впервые обеспечила полное представление о прибрежной зоне США, объединяя батиметрию открытых акваторий с топографией суши. Данные этого массива в настоящее время успешно используются в системе оперативного моделирования цунами, разработанной Тихоокеанской морской лабораторией (PMEL — Pacific Marine Environmental Laboratory) [48].
Перечисленные массивы батиметрии являются открытыми и могут свободно использоваться для научных и образовательных целей. Вместе с тем многие страны имеют в своем распоряжении существенно более детальные данные о рельефе морского дна с пространственным разрешением до 1 с и менее, получаемые в результате геофизических и инженерных исследований прибрежных территорий, входящих в их зону экономической ответственности. Еще более детальными массивами цифровой батиметрии (с разрешением до 1 м и менее) располагают государственные и частные компании, ведущие разведку и добычу нефти и газа на шельфе. Однако такая информация обычно не выкладывается в открытый доступ; более того, в ряде стран она вообще относится к категории закрытой, доступ к таким данным может быть получен лишь на основе специальных договоров и соглашений.
3. Оценка качества батиметрических данных в северо-восточной части сахалинского шельфа
Батиметрическая съемка на шельфе и склоне северо-восточной части о. Сахалин (рис. 1, 2, а) проводилась в рамках трех международных проектов: Kurile-Okhotsk Marine Eхperiment (KOMEX), Hydro-Carbon Hydrate Accumulations in the Okhotsk Sea (CHAOS), Sakhalin Slope Gas Hydrates (SSGH) [49, 50, 51]. В дальнейшем по ходу изложения мы будем использовать сокращение KCS при упоминании этих трех проектов. Полученные данные промеров были использованы для составления батиметрической карты с разрешением 0.5 угловых минут (шаг по долготе Дх « 600 м, шаг по меридиану
Рис. 1. Карта северо-восточной части Тихого океана, построенная на основе данных GEBC0_2014, версия 20141103. Прямоугольником отмечен район, где выполнялась батиметрическая съемка в рамках проектов KCS. Нанесены изолинии косейсмической деформации дна, вызванной землетрясением 15 ноября 2006 г. (M = 8.3). Изолинии проведены через 0.2 м: красным цветом обозначено поднятие, синим — опускание дна. На карте красным треугольником отмечено положение глубоководного регистратора цунами DART 21414
Ду « 920 м), представленной на рис. 2, б, а также для сопоставления с картой, построенной по данным ОББСО (рис. 2, в).
Следует отметить, что карта на рис. 2, б (КСБ) выглядит «естественнее», чем карта ОББСО (рис. 2, в): на ней отчетливо выделяется граница сахалинского шельфа, а склон отмечен сгущением параллельных изобат. Различие между глубинами на картах колеблется в диапазоне от — 222 до 661 м. Среднеквадратическое значение разницы глубин составляет 126 м (рис. 2, г), при том что среднее значение глубины в пределах полигона около 800 м. Такие значительные отличия в батиметрии должны приводить к значимым расхождениям высот цунами на побережье. К сожалению, нам недоступна более точная (чем ОББСО) батиметрическая карта расчетной области. Тем не менее можно попытаться оценить неопределенности, вносимые неточным знанием рельефа дна, сравнив волновое поле цунами, полученное в расчетах только по данным ОББСО и по данным со «вставкой» данных батиметрии КСБ.
4. Численный эксперимент по оценке влияния качества батиметрии на распространение волн цунами
Оценить влияние ошибок в задании рельефа дна на результаты модельных расчетов довольно сложно. Общепринятой методики для таких оценок не существует. Фактически искажение волнового поля, рассчитанного с использованием неточных батиметрических данных, можно интерпретировать как эффект рассеяния волн цунами на неоднородностях дна, заданных в виде массива разниц глубин модельного и истинного рельефа. Очевидно, что влияние рассеяния носит интегральный («накопительный») характер — чем дальше от источника расположен наблюдатель, тем сильнее искажение волнового поля. Наибольший интерес представляет оценка точности расчетов максимальных высот цунами вблизи берега. Именно эти величины обычно используются при построении карт цунамирайонирования.
Восточное побережье о. Сахалин относится к зонам умеренной цунамиопасности. Согласно исследованию [52], ожидаемая высота волн может дости-
Рис. 2. а — карта эхолотных промеров, выполненная в рамках проектов КСБ; б — карта рельефа шельфа и склона северо-восточной части о. Сахалин по данным батиметрической съемки КСБ; в — карта рельефа дна, построенная по данным ОЕБСО; г — разница глубин между картами б и в. Все значения глубины
приведены в метрах
гать 1.5 м. Известно [53], что для побережья Охотского моря наибольшую опасность представляют удаленные очаги цунами, расположенные в наиболее сейсмоактивных зонах Тихого океана. Несмотря на «экранирующий» эффект Курильской гряды, существенная часть волновой энергии цунами может проникать в Охотское море, вызывая значительные колебания уровня вдоль побережья Сахалина.
В настоящем исследовании мы попробуем оценить влияние качества батиметрических данных на оценки максимальных высот цунами на восточном побережье Сахалина на примере события 15 ноября 2006 г. (М = 8.3) в районе центральных Курильских островов. Этот выбор объясняется тем, что параметры этого очага хорошо известны по литературным источникам, при этом инструментально подтверждено, что волны цунами проникли в Охотское море и достигли побережья Сахалина [54]. Фактически это дает нам возможность сравнить модельные расчеты с имеющимися данными измерений, при этом оценить разброс оценок высоты цунами на восточном побережье Сахалина при учете и без участка батиметрии КСБ.
Численный эксперимент выполнялся на сетке с разрешением 0.5', 4801 х 2761 узлов, покрывающих северо-западную часть Тихого океана, включая Охотское море (рис. 1). Цифровой массив глубин был сформирован на основе данных батиметрии
0ЕБС0_2014. Для расчетов использовалась модель, являющаяся разновидностью хорошо известной программы ТиЫАМ1 для численного расчета распространения волн цунами [6], в которой реализована конечно-разностная аппроксимация линейных уравнений мелкой воды в сферических координатах [55].
Для расчета параметров источника цунами 2006 г. использовалась модель иБОБ [56], в которой были рассчитаны смещения вдоль плоскости разрыва прямоугольной формы. Эти смещения в дальнейшем были пересчитаны в вертикальные деформации морского дна с помощью модели [24]. В рамках приближения мелкой воды, используемого обычно для расчета цунами, начальные отклонения поверхности океана в области источника полностью совпадают с подвижками морского дна. Однако если характерные горизонтальные масштабы косейсми-ческих деформаций сопоставимы с глубиной океана, гидростатическое приближение может привести к заметным ошибкам. Поэтому в качестве начального возмущения поверхности океана мы использовали «негидростатическое» решение уравнения Лапласа [57]. По сравнению с деформацией дна оно оказывается сглаженным. В частности, максимальное отклонение свободной поверхности для очага 2006 г. уменьшилось с 2.7 до 1.9 м [58, 59].
На рис. 3 показано сравнение наблюдаемых и рассчитанных по модели колебаний уровня океана
DART21414 — Модель—"
- Г / / / / Л
1 J (V \\ \ \ \ \ Л/'тГ
о \ \ J 1 J \ / •ьЛ
1 \
сЗ К
О о Л
о
Z
Время после толчка, ч
Рис. 3. Запись колебаний уровня моря на станции DART 21414, вызванных цунами 15 ноября 2006 г. (черная сплошная линия). На графике нанесена кривая, соответствующая модельным расчетам волнения (красная штриховая линия)
в точке расположения глубоководного датчика DART 21414.
Видно, что приход волны цунами на записи глубоководного датчика, расположенного на расстоянии от источника примерно 1600 км, отмечается через
2 ч после землетрясения. Модель неплохо воспроизводит как время вступления сигнала, так и максимальную высоту первого гребня. Прибор расположен на глубине 5375 м приблизительно в 270 км к югу от о. Амчитка (Алеутские о-ва). Такое соответствие между зарегистрированной и модельной волнами говорит о хорошем качестве заданных в модели начальных условиях (очага цунами).
В настоящем исследовании оценка влияния качества данных батиметрии при определении цуна-миопасности основываются на сравнении значений максимальных высот волн вдоль северо-восточного побережья Сахалина, полученных в результате численного моделирования события 15 ноября 2006 г. с использованием данных батиметрии 0БВС0_2014, а также с учетом «вставки» рельефа, построенного по данным промеров КСБ в прямоугольной области, показанной на рис. 1. Чтобы оценить искажение волнового поля цунами из-за рассогласования локальной батиметрии в этом прямоугольнике, при-
Рис. 4. Записи колебаний уровня моря, рассчитанные для трех пунктов побережья северо-восточной части Сахалина по модели батиметрии 0БВС0_2014 модели батиметрии 0БВС0_2014 (красная штриховая линия) и комбинированной батиметрии с учетом промеров КСБ (синяя сплошная линия)
летающем к северо-восточному побережью Сахалина, была сгенерирована комбинированная сетка глубин, включающая данные GEBCO_2014 и промеры KCS. Чтобы избежать появления скачков глубины, была выполнена интерполяция значений глубин, обеспечивающая плавный переход к батиметрии GEBCO в полосе шириной 0.5°, охватывающей полигон KCS.
На рис. 4 приведены графики колебаний уровня моря, полученные в результате расчета волнового поля цунами 15 ноября 2006 г. за 5-часовой промежуток времени в трех точках на побережье Сахалина (рис. 5) для двух моделей батиметрии — GEBCO_2014 и комбинированного рельефа, включающего данные GEBCO_2014 и KCS. Записи колебаний уровня в самой южной точке №1 для этих двух моделей оказываются весьма близки. Отличие становится заметным для точек 2 и 3, расположенных непосредственно напротив полигона KCS. Наиболее очевидно расхождение записей колебаний уровня для точки 3. Обращает внимание, что, если для комбинированного рельефа дна (включающего батиметрию KCS) максимальная высота цунами проявляется для первой волны, для модели GEBCO_2014
максимум достигается примерно через 1.5 ч после прихода первой волны.
На рис. 5, а представлена карта распределения максимальных значений высоты цунами. Хорошо видно, что наибольших значений высота волн достигает в окрестности северо-восточной части сахалинского побережья. Следует подчеркнуть, что именно в этом районе сосредоточены объекты топливно-энергетического комплекса (ТЭК). Справа на рис. 5, а показан профиль распределения высот цунами вдоль побережья Сахалина (на изобате 10 м). Характерные высоты цунами колеблются в пределах от 0.5 до 1 м, достигая минимальных значений на восточном побережье п-ва Шмидта.
На рис. 5,б показано распределение модуля разницы высот цунами, полученных при расчете волнового поля по модели GEBCO и комбинированной модели батиметрии. Справа на рис. 5, б показан профиль распределения модуля разницы высот цунами вдоль побережья Сахалина (на изобате 10 м), полученных по упомянутым моделям. Среднеквад-ратическая величина «рассогласования» составляет более 10 см для участка побережья от 51.5 до 54° с. ш. и около 5 см для более южного участка
Рис. 5. а — карта распределения максимальных высот волны цунами вблизи северо-восточного побережья о. Сахалин, рассчитанных по модели батиметрии GEBCO_2014. Справа (сплошная фиолетовая линия) — распределение максимальных высот цунами вдоль побережья на глубине 10 м; б — карта распределения модуля разницы максимальных высот волны цунами вблизи северо-восточного побережья о. Сахалин рассчитанных по модели батиметрии GEBCO_2014 и комбинированной батиметрии с учетом данных KCS. Справа (сплошная красная линия) — распределение модуля разницы максимальных высот волны цунами вдоль побережья на глубине 10 м. Прямоугольником отмечен полигон KCS, кружки — положения точек, для
которых рассчитывались модельные мареограммы
от 54 до 50° с. ш. Причем отдельные выбросы могут достигать значений 0.4 м. То есть непосредственно напротив полигона КСБ ошибка в задании глубин внутри прямоугольника КСБ приводят к расхождению оценок высоты цунами в 15-20%, а в отдельных точках до 40-50%.
Заключение
Прогресс в развитии численного моделирования процесса распространения цунами в глубоком океане и на шельфе в настоящее время достиг того уровня, когда отдельные научные вопросы становятся в большей степени проблемами технологии и инженерии. Современные возможности вычислительной техники (объемы памяти, быстродействие) позволяют использовать в моделях сеточные массивы батиметрии с разрешением, заведомо превосходящем доступные картографические данные. И именно поэтому в последние годы все большее внимание при изучении цунами стало уделяться качеству данных батиметрии. В настоящем исследовании мы ограничились оценкой влияния ошибок батиметрии на качество моделирования для относительно небольшого участка Охотского моря размером 334 х 130 км. Этот прямоугольник охватывает область северо-восточного шельфа и склона о. Сахалин, и соответственно наибольшему влиянию ошибок в задании глубин при расчете высот цунами подвержена часть восточного побережья от 52 до 54° с. ш. Выполненные численные эксперименты (предполагалось, что массив глубин КСБ близок к истинному) показали, что качество сеточного массива, заданного на основе данных 0БВС0_2014, может приводить к ошибкам в максимальной высоте цунами примерно 15-20% с отдельными выбросами до 50%.
Необходимо отметить, что полученные выводы носят локальный характер, связанный с качеством батиметрии непосредственно в области, прилегающей к северо-восточной части Сахалина. Очевидно, что неточности в задании рельефа дна на всей трассе распространения волн приведут к еще большим ошибкам.
Результаты этого исследования наиболее важны при проведении численных расчетов при цуна-мирайонировании участков побережья. Завышение или занижение значений максимальных высот цунами при оценке цунамиопасности могут привести к негативным последствиям: излишним расходам при строительстве прибрежных инфраструктурных сооружений в случае преувеличенных оценок или к трагическим людским потерям и экономическому ущербу в случае занижения уровня опасности. Разработка общедоступных батиметрических массивов данных высокого качества особенно важна для прибрежных районов Дальневосточного региона России, наиболее подверженных воздействию волн цунами: побережья Камчатки, Курильских островов, Сахалина, Приморья и Магаданской области.
Работа выполнена при финансовой поддержке грантов РНФ (проекты 14-50-00095 и 14-17-00219).
Список литературы
1. Rabinovich A.B., Eble M.C. // Pure and Appl. Geophys. 2015.
2. Gonzalez F., Geist E., Jaffe B., Kanoglu U. et al. // J. Geophys. Res. 2009. 114. P. C11023.
3. Ле Меоте Б. Введение в гидродинамику и теорию волн на воде. Л.: Гидрометеоиздат. 1976. (Le Me-haute B. An Introduction to Hydrodynamics and Water Waves. Springer Science & Business Media, 2013. Dec 18.)
4. Levin B.W., Nosov M.A. Physics of Tsunamis. Second Edition. Springer, 2016.
5. Shuto N. // Tsunami Hazard: A Practical Guide for Tsunami Hazard Reduction / Ed. by E. N. Bernard. Kluwer Academic Publisher, 1991. P. 171.
6. Imamura F. // Long-Wave Runup Models / Ed. by H. Yeh, P. Liu, and C. Synolakis. River Edge, NJ: World Scientific, 1995. P. 43.
7. Isozaki I., Unoki S. // Studies on Oceanography / Dedicated to Prof. Hidaka in Commemoration of His Sixtieth Birthday. 1964. 389.
8. Ueno T. // Oceanographical Magazine. 1965. 17. P. 87.
9. Aida I. // Bull. Earthq. Res. Inst. 1969. 47. P. 673.
10. Mogi А., Kawamura В., Iwabuchi Y. // J. Ged. Soc. Japan. 1964. 10. P. 180.
11. Hatori T. // Bull. Earthquake Res. Inst. Univ. Tokyo. 1965. 43, N 1. P. 129.
12. TakahasiR., Hatori T. // Bull. Earthq. Res. Inst. 1962. 40. P. 873 (in Japanese).
13. Hwang L.-S., Divoky D. // J. Geophys. Res. 1970. 75, N 33. P. 6802.
14. Hwang L.-S., Buttler H., Divoky D. Rat island tsunami model generation and open sea characteristics: Tetra Tech Inc. Report. Pasadena, USA, 1971.
15. Алексеев А.С., Гусяков В.К., Чубаров Л.Б., Шо-кин Ю.И. // Изучение цунами в открытом океане. М.: Наука, 1978. С. 5.
16. Гусяков В.К., Чубаров Л.Б. // Эволюция цунами от очага до выхода на берег. М.: Радио и связь, 1982. С. 16.
17. Chubarov L.B., Shokin Yu.I., Gusiakov V.K. // Computers and Fluids. 1984. 12, N 2. P. 123.
18. Liu P.L.-F, Synolakis C.E., Yeh H.H. // J. Fluid Mech. 1991. 229. P. 675.
19. Long-Wave Runup Models: Proc. of the Second International Workshop on Long-Wave Runup Models / Ed. by H. Yeh, P. Liu, C. Synolakis. Friday Harbor, USA, 12-17 September 1995. River Edge, NJ: World Scientific, 1997.
20. Dziewonski A., Chou Т., Woodhouse J. // J. Geophys. Res. 1981. 86. P. 2825.
21. Ekstrom G., Nettles M. // Phys. Earth Planet. Inter. 1997. 101. P. 219.
22. Aida I. // J. Phys. Earth. 1978. 26. P. 57.
23. Mansinha L., Smylie D. // Bull. Seism. Soc. Amer. 1971. 61. P. 1433.
24. Okada Y. // Bull. Seism. Soc. Amer. 1985. 75. P. 1135.
25. Chinnery M. // Bull. Seis. Soc. Am. 1951. 51. P. 355.
26. Maruyama T. // Bull. Earthq. Res. Inst. 1964. 42, N 2. P. 289.
27. Press F. // J. Geophys. Res. 1965. 70, N 10. P. 2395.
28. Savage J, Hastie L. //J. Geophys. Res. 1966. 11, N 20. P. 4897.
29. Sato R, Matsu'ura M. // J. Geophys. Earth. 1974. 22, N 2. P. 213.
30. Гусяков В.К. // Условно-корректные задачи математической физики в интерпретации геофизических наблюдений. Новосибирск: ВЦ СО РАН, 1978. C. 23.
31. Ji C, Wald D.J., Helmberger D.V. // Bull. Seism. Soc. Amer. 2002. 92, N 4. P. 1192.
32. Synolakis C.E., Bernard E.N. // Phil. Trans. Roy. Soc. A. 2006. 364. P. 2231.
33. Носов М.А. // Изв. РАН. Физика атмосферы и океана. 2014. 50, № 5. С. 540.
34. Bolshakova A.V., Nosov M.A. // Pure and Appl. Geophys. 2011. 168. P. 2023.
35. Nosov M.A., Bolshakova A.V., Kolesov S.V. // Pure and Appl. Geophys. 2014. 171. P. 3515.
36. Большакова А.В., Носов М.А., Колесов С.В. // Вестн. Моск. ун-та. Физ. Астрон. 2015. № 1. С. 61. (Bolshakova A.V., Nosov M.A., Kolesov S.V. // Moscow University Phys. Bull. 2015. 70, N 1. P. 62.)
37. Gonzalez F., Ernard E.N., Meinig C. et al. // Nat. Hazards. 2005. 35(1). P. 25.
38. Titov V.V., Gonzalez F.I., Bernard E.N. et al. // Nat. Hazards. 2005. 35. P. 41.
39. Wei Y., Bernard E.N., Tang L. et al. // Geophys. Res. Lett. 2008. 35, N 4.
40. Titov V.V. // The Sea. 15. Tsunamis / Ed. by A. Robinson, E. Bernard. Cambridge, USA: Harvard University Press, 2009. P. 371.
41. Data Announcement 88-MGG-02. Digital relief of the Surface of the Earth. Boulder, Colorado: NOAA, National Geophysical Data Center, 1988. http://www.ngdc. noaa.gov/mgg/global/etopo5.html
42. Jones M.T., Tabor A.R. Weatherall P. CD-ROM and Supporting Volume to the GEBCO Digital Atlas. Birkenhead, UK: British Oceanographic Data Center, 1994.
43. Smith W.H.F., Sandwell D. // Science. 1997. 277. P. 1956.
44. Centenary Edition of the GEBCO Digital Atlas. Published on CD-ROM on Behalf of the Intergovernmental Oceanographic Commission and the International Hy-drographic Organization as Part of the General Bathy-metric Chart of the Oceans. Liverpool, UK: British Oceanographic Data Centre, 2003. http://www.gebco. net/data_and_products/gridded_bathymetry_data/gebco_ one_minute_grid
45. British Oceanographic Data Center. The GEBCO_2014 Grid. 2016. https://www.bodc.ac.uk/data/documents/ nodb/301801/
46. Amante C., Eakins B.W. ETOPO1 1 Arc-min Global Relief Model: Procedures, Data sources and analysis.
NOAA Technical Memorandum NESDIS NGDC-24.
Boulder, Colorado: National Geophysical Data Center, Marine Geology and Geophysics Division, March 2009. http://www.ngdc.noaa.gov/mgg/global/relief/ETOPO1/ docs/ETOPO1.pdf
47. Jakobsson M., Mayer L.A., Coakley J.A. // Geophys. Res. Lett. 2012. http://www.ngdc.noaa.gov/mgg/bathy-metry/arctic/2012GL052219.pdf
48. NOAA National Centers for Environmental Information. U.S. Coastal Relief Model. http://www.ngdc.noaa.gov/ mgg/coastal/crm.html
49. Matveeva T., Soloviev V., Shoji H., Obzhirov A. Hydro-Carbon Hydrate Accumulations in the Okhotsk Sea (CHAOS Project Leg I and Leg II). Report of R/V Akademik M. A. Lavrentyev Cruise 31 and 32. St. Petersburg: VNIIOkeangeologia, 2005.
50. Shoji H., Jin Y.K., Obzhirov A., Baranov B. Operation Report of Sakhalin Slope Gas Hydrate Project 2011. R/V Akademik M. A. Lavrentyev Cruise 56. Kitami: New Energy Resources Research Center, Kitami Institute of Technology, 2012.
51. Cruise Report: Komex (Kurile Okhotsk Sea Marine Experiment) RV Akademik M. A. Lavrentyev Cruise 29, Leg 1 and Leg 2. Vladivostok — Pusan — Okhotsk Sea — Pusan — Okhotsk Sea — Pusan — Vladivostok. May 25 — August 05 2002 / Ed. By N. Biebow, R, Kulinich, B. Baranov. 2003. 190.
52. Шевченко Г.В., Файн А.В., Рабинович А.Б., Мансуров Р.Н. // Природные катастрофы и стихийные бедствия в Дальневосточном регионе. Владивосток, 1990. Т. 1. C. 253.
53. Гусяков В.К., Бейзель С.А., Чубаров Л.Б. // Вулканология и сейсмология. 2015. № 4. C. 59.
54. Лобковский Л.И., Куликов Е.А., Рабинович А.Б. и др. // Доклады РАН. 2008. 418, № 6. C. 829. (Lob-kovskii L.I., Kulikov E.A., Rabinovich A.B. et al. // Doklady Earth Sciences. 2008. 419, N 2. P. 320.)
55. Fine I.V., Kulikov E.A., Cherniawsky J.Y. // Pure and Appl. Geophys. 2013. 170, N 6-8. P. 1295.
56. Hayes G. Preliminary Finite Fault Results for the Nov 15, 2006 Mw 8.3 Kuril Islands Earthquake (Version 1). 2015. http://earthquake.usgs.gov/earthquakes/ eventpage/usp000exfn#finite-fault
57. Файн И.В., Куликов Е.А. // Вычислительные технологии. 2011. 16, № 2. C. 111.
58. Носов М.А., Колесов С.В. // Вестн. Моск. ун-та. Сер. 3. Физ. Астрон. 2009. № 2. C. 96. (Nosov M.A., Kolesov S.V. // Moscow University Phys. Bull. 2009. 64. P. 208.)
59. Nosov M.A., Kolesov S.V. // Pure and Appl. Geophys. 2011. 168, N 6-7. P. 1223.
Numerical tsunami modeling and the bottom relief
E.A. Kulikov 1,a, V. K. Gusiakov2,3, A. A. Ivanova1, B.V. Baranov1
1 Shirshov Institute of Oceanology, Russian Academy of Sciences. Moscow 117997, Russia.
2 Institute of Computational Technologies, Siberian Branch of Russian Academy of Sciences. Novosibirsk 630090, Russia.
3 Institute of Computational Mathematics and Mathematical Geophysics, Siberian Branch of Russian Academy of Sciences. Novosibirsk 630090, Russia.
E-mail: a [email protected].
The effect of the quality of bathymetric data on the accuracy of tsunami-wave field calculation is considered. A review of the history of the numerical tsunami modeling development is presented. Particular emphasis
is made on the World Ocean bottom models. It is shown that the modern digital bathymetry maps, for example, GEBCO, do not adequately simulate the sea bottom in numerical models of wave propagation, leading to considerable errors in estimating the maximum tsunami run-ups on the coast.
Keywords: tsunami, numerical model, bottom relief, bathymetric data. PACS: 91.30.Nw. Received 31 May 2016.
English version: Moscow University Physics Bulletin. 2016. 71, No. 6. Pp. 527-536.
Сведения об авторах
1. Куликов Евгений Аркадьевич — доктор физ.-мат. наук, зав. лабораторией; тел.: (499)-124-87-13, e-mail: [email protected].
2. Гусяков Вячеслав Константинович — доктор физ.-мат. наук, зав. лабораторией; тел.: (383)-330-70-70, e-mail: [email protected].
3. Иванова Анастасия Алексеевна — мл. науч. сотрудник; тел.: (499)-124-87-13, e-mail: [email protected].
4. Баранов Борис Викторович — канд. физ.-мат. наук, зав. лабораторией; тел.: (499)-124-79-42, e-mail: [email protected].
