Водные ресурсы, 2023, T. 50, № 6, стр. 751-766

ВВЕДЕНИЕ

Несмотря на то, что вопросам изменений в климатической системе в последнее время уделяется повышенное внимание (что отражено, например, в обширной библиографии в последнем отчете Международной группы экспертов по изменениям климата IPCC (Intergovernmental Panel on Climate Change)) [13] – колебания составляющих водообмена в глобальном пространственном масштабе и в климатических временны́х масштабах изучены совершенно недостаточно. Так, остается неисследованной следующая важнейшая проблема. Почему при несомненном росте глобальной температуры воздуха за последние 100–150 лет и кажущемся неизбежным соответствующем росте испарения с океана и интенсификации всего глобального водообмена – значительного увеличения стока большинства рек мира не наблюдается [2]. Последний вывод сформулирован в работах [2, 3] на основе анализа многолетних изменений годового речного стока в отдельных створах и суммарного стока в океан. Так, исследование всех доступных длинных (более 40 лет) около 1000 рядов стока рек основного типа (не озерных и не зарегулированных) показало, что только в 4% створов наблюдаются статистически значимые монотонные многолетние тренды [3]. Изучение ряда величин глобально суммированного стока рек в океан длительностью 87 лет [2] приводит к выводу о стационарности и этих изменений. Аналогичные выводы делаются во многих посвященных этому вопросу работах, их анализ приведен в книге [2]; например, авторы фундаментальной монографии [1] (с. 131) делают вывод о том, что “значимых трендов водных ресурсов рассматриваемых речных бассейнов и бассейнов морей не существует”.

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

Определенные общие закономерности изменений компонентов глобального водообмена могут быть получены при помощи трактовки результатов расчетов на детальных гидродинамических моделях климатической системы – “климатических моделях” (КМ), описывающих большое количество процессов в океане, атмосфере и на поверхности суши. Сравнительные результаты расчетов по более чем 50 КМ представлены в рамках проектов CMIP-5 и CMIP-6 и имеются в открытом доступе. Особенный интерес представляют так называемые “исторические” эксперименты на указанных моделях, воспроизводящие изменения различных модельных гидрометеорологических характеристик на протяжении периода времени с середины XIX в. вплоть до второго десятилетия XXI в. В качестве “внешних” сигналов (форсинга) в таких моделях учитываются два основных типа процессов: медленное накопление парниковых газов в атмосфере и быстрые (“взрывные”) вулканические извержения. Процессы других временны́х масштабов, как правило, рассчитываются внутри моделей, причем все модели приблизительно сходно реагируют на внешние воздействия указанных двух типов, демонстрируя сходные хронологические ответы на эти сигналы. В свою очередь, основным “внутренним” процессом, влияющим на поведение климатических характеристик моделей, является так называемое “погодное возбуждение” климатической системы (по К. Хассельманну [12]) – результат воздействия сравнительно быстрых атмосферных процессов синоптического масштаба на более инерционные компоненты системы: океан, поверхность суши и покровное оледенение. В работах [12, 14] и в ряде других работ сотрудников и соавторов К. Хассельманна на основе данных наблюдений, палеореконструкций, численных расчетов с помощью стохастических моделей, а также в результате теоретического анализа было убедительно показано, что в глобальном пространственном масштабе и в широком диапазоне климатических масштабов времени “погодное возбуждение” приводит к вариациям средней глобальной температуры по типу случайных блужданий (или винеровских процессов). С учетом этой концепции различные модели экспериментов CMIP-6 должны давать различные хронологические отклики на “погодное возбуждение”– в отличие от указанных выше двух типов “внешних” сигналов. Причиной этого является то, что модели в “исторических экспериментах” не воспроизводят реальную хронологию синоптических изменений.

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

