РЕГИОНАЛЬНАЯ ГИДРОДИНАМИЧЕСКАЯ МОДЕЛЬ
краткосрочного прогноза полей метеорологических и осадков на срок до 48 часов


Общие сведения
В Гидрометцентре России разработана региональная модель гидродинамического прогноза полей основных метеорологических величин высокого пространственного разрешения (шаг по горизонтали не более 50 км) по территории Евроазиатского континента в s - системе координат и метод прогноза осадков синоптического масштаба, включая фронтальные большой интенсивности, на срок до 48 часов. Модель включена в оперативную технологию регионального прогноза Гидрометцентра России.

Новая региональная неадиабатическая модель имеет 30 расчетных уровней по вертикали, что позволяет усваивать исходные данные о ветре и геопотенциале изобарических поверхностей 30, 50, 70 гПа в нижней стратосфере и более детально, за счет улучшения аппроксимации по вертикали, описывать пограничный слой атмосферы.

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

При реализации прогностической модели в s - системе координат решены и другие возникшие проблемы. Это выбор рельефа для данной территории, переход из p - системы координат в s - систему и обратно, разложение на вертикальные моды и инициализация начальных данных, нелинейная вычислительная неустойчивость, аппроксимация барического градиента в уравнениях движения (особенно в районах с большим наклоном рельефа). Их решение предшествовало учету параметризации радиации, пограничного слоя и других неадиабатических процессов.

В модели интегрирование уравнений по времени осуществляется методом центральных разностей на сетке — "С" Аракавы для территории, имеющей вид прямоугольника на карте стереографической проекции с количеством узлов 137Х209 и шагом по горизонтали 75 км. Область расчета охватывает Европу, большую часть Азии, включая территории всех стран СНГ. Расчетная схема сделана достаточно гибкой, т.е. заданием входных параметров изменяется количество расчетных s - уровней по вертикали, шаг по горизонтали, область прогноза. Реализованный вариант модели требует для расчета на ЭВМ CRAY 30 мин.

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

Начальные поля геопотенциала изобарических поверхностей, зональной и меридиональной составляющих скорости и ветра, влажности воздуха задаются на стандартных изобарических поверхностях 10, 30, 50, 70, 100, 150, 200, 250, 300, 400, 500, 700, 850, 925 и 1000 гПа. Расчет по модели ведется для 30 s - уровней, на которые интерполируются исходные данные. По окончании прогноза осуществляется переход к Р - уровням. Проблема прямой и обратной интерполяции непроста и во многом от ее успешного решения зависит адекватность учета рельефа.

Система уравнений
Переход к s - системе координат означает, что интегрирование уравнений модели будет проводиться на спрямленной относительно рельефа поверхности. Высота рельефа задается как функция x,y(Г=Г(x,y)). Вертикальная координата s=p/ps , где ps - давление на земной поверхности, которое в системе уравнений модели является эволюционной величиной. Очевидно, что уровень s=1 совпадает с поверхностью Земли, на этом уровнем значение геопотенциальной высоты H  всегда равно высоте рельефа.

Исходными уравнениями в s - системе координат на плоскости стереографической проекции являются два уравнения движения, записанные в форме Громеко-Лэмба, которые описывают эволюционные изменения компонентов скорости ветра u и v, уравнение притока тепла, уравнение переноса влажности, уравнение неразрывности, и уравнение гидростатики в стереографической системе координат. Кроме этого в модели имеются блоки параметризации радиации, пограничного слоя атмосферы и расчета осадков.

Область интегрирования по горизонтали — прямоугольник на карте стереографической проекции, по вертикали — от s=0 до s=1.







Чтобы алгоритм не стал громоздким, здесь не вошли блоки параметризации радиации и пограничного слоя атмосферы.

Для заблаговременностей прогнозов 12, 18, 24, 36, 42 и 48 ч проводится обратная интерполяция значений метеорологических величин в о - точках расчетной сетки с s - уровней на стандартные Р–уровни. Прогностические результаты оперативной версии записываются в базу данных «Прогноз».


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


Подготовка начальных данных и обработка результатов расчета
Начальными данными для прогноза являются результаты оперативной схемы объективного анализа (ОА) полей метеорологических величин на полусфере, которые берутся из базы данных «Прогноз». Из сферической сетки с шагом по широте и долготе 2,5° проводится интерполяция на сетку области регионального прогноза в стереографической проекции. Для каждой вертикали в о - точках (рис.1а) значения метеорологических величин на поверхностях 30, 50, 70, 100, 150, 200, 250, 300, 400, 500, 700, 850, 925 гПа и p0 интерполируются на σ-уровни: 0.05,0.07,0,085,0.1,0.15,0.2, 0.25, 0.3, ......1. Всего 30 уровней, которые сгущаются в пограничном слое атмосферы, то есть около s=1.

Функция рельефа Г(x, y) – задана. Она должна быть достаточно гладкой, поскольку в s-системе рельеф неявно входит в Н и Т ln ps и при дифференцировании этих величин по х берутся производные от рельефа. Переход из изобарической в s-систему координат и обратно осуществляется с помощью сплайн-функций.

Нелинейная неустойчивость
Весь описанный алгоритм может быть неработоспособен при наличии нелинейной неустойчивости. В баротропном варианте модели нелинейной неустойчивости нет, в ней достаточно применить фильтр Асселини и учесть критерий Куранта, чтобы не возникло вычислительных проблем. В бароклинной модели этого недостаточно.



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


Проблема барического градиента



Аналитически последнее выражение тождественно равно нулю. В схеме и первый и второй члены дифференцируется в конечных разностях. Как показано в работе /5/ для сетки "С" Аракавы остаточный член аппроксимации выражения в скобках равен   

При шаге по горизонтали величиной менее 150 км и реальных f' остаточный член невелик. Он может казаться только над горными районами. Тем не менее, можно записать барический градиент в форме, при которой аналогичная ошибка отсутствует, а именно:



Аналогично записывается барический градиент по у.
Очевидно, что при s=1 первые два члена нового выражения взаимно уничтожаются. Установлен факт: если рельеф не сглаживать, то при обычном расписывании барического градиента возникает нелинейная неустойчивость, а при применении последней формулы ее нет.


Инициализация начальных полей геопотенциала изобарических поверхностей, ветра и приземного давления
Основная идея инициализации принадлежит Махенхауеру /9/. Решение бароклинной системы уравнений атмосферной гидро- термодинамики может быть разложено на вертикальные моды — собственные решения соответствующего линейного дифференциального оператора. Коэффициенты разложения по вертикальным модам описываются системой трех эволюционных уравнений, аналогичной той, которой удовлетворяют геопотенциал изобарических поверхностей и два компонента скорости ветра в баротропной модели атмосферы. «Баротропные» системы, отвечающие разным модам, различаются одним параметром — скоростью распространения гравитационных волн, поэтому метод нелинейной инициализации сначала разрабатывается для баротропной модели, а в бароклинной модели инициализации подвергаются указанные коэффициенты разложения.

Решение линейного оператора баротропной задачи представлено в виде суммы решений, два из которых (быстрые моды) отвечают гравитационным волнам и одно (медленная мода) — волне Россби. Удовлетворение требованиям стационарности операторов, найденных для быстрых мод, и неизменности в процессе инициализации оператора, отвечающего медленной моде, позволяет так изменить начальные поля ветра и геопотенциала изобарических поверхностей, что при интегрировании по времени амплитуды быстрых гравитационных мод не будут расти. Этим требованиям могут быть удовлетворены при организации специальной итерационной процедуры. Каждая итерация состоит из двух этапов: нахождение производных по времени по начальным для итерации полям скорости ветра и геопотенциала изобарических поверхностей (делается шаг по времени численной прогностической схемы), которые входят в правую часть некоторой линейной системы уравнений, и решение этой системы. Определяющие изменения указанных полей производятся на этой итерации. Этот алгоритм описан в работах /1/ и /2/. Он оказался быстро сходящимся — изменения ветра и геопотенциала изобарических поверхностей уже на десятой итерации были малы (в сравнении с первой итерацией).

Инициализация позволяет избежать ложных колебаний во временных изменениях метеорологических полей с периодом порядка 6 ч, особенно заметных в первые 12 ч прогноза. Например, геопотенциал на уровне 500 гПа в одной из центральных точек расчетной области без инициализации начальных данных имеет амплитуду ложных колебаний порядка 4-6 дам при характерной изменчивости 4 дам/сут. После инициализации такие колебания практически отсутствуют. Плавность расчетных тенденций позволит более эффективно подойти к проблеме усвоения данных и параметризации процессов фазовых переходов атмосферной влаги.


ПАРАМЕТРИЗАЦИЯ НЕАДИАБАТИЧЕСКИХ ФАКТОРОВ В РЕГИОНАЛЬНОЙ МОДЕЛИ
Одной из важнейших величин, которую можно уточнить, привлекая параметризацию неадиабатических факторов, является температура поверхности суши, моря и ледяного покрова Ts . Ее расчет приводится с помощью уравнения теплового баланса:


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

Физика радиационных процессов реализована в виде процедур для расчета радиационных притоков к атмосферным слоям и для расчета компонент радиационного баланса

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

Предварительно для каждого из слоев (промежуток между расчетными уровнями) рассчитывается водность облаков (если она есть), а также эффективные поглощающие массы H2O, CO2 и O3 в каждом из спектральных участков с учетом температуры и давления в слое.

Далее, с помощью интегральных функций пропускания атмосферных газов H2O, CO2 и O3, вычисляется эффективный поток длинноволновой радиации F(z). Расчет этой функции зависит также от наличия облачности. После вычисления F(zi) определяем длинноволновый радиационный баланс подстилающей поверхности

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


Параметризация радиационных процессов
При увеличении числа уровней по вертикали до 30 в региональной модели расчет радиационных потоков для коротковолнового (видимого и инфракрасного (ИК)) диапазона спектра, содержащего более 90 % всей лучистой энергии Солнца, становится настолько громоздким и сложным, что прогноз с учетом неадиабатических факторов в реальном масштабе времени становится проблематичным даже на супер-ЭВМ CRAY. Приходится делать упрощения, в результате которых понижается точность расчетов.

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

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

Расчет радиационных потоков видимой части спектра
Наиболее наглядно рекуррентный метод просматривается для интервала длин волн видимого излучения (0,3 мкм ≤ λ≤ 0,75 мкм), составляющего примерно 47% от суммарного потока солнечной радиации. Для этого диапазона каждый слой частично пропускает, частично отражает поток солнечной радиации, но нет источников и стоков. Поэтому важно рассчитать какая часть падающего сверху потока дойдет до поверхности Земли.

Расчет потоков инфракрасного излучения
Более сложным является расчет потоков инфракрасного (ИК) излучения (0,75 мкм ≤ λ≤ 4 мкм), составляющего 44% приходящей от Солнца радиации. В этом участке спектра учитываются поглощение в облачных и безоблачных слоях, а также отражение, но только от облачных слоев и от поверхности Земли. Для этого диапазона длин волн не рассматриваем вопрос определения альбедо, считая его значение известным.

Расчет потоков длинноволнового излучения
Расчет потоков теплового излучения (4 мкм ≤ λ≤ 50 мкм) на данном этапе по существу не содержит новых разработок. При определении потоков этого диапазона длин волн учитывается только поглощение. Принимаем, что для этого участка спектра выполнены известные соотношения, их вид громоздок и здесь не приводится. Отметим только, что длинноволновые потоки существенно зависят от удельных плотностей водяного пара, водяного конденсата (в наибольшей степени), углекислого газа и озона.


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

где f - любая из функций u,v,T,q, соответственно компоненты скорости ветра, температура и влажность. Коэффициент турбулентности (Kf) является функцией координат и времени. Его нахождение — основная задача теории пограничного слоя атмосферы. В основе алгоритма лежит метод решения эволюционного уравнения для функции Е - энергии турбулентных пульсаций, являющейся при учете ПСА такой же переменной, как u,v,T,q. Это уравнение имеет второй порядок по вертикальной координате z, поэтому требует два граничных условия: сверху и снизу.





Блоки данного алгоритма не представляют особой сложности кроме громоздкости выражений. В целом схема прогноза полей метеорологических величин с учетом параметризации радиации и пограничного слоя имеет вид:
1. На уровне z1 задаются все искомые функции.
2. Рассчитываются радиационные потоки видимого, инфракрасного и длинноволнового излучений, после чего из уравнения теплового баланса находятся приземные значения температуры и влажности.
3. Решается задача для приземного подслоя. Находится масштаб турбулентности и через него динамическая скорость u*, потоки тепла, импульса и влаги на уровне z1 и граничные значение для Е - энергии турбулентных пульсаций.
4. Делается шаг по времени для Е, по Е находятся Kf (f - u,v,θ,q).
5. Зная Kf, делается шаг по времени для u,v,θ,q.
6. Повторяются пп. 2-5, пока не завершится прогноз на заданный интервал времени.


РАСЧЕТ ОСАДКОВ
При расчете осадков реализована идея Санквиста, согласно которой конденсация пара и выпадение осадков начинается после достижения относительной влажностью U критического значения U0 и продолжается с ростом U0 до значения Us=1. При этом предполагается, что если среднее по расчетной ячейке значение U удовлетворяет условию U0<U<1, то часть ячейки размером b является облачной, в ней U=1 и возможно наличие водности . Другая часть ячейки размером 1-b является безоблачной, в ней пар не насыщен, причем U=U0. Водности в безоблачной части ячейки нет. Величина b является баллам облачности. Общая водность ячейки


Здесь AT, Aq, Am... &3151; изменения соответствующих метеорологических величин вследствие влияния всех факторов кроме гидрологического. В них входят адвективные, турбулентные и радиационные факторы. Остальные величины имеют следующий смысл: Ec - испарение капель из облачной зоны в безоблачную, Er - испарение дождевых капель, падающих сверху, которое происходит лишь в безоблачной зоне, Q - конденсации пара, Р - осадки, выпадающие лишь из облачной зоны. Параметры L и cp - соответственно удельная теплота испарения и теплоемкость воздуха при постоянном давлении. Предполагается также, что адвекция водности происходит адвектируется лишь в облачную зону.

Далее принимается, что интенсивность испарения водности и падающих сверху осадков, а также образование осадков параметризованы через значения основных метеоэлементов: Ec - через водность m и дефицит влажности 1 - U, Er - через интенсивность осадков сверху Pr и ту же разность 1 - U, Р - через m. Интенсивность осадков сверху есть



Предполагается также определенным алгоритм расчета в случае, если U<U0, а m≠0 (например, уменьшение водности вследствие испарения, и повышение U вплоть до 1).

Система уравнений (1) - (5) будет полностью замкнутой, если определить интенсивность конденсации Q. Надо ясно понимать, что для замыкания необходимо сделать некоторое физически понятное предположение (гипотезу), справедливость которого можно выяснить только из практики. В работе /7/ формула для Q была получена сложным способом. Помимо этого создается иллюзия, что никакого дополнительного предположения не делалось, а некоторые рассуждения (именно с их помощью замыкается система) остаются физически непонятными.

Формула для Q можно вывести проще, исходя из предположения, что пар, поступающий за счет адвекции в облачную долю b ячейки, идет на поддержание насыщения в этой доле. За счет этого не меняется значение величины b, при этом излишки конденсируются, а нехватка восполняется за счет испарения (конденсация с обратным знаком) части водности.

Записать это можно следующим образом:



Теперь система (1) - (5), (8) является замкнутой. Определение конденсации по формуле (8) физически понятно, хотя справедливость самого предположения надо еще подтвердить численными экспериментами. Эта формула используется при расчете осадков в модели регионального прогноза Гидрометцентра России.


Выходная продукция
Региональной модели являются прогностические поля следующих метеорологических величин:
–  приземного давления и геопотенциала изобарических поверхностей 925, 850, 700, 500, 400, 300, 250, 200, 100 гПа,
–  температуры воздуха у поверхности Земли и на изобарических поверхностях 925, 850, 700, 500, 400, 300, 250 и 200 гПа,
–  влажности воздуха на изобарических поверхностях 925, 850, 700, 500 и 400 гПа,
–  скорости ветра на уровнях изобарических поверхностей 925, 500, 300, 250 и 200 гПа,
–  обложных осадков.
Заблаговременность прогнозов 12, 24, 36 и 48 ч.


Иллюстрации прогнозов (рисунки)


Результаты апробации и сведения о внедрении региональной модели


Вариант региональной 11- уровенной Р-модели атмосферы прогноза полей геопотенциала, ветра, температуры воздуха и метода прогноза осадков с пространственным разрешением 100 км и заблаговременность 48 ч адаптирован для территории Дальнего Востока России в отделе гидрометеорологических исследований и прогнозов (ОГМИП) ДВНИГМИ совместно с автором модели. Результаты испытания данного варианта модели в Дальневосточном УГМС вполне удовлетворительные. На основании этих результатов Центральная методическая комиссия по гидрометеорологическим и гелиогеофизическим прогнозам рекомендовала внедрение данного варианта модели в ДВ РВЦ в качестве основной технологии численного краткосрочного прогноза метеорологических полей и использование прогностической продукции в оперативной практике Дальневосточного УГМС.

Автор модели: Владимир Маркович Лосев, Гидрометцентр России, лаборатория региональных краткосрочных прогнозов погоды
E-mail: losev@rhmc.mecom.ru

Ответственный исполнитель работ: ОГМИП ДВНИГМИ в г.Хабаровске,
Евгения Митрофановна Вербицкая
E-mail: werba@mail.kht.ru


© Методический кабинет Гидрометцентра России