CC BY
34
Известия ТИНРО
2005 Том 143
УСЛОВИЯ ОБИТАНИЯ ПРОМЫСЛОВЫХ ОБЪЕКТОВ
УДК 551.467(265.54)
А.Н.Четырбоцкий (ДВГИ ДВО РАН, г. Владивосток)
ПАРАМЕТРИЧЕСКАЯ ИДЕНТИФИКАЦИЯ МАТЕМАТИЧЕСКОЙ МОДЕЛИ ЭВОЛЮЦИИ ЛЕДЯНОГО ПОКРОВА ЯПОНСКОГО МОРЯ
Выполнен анализ данных параметров декадных состояний ледяного покрова Японского моря. Качественные результаты анализа используются для математической формализации модели термической динамики толщины покрова и модели термической динамики площади отдельной льдины. Далее формулируется модель эволюции распределения толщин льда. Выполнено аналитическое исследование модели. Разработана методика параметрической идентификации модели и оценка ее адекватности наблюдаемым распределениям. Предлагаемая модель и методика оценки параметров могут быть использованы для прогноза состояния ледяного покрова Японского моря.
Chetyrbotsky A.N. Parametrical identification of the mathematical model of ice cover evolution in the Japan Sea // Izv. TINRO. — 2005. — Vol. 143. — P. 149-178.
Distribution of 10-days parameters of ice cover in the Japan Sea for the period 1961-1989 is analyzed together with daily air temperature above the ice and wind rate in 1960-2001. The air temperature values over the ice of initial formation and initial thermal destruction are found to be equal. Phenomenology of the ice cover is formalized mathematically, and the models of thermal dynamics of ice thickness and thermal dynamics of a separate ice floe are created. The model of ice thickness is based on assumption on physical limit of marine environment and, hence, a limit of the ice thickness increasing. The model of an ice floe is based on conception of "competition" of ice floes for open water. Within the framework of the accepted positions and assumptions, both models have identical form of mathematical representation. Further, a kinetic model of ice cover evolution is formulated. Integration of the model equations over ice covered square results in the model of existential distribution of ice thickness. In the models, fast ice and drift ice are separated, but the drift ice in coastal areas can transform to the fast one during ice cover formation, and an opposite situation takes place in the time of the fast ice thawing and mechanical crushing. Some special cases of the models are analyzed. Technique of parametrical identification of the ice cover model is developed , and the model adequacy to observed distribution is estimated. Initial values of the parameters are set on the basis of selective distributions of drifting ice covering and fast ice thickness. The suggested model and technique of its parameters estimation can be used for forecast of ice cover in the Japan Sea.
Японское море является важным стратегическим природным объектом России. В его водах присутствуют богатейшие биоресурсы, а на шельфе выявлены перспективные для разработки месторождения нефти и газа. Кроме того, на акватории моря расположены интенсивные транспортные морские пути. В осенне-
весенний период года акватории определенной части его прибрежных районов покрыты толщей льда (Плотников, 2002). Данное обстоятельство существенным образом влияет на деятельность рыбопромысловых организаций и разработку месторождений на шельфе. Кроме того, весной, когда происходит взлом припая, его обломки вместе с рыбаками и техникой выносятся в открытое море. Поэтому для обеспечения безопасности жизнедеятельности населения актуальной является проблема научно обоснованного прогноза состояния ледяного покрова моря.
Построение указанного прогноза возможно только при наличии соответствующей математической модели эволюции ледяного покрова, где учитывается воздействие на эволюцию покрова реально измеряемых на практике динамических факторов (температуры атмосферы, скорости ветра и т.д.). В этом случае специфика условий наблюдений должна учитываться в уравнениях модели, что делает эти условия неотъемлемой частью объективного описания явления (При-гожин, Стенгерс, 2003). Как правило, параметрами натурных исследований морского ледяного покрова являются суточные измерения температуры воздуха и составляющих скоростей ветра, а также усредненные за декаду измерения геометрических характеристик покрова. Здесь используются статистическая выборка температур на 2-метровом горизонте слоя воздуха и статистическая выборка скоростей ветра на 10-метровом горизонте для каждого района акватории моря, где отмечается присутствие льда (материал за период 1960-2001 гг. предоставлен проф. В.В.Тунеголовцем). Статистическая выборка наблюдений параметров состояния ледяного покрова для этих районов за период 1961-1989 гг. предоставлена проф. В.В.Плотниковым.
Актуальность проведения статистических исследований обусловлена тем, что их результаты позволяют выявить функциональную структуру процесса эволюции покрова Японского моря, что является предварительным этапом построения математической модели эволюции покрова. Экспериментальный материал здесь также используется для параметрической идентификации моделей и оценки их адекватности наблюдаемым распределениям параметров состояния покрова.
Статистические распределения параметров состояний ледяного покрова
В практике исследования ледяного покрова проведение измерений обычно выполняется по следующей схеме. Акватория моря разбивается на ряд сравнительно однородных районов. Далее по картам ледовых условий в каждом районе определяются средняя за декаду сплоченность ледяного покрова 5 (отношение площади льда к площади района), его преобладающая толщина льда (возраст) к и преобладающий линейный размер льдин (форма) И. Их осредненные в рамках района значения приписываются его центру. Степень детализации сведений о ледовых условиях в данном случае напрямую зависит от размеров элементарных ячеек сетки, которыми покрывается акватория моря. В прибрежных зонах, наиболее интересных с точки зрения хозяйственного освоения, площади районов для снятия ледовых характеристик обычно уменьшаются. При этом оценивание численного значения 5 и И выполняется непосредственно с ледовых карт. Оценивание к выполняется на основании доступного для выполнения измерений материала (карт ледовых авиаразведок, спутниковых снимков и т.д.). Здесь и везде в дальнейшем используются обозначения и терминология из "Международной символики для морских ледовых карт и номенклатуры морских льдов" (1984).
Для выполнения статистического анализа используется выборка усредненных за декаду измерений к и И за 1961-1989 гг. Каждое наблюдение выборки определяется индексным набором (й, г, У), где й — номер декады (начинается с первой декады декабря); г — номер района (нумерация районов представлена на рис. 1, который заимствован из работы В.В.Плотникова
(2002)); У — год наблюдения. В выборке оценки к и Р представлены в баллах порядковых шкал их измерений. Согласно существующей методике (Международная символика ..., 1984), границы соответствующих диапазонов
определяются фиксированными наборами чисел {£;(т) :Б(т) = г/10}г=0Ч0, {й/т)};=0^6,
(Г) };=0^7 и = } = ^0(Г) = 0. Восстановление пропущенных значений выполнялось методами непараметрической регрессии (Четырбоцкий, Плотников, 2003).
130 131 132 133 134 133 133 137 138 139 140 141 142 143
Рис. 1. Порядок и нумерация районов акватории Японского моря Fig. 1. The order and numbering of areas of water area of the Japan Sea
Вычисление элементарных статистик для ¿-той декады г-того района следует выполнить на основании методики формирования исходных наблюдений: оцифрованное балльное значение каждого параметра в рамках его определенной градации (оно обозначается символом измеряемого параметра) есть равномерно распределенная случайная величина (с.в.). Предполагается статистическая состоятельность оценок и непрерывность оценок как функций случайных аргументов (Рао, 1968). Обоснованность принятия этих положений следует из анализа эмпирических распределений каждого из параметров для каждой конкретной декады отдельных районов. Так, декадные коэффициенты вариации рассматриваемых параметров состояния покрова для зал. Петра Великого (центр района — 131,75° в.д. 43,25° с.ш.) не превышают 0,4. Результаты анализа показывают допустимость линеаризация оценок в окрестности математического ожидания (м.о.) аргументов. В этом случае выборочное среднее и выборочная дисперсия некоторой функции X = (р(хг х2,..., хп) определяются соотношениями
п п
X - ф( XI, х2,..., хп) и а2 (х) - £ [(Эф / Эх,- )2 а2 ( х,) + Эф / Эх,- £ (Эф / Эх.) , х.)] ,(1)
где частные производные вычисляются в точке м.о.; соу(х4, х) — выборочная ковариация между хк и х. В случае, когда переменная хк является равномерно
распределенной на полуинтервале (хк-),х^)] с.в., ее м.о. и дисперсия определяются выражениями Е(хк) = (хк-) + х^у)/2 и а2(хк) = (х[г) -хк-))2 /12. Для оценки выборочной сплоченности Эаг в качестве аргументов выступают набор (£г};=1+10
и набор {р^},=Ь10 частот встречаемости площади льда с г-тым баллом сплоченности. Тогда
Е^) -£Е(Р^Еф),
1=1
10 10 а2^) -£[Е2(Sг.)а2(Р£) + Е2(Р^а2^)-Е)Е(Р^)£E(S,)Е(Р£)/(л,-1)],
¿=1 ¡Ф1
где Е (р^) — выборочная оценка средней частоты встречаемости, выборочная
дисперсия которой а2 (р^); П — общее число наблюдений для ¿-той декады (для наблюдений осеннего сезона это число равно 28, для наблюдений остальных сезонов — 29). Действительно, с.в. Бг(1 = 1 ■ 10) не зависит от частоты ее
встречаемости. Поэтому cov( р^, Si) = 0. Кроме того, согласно методике проведения наблюдений, 10-элементный набор появления льда различных сплоченностей в каждом наблюдении имеет только один ненулевой элемент. Тогда для г Ф j
cov(p(fJ,р^) = -Е(р^)Е(р./(пй -1). Оценки выборочных средних и выборочных дисперсий для толщин льда Нйг и преобладающих размеров льдин Гйг записываются аналогичным образом:
Е(кЛг) -£Е(р^)Е(к,),
¿=1
а2(кЛг) - £[Е2(к,)а2(р%) + Е2(Р£)а2(к,)-Е(к,)Е(Р£)£Е(к.)Е(р. /(пй -1)],
¿=1
6
е (^)-£ Е ()) Е (*;),
=1
а2(ЕЛг)-£[Е2&)а2(Р?) + Е2(Р?)а2(^)-Е(^)Е(РГ?)£Е(^.)Е(Р£)-1)] .
152
Здесь Е( — выборочная оценка средней частоты встречаемости отдельных полей льда с г-тым баллом толщины, а о2 (р— ее дисперсия; Е()) — выборочная оценка средней частоты встречаемости отдельных полей льда с г-
тым баллом преобладающего размера льдин, а о 2( — ее дисперсия. Вследствие существенного превышения масштабом Р масштабов 5 и к при анализе используется логарифмическая шкала ее измерений. Обоснованность подобных преобразований обусловлена тем обстоятельством, что для реальных ситуаций важен поиск зависимостей не только между исходными параметрами, но и между их некоторыми монотонными функциями.
Особенности изменения геометрических параметров ледяного покрова Японского моря для рассматриваемого периода наблюдений представлены на рис. 2. Для различия гистограмм их оси абсцисс отмечены маркерами соответствующих переменных. Статистические оценки отдельных параметров указывают на следующее. Гистограмма сплоченности на рис. 2 (А) соответствует и-образному распределению. Первая мода соответствует льдам низкой сплоченности (осеннее формирование и весеннее разрушение покрова), а вторая — покрову тех районов акватории, где льдом занято от 0,6 и более частей площади. Значение моды приходится на значение сплоченности, равное 0,8. Интересно отметить, что у В.О.Ивченко, Д.Е.Хейсина (1974) и С.Н.Овсиенко (1976) это значение фигурирует в качестве верхней оценки сплоченности для ненулевого значения аналога гидростатического давления. Наблюдения в указанной зоне соответствуют ледяным полям большой протяженности и зрелым льдам припая, что имеет очевидное объяснение: на акватории моря есть районы, где льды покрывают большую часть до момента их весеннего термического разрушения. Несмотря на очевидность этого факта, и-образное распределение сплоченности (рис. 2, А) предопределяет требования для математической модели: должна быть предусмотрена ситуация, когда при максимальной сплоченности дальнейшее понижение температуры не изменяет ее величину или ее полная производная в этом случае равна нулю.
Гистограмма толщин на рис. 2 (В) соответствует бимодальному распределению. Согласно ей, первая мода соответствует льдам осеннего и зимнего этапов, а вторая — зрелому состоянию припая в северных районах моря.
Гистограмма логарифма размеров отдельных льдин на рис. 2 (С) также соответствует бимодальному распределению, что достаточно очевидно: при изменении сплоченности соответствующим образом изменяются размеры отдельных льдин. Достаточно отметить, что выборочный коэффициент корреляции между Б^) и 1п Б(Р^) равен 0,658.
Распределение выборочных точек на рис. 2 Ш) интуитивно предсказуемо: чем ближе сплоченность к 1, тем ниже разброс ее значений. В качестве меры разброса выступает величина среднеквадратического отклонения (с.к.о.) о-пара-метра. В правой части конфигурации сосредоточены выборочные точки в основном зимнего периода эволюции. При этом общее изменение сплоченности за весь цикл эволюции в определенной степени удовлетворяет требованиям обратимого процесса, когда осеннее формирование покрова можно описать в терминах его весеннего разрушения и наоборот. Небольшие отличия между осенними и весенними этапами обусловлены следующим: весной, кроме таяния льда, имеет место его дробление, что приводит к более быстрому его разрушению.
Распределение выборочных точек на рис. 2 (Е) указывает на существенное различие осеннего и весеннего этапов эволюции покрова. При осеннем формировании покрова значение о(ка) следует почти линейному закону Б(ка ). В отличие от предыдущего случая, здесь отсутствует период падения о(к) для больших значений аргумента. Указанная ситуация обусловлена непрерывным характером
изменения толщины и неоднородной структурой покрова всей акватории моря. Отличительной особенностью весеннего разрушения является снижение Б(кл) при фактически неизменном с(ка ). Данное обстоятельство обусловлено влиянием теплопроводности льда и проникновением коротковолновой радиации в толщу льда (Аппель, Гудкович, 1992). Повышение температуры атмосферы сопровождается таянием ледяной толщи как изнутри, так и с обеих сторон ее поверхности. При этом рост свободной ото льда площади воды, отражательная способность которой существенно ниже поверхности льда, способствует интенсивному прогреву нижележащих слоев воды. Поэтому первыми выбывают из рассмотрения тонкие льды, а достаточно мощные ледяные массивы не претерпевают столь существенных изменений. Сообразно вышеизложенному сначала уменьшается толщина покрова, а затем снижается его неоднородность с(ка ). Различие осеннего и весеннего этапов эволюции покрова обусловлено также разницей физико-химических свойств подстилающей морской воды. Действительно, осенью к моменту появления льда в толще воды формируется верхний опресненный квазиоднородный слой, что обусловлено опусканием охлажденных и, следовательно, более тяжелых масс морской воды и поднятием более теплых легких глубинных масс. При этом охлажденные воды имеют повышенную соленость, а воздействие ветра существенно увеличивает интенсивность перемешивания. Поэтому начальное формирование покрова протекает в достаточно однородной среде, а весной его разрушение происходит уже в неоднородной среде.
W JV W
Рис. 2. Гистограммы S, h, F (A-C); совместные распределения выборочных средних и их среднеквадратических отклонений (D-F); парные распределения выборочных средних значений (G-I). Сезоны: 1 — осенний, 2 — зимний, 3 — весенний
Fig. 2. Histograms of S, h, F (A-C); joint distributions selective average and them mean-scuare deviations (D-F); pair distributions of selective average values (G-I). Seasons: 1 — autumn, 2 — winter, 3 — spring
Размещение выборочных точек на рис. 2 (F) показывает, что распределение ln E(Fdr) идентично распределению E(Sd). Отличия касаются наклонов огибающих и наличия участка, где значения с.к.о. принимают нулевые значения. Первая
часть утверждения обусловлена различием удельных темпов роста /уменьшения значений этих переменных. Другая часть наблюдаемого эффекта объясняется следующим. Как правило, крупные льдины (их размеры превышают 104 м) присутствуют в береговой зоне, представляя собой припай, который в осенне-зимний период устойчив к изменению своих размеров. Скачок о(1п Р) соответствует формированию, а также таянию и разрушению крупных льдин припая.
Парные распределения параметров представлены на рис. 2 ^—1). Конфигурации выборочных точек для парных распределений (рис. 2, G, I) указывают на отсутствие значимых корреляций между Б(Эаг) и Б(кл ), а также между 1п Б(Ра) и Б(ка ). При этом отмечается широкий разброс Б(Эаг) для тонкого льда (0 < к < 0,5), что может быть обусловлено наличием точки перегиба кривой толщины. Действительно, при понижении температуры на акватории сначала формируется тонкая пленка льда, и только затем следует заметный рост Б(ка ). Случай Н указывает на взаимно-однозначное соответствие Б(Эаг) и измеренного в логарифмическом масштабе значения Б(Ра ) (коэффициент корреляции равен 0,658). Поэтому можно допустить, что в рамках принятой пространственной дискретности на образование торосов при агрегации отдельных льдин расходуются их незначительные площади. Действительно, если бы наблюдалось обратное утверждение, то оно относилось бы и к преобладающему размеру льдин, а это не соответствует реальной ситуации.
Анализ распределений температуры показывает, что температура на стандартном при исследованиях покрова 2-метровом горизонте слоя воздуха в начальный период формирования покрова и в начальный период его таяния совпадает. Так, средняя температура атмосферы за предшествующую декаду первичного присутствия льда на акватории Японского моря Т0 = -(8,4 ± 4,2) оС, а в момент первичного разрушения сплоченности (переход через максимум ее значений) — Т51 = -(9,1 ± 4,7) оС. Если за момент начального таяния покрова принять декаду перехода через максимум толщины покрова, то Тк1 = -(7,8 ± 4,3) оС (для различия величин им приписан соответствующий подстрочный индекс). Согласно критерию Стьюдента, между этими величинами нет статистически значимых различий. Таким образом, согласно результатам анализа декадных распределений сплоченности и толщины покрова Японского моря справедливым является допущение о совпадении температуры атмосферы в моменты первичного формирования покрова и первичного разрушения.
Кинетическая модель эволюции ледяного покрова
Исследование ледяного покрова обычно включает изучение сплоченностей, толщины и массы льда. Если в качестве динамической переменной системы используется масса образований льда, то для расчета скорости дрейфа требуется определить их площадь (дрейф льда обусловлен касательными напряжениями на обеих сторонах образований). Здесь независимыми динамическими переменными системы являются толщина к и площадь отдельной льдины а. Их независимость проявляется из анализа гистограммы (рис. 2, G), где отсутствует статистически значимое согласованное изменение толщины и площади покрова. Поэтому в дальнейшем полагается, что в произвольные моменты времени совокупности льдин можно упорядочить согласно сформированной на основании исходной выборки двумерной таблице размеров площади отдельных льдин и их толщины. Тогда эволюции покрова соответствует определенный переход величин из одной ячейки этой таблицы в другую. На основании результатов обработки данных и простых качественных соображений эволюцию покрова можно представить сочетанием процессов: 1 — дрейф льда; 2 — непосредственное образование изолированных первичных льдин; 3 — последовательный переход льдин низких площадей в льдины более значительных размеров в результате их индивидуально-
го роста; 4 — формирование площадей k-того размера в результате агрегации площадей m-того и j-того размера ak = am + a , где m + j = k и в рамках принятых на практике пространственных масштабов потерями площадей льда на формирование торосов можно пренебречь; 5 — выбывание определенного числа льдин из k-того размера вследствие их агрегации с льдинами произвольной площади. Достоверность 4-го допущения обусловлена тем обстоятельством, что потери площадей отдельных льдин при их столкновениях незначительны (достоверность этого факта косвенно следует из рис. 2, H). Площадь льдины является аддитивным динамическим параметром эволюции покрова: при агрегации льдин формируется льдина суммарной их площади. Более того, непосредственные наблюдения указывают на преимущественно парный характер столкновений льдин, столкновения более высоких порядков крайне редки. Поэтому в дальнейшем рассматривается только случай бинарной агрегации льдин.
Временная дискретность используемых здесь распределений температуры воздуха и скоростей ветра составляет одни сутки, поэтому в дальнейшем рассматривается суточный временной шаг. Анализ соответствующих уравнений Эйлера для скорости дрейфа показывает, что для этого временного шага она имеет квазистационарный характер и определяется простыми соотношениями (Масловский, 1982; Тимохов, Хейсин, 1987). При этом предполагается массовый дрейф льда. Достоверность этого положения следует из анализа спутниковых снимков дрейфа льда: для массивов льда более 50 км отклонения в движениях отдельных льдин не превышают 5 % от осредненного непрерывного движения ледяного покрова (Легеньков, 1992).
При выполнении указанных допущений динамика плотности числа льдин nah = n(x, y, t, a, h) в открытых районах моря определяется результатом интегрирования уравнения Больцмана по импульсам отдельных льдин и положениями уравнения Смолуховского о характере парных столкновений:
dnah /dt + дщпак /+ diinah / dh + danah /да = fah + Dahd2nah / dh2 + Qah + Rah, (2)
где u = (u u2) — скорость дрейфа льда; x = (xv x2) — пространственные координаты; h = dh / dt — скорость термического роста толщины; а = da / dt — скорость термического роста площади изолированной льдины; fah = f(a, hv T) — скорость образования/выбывания из системы изолированных льдин начальной градации толщины h{ и площади a; Dah — коэффициент диффузии (для простоты изложения принимается const); Qah = Q(x, y, a, h, T), Rah = R(x, y, a, h, T) — члены, которые представляют собой статистическое описание динамики агрегации и дробления частиц-льдин (Лушников, Пискунов, 1976; Волощук, 1984). Дополнительно полагается, что в отдельном районе происходит выравнивание толщин и этот процесс количественно описывается диффузией.
Представляется целесообразным выполнить параметризацию h и a на основании положения о пространственной ограниченности, вмещающей покров морской среды (Четырбоцкий, 1999, 2001). Понятно, что h пропорциональна как текущему значению h, так и текущему значению доступного для нее ресурса hmax — h, где hmax — максимальная толщина льда за многолетний период наблюдений. Действительно, на этапе первичного формирования покрова, когда h мало, скорость роста является ее линейной функцией. По мере ее приближения к зрелому состоянию покрова h ~ hmax, когда весь ресурс hmax — h уже "израсходован", толщина стабилизируется, и поэтому на этапе зрелого состояния h равна нулю. На этапе же таяния и механического разрушения возникает обратная ситуация. При этом в момент окончательного разрушения льда эта функция обращается в ноль. Интенсивность процессов определяется температурой атмосферы, толщиной снежного покрова hs и рядом динамических факторов b = (bp b2,...bk) (альбедо по-
верхности снега и т.д.). Поскольку результаты статистического анализа распределений температуры воздуха показывают совпадение температур начального формирования и начального таяния покрова (в дальнейшем обозначается через Т*), то
к = /(Г,Г*,к,¿)(ктах -к)к.
Аппроксимацией / (Г, Г , ку, 6) может служить ее линейное представление
к
а(Г - Г *) + +£у, (6, - 6*), где 6* — значение факторов внешней среды, при
I=1
которых начинается первичное формирование покрова. Значения оценок соответствующих коэффициентов могут быть найдены на основании многолетних рядов наблюдений. Однако для покрова Японского моря подобные наблюдения имеют фрагментарный характер или отсутствуют. Поэтому далее рассматривается частный случай, когда воздействие внешней среды на толщину покрова характеризуется температурой воздуха и не учитывается толщина снежного покрова, а значения остальных внешних динамических переменных фиксированы. Таким образом, термическая эволюция толщины покрова записывается в виде
Н = (Т - Т * )[а Л©1 (Т) + а;© 2 (Т )](Атах - Н)Н,
©1(7) = ©(Т* - Т), (3)
© 2 (Т) = ©(Т - Т *),
где ©(2) — функция Хевисайда, равная 1 при г > 0, в ином случае — 0; а;, а'; — неотрицательные коэффициенты пропорциональности. Их размерность — (м ■ 0С ■ сут)-1, т.е. их численное значение определяет суточное изменение толщины покрова при изменении температуры атмосферы на 1 оС. Перемена
знака к происходит в тот момент, когда температура слоя воздуха превышает значение Т*. В настоящем случае Т характеризует определенное сочетание внешних условий, при которых начинается первичное льдообразование. Очевидно
также, что ее значения для припая ГВ и льда открытой части прибрежных районов моря Т* различаются. Действительно, воды припая являются относительно неглубокими и более пресными вследствие опреснения их речными водами и стоками промышленных предприятий, и поэтому льдообразование в открытых районах моря начинается при более низких температурах атмосферы, чем в прибрежных районах. Термическое разрушение льда в этих районах также начинается раньше разрушения припая, при более низких температурах, — в этот период активно проявляет себя сочетание воздействия солнечной радиации и теплых течений Японского моря (Якунин, 1986). Численные значения оценок аак, Г * и ТВ определяются при проведении серии вычислительных экспериментов на основании выборки распределений температуры и толщины покрова.
Если Т — периодическая функция, то результат ее интегрирования сдвинут во времени по отношению к периодичности Т. Поэтому наблюдается временное запаздывание между минимумом температуры атмосферы Т и максимумом толщины покрова. Указанная закономерность отмечается в многочисленных исследованиях (Тимохов, Хейсин, 1987), а для покрова Японского моря, в частности, результатами настоящего статистического анализа выборки. Независимо от темпа охлаждения температуры атмосферы Г кривая Н(0 в момент ^) (корень уравнения Т = -а;Ншах(Т* -Т)2) имеет точку перегиба
й(4А)) = {¿тах -Т(4А))а-1[Т* -Т^®)]"2}/2. В этом случае существует определен-
157
ный "инкубационный" период в течение которого при понижении температуры атмосферы собственно и происходит начальное формирование покрова. В реальных условиях наличие этого периода обусловлено следующим. В начальный этап эволюции образуются отдельные изолированные скопления зародышей льда. Затем площади зародышей растут, происходит их агрегация, а далее акватория района покрывается тонкой пленкой льда. Толщина этой пленки практически не изменяется в течение некоторого промежутка времени, после чего следует ускоренный рост толщины до ее предельного значения ктах. Исключение представляет ситуация, когда за похолоданием Т* - Т > 0 следует резкое потепление Т* - Т < 0.
При рассмотрении термической эволюции площади отдельной льдины
Ашах А*
ее текущий ресурс определяется площадью А* - | |п(а, h)adadh = А* - А,, где
0 0
А — площадь района акватории. Сообразно (3) тогда для а следует:
а = (Т - Т* )[аа©! (Т) + а'а02 (Т)](А* - А,)а, (4)
где а а, а'а — неотрицательные коэффициенты пропорциональности.
Представляется естественным определить / посредством быстро убывающей по а кривой. Например, гиперболой С(Т, Л)(а1 + а)-к, где первый сомножитель отличен от нуля только для к = Н1, а к > 2 — некоторое число. Функция С(Т,
hmax А
к) определяется из условия 1 = | | fahи представляет собой плотность
0 0
числа образованных при формировании покрова льдин. При таянии льда этот интеграл определяет плотность числа растаявших льдин. В первом случае разумно предполагать, что I пропорциональна доступному для льда ресурсу А* - А, , а
во втором случае — текущей площади тонкого льда А1. В обоих случаях интенсивность процессов определяется переохлаждением Т* - Т надледного слоя воздуха, поэтому I = (Т* -Т)[С(А* - А,)01 (Т) + С'А102(Т)], где С и С' — неотрицательные коэффициенты, характеризующие скорость образования/удаления из системы числа первичных льдин единицей, доступной для льда морской воды/ единицей площади тонкого льда. Поскольку а1 << А*, то
С, (Т, И) - (Т* - Т)[С(А* - А, )0 (Т) + С'А102 (Т)](к - 1)а,к-15М1,
где 5Й и — символ Кронекера. При записи (2) полагается, что в произвольный момент времени совокупность отдельных образований льда может быть однозначно упорядочена согласно сформированной на основании выборки двумерной
таблице размеров площадей и толщин 0 = {(а,к): 0 < а < А*,0 < к < ктх}. При условии, что площадь результата агрегации не превышает площадь района А*, результат представляет собой элемент таблицы. Поэтому здесь так называемый коагуляционный член (Волощук, 1984) принимает вид
а*=(а* - А, )
2 Яв(а* - а'*'' а'*')Па-а'м„-а'„')/(а-а')Па',*'^'- ПаН ДОРИ. а'*')Па',*'dа'dh' 2 О
где Оа* = {(а', *'): 0 < а' < а,0 < < *}; в(2, у) — ядро кинетического уравнения, которое представляет собой симметричную функцию. Введение дополнительного
сомножителя (А - А,) позволяет учесть то обстоятельство, что агрегация воз-
158
можна только при наличии свободной ото льда акватории района. При записи в( 2, у) учитывается тот факт, что обычно агрегирование формализуется в терминах объемов или масс, а не в терминах площадей отдельных льдин. Поэтому предлагаемый механизм агрегации состоит в следующем: площадь результата агрегации равна сумме площадей льдин, а толщина пересчитывается — сумма объемов исходных льдин делится на суммарную площадь. Присутствие коэффициента / обусловлено двойным пересчетом пар частиц-льдин.
Здесь при записи Яак предполагается, что результат дробления льда имеет толщину исходной льдины. Модификация соответствующего члена кинетического уравнения Мелзака (Лушников, Пискунов, 1976) принимает вид
А а
Rah = |у(а', а, й)па^а' - а~хпаЛ |у(а, а , й^а' ,
а 0
где у(а', а, к) — вероятность образования льдин площади а при дроблении льдин площади а' > а. Нормирование у(а',а,к) выполняется таким образом, что интег-
а
рал Рак = а-1 |у(а,а',к^а' равен вероятности распада льдины площади а в еди-0
ницу времени.
Поскольку рассматривается полный цикл эволюции покрова, то начальное распределение полагается равным нулю. Граничные условия естественным образом следуют из (2) и (3): отсутствие потоков числа льдин на крайних границах градаций размеров льдин и толщин. Таким образом, цикл эволюции покрова рассматривается в рамках закрытой для потоков вещества системы.
Распределение толщины ледяного покрова
Уравнение для распределения толщины льда Ак = А(х, у, I, к) получается умножением (2) на а и последующим интегрированием по этой переменной полученного результата. Для районов открытого моря получаем
дАк /3* + дМк /дх, +ЭИ /дк = (Т* -Г)/Ак + Якд2Л /дк2,
^ = (а а Ак +а аА)(А* - А, ^ (Т) + [< (А* - А,) Ак +а;Д А^ (Т), (5)
кш„
А,= | Ак^к .
0
В прибрежных районах различают льды открытой части Ак и льды припая Ак(В), обломки которых при весеннем термическом разрушении и механическом дроблении выносятся в открытую часть и переходят в категорию льдов открытого моря. Для учета этого обстоятельства в правую часть эволюционного уравнения из (5) следует добавить член, который отражает прирост площадей различных толщин льда открытого моря за счет выноса обломков припая. Соответствующим образом площади припая убывают на эту величину. Понятно, что этот член пропорционален текущей площади припая Ак(В), т.е. его можно записать в виде Ь(Т, к)Ак(В). Таким образом, в прибрежных районах динамика площадей льда может быть записана в виде
дАк / дг + дм;Ак / дх + дкАк / дк = (Т* - Т)fЛh + Бакд2Ак / дк2 + , к)АкВ),
дАкВ) / дг + дкАкВ) / дк = (Т* - Т)fAh + ^кд2АкВ) / дк2 - Ь(Т, к)Акв), (6)
ктах
А, = | (Ак + Акв Ук.
о
Начальные и граничные условия для (5)-(6) принимают вид
Л (х, у,0) = 0 и НЛ,1 = 0, (7)
лк*чх, у,о) = 0 и н4В)
= 0.
к=0,кт
тах
Адекватность моделей (5)-(7) реальным наблюдениям и оценка их параметров выполняется на основании серии вычислительных экспериментов.
Понятно, что этот коэффициент пропорциональности зависит от температуры атмосферы и толщины льда: интенсивность дробления тонкого льда выше, чем толстого. Представляется, что для описания отмеченных особенностей дробления достаточно ограничиться линейным представлением
Ь(Т, к) = [Ьо + Ьт (Т - Т*) - Ькк]&2 (Т), (8)
где Ь0, Ьт, Ьк — неотрицательные параметры дробления припая. Их предметная интерпретация может заключаться в следующем. Параметр Ь0 отражает интенсивность усталостного разрушения припая как твердого тела, Ьт — интенсивность разрушения припая вследствие тепловых процессов и Ьк — интенсивность инерционных процессов разрушения.
В частном случае, когда припай на этапе формирования и зрелого состояния составляет основную часть ледяного покрова района (относительно закрытые прибрежные районы и небольшие озера), интегрирование для этого этапа второго уравнения в (6) по к при условии (7) приводит к уравнению для суммарной площади припая:
4В) (*) = А* 11--*-а а А* + а акк-[, (9)
[ ааА +аа/Д ехр[(ааА +аакк!)(ТВ -Т)*]\
— 1'
где Т = I Tdt — текущая средняя температура воздуха в прибрежном районе
t ^
1 о
акватории. В отличие от к($) кривая А^^) при низком темпе понижения температуры -Т(^А)) < (Р^-ааА*)[Т* -Т(^А))]2 может и не иметь точки перегиба
(^А) корень уравнения - Т(^Л)) = (Рк1 - ааЛ*)[Т* - Т(^Л))]2). В этом случае акватория равномерно покрывается льдом. При нарушении неравенства значение
А= {А* -а-1р^ -ТХ^А))а-1[Т* -Т(^А))]"2}/2. Данная ситуация соответствует постепенному росту площадей начальных льдин, а затем следует их скачкообразный рост.
Стратегия оценивания параметров
Центральным объектом теории оценивания параметров для математического представления рассматриваемого явления или процесса на основании статистической выборки является матрица моментов остатков (Кендалл, Стьюарт, 1973):
М(0) = ]Г ^(0)<(Э),
Ц=1
где п — число проведенных экспериментов; е^ = у^ - fg (хц, 0) — остатки или
разность между наблюдаемым в ц-том эксперименте значением £-того элемента вектора зависимых переменных (^ = 1 ■ т) и его соответствующим модельным
значением; хц = {х }е=1+/ — независимые переменные модели; 0 = {0а }а=м —
вектор оцениваемых параметров. В этих обозначениях общая запись целевой функции принимает вид
Ф(0) = ¥(M (0)),
где ¥ — заданная функция. В большинстве практических случаев при ее построении полагают, что остатки имеют нормальное распределение с одной и той же ковариационной матрицей. Принятие допущения о нормальности в реальных ситуациях является типичным: если при большом объеме выборки наблюдаемое значение случайной переменной формируется в результате воздействий многих аддитивных и независимых факторов, то выполняются условия центральной предельной теоремы и распределение случайной величины подчиняется нормальному закону (Айвазян и др., 1985).
При использовании для нахождения экстремума функционала ф (0) метода Гаусса элементы вектора градиента q = {qa } ^ и приближения матрицы Гессе N = {Naß }aß=1ik определяются матричными соотношениями (Бард, 1979)
q = -2¿ < и N = 2¿ BT Г£ц,
где Bw>a = /Э0а = Э/№ /Э0а и rgb = /ЭМгЬ.
Если остатки удовлетворяют нормальному распределению, то логарифм функции правдоподобия принимает вид
ln L(0, V) = -nk ln 2n- (n / 2) ln det V - Tr[ V -1M (0)] / 2,
где V — ковариационная матрица остатков, одинаковая для всех экспериментов. При записи ln L(0, V) принято допущение о том, что в каждом эксперименте остатки распределены с одной и той же ковариационной матрицей. При неизвестной ковариационной матрице V справедливы аппроксимации
V = M (0) / n , Г = (n /2)M-1 (0) и ф (0) = (n / 2) lndet M (0).
Следовательно, максимизация логарифма функции правдоподобия ln L(0, V) эквивалентна задаче минимизации функционала Ф(0). Для так называемого од-нооткликового случая нелинейного метода наименьших квадратов m = 1
B,a = / / Э0а, Г = 1, Naß = 2^ (Э/ц / Э0а)(Э/Ц / Э0р), V0 = 2С2N-1,
ц=1
где V0 — оценка ковариационной матрицы параметров. При неизвестной дисперсии a2 она заменяется ее оценкой e^ (0)бц (0) /(n - k).
Предметная интерпретация вышеприведенных положений для параметрической идентификации модели эволюции ледяного покрова состоит в следующем. Эволюция ледяного покрова характеризуется выборкой декадных наблюдений за 1961-1989 гг. В настоящем случае экспериментами являются выполненные за отмеченный период в каждом из 114 районов расчетной карты акватории 22 съемки декадных значений сплоченности, толщины и преобладающих размеров льдин. При этом съемки с 31-й по 36-ю декаду выполнялись за период 19611988 гг., а с 1-й по 16-ю декаду — за период 1961-1989 гг. В прибрежных районах (их общее число на рис. 1 равно 51) съемки выполнялись отдельно как для припая, так и для покрова в открытой части акватории этого района. Фактически выборка характеризует покров 165 районов. Для статистического анализа использовалась выборка, которая насчитывает 11130 наблюдений. На основании этой выборки была сформирована выборка усредненных за период 1961-1989 гг. наблюдений, каждое из которых характеризует площадь льда конкретной g-той
градации толщины (число исходных градаций толщины равно 6; Международная символика ..., 1984) А^ в отдельной декаде отдельного района акватории моря. Полученная выборка насчитывает 8022 наблюдений. Так как оценка параметров модели выполняется для всей совокупности районов, то номер наблюдения / оп-
г-1
ределяется выражением [ё1(г') - ё0 (г')] + ё - ё0 (г) + г -1} + g. Здесь г —
г'=1
номер района, г = 1 ■ 165; й — номер декады; й0(г) и й1(г) — начальная и конечная декады эволюции в г-том районе. Для каждого / набор независимых переменных составляет: текущие сутки года, среднесуточная температура на 2-метровом горизонте слоя воздуха, среднесуточные составляющие и скорости ветра на 10-метровом горизонте.
В качестве зависимых переменных выступают модельные площади льда отдельных толщин А, , . для текущего /-того дня года из равномерной системы разбиения диапазона (0, Нтах) на полуинтервалы. Введение этого разбиения обусловлено тем, что исходное разбиение является неравномерным, а для численного моделирования требуется как раз равномерное разбиение с некоторым шагом АН. Вхождение полуинтервалов модифицированного разбиения в полуинтервалы исходного разбиения можно указать посредством совокупности {Jg }^1+6. Каждый элемент ]представляет собой перечисление номеров полуинтервалов модифицированного разбиения, которые содержатся в исходном разбиении. В настоящем рассмотрении АН = 0,15 м, а первичная градация толщины Н1 = 0,075 м.
Набор искомых параметров 0 = (аh, а'А, аа, а'а, ааках, а ак, Т*в, Т*, Dah, Ь0, Ьт, Ьк) насчитывает 12 элементов. Тогда удобные для проведения вычислений выражения принимают вид
г-1
Л = 6{ХМ(г') - й>(г')] + й - ¿о(г) + г -1} + g,
г' =1
1о
г ,10(4-1)+й', у (0)
й'=1
е»(0) = -/0
ф(0) = ет (0)е(0), (10)
Мар= (Э/ /Э0„)(Э/ /Э0р )
л=1
о\0) = Ф(0) /(п -12), V(0) = 2а2 (0)N-.
При записи выражения для остатков учитывается еще тот факт, что А^ и X ,у (0) имеют различные временные масштабы измерений: первая характе-
УЕ Зг
ризует декадное значение площади £-той градации толщины, а вторая — ее суточное значение.
При оценивании параметров могут понадобиться конкретные практические рекомендации относительно специфики выбора определенной вычислительной процедуры минимизации и аппроксимации частных производных конечными разностями. Ниже приведена методология решения подобных вопросов, которые используются в настоящем исследовании. Отдельные методики их реализации рассматривают-
ся в многочисленных публикациях по различным тематикам. Представляется, однако, что заслуживает особого внимания их подробное изложение в некоторых работах (Стронгин, 1978; Бард, 1979; Айвазян и др., 1985; Ватутин и др., 1985), где на примерах из различных областей приводится методология решения подобных задач.
Значения искомых параметров имеют разные диапазоны изменения. Так, порядок чисел Тв и T * находится в первом десятке градусов Цельсия, а других параметров — в сотых долях соответствующих им масштабов. Для нахождения экстремума функционала Ф(0) целесообразно использование такого метода, где должным образом учитывается эта особенность. Подобными характеристиками обладает простой в численной реализации метод Марквардта. В численной реализации метода на его ¿-той итерации выбор вектора направления поиска V. определяется соотношением Vг- =-[N(0,) + ХгСг2]-1 qi, где X1 — параметр метода (он определяется в процессе линейного поиска экстремума вдоль выбранного направления); С — вспомогательная диагональная матрица, которая вводится для учета различия масштабов искомых параметров. Значения элементов этой матрицы вычисляются на каждом шаге процедуры минимизации. Тогда 0 = 0. + pvi, где р — размер шага. Останов вычислительной процедуры происходит как в случае превышения заранее заданного количества итераций, так и при низком значении размера шага
р< 10-2 шт [(0, а+10-3)/V,„]. (11)
а
Для аппроксимации градиента с помощью конечных разностей обычно используются односторонние и центральные разности, где приращение аргумента 30а целевой функции удовлетворяет условию
10-5 | 0а |<| 50„ |< 10-2 | 0а |.
Вопрос о применении односторонних или центральных разностей решается в каждом случае отдельно. Обычно первый тип используется на нескольких итерациях до момента, пока не нарушится сходимость итерационного процесса. Если есть предпосылки ее возникновения, то следует перейти ко второму типу. В соответствующей вычислительной процедуре можно автоматически предусмотреть подобную ситуацию и выработать способы ее разрешения.
В рассматриваемом классе многоэкстремальных задач оценивания требуется задание начального приближения искомых параметров. Его успешный выбор способствует эффективному поиску нахождения экстремума и более быстрому схождению к решению задачи поиска оптимума. При его поиске в значительной степени полагаются на априорные сведения о предмете изучения, выполнение определенных свойств решений математической модели и интуицию исследователя. Для нахождения начального приближения целесообразно использовать все доступные методы, и только в случае окончательной неудачи этой стратегии прибегать к более изощренным процедурам. Таким довольно часто применяемым на практике приемом решения многоэкстремальных задач является многоступенчатое оценивание, или редукция первоначальной задачи в связанную цепочку менее сложных. Данный прием целесообразен еще и тем, что простейшие частные случаи способствуют детальному пониманию существа рассматриваемых проблем. Действенным способом также оказывается упрощенный метод оценивания, когда целевая функция заменяется относительно простым функционалом.
Оценка начальных приближений параметров
Надлежащие рекомендации в полной мере применимы для построения начального приближения для параметров ак ,аь,Т*, аакН1,аа ,аа ,аакН1,Ь0,Ьт и на
основании временного изменения средних значений толщины и полных площадей припая. В этом случае объектом изучения является модель термической эволюции толщины (2) и модель эволюции общей площади припая. При этом на этапе формирование — зрелое состояние используется модель (8), а для этапа весеннего таяния и дробления — ее модификация. Целесообразность их применения обусловлена тем, что для этих моделей отсутствует необходимость вычисления частных производных по пространственным и фазовым переменным исследуемой системы. Подобное обстоятельство существенным образом упрощает вычисление оценок для этих параметров эволюции ледяного покрова.
Модель термической эволюции толщины покрова
Для оценивания начальных значений аъ ,а'ь и Т* модели (2) следует предварительно исходные декадные гистограммы средних показателей толщины припая трансформировать соответствующим образом в распределения толщины для каждого дня эволюции. Необходимо отметить, что в данной ситуации форма кривой конечного распределения интуитивно предсказуема: для большинства прибрежных районов акватории она имеет вид кривой на рис. 2 (В). Простые функциональные формы подобного изменения могут быть заданы, например, выражениями типа разложения по нормальным ядрам (Фукунага, 1979) или простой квадратичной зависимости от текущего времени эволюции:
К &*, *) = Кд ехр[-Кг>2 (* - 1Г* )с ], (12а)
К (*) = ~г, 1 + ~г, 2* + ~г,3*2, (12б)
где г — номер прибрежного района; /гВ — среднее первых дней начальной и конечной декад эволюции припая в г-том районе; / — текущее время;
Кг 15 Кг 2, Кг 15 Кг 2 и К 3 — искомые коэффициенты; с определяется методом перебора в процессе вычислений. Начальной точкой временного отсчета принят 310-й день текущего года, что является следствием методики проведения непосредственных наблюдений.
При заданных значениях /гВ и с логарифм от левой части (12а) линеен относительно 1п Н1 и Нг Поэтому для вычисления их значений можно использовать стандартный метод наименьших квадратов (МНК). В этом случае
Ка = X (Уг(г)-Уг)(Х(г)-хг)/ Х(Х(г)-хг)2,
г=г„ (Г) г=г„ (г)
К,1 = еХР(7г - К 2Хг ), где t0 (г) = 10 (г) +1 — первый день декады, когда на акватории г-того района отмечается присутствие льда; г1(г) = 10ё1(г) — последний день декады присутствия льда на акватории; Уг (/) = 1п Нг /10]), Ьг (ё) — выборочные значения средней толщины припая в й-той декаде; Хг (г) = (г - гг в)с. В этих соотношениях [//10] — целая часть аргумента, а знак черты — среднее значение величины за период присутствия льда.
Для нахождения коэффициентов представления (12б) требуется сформировать вектор зависимых переменных и соответствующую матрицу плана МНК. Элементы этого вектора определяются значениями средней толщины припая, а /тая строка — (1, /, /2). На основании этих массивов вычисляются Кгд,Кг 2 и Кг3. Процедура оценки выполняется только для тех прибрежных районов, выборки которых содержат необходимое в плане достоверности результатов количество
164
наблюдений. Таких районов всего 4, а остальные не обеспечены необходимым для вычислений объемом наблюдений. Оказывается, что наибольшее соответствие между толщиной h(tв, ^ из представления (12а) и ее декадным прототипом достигается при с = 4. В этом случае коэффициент корреляции между ними для каждого из прибрежных районов изменяется в диапазоне 0,8160,923, что показывает приемлемость этого представления в качестве аппроксимации декадных распределений толщины припая. Аналогичная ситуация отмечается и для (12б). Корреляция в этом случае несколько выше — от 0,869 до 0,964.
Согласно аппроксимации (12а), для этапа формирование — зрелое состояние максимальная скорость прироста толщины припая для Японского моря составляет 2,265 ' 10-2 м/сут при среднем значении прироста к = 4,651 ' 10-3 м/сут и среднеквадратичном ее отклонении о(к) = 6,876 ' 10-3 м/сут, а для (11б) —
2,113 ' 10-2, к = 7,818 ' 10-3 и с(к ) = 6,253 ' 10-3. Для припая сопоставляемых районов эти величины принимают значения: для зал. Петра Великого (центр
района — 131,75° в.д. 43,25° с.ш.) = 0,788 м, = 2,326 ' 10-7 м, тах 4 =
= 2,645 ' 10-2 м/сут при t = 58 сут, к8 = 7,909 ' 10-3 м/сут, о(к8) = 9,405 ' 10-3 м/сут; для района Татарского пролива (центр района — 140,75о в.д. 52,25о с.ш.) ^141 =
= 1,198, hlU2 = 2,998 ' 10-8, тах И114 = 2,122 ' 10-2 при t = 44, к,14 = 9,224 ' 10-3,
о(к1 14) = 7,539 ' 10 3. Выбор именно этих районов обусловлен их крайним положением на акватории моря, что способствует исследованию эволюции покрова в крайних ее положениях. При этом припай заполняет их основную часть, что позволяет рассматривать эти районы как замкнутые акватории, в которых влияние льдов открытого моря не столь заметно (по крайней мере, на этапе оценки начальных приближений параметров). Поэтому покров в этих районах является вычислительным полигоном для тестирования рассматриваемых здесь моделей.
Необходимо отметить определенное совпадение оценок и большую, относительно средних значений, величину среднеквадратичных отклонений (в первом случае коэффициент вариации превышает 1, а во втором — близок к ней). Последнее обстоятельство явно указывает на непостоянство скорости прироста толщины и ошибочность использования ее среднего значения в качестве показателя динамики толщины, что свидетельствует о приемлемости подобных (5)-(7) моделей эволюции морского ледяного покрова.
Приведенные выше численные оценки Т и Т показывают, что диапазон их изменения определяется интервалом минус 10 оС ... минус 5 оС средних за декаду температур 2-метрового горизонта надледного слоя воздуха. Поэтому при оценке
ак и а'к температура Тв считается фиксированной. Тогда для их вычисления можно организовать следующую вычислительную процедуру. При фиксированной температуре 7в оценки искомых величин, согласно схеме МНК, определяются значениями выражений
аА = ХХ4(0Х(0/XX*2(0 и а = ХХ4(О*Ц)/XX*2(о,
г ШгВ г ШгВ г ЫгВ г ЫгВ
где Хг (t) = [Тв* - Тг ^)][ктах - кг (t)]кг ^). В первом случае суммирование выполняется по районам и только для периода t < trB формирования и зрелого состояния покрова, а во втором — по районам и для периода его весеннего таяния
t > trB. Далее с определенным шагом ДТВ* изменяется Т*, а затем вычисления повторяются. Критерий останова вычислительной процедуры — стабилизация значений оценок. При выполнении процедуры допускается комбинирование использования аппроксимаций: (12а) для оценивания ак, а (12б) — а'к. Результаты использования этой процедуры приводят к следующим оценкам средних значений этих параметров и границ их доверительных интервалов:
ак = (2,597 ± 0,768) • 10-3,
ак = (8,243 ± 2,151) • 10-3,
Т* = -(6,9 ± 1,5)оС.
Границы доверительных интервалов определяются стандартным образом на основании среднеквадратичного отклонения и ^статистики распределения Стьюдента (Фукунага, 1979). Чтобы не загромождать рассуждения излишними подробностями, в дальнейшем единицы измерений оценок параметров опускаются. Поэтому далее приводятся именно численные значения оценок параметров.
Приемлемость указанных оценок в качестве достаточно близких начальных приближений истинных значений а К и а'К обусловлена значимой корреляцией между наблюдаемым распределением толщин и его модельным представлением (2). Соответствующие коэффициенты корреляции для различных районов не опускаются ниже 0,682. Кроме того, вычисленные начальные приближения имеют узкие рамки доверительных интервалов. В настоящей статье оценки вычисляются на основании объединенной выборки наблюдений для отдельных районов, что позволяет увеличить объем выборки. Подобные действия приводят к увеличению статистической состоятельности и статистической эффективности искомых оценок. В дальнейшем эти значения уточняются при оценке всей совокупности параметров системы на основании выборки декадных наблюдений площадей льда различных толщин.
Модель эволюции суммарной площади припая
Для оценивания параметров этой модели следует выполнить аналогичную предыдущему рассмотрению процедуру аппроксимации декадных гистограмм площадей припая и льда открытого моря для каждого дня эволюции покрова. Поскольку вычисления показывают практически полное совпадение обоих представлений для исходных декадных распределений суммарной площади припая (коэффициенты корреляции между различными формами аппроксимаций превышают 0,984), то используется только подобная (12а) форма:
А( в)(^, t) = «<?ехр[-в< > - Кв )С ]
А (^в ,t) = аг,1 еХР[-аг,2 (t - Кв )С ] ,
где г — номер прибрежного района; trB — среднее первых дней начальной и конечной декад эволюции припая в г-том районе; t — текущие сутки года;
(В) (В)
аг1 ,аг 2, аг 1 и аг 2 — искомые коэффициенты; с — определяется методом перебора в процессе вычислений. Понятно, что значение trB совпадает с соответствующей величиной из (12а). Вычисления показывают аналогичное совпадение и для коэффициента с. По сравнению с (12б) представление декадных распределений общих площадей припая в виде нормальных ядер предпочтительней, поскольку оно отражает существо процесса: в начальный и конечный моменты эволюции припая кривая площадей касается оси времени. Методика их оценки такая же, как и при рассмотрении динамики средних значений толщины припая.
Согласно результатам аппроксимации, максимальная относительная скорость прироста площади припая max Ar( B) / A* для Японского моря составляет 2,881 %
r
в сутки (A* — площадь акватории района). Среди районов, выборки которых содержат достаточный для вычислений объем, эта скорость отмечается в зал. Петра Великого. Здесь она наблюдается на 37-е сут от начала отсчета времени, а коэффициенты аппроксимации и элементарные статистики для выборки принимают значения: a^ = 0,533, а^ = 1,106 • 10-7, A8(B) = 1,168 %, = 1,033 % и коэффициент вариации 0,855. В районе Татарского пролива эти величины
принимают следующие значения: max Ajj^t) / Aj*14 = 1,875 % достигается на
t
29-е сут; a<?i = 1,015, а^ = 9,831 • 10-9, A/^ = 0,945 %, o(A(f4)) = 0,668 % и
коэффициент вариации 0,707. Сопоставление коэффициентов вариаций указывает на более плавный характер эволюции площади припая в Татарском проливе, чем в зал. Петра Великого. Наличие временного запаздывания между изменением толщины и площади следует из сравнения времени перехода соответствующих кривых через точку их перегиба: для зал. Петра Великого толщина имеет точку перегиба на 58-е сут ее эволюции, а площадь — на 37-е сут; для района Татарского пролива это происходит на 44 и 29-е сут. В первом случае величина запаздывания составляет 21 сут, а во втором — 15 сут.
Модель эволюции площадей припая следует из второго уравнения (6), где потоки и условия на крайних границах диапазона толщины равны нулю. Если выполнить интегрирование (7) по всему интервалу толщины, то для этапа формирование — зрелое состояние припая получается уравнение
) = (Г* - Г)(A* - Д. )(аж + а„A^
где — общая площадь припая, Ar = Ar( + Ar(B) — общая площадь покрова в
прибрежном районе, A( S) — площадь льда в открытой части r-того прибрежного района. Это уравнение насчитывает два линейных параметра, aah и aa (параметр полагается фиксированным). Для их оценки можно воспользоваться его эквивалентной формой записи
Y (t„) = a ah А +а a X r (tj , где i — номер наблюдения для периода формирования и зрелого состояния покрова; Yr (tri) = Ar(B) (tri) /[(T* - Г(tri))(A* - A(tri)]; tri — текущий день эволюции покрова в r-том районе (изменяется от начального дня первой декады присутствия льда на акватории до trB); Xr (tr;) = Ar(B)(tr;). Тогда значения искомых оценок определяются соотношениями
a a = £ [( I Yr (t) - Y )( £ Xr (t) - X)] / £ £ (Xr (t) - X) 2 и a ah = (Y - a aX) / V
r t <t,B t <t,B r t <t,B
где Y = HY (t) / K; X = I£Xr (t) / K; K = £ trB. Непосредственные вычис-
r t<trB r tStrB r
ления приводят к оценкам численных значений параметров afl = (2,693 ± 0,659) • 10-2, aahkl = (3,257 ± 0,581) • 10-3.
167
Эволюция общей площади припая на этапе его весеннего таяния и распада определяется уравнением
A(B ) = a T - T)(A? - A(S) - A(B )) )+а;л(г; - T) A^ - b(T , h) Ar(B),
где Ar(f) — площадь тонкого льда припая в r-том прибрежном районе; b(T, h) — определяется из (8) и характеризует интенсивность распада единицы площади припая толщины h; h(B} — средняя толщина припая. Поскольку в исходной выборке отсутствуют наблюдения временного распределения Ar(*), то параметр a'ah можно оценить только косвенным образом на стадии поиска экстремума функционала Ф(0) из (10). Поэтому для оценки начального приближения а'а, ba и bh параметр a'ah полагается равным нулю. Тогда для нахождения искомых значений параметров используется уравнение
AB) / AB) = a'a T - T)(AP - A(S) - A(B)) - ba (TB - T) + bhh(B>,
h
max
где h(B) = J hA(h B) dh / Ar(B) — текущая средняя толщина припая. Чтобы воспользо-
0
ваться стандартной вычислительной схемой МНК, следует привести это уравнение к матричному виду
Y=a'aX<? -baX+ bhXS>,
где ты = £B\tr;)/A*; X((J = [T* -Tr(tr4)][A; -)- A(BB(tri)];
X® = T** -Tr(tri); Xr(3) = h (B)(tri). При этом учитываются только те наблюдения, для которых tr • > tr,B. Так же как и в предыдущих случаях, увеличение количества наблюдений способствует повышению статистической состоятельности и статистической эффективности первоначальных оценок. Результаты применения вычислительной процедуры приводят к значениям
a'a = (5,138 ± 1,164) • 10-3, b0 = (6,310 ± 8,471) • 10-7, ba = (9,694 ± 3,543) • 10-3, bh = (2,056 ± 0,982) • 10-3.
Приемлемость этих оценок в качестве начального приближения можно оценить на основании определенной меры близости между результатами моделирования и наблюдаемым распределением общих площадей припая. Например, коэффициента корреляции между ними. В данном случае он равен 0,808, что свидетельствует о достаточно приемлемой степени адекватности модельного представления.
Результаты вычислительной процедуры оценки параметров
На основании результатов предыдущих рассуждений получены начальные оценки 9 из 12 искомых начальных значений параметров. Начальное приближение для Т* = -7,5 °С известно из литературных источников (Болч, Хуань, 1979). Начальное значение для коэффициента диффузии Dah выбираем из соображений его размерности. По своему физическому смыслу она (размерность) имеет размерность квадрата приращения толщины, деленного на размерность времени. В рассматриваемой ситуации Ah = 0,15 м, поэтому начальное значение коэффициента диффузии полагается равным 2,25 • 10-2 м2/сут.
Начальное значение а'а1г задавалось равным ааИа1. Таким образом, определены начальные значения для всех элементов искомого вектора параметров
0 = (а h , a'h , а a , a'a , а ahK a'ahA. TB , T ' , D ah , ЬТ , bA ', ^^ Далее выполняется про-
цедура поиска минимума функционала Ф(0) из (10).
Критерием останова вычислительной процедуры является динамическое условие минимальности шага (9). На первых пробных ее шагах целесообразно использование ограничений на число итераций вычисления функционала. Допускается также использование и других условий останова. Результаты приводят к следующим численным оценкам параметров
Dah = (1,310 ± 0,021) • 10-3; b0 = (6,324 ± 6,431) • 10-7, bT = (8,740 ± 1,096) • 10-3 и bh = (2,185 ± 0,046) • 10-3.
Анализ выборочных коэффициентов корреляции показывает высокую значимую корреляцию между оценками b0 и br (коэффициент равен 0,913). Данный факт означает, что при изменении bT происходит соответствующее изменение b0, а это свидетельствует о параметрической вырожденности модели. Поэтому в дальнейшем свободный член для функции b(t, h) опускается из рассмотрения. Значения остальных коэффициентов не превышают 0,344, поэтому можно сделать вывод о независимости параметров. После исключения b0 параметры модели были пересчитаны заново. Оказывается, что значение Ф(0) равно 19,724 при его ненулевом значении и 19,722 в противном случае. Таким образом, значение функционала практически не изменяется.
Анализ значений оценок параметров как показателей значимостей отдельных составляющих процесса эволюции показывает справедливость принятых при построении моделей допущений. Согласно результатам вычислений, мера интенсивности процессов весеннего таяния и разрушения существенно превосходит соответствующую меру для процесса формирование — зрелое состояние покрова. При этом вклад в эволюцию покрова диффузионного процесса, при котором происходит выравнивание площадей различных толщин, не столь заметен. Узкие границы доверительных интервалов свидетельствуют также о достаточно высокой степени надежности их определения.
Процедура оценки адекватности модели насчитывает ряд этапов, которые выполняются в определенной последовательности. В первую очередь необходимо выяснить, насколько верны представления о структуре информационных потоков в системе. Если они соответствуют реальности, то требуется также указание границ применимости модели для различных частных ситуаций. Действенным способом оценки адекватности является изучение статистических свойств искомых параметров. В рассматриваемом здесь случае узкие границы доверительных интервалов свидетельствуют об их значимом отличии от нуля. Особый интерес представляет рассмотрение их статистической состоятельности или решение вопроса, насколько изменяются значения оценок, когда они вычисляются на основании отдельных частей исходной выборки. Если имеет место существенное изменение значений, то эти оценки не имеют практического смысла, поскольку увеличение объема исходных наблюдений не будет прибли-
ah = (2,842 ± 0,209) • 10-3 аа = (2,486 ± 0,175) • 10-2 aahk, = (1,843 ± 0,162) • 10-3 T* = -6,3 ± 0,8
в
и ah = (8,051 ± 1,137) • 10-3;
и а'а = (5,980 ± 0,186) • 10-2;
и a'ahhl = (1,473 ± 0,125) • 10-3;
и т * = -7,6 ± 0,9;
Оценка адекватности модели
жать к истинному значению. Поэтому свойство состоятельности должно проверяться в первую очередь.
Для исследования состоятельности оценок исходная выборка разбивается на отдельные совокупности наблюдений, далее выполняется оценка параметров на основании каждой части, затем полученные результаты сравниваются между собой. В настоящем случае выборка разбита на две приблизительно равные по объему наблюдений части: одна содержит наблюдения за 1961-1975 гг., а другая — с 1976 по 1989 г. В качестве начальных значений выступает один и тот же фиксированный набор. Результаты приведены в табл. 1. Анализ ее элементов показывает, что оценки параметров при разбиении выборки меняются незначительно. Данное обстоятельство можно интерпретировать как факт статистической состоятельности вычисленных оценок. При этом границы доверительных интервалов оцениваемых параметров для отдельных выборок исходной совокупности данных существенно выше, поскольку при усреднении величин их дисперсии снижаются.
Таблица 1
Оценка статистической состоятельности параметров
Table 1
Estimation of a statistical solvency of parameters
Параметр 1961-1976 1977-1989 1961-1989
Ф(в) 26,339 25,465 19,724
ah х103 3,016 ± 0,473 2,731 ± 0,484 2,842 ± 0,209
ah х103 8,962 ± 2,626 7,488 ± 2,510 8,051 ± 1,137
aa х102 2,336 ± 0,315 2,401 ± 0,287 2,486 ± 0,175
<х102 5,811 ± 0,283 5,997 ± 0,301 5,980 ± 0,186
a ahh1 Х 103 1,769 ± 0,374 1,926 ± 0,418 1,843 ± 0,162
a'ah х 103 1,215 ± 0,410 1,384 ± 0,376 1,148 ± 0,125
Т * - 7,1 ± 1,4 - 6,8 ± 1,3 - 6,3 ± 0,8
Т * - 7,7 ± 1,6 - 7,8 ± 1,4 - 7,6 ± 0,9
Dah х 103 1,440 ± 0,156 1,299 ± 0,183 1,310 ± 0,021
bT х103 7,231 ± 2,168 9,785 ± 2,003 8,740 ± 1,096
bh х 103 1,989 ± 0,226 2,246 ± 0,209 2,185 ± 0,046
Примечание. После знака "±" следует среднеквадратичное отклонение параметра.
В качестве меры адекватности модельного и наблюдаемого значения суммарной площади льда в ¿-той декаде У-того года выступало принятое в практике исследований морского льда соотношение
р.= К>г /па>у) • 100 %,
где mdY — число оправдавшихся прогнозов суммарной площади льда; ndY — общее число прогнозов, которое можно оценить при наличии фактического материала. Значение PdY не ниже 68 %, что указывает на соответствие модели рассматриваемому процессу. Чтобы воспользоваться этой метрикой близости для сопоставления площадей льда в прибрежных районах, следует площади льда перевести в соответствующие сплоченности льда (относительных площадей льда). Для этого надо площадь льда поделить на площадь данного района. Результаты этого тестирования приведены в табл. 2.
Таблица 2
Оценка достоверности численных прогнозов сплоченности Р, %
Table 2
Estimation of reliability of numerical forecasts of unity Р, %
Номер декады года 1987 1988 1989 Среднее
S s S * S S s S * S S s S * S S s S * S
31 72 100 86 68 91 80 56 92 74 65 94 80
32 51 70 61 56 92 74 62 81 72 56 81 69
33 64 76 70 62 84 73 70 79 75 65 80 73
34 66 82 74 70 86 78 72 81 77 69 83 76
35 75 86 81 75 83 79 65 78 72 72 82 77
36 73 88 81 81 87 84 87 94 91 80 90 85
1 76 82 79 77 81 79 75 82 79 76 82 79
2 66 73 70 71 74 73 81 76 79 73 74 74
3 68 75 72 72 75 74 72 74 73 71 75 73
4 72 76 74 77 73 75 69 82 76 73 79 76
5 75 80 78 73 78 76 71 80 76 73 79 76
6 74 73 74 70 73 72 70 75 73 71 74 73
7 77 86 82 75 80 78 75 76 76 76 81 78
8 81 88 85 82 71 77 82 82 82 82 80 81
9 79 82 81 74 78 76 70 76 73 74 79 77
10 74 92 83 85 90 88 72 87 80 77 90 83
11 82 96 89 93 78 86 89 86 88 88 87 87
12 79 76 78 78 89 84 81 79 80 79 81 80
13 80 69 75 81 79 80 88 78 83 83 75 79
14 76 90 83 75 84 80 74 82 78 75 85 80
15 87 84 86 82 83 83 78 86 82 82 84 83
16 93 100 97 79 91 85 89 90 90 87 94 90
Среднее 79 87 83 79 87 83 79 86 82 79 86 83
Примечание. — сплоченность льда открытого моря, — сплоченность припая, 5 — общая сплоченность района.
При построении табл. 2 использовались наблюдения всех 114 районов расчетной области, из которых в 51 районе имеется еще и припай. Для каждого декадного наблюдения отдельного района вычислялась выбранная метрика при допустимой погрешности прогнозов ±1 балл, а затем значения суммировались.
171
Анализ данных табл. 2 показывает: модель эволюции площадей припая в большей степени, чем модель эволюции площадей открытого моря, соответствует исходным распределениям. Объяснением этого обстоятельства является специфичность припая как природного объекта: он более устойчив для различного рода внешних воздействий, что обусловлено физической средой его непосредственного расположения.
Особое значение в оценке адекватности модели выполняет сопоставление выборочных распределений в отдельных районах моря и соответствующих результатов численного моделирования. На рис. 3 представлены частные случаи такого сопоставления для общей площади и средних толщин льда на акватории. Их анализ показывает хорошую степень соответствия исходных и модельных распределений площадей покрова. Для каждого из них (рис. 3) коэффициенты корреляций между наблюдаемыми и модельными динамическими переменными определяются значениями соответственно 0,933, 0,945, 0,872 и 0,884. Необходимо отметить, что между модельными распределениями для зал. Петра Великого и Татарского пролива имеется определенная доля несоответствия в начальный период эволюции. Поскольку в модели максимальная толщина задавалась равной 1,5 м, то имеет место превышение толщин. Дальнейшее рассмотрение случаев отдельных районов также позволяет выявить подобные особенности. Однако оценку адекватности модели и границы ее применимости следует выполнить для всей совокупности выборочных распределений в отдельных районах.
Рис. 3. Наблюдаемые (1) и модельные (2) распределения площадей льда (А, С), средних толщин (B, D) в зал. Петра Великого (А, B) и в районе Татарского пролива (C, D) (распределения площадей льда приводятся в единицах площади открытого района моря; ось абсцисс day является временной осью, а ее начальная точка отсчета соответствует моменту начала эволюции ледяного покрова на акватории моря)
Fig. 3. Observable (1) and modelling (2) distributions of the areas of ice (А, С), average thickness (B, D, in meters) in Peter the Great Bay (A, B) and in area of Tatar Strait (C, D) (distributions of the areas of ice it is resulted in terms of the area of open area of the sea; the axis absciss day is a time base, and her index point of readout corresponds to the moment of the beginning of evolution of an ice cover on water area of the sea)
Границы применимости предлагаемых моделей непосредственно следуют из анализа карт изолиний корреляций между выборочными и модельными распределениями в отдельных районах моря (рис. 4).
131 136 141 131 136 141
Рис. 4. Карты изолиний корреляций между исходными и модельными распределениями: А — площади, B — толщины. Условные обозначения в тексте
Fig. 4. Isocline's maps for correlations between initial and modeling distributions (A — the areas, B — thickness). Signs are in the text
Временная дискретность исходных распределений составляет одну декаду. Поэтому для вычисления корреляций результаты моделирования были представлены декадными распределениями. Маркировка областей на рис. 4 соответствует трем градациям значений коэффициента корреляции. Так, области, где отмечается наибольшее соответствие (корреляция r от 1,0 до 0,8), помечены цифрой 3; если r изменяется от 0,6 до 0,8, то области помечены 2; цифрой 1 помечены области, где r ниже этих значений. Анализ этих изолиний указывает на более высокое соответствие средних модельных распределений площадей и их выборочных значений, нежели толщин покрова. Данная ситуация следует из того, что вариабельность толщин существенно выше вариабельности площадей (см. рис. 2). Действительно, в период уже стабильного состояния площадей покрова все еще наблюдается изменение его толщин. А разрушение толщи начинается ранее заметного разрушения площадей.
Выявленное хорошее соответствие указывает на отсутствие каких-либо причин для неприемлемости исследуемой модели на основании используемой выборки данных. В действительности никакое количество данных не может никогда доказать правильность модели: наибольшее, на что можно положиться, — эти данные не опровергают модель. Весь вопрос заключается в том, насколько привлеченная для исследований выборка отражает реальное состояние процессов. В настоящем рассмотрении представленные численные расчеты и эвристические
доказательства основаны на многолетнем фактурном материале, выполненный статистический анализ которого показывает, что он на самом деле является статистически представительным. Поэтому в данном случае не опровергается допущение о приемлемости предлагаемой модели и ее адекватности процессам эволюции площадей льда.
Вычислительные эксперименты
Вычислительные эксперименты условно подразделяются на два типа. Первые из них ориентированы на оценку значений параметров модели и установление их адекватности экспериментальным наблюдениям. Основное предназначение экспериментов второго класса заключается в оценке различного рода прогнозов долгосрочных изменений, которые могут вызываться естественными или искусственными обстоятельствами.
На рис. 5 и 6 представлены нормированные на площадь открытого района моря распределения толщины покрова в важном южном районе акватории (зал. Петра Великого) и самом северном районе моря (устье Татарского пролива).
A J А*/
Л
зю зет 50 100 150 day
А
А®«/Г
Л
310 360 50 100 150 с/эу
В
h .'■
Рис. 5. Распределения толщины покрова на акватории зал. Петра Великого: А — в открытой части, B — припая, С — совокупное распределение
Fig. 5. Distributions of thickness to water areas of Peter the Great Bay: A — in an open part, B — fast ice, C — cumulative distribution
0 31Г
О 1583
1.5Q
310
360
Я
С
100
™MaW 0.7S ь 0.0Ü
150 day
Рис. 6. Распределения толщины покрова в устье акватории Татарского пролива: А — в открытой части, B — припая, С — совокупное распределение
Fig. 6. Distributions of thickness in a mouth of water area of Tatar Strait: A — in an open part, B — fast ice, C — cumulative distribution
Анализ распределений толщины покрова в открытой части районов (рис. 5, А; 6, А) показывает, что при таянии покрова отмечается рост площадей льда в открытых участках акватории районов. Этот рост обусловлен выносом сюда обломков льдов припая.
Поскольку установлена адекватность модели, то ее (модель) можно использовать для проведения различного рода прогностических экспериментов. В частности, экспериментов, направленных на изучение последствий воздействия на ледяной покров Японского моря возможных климатических изменений.
При современном темпе выброса в атмосферу диоксида углерода С02 антропогенного происхождения вероятность таких изменений достаточно высока. Весьма вероятное повышение температуры обусловлено существованием так называемого "парникового эффекта" Земли: коротковолновое излучение Солнца нагревает Землю, а ее излучение лежит в области более длинных волн. Поэтому увеличение количества С02 значительно ограничивает выход собственного излучения Земли в космос и способствует повышению температуры. Так, согласно результатам математического моделирования системы "атмосфера—океан" (Наставление ..., 1982), при удвоении концентрации С02 в атмосфере ее современная средняя температура повысится с минус 19,2 оС до минус 17,54 оС. При этом уровень атмосферных осадков с 2,04 мм/сут увеличится до 2,15 мм/сут. Анализ данных измерений толщины льда с американских и английских подводных лодок в 29 точках Северного Ледовитого океана показал уменьшение средней толщины примерно на 43 % за период с 1958 по 1997 г. (Rothrock et al., 1999; Самарский, Михайлов, 2002).
Средняя толщина льда уменьшилась за этот период примерно с 3,5 до 2,0 м. По данным спутниковых снимков, площадь арктических льдов уменьшилась менее чем на 3 % за период с 1979 по 1999 г. (Vinnikov et al., 1999; Wadhams, Davis, 2000). Для ледяного покрова Японского моря подобные наблюдения отсутствуют. Интерес к данной проблеме очевиден по многим обстоятельствам. Прежде всего для проведения экологических экспертиз требуются различного рода количественные оценки антропогенного воздействия промышленных и производственных предприятий на окружающую среду. Кроме того, эти прогнозы полезны для комплексной оценки возможного состояния климатической системы.
Расчеты выполнялись при том же самом ветровом режиме 2-метрового надледного слоя воздуха, что был использован для оценки адекватности модели.
Температура атмосферы изменялась в сторону ее увеличения с шагом 0,5 оС. Результаты расчетов в форме графиков изменения общих объемов льда приведены на рис. 7. Согласно этим результатам, при повышении температуры на 1 оС общие потери объемов льда составляют 1,695 ' 1012 м3 (19 % от современного уровня), а при повышении на 5 оС — 6,092 ' 1012 м3 (69 % от современного уровня). Чтобы реально представить порядок этих огромных чисел, стоит перевести их в площадь открытого района моря. В первом случае суммарные потери составляют ледяной покров почти 249 таких районов, толщина льда в которых равна 1 м; во втором случае эта цифра соответствует покрову с площадью 826 районов.
V X Ю'10 т3
310 360 50 100 150
day
Рис. 7. Динамика объемов льда Японского моря при изменении температуры атмосферы (масштаб измерений объемов льда — 1010 м3): 1 — наблюдаемая; 2-4 — модельная (2 — при заданной температуре, 3 — при ее повышении на 1 °С, 4 — при ее повышении на 5 оС)
Fig. 7. Dynamics of volumes of ice of the Japan Sea at change of temperature of an atmosphere (scale of measurements of volumes of ice — 1010 m3): 1 — observed; 2-4 — model (2 — for fixed temperature, 3 — for increased in 1 оС, 4 — for increased in 5 оС)
Понятно, что полученные результаты дают грубое представление о возможных последствиях климатических изменений. Ведь при глобальном потеплении увеличивается уровень осадков, вследствие чего изменяется солевой режим моря, а стало быть, и условия и температура атмосферы, при которой первоначально формируется покров. Кроме того, возможны не учитываемые в модели обстоятельства. Однако даже и на данном уровне представлений может быть выполнена оценка этих последствий.
В заключение отметим, что полученные результаты основываются на том, что при адаптации морской воды к термическому и ветровому воздействию формируется морской лед, динамика параметров состояний которого определяется набором лимитирующих факторов. В частности, рост за счет теплопроводности толщины льда и снега приводит к уменьшению потока тепла через лед. В дальнейшем этот поток выходит на уровень величины потока тепла из воды в лед. Миграция изо льда в воду рассола (жидкой фазы с высокой концентрацией солей) также приводит к снижению потока тепла через лед и к осолонению подледных вод и повышению плотности морской воды, поэтому они трансформируются в лед при более низких температурах, чем температура уже сформированного слоя льда. При этом осолонение воды приводит к ее конвективному перемешиванию: на смену уже охлажденным соленым водам поднимаются более теплые глубинные, увеличивающие поток тепла к нижней поверхности льда. Сочетание этих и ряда других процессов останавливают нарастание толщины льда.
Построенная кинетическая модель эволюции ледяного покрова здесь используется как промежуточное звено для построения модели распределения толщины, однако представляется перспективным ее использование для детального изучения эволюции покрова замерзающих морей.
Литература
Айвазян С.А., Енюков И.С., Мешалкин Л.Д. Прикладная статистика. Исследование зависимостей. — М.: Финансы и статистика, 1985. — 487 с.
Аппель И.Л., Гудкович З.М. Численное моделирование и прогноз эволюции ледяного покрова арктических морей в период таяния. — СПб.: Гидрометеоиздат, 1992. — 144 с.
Бард И. Нелинейное оценивание параметров. — М.: Статистика, 1979. — 349 с.
Болч Б., Хуань К.Дж. Многомерные статистические методы для экономики. — М.: Статистика, 1979. — 317 с.
Ватутин В.А., Телевинова Т.М., Чистяков В.П. Вероятностные методы в физических исследованиях. — М.: Наука, 1985. — 208 с.
Волощук В.М. Кинетическая теория коагуляции. — Л.: Гидрометеоиздат, 1984. — 283 с.
Ивченко В.О., Хейсин Д.Е. Определение внутренних напряжений в ледяном покрове, возникающих при дрейфе льда // Проблемы Арктики и Антарктики. — 1974. — Вып. 43. — С. 84-89.
Кендалл М.Дж., Стьюарт А. Статистические выводы и связи. — М.: Наука, 1973. — 890 с.
Легеньков А.П. Деформации дрейфующего льда в Северном Ледовитом океане. — СПб.: Гидрометеоиздат, 1992. — 104 с.
Лушников А.А., Пискунов В.Н. Коагуляция в присутствии внешних источников // ДАН СССР. — 1976. — Т. 231. — С. 1403-1406.
Масловский М.И. Математическое моделирование короткопериодного ветрового дрейфа и перераспределения морского льда различной сплоченности (на примере Южного океана) // Тр. ААНИИ. — 1982. — Т. 387. — С. 116-136.
Международная символика для морских ледовых карт и номенклатура морских льдов. — Л.: Гидрометеоиздат, 1984. — 56 с.
Наставление по службе прогнозов. Раздел 3. Часть III. — Л.: Гидрометеоиздат, 1982. — 144 с.
Овсиенко С.Н. О численном моделировании дрейфа льда // Изв. АН СССР. Сер. физика атмосферы и океана. — 1976. — Вып. 12, № 11.
Плотников B.B. Изменчивость ледовых условий дальневосточных морей России и их прогноз. — Владивосток: Дальнаука, 2002. — 172 с.
Пригожин И., Стенгерс И. Порядок из хаоса. Новый диалог человека с природой. — М.: Едиториал, УРСС, 2003. — 312 с.
Рао С.Р. Линейные статистические методы и их применения. — М.: Наука, 1968. — 547 с.
Самарский А.А., Михайлов А.П. Математическое моделирование: Идеи. Методы. Примеры: 2-е изд., испр. — М.: Физматлит, 2002. — 320 с.
Стронгин Р.Г. Численные методы в многоэкстремальных задачах (информационно-статистические алгоритмы). — М.: Наука, 1978. — 240 с.
Тимохов Л.А., Хейсин Д.Е. Динамика морских льдов (математические модели). — Л.: Гидрометеоиздат, 1987. — 272 с.
Фукунага К. Введение в статистическую теорию распознавания образов. — М.: Наука, 1979. — 368 с.
Четырбоцкий А.Н. Идентификация моделей вертикального распределения плотности океанских водных масс // Информатика и моделирование в океанологических исследованиях. — Владивосток: Дальнаука, 1999. — С. 131-143.
Четырбоцкий А.Н. Локальная эволюция толщины ледяного покрова водных поверхностей // Тр. ТОВМИ им. адм. Макарова. — 2001. — Вып. 23. — С. 117-123.
Четырбоцкий А.Н., Плотников B.B. Ледяной покров Японского моря: исходные данные и процедуры восстановления пропущенных значений // Электронный журнал "Исследовано в России". — 2003. — № 7. — С. 88-93. http://zhurnal.ape.relarn.ru/ articles/2003/007.pdf
Якунин Л.П. Количество льда и затраты тепла на его таяние в Дальневосточных морях // Проблемы Арктики и Антарктики. — 1986. — Вып. 62. — С. 93-96.
Rothrock D.A., Yu Y., Maykut G.A. Thining of the sea-ice cover // Geoph. Res. Letters. — 1999. — Vol. 26. — P. 3469-3472.
Vinnikov K., Robock F., Stouffer R.J. Global warming and northern hemisphere sea ice extent // Science. — 1999. — Vol. 286. — P. 1934-1937.
Wadhams P., Davis N.R. Futher evidence of ice thining in the Arctic Ocean // Geoph. Res. Letters. — 2000. — Vol. 27. — P. 3973-3975.
Поступила в редакцию 22.07.05 г.