Важным преимуществом детерминистических гидродинамических климатических моделей экспериментов CMIP-5 и CMIP-6 является то, что ансамбль моделей – в отличие от инструментальных наблюдений и данных реанализа – позволяет воспроизводить большое количество различных “реализаций” эволюции климатической системы, т. е. большее количество теоретически возможных вариантов поведения системы в предположении о том, что внешние сигналы для прогона каждой модели остаются одними и теми же, а эффекты “погодного возбуждения” различны. Таким образом, всю совокупность прогонов разных детерминистических моделей из общего ансамбля моделей на протяжении исторических экспериментов можно рассматривать как совокупность выборочных траекторий случайного процесса. Как следует из теории стохастических моделей климата К. Хассельманна, осреднение бесконечного числа реализаций тех изменений (винеровского характера) глобально осредненных параметров климатической системы, причина которых – погодное возбуждение, в масштабах времени порядка десятков и сотен лет должно дать отсутствие изменений за счет этого фактора. Так как число используемых авторами статьи моделей ограничено, трудно ожидать полного подавления глобальных изменений вследствие “погодного возбуждения” при осреднении по всем моделям, однако такое осреднение способно сильно ослабить указанные эффекты – и таким образом более четко выявить роль парникового эффекта.

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

Цели настоящей работы можно сформулировать следующим образом:

1) изучение при помощи ансамблевых оценок на моделях климатической системы особенностей стационарного/нестационарного поведения различных составляющих глобального водообмена в климатических масштабах времени – от межгодовых колебаний до изменений на протяжении полутора столетий;

2) выявление причины на этой основе отмеченного выше парадокса – отсутствия значимого монотонного тренда в глобальном речном стоке при несомненном росте средней глобальной температуры и тренда в среднем испарении с поверхности океана;

3) выявление также механизма, за счет которого спектральные плотности изменений глобального речного стока имеют в основном “белошумный” характер – неизменны по частоте;

4) в целом – предложение непротиворечивой гипотезы относительно основных особенностей и закономерностей многолетних изменений компонентов глобального водообмена на основе “исторических” модельных экспериментов.

ДАННЫЕ МОДЕЛЬНЫХ РАСЧЕТОВ

Детали использования результатов моделей проекта CMIP-6 заключаются в следующем. Проект CMIP-6 [19] объединяет множество как глобальных, так и региональных (LAM) гидродинамических моделей. Из этого ансамбля моделей рассмотрены лишь те глобальные модели, результаты (поля) которых представлены для общего пользования, т. е. имеющие на сайтах (нодах) проекта ESGF (Earth System Grid Federation) [11] данные о глобальном поле испарения/эвапотранспирации, осадках (общего типа) и общем стоке.

Так как сам проект CMIP-6 не задает жестких рамок и правил, лишь 48 глобальных моделей предоставили эти поля на сегодняшний день. В дальнейшем возможно расширение этих данных и поступление сведений от других моделей. Также в настоящее время еще продолжается работа в рамках предыдущего проекта, CMIP-5, и эти данные (поля) также могут быть использованы. Результаты обоих проектов в целом аналогичны, но точки доступа к информации находятся в разных базах данных.

В рамках представляемой работы для сопоставления используемых моделей между собой была осуществлена интерполяция различных модельных глобальных полей на единую гауссову сетку T62 с треугольным усечением. Ее характерное разрешение – 2 × 1.7 градуса. Более плотные сетки для интерполяции были использованы авторами статьи ранее в других работах, например [5].

“Маска” суши и океана, использованная в проекте, показана на рис. 1. Результаты тестовых экспериментов с интерполяцией данных на маску T255 или T254 (с разрешением ~0.5 градуса) не изменяют существенно полученные результаты (средние глобальные значения), но намного увеличивают время множественных расчетов. Для интерполяции, выборки, интегрирования по поверхностям суши и океана использовался пакет программ CDO Немецкого метеорологического института им. Макса Планка [15], весьма удобный для быстрой пакетной обработки результатов численного моделирования и рекомендованный американской университетской корпорацией атмосферных исследований [17].

В конечном итоге для анализа использовались данные всех доступных моделей климата проекта CMIP-6, для которых удалось обобщить информацию, полученную в результате “исторических экспериментов” по составляющим водообмена. Наборы параметров, для которых были доступны данные исторических экспериментов, несколько различались. Кроме того, несколько моделей были исключены из рассмотрения после анализа результатов, так как они давали кардинально отличающиеся от остальных моделей оценки. Наконец, при сопоставлении ансамблевых оценок по различным параметрам, а также для оценки расчетных параметров необходимы были наборы моделей, по которым имелись результаты в отношении всех требуемых комбинаций параметров. Например, для оценки изменений годовых величин испарения с поверхности океана использовались результаты 41 модели, для оценки изменений осадков над океаном и “эффективного” испарения с океана (разности “испарение минус осадки”) – 40 моделей, годового слоя стока с суши – 34 модели.

В итоге, в работе использовались результаты экспериментов на 41 модели проекта CMIP-6. Информация о моделях приведена в табл. 1.

Почти во всех случаях (за исключением 1–3 моделей) использовались данные “исторических экспериментов” за одинаковый промежуток времени: с 1850 по 2014 г., всего 165 лет. Для всех моделей скачивались также месячные данные по каждому из параметров, однако в связи с ограниченностью объема статьи в нее не включены результаты их анализа.

По унифицированной сетке, соответствующей “маске” на рис. 1, результаты численного моделирования приходилось интерполировать, поскольку разные модели используют разные расчетные сетки: “гауссову”, “географическую” (“lonlat”) или “криволинейную” (“curvilinear”). Интерполяция и выборка данных, суммирование стока по площади проводились с помощью пакета программ CDO [9].

Для расчетов среднего, максимального и минимального стока использовались результаты моделирования с месячным осреднением, которые загружались с сайта PCMDI (Program for Climate Model Diagnosis and Interpretation) [16].

МЕТОДЫ АНАЛИЗА РЯДОВ ДАННЫХ

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

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

Авторами предложена новая система статистического и стохастического анализа сегментов реализаций стационарных случайных последовательностей, приводящая в том числе к решению вопроса о степени нестационарности изучаемых рядов и основанная на принципах теории случайных функций (например, [7, 8, 20, 21]). Так как имеющиеся ряды, тем более их сегменты, сравнительно коротки, возникает задача получения оценок параметров стохастических моделей, превосходящих по точности обычно используемые методы.

Использованные варианты методов анализа временны́х рядов подробно описаны в работах [2–4]. Вкратце напомним их суть. В связи с тем, что теория случайных функций и ее приложения дают наиболее адекватные результаты в случае гауссовского характера изучаемых сегментов, целесообразно преобразовывать исходные сегменты рядов в сегменты выборочных значений нормально распределенных величин (такого рода анализ дублировался авторами анализом рядов в их исходном виде). С этой целью авторами был предложен новый алгоритм соответствующего пересчета, который оказался намного более эффективным, чем широко распространенные алгоритмы такого плана, например – алгоритм Вычислительного центра (ВЦ) РАН. В частности, проверка на 1 млн сгенерированных равномерно распределенных (т. е. максимально отличающихся от гауссовских) псевдослучайных чисел показывает, что среднее значение (т. е. смещенность) полученных этим методом чисел равно +0.001734, что по модулю на порядок меньше, чем при использовании алгоритма ВЦ РАН. Стандарт оказался равен единице с точностью до шестого знака после запятой (ошибка по крайней мере на 4.5 порядка меньше), асимметрия –0.029977 (в 5 с лишним раз лучше по модулю, чем с применением алгоритма ВЦ РАН).

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

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

смещенность среднего можно устранить практически полностью – в указанном эксперименте ее величина <0.000001;

стандарт величин можно привести к единичному – его отличие от единицы в эксперименте также <0.000001;

коэффициент асимметрии уменьшается по модулю в 390 раз по сравнению с алгоритмом ВЦ РАН и составляет в данном случае –0.000391;

эксцесс относительно нормального уменьшается в >2 раза по сравнению с алгоритмом ВЦ РАН и составляет по модулю <1/30 от эксцесса гауссовского распределения (–0.098242).

В связи с тем, что сегменты изучаемых рядов рассматриваются как сегменты реализаций стационарных случайных процессов, для их анализа в работе использован метод наибольшей энтропии (МНЭ). В классической работе [18] показано, что максимизация энтропии, порождаемой отсутствующей информацией (отсутствующими частями анализируемых реализаций, предполагаемых бесконечными в теории случайных функций), и минимизация энтропии, порождаемой имеющейся информацией, приводят к оптимальному представлению рядов в виде отрезков реализаций процессов авторегрессии конечного порядка (АР):

где X – моделирующая изучаемый процесс случайная последовательность; t – целочисленное безразмерное время (в данном случае – безразмерные годы); a(t) – некоррелированная по времени одинаково распределенная случайная величина. Ключевым моментом является выбор порядка авторегрессии М и оценка коэффициентов авторегрессии К_i – прежде всего первого – близкого по смыслу к коэффициенту автокорреляции для единичного сдвига по времени. Решение этих вопросов заключалось в следующем.

На основе предложенного авторами нового алгоритма генерирования гауссовских псевдослучайных величин и нового алгоритма пересчета исходных рядов в ряды выборочных значений нормально распределенных величин [2, 3] методом Монте-Карло получены новые формулы для оценок статистических и стохастических параметров рядов. Пробные параметры процессов АР оценивались авторами при помощи алгоритма Юла–Уокера с использованием рекурсивных формул Дарбина–Левинсона. Кроме того, рассмотрение изменений многочисленных гидрометеорологических параметров в климатических масштабах времени показало, что преобладающую роль играют стохастические модели нулевого (“белый шум”) или первого (“простая марковская цепь”) порядка. И в том, и в другом случае решающей является оценка первого пробного коэффициента авторегрессии, на практике вычисляемого как коррелированность смежных величин ряда. Авторами предложена новая формула для расчета этого ключевого при идентификации модели параметра – выборочной оценки (r_1выб) величины коррелированности между смежными значениями ряда (N – длина ряда в безразмерных единицах времени). Соответствующее выражение имеет вид поправки к обычной, “несмещенной” оценке:

Сравнение оценок (2) с широко используемыми в гидрологии оценками Резниковского–Костиной показывает, что первые обладают существенно меньшей смещенностью. Оценки Резниковского–Костиной предпочтительнее только для значений коэффициента автокорреляции ~ 0.9, которые почти не встречаются на практике изучения предполагаемых стационарными сегментов рядов. Во всех остальных случаях – при значениях коэффициента от –0.9 до 0.8 – предлагаемые в настоящей работе оценки обладают преимуществом: для рядов длиной 20 ординат средняя смещенность оценок авторов статьи по всем коэффициентам меньше в 27 раз, а для рядов длиной 40 ординат – в 55 раз.

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

Также предложена новая модификация метода Акаике оценки порядков стохастических (авторегрессионных) моделей для описания изучаемых рядов. Эта модификация устраняет основной недостаток существующих методов: сильную зависимость оценок от длины рядов наблюдений. В предлагаемом подходе все длины рядов условно приравниваются к единой (стандартной), технически такое “уравнивание” достигается умножением критерия Акаике при нулевом порядке на специально подобранный корректирующий множитель (r₁ – аналог r_1выб):

Примеры значений коэффициентов а для различных областей значений N и r₁ приведены в монографии [3] (с. 25).

В качестве примера эффективности предложенного метода оценки порядков процесса авторегрессии в упомянутой выше монографии [3] приведено сравнение порядков, рассчитанных при помощи критерия Акаике и при помощи предлагаемого варианта критерия. Для сравнения использованы все (764) имевшиеся в распоряжении длинные (не менее 20 лет) ряды годового стока на территории Российской Федерации. Сравнение показало, что использование традиционного критерия Акаике давало искусственное увеличение порядка модели авторегрессии на очень большую, критичную величину (0.7) при переходе от наименьших длин к наибольшим рядов, в то время как применение предложенного авторами статьи критерия демонстрировало неизменность его оценок.

Для изучения вопроса о стационарности изменений стока были предложены новые двусторонние критерии – раздельно для стационарности/нестационарности по среднему значению (математическому ожиданию М_выб), стандарту (среднеквадратичному отклонению σ_выб), автокоррелированности r_выб. Указанные критерии основывались на оценке стандартов разностей между выборочными средними значениями указанных параметров для первой и второй половин рядов (длиной не менее 40 лет) соответственно, каждая из которых предполагается отрезком реализации стационарного процесса. Формулы для них были получены путем построения несколькими методами трехмерных поверхностей зависимостей искомых величин от специально подобранных функционалов ключевых параметров и последующей аппроксимации этих поверхностей многочленами второго порядка. Так, формула для М_выб выглядит следующим образом:

На основе формулы (4) и аналогичных ей рассчитывались значения соответствующих критериев стационарности/нестационарности рядов, пересчитанных в выборочные значения гауссовских случайных величин: по математическому ожиданию I_SM, по стандартам (дисперсиям) I_SS, по автокоррелированности I_sr. Так, формула для I_SM выглядит следующим образом:

где $М_{{{\text{выб}}}}^{2} - М_{{{\text{выб}}}}^{1}$ – разность между выборочными средними значениями второй и первой половин ряда стока, вычисляемыми по обычным формулам; $\sigma {\text{(М}}_{{{\text{выб}}}}^{{{\text{2,МК}}}} - {\text{М}}_{{{\text{выб}}}}^{{{\text{1,МК}}}}{\text{)}}$ – стандарт аналогичной разности, оцененный методом Монте-Карло для стационарной последовательности соответствующей длины и временнóй коррелированности – и аппроксимированный аналитическим выражением (4).

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

ОСНОВНЫЕ РЕЗУЛЬТАТЫ

ЗАКЛЮЧЕНИЕ

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

Анализ 41 модельного ряда годовых величин суммарного испарения с поверхности океана в общем подтверждает предположения о росте испарения за последние полтора века вследствие роста глобальной температуры: осредненное по всем моделям испарение обнаруживает детерминированный монотонный положительный тренд. Если рассматривать индивидуальные модельные реализации, то априорная вероятность повышения испарения составляет ~2/3, вероятность уменьшения испарения составила 1/3 – очевидно, в основном вследствие внутренних процессов типа “погодного возбуждения”. Указанную выше вероятность повышения можно объяснить в равных пропорциях воздействием глобального потепления, с одной стороны, и “повышательным” воздействием внутренних климатических процессов винеровского типа – с другой.

Привлеченные к анализу результаты исторических экспериментов позволили авторам ответить на важный вопрос: почему вопреки современному росту испарения с поверхности океана не наблюдается значительного увеличения стока большинства рек мира? Причиной, по-видимому, является механизм преобразования испарения в ‘эффективное испарение” – разность между испарением с океана и осадками над океаном.

Разность “осадки минус испарение” над океаном неизбежно попадает на сушу. Метод оценки суммарного стока с суши при помощи расчета эффективного испарения с океана, таким образом, можно рассматривать как обратный так называемому “аэрологическому методу в гидрологии” (АМГ). Если АМГ заключается в расчете баланса влаги в пределах определенного бассейна на суше при помощи расчета результирующего горизонтального потока влаги в атмосфере через контур этого бассейна, то в исследуемом случае баланс влаги на поверхности океана используется для расчета результирующего потока через береговую линию материков. Сравнение модельных рядов эффективного испарения с океана и рядов модельного суммарного стока с суши показывает, что коэффициенты парной корреляции чрезвычайно велики – в среднем 0.85 для 34 моделей с сопоставимыми данными. Выводы, сделанные выше об эволюции эффективного испарения с океана, следовательно, применимы и к суммарному стоку с суши. По-видимому, именно осадки над океаном оказывают стационаризирующую роль в отношении испарения с океана и в очень значительной степени отфильтровывают монотонный тренд вследствие глобального потепления в том, что касается глобальных изменений речного стока. Осадки над океанами могут также играть существенную роль в белошумном (не меняющемся вдоль оси частот) характере спектральной плотности изменений годовых величин глобально суммированного стока [2].

Выводы относятся к глобально и ежегодно осредненным величинам. Можно сформулировать следующую гипотезу: региональные и внутрисезонные особенности изменений водообмена имеют более сложную и более мозаичную структуру, которая впоследствии также может быть детально исследована при помощи “исторических” экспериментов на глобальных гидродинамических климатических моделях.