Comparison of mesoscale and large-eddy simulation results with observational data in the atmospheric boundary layer

全文:

1. Введение

Интенсивное развитие численных моделей прогноза погоды и климата, обусловленное экспоненциальным ростом производительности суперкомпьютеров и расширением сети наблюдений, происходит в двух направлениях. В одном из подходов основные усилия сосредоточены на развитии моделей глобальной циркуляции, решающих упрощённые уравнения геофизической гидродинамики, и направлены на более детальное описание ключевых физических процессов, измельчение расчётных сеток до субкилометрового масштаба и развитие методов усвоения данных. Альтернативный подход сосредоточен на развитии микрометеорологических моделей, явно разрешающих длинноволновую часть и часть инерционного интервала спектра турбулентных вихрей в атмосферном пограничном слое (АПС), направлен на увеличение расчётной области и подключение дополнительных модулей, описывающих физические процессы, при практически постоянном пространственном разрешении [Schalkwijk et al., 2015].

В мезомасштабных моделях, которые занимают промежуточное положение между глобальными и микрометеорологическими моделями, обеспечивают ограниченное покрытие территории и активно применяются для краткосрочного прогноза погоды и регионального моделирования климата, используются уравнения гидродинамики для сжимаемой атмосферы во вращающейся системе координат. Для описания АПС используются рассматривающие только вертикальный обмен одномерные схемы, которые в условиях неоднородной земной поверхности перестают быть применимыми [Goger et al., 2019; Liu et al., 2020]. Кроме этого, относительно грубое горизонтальное разрешение (1–10 км) не позволяет должным образом учитывать в мезомасштабных моделях влияние микрометеорологической турбулентности внутри АПС на его динамику, обмен между атмосферой и поверхностью, потоки и распределение тепла, водяного пара и примесей. Уменьшение размера ячеек мезомасштабных моделей, применяемых в условиях дневной конвекции, приводит к тому, что модель начинает разрешать конвективные движения в АПС, которые параметризуются, тем самым дублируя описание одного и того же процесса [Honnert et al., 2011]. Это влечёт риски получения вычислительных артефактов в виде структур роликовой циркуляции с горизонтальной осью и осцилляций в распределении турбулентной кинетической энергии, плотностей потоков импульса и тепла [Ching et al., 2014; Haupt et al., 2019; Rai et al., 2019; Kadash et al., 2021].

Аккуратное воспроизведение турбулентной динамики АПС должно охватывать диапазон масштабов от синоптической квазидвумерной турбулентности до вихрей, характерные размеры которых в несколько раз меньше высоты АПС. Поэтому тенденция к увеличению пространственного разрешения в многомасштабном моделировании атмосферных течений реализуется за счёт применения последовательности «телескопически» вложенных доменов с размером ячеек от 9–12 км для внешней области с горизонтальными размерами ~ 10³ км до 10–10² м для внутренней области с размерами 1–4 км [Moeng et al., 2007; Gibbs et al., 2011; Talbot et al., 2012; Muños–Esparza et al., 2014; Rai et al., 2017; Udina et al. 2020; Kim et al., 2023]. В более простом и часто применяемом подходе одностороннего сопряжения решение во внешнем более крупном домене обеспечивает входными данными решение в дочернем вложенном домене, которое не оказывает обратного влияния на родительский домен. В переходном между разрешением мезомасштабных и микрометеорологических моделей диапазоне 10²–10³ м выделяют, так называемую, «серую зону», представляющую область, в которой размер ячейки модели сопоставим с характерными размерами вихрей энергонесущего интервала турбулентности в АПС [Wingaard, 2004]. Для расширения применимости мезомасштабных моделей на эту область требуются более общие трёхмерные схемы подсеточного переноса турбулентной кинетической энергии (ТКЭ) [Zhang et al., 2018].

Вихреразрешающее моделирование LES (large-eddy simulations) широко используется в качестве инструмента исследования турбулентных течений в атмосфере с числами Рейнольдса Re ~ 10⁸–10⁹, которые соответствуют горизонтальным масштабам в несколько десятков километров [Stoll et al., 2020]. Одной из важных задач, решаемых с помощью LES, является разработка параметризаций для мезомасштабных региональных моделей [Ayotte et al., 1996]. Основное упрощение, используемое в LES для того, чтобы сделать возможным численное решение уравнений динамики атмосферных течений, переноса включённых в них примесей и других связанных процессов, состоит в эффективном уменьшении Re посредством фильтрации мелкомасштабных флуктуаций компонент скорости в исходных уравнениях гидродинамики и использовании анэластического приближения или приближения Буссинеска с учетом ускорения Кориолиса [Spiegel and Veronis 1960; Mihaljan, 1962; Bannon, 1996; Lilly, 1996]. При достижении нижней границы горизонтального разрешения мезомасштабной модели для дальнейшей детализации возникает необходимость подключать к моделированию LES, что даёт возможность явно разрешать трёхмерные вихри в АПС, учитывать формы рельефа, горизонтальную неоднородность напряжений Рейнольдса и потоков у земной поверхности.

Последовательное уменьшение масштабов, соответствующее увеличению пространственно-временного разрешения, сталкивается с несколькими проблемами, одна из которых связана с необходимостью возбуждения мелкомасштабных компонент течения в процессе переноса крупномасштабных компонент во вложенном домене на сетку с меньшим размером ячеек. Существует несколько методов инициализации турбулентных флуктуаций, отсутствующих в решении на более крупной сетке, вычислительная эффективность которых зависит от возмущаемой величины (температуры или компонент скорости), стратификации АПС и соотношения размеров ячеек внешнего и внутреннего доменов [Lund et al., 1998; Mayor et al., 2002; Moeng et al., 2007; Muños–Esparza et al. 2014; Muños–Esparza and Kosović, 2018; Haupt et al., 2019]. Вторая проблема связана с адаптацией решения мезомасштабной модели для сжимаемой атмосферы с отличной от нуля дивергенцией поля скорости и решения в LES-домене, в котором, как правило, ищется бездивергентное поле скорости, соответствующее условию несжимаемости, если используется приближение Буссинеска или слабой стационарной сжимаемости в анэластическом приближении [Moeng et al., 2007; Nakayama et al., 2012; Wyszogrodzki et al., 2012]. Аргументация отказа от полноценного учёта сжимаемости в LES основана на малости вертикальной протяженности LES-домена по сравнению с характерным масштабом изменения плотности с высотой, малости чисел Маха (M < 0.1) и малого влияния акустических волн на турбулентную динамику метеорологических переменных. Так, отношение кинетических энергий потенциальной и соленоидальной компонент поля скорости является малым параметром даже при сверхзвуковых течениях (M < 10), и при M < 0.2 может быть оценено как ⟨u_p²⟩/⟨u_s²⟩ < 10^–2 [Dobler et al., 2003]. В анэластическом приближении для атмосферы divu ~ 10^–4 с^–1, где u – вектор скорости течения [Ogura and Philips, 1962], что соответствует величине связанных с эффектами сжимаемости флуктуаций скорости ветра ~ 1 м/с на масштабе ~ 10 км, что обычно не принимается в расчёт в микрометеорологических моделях, если специально не исследуются процессы распространения акустических волн. На оптимизацию соотношения размеров ячеек внешнего домена мезомасштабной модели и вложенного LES-домена также следует обратить внимание, поскольку это соотношение, с одной стороны, не должно быть слишком большим, так как в этом случае усугубляется проблема развития мелкомасштабной турбулентности в достаточно гладком мезомасштабном течении. С другой стороны, нельзя допускать, чтобы размер ячеек в LES-домене оказался в «серой зоне», так как при столь грубом разрешении, соизмеримом с высотой АПС, существенные детали турбулентного течения окажутся потерянными [Wingaard, 2004].

Одна из наиболее широко используемых в мире численная модель прогноза погоды и система моделирования атмосферных процессов WRF (Weather Research and Forecasting) предлагает спектр возможностей для широкого круга исследовательских задач и оперативных приложений [Powers et al., 2017; Skamarock et al., 2019]. Цель настоящей работы состоит в оценке прогностических возможностей комбинированного многомасштабного моделирования дневного АПС над неоднородной поверхностью с использованием сопряжения WRF c моделью PALM [Maronga et al., 2020], представляющей разрабатываемую более 20 лет систему для вихреразрешающего моделирования пограничных слоёв атмосферы и океана. Мотивация сопряжения мезомасштабной и вихреразрешающей моделей оправдана более высокой специализацией каждой из указанных моделей в своём диапазоне масштабов, поэтому реализуемый подход предоставляет больший выбор опций для моделирования внутренней области и подстилающей поверхности, чем использование доступной в WRF опции вихреразрешающего моделирования. Использование различных инструментов акустического и микроволнового зондирования атмосферы, а также средств локальных наблюдений метеорологических параметров на Геофизической обсерватории «Борок» ИФЗ РАН позволяет сравнить результаты моделирования с таким набором данных наблюдений, который не рассматривался в предыдущих аналогичных работах.

2. Инструменты и методы измерений

Наблюдения метеорологических переменных, параметров турбулентности, плотности потока коротковолновой солнечной радиации и облачности проводились с помощью измерительного комплекса Геофизической обсерватории «Борок» ИФЗ РАН (58° 04ʹ с.ш., 38° 14ʹ в.д.). Измерительный полигон обсерватории расположен вдали от возможных источников антропогенных возмущений естественного состояния природной среды и представляет собой ровное поле площадью около 40000 м², окружённое редким лесом и болотистой местностью. К востоку от измерительного полигона на расстоянии 1.5 км расположено Рыбинское водохранилище, побережье которого представляет собой заболоченную местность и редкий лиственный лес. С юга на расстоянии более 500 м расположены низкоэтажные строения. На север и запад от измерительного полигона находятся поля с травянистой растительностью протяженностью более 2 км.

Наземные измерения набора метеорологических и турбулентных параметров атмосферы проводились с помощью двух ультразвуковых метеостанций «МЕТЕО-2Н», установленных в поле на высоте 2 м от поверхности, и ультразвуковой метеостанции «АМК-03П/1», установленной на крыше регистрационного корпуса на высоте 10 м. Метеостанции обладают высоким быстродействием, выполняя измерения температуры и трёх компонент скорости ветра с частотой 10 Гц, а также относительной влажности и статистических параметров турбулентности, включая вертикальный турбулентный поток явного тепла, один раз в минуту, передавая результаты в цифровом виде в систему сбора информации.

Измерения температуры воздуха до высоты 1 км проводились с помощью температурного профилемера MTP-5, представляющего собой прибор дистанционного зондирования атмосферы. Принцип работы прибора основан на регистрации микроволнового излучения кислорода на частоте 56.6 ГГц. Периодичность измерения 5 мин. Трёхкомпонентный акустический доплеровский локатор (SODAR – Sound Detecting And Ranging ) «Волна-3» использовался для измерения скорости и направления ветра на высотах от 50 м до 800 м с разрешением по высоте 5 м. Рабочая частота локатора «Волна-3» – 1700 Гц, период цикла зондирования – 16 с, диапазон измерения средней скорости ветра – от 0.3 м/с до 20 м/с, приборная погрешность – 0.3 м/c. Температурный профилемер MTP-5 и СОДАР «Волна-3» полностью автоматизированы и интегрированы в систему сбора информации ГО «Борок» ИФЗ РАН.

Непосредственные высотные измерения проводились с помощью датчиков, расположенных на аппаратной платформе [Анисимов и др., 2023; Anisimov et al., 2021a], которая поднималась заполненным гелием аэродинамическим привязным аэростатом Au-6 (ООО «Авгур – аэростатные системы», Россия). Аэростатное оборудование позволяет длительное время удерживать полезную нагрузку на требуемой высоте над поверхностью земли. Высота удержания полезной нагрузки мало меняется в рабочем диапазоне скоростей ветра аэростата. В составе аэростатной аппаратной платформы присутствует метеорологический модуль, регистрирующий температуру и относительную влажность воздуха; барометрический альтиметр, регистрирующий высоту подъёма платформы с точностью до 0.1 м; акселерометр-гироскоп и цифровой микроэлектромеханический компас, предназначенные для регистрации телеметрической информации. Для определения температуры и относительной влажности воздуха на аэростатной платформе применяется датчик DHT11, установленный снаружи аэростатной платформы в защитном вентилируемом экране. Датчик DHT11 позволяет измерять температуру окружающего воздуха с точностью ±0.2 °C в диапазоне от 0 до +50 °C. Абсолютная погрешность измерения относительной влажности не превышает ±5% в диапазоне от 20% до 95%. Временное разрешение датчика DHT11 – одно измерение в секунду. С целью оценки корректности работы датчика производились синхронные наблюдения температуры и относительной влажности воздуха с помощью датчика DHT11 и ультразвуковой метеостанции «AMK-18424». В ходе наблюдений выявлено, что отклонения показаний датчика DHT11 от показаний метеостанции «AMK-18424» не превышают паспортных значений. Пространственное положение платформы регистрируется с помощью GPS/ГЛОНАСС оборудования. Автоматизированная автономная система сбора данных в составе аэростатной аппаратной платформы, выполнена на основе одноплатного микроЭВМ RaspberryPI и работает под управлением Debian GNU/Linux. Системное время микроЭВМ также синхронизируется с помощью GPS/ГЛОНАСС оборудования, и весь сбор данных осуществляется с учетом этого времени. Аэростатные высотные измерения проводились в хорошую погоду при следующих условиях: скорость горизонтального ветра у поверхности земли не более 5 м/с (на высоте свыше 100 м – не более 10 м/с); облачность – менее 10%; отсутствие осадков или иных ухудшающих видимость явлений. Плотность потока приходящего коротковолнового солнечного излучения измерялась с помощью пиранометра CMP3 Kipp & Zonen.

3. Численное моделирование

Ранее были разработаны интерфейсы вложения PALM (динамический драйвер) для сопряжения с региональной моделью COSMO (Consortium for Small-scale Modelling) [Kadash et al., 2021], которая не является моделью с открытым исходным кодом, и с более распространённой мезомасштабной моделью WRF, имеющей открытый исходный код [Lin et al., 2021]. Поэтому многомасштабное моделирование динамики АПС выполнялось с помощью численных моделей WRF-ARW v4.5 [Skamarock et al., 2019] и PALM v. 6.0 [Maronga et al., 2020] с четырёхкратным вложением расчётных областей (рис. 1) и однонаправленным нисходящим сопряжением мезомасштабной и вихреразрешающей моделей. Для мезомасштабного моделирования WRF выбраны три области (рис. 1). Во внешнюю область D₁ с размерами 1800×1800 км², покрывающую европейскую часть России и имеющую шаг сетки 9 км вложена промежуточная область D₂ с размерами 360×360 км² и шагом сетки 3 км. Вложенная в D₂ внутренняя область D₃ с размерами 90×90 км² имеет шаг сетки 1 км. Все перечисленные сетки имеют совпадающие 50 уровней в вертикальном направлении до давления в 5 кПа. Центры расчётных областей совмещены и соответствуют расположению измерительного комплекса ГО «Борок» ИФЗ РАН (58.070531 с.ш., 38.232595 в.д.). Такая конфигурация детально воспроизводит дневную конвекцию во всей расчётной области с шагом внутренней сетки достаточным для анализа мезомасштабной динамики АПС. Во избежание попадания в «серую зону» с последствиями в виде получения вычислительных артефактов разрешение мезомасштабной модели не следует увеличивать, поскольку оно соответствует характерной высоте АПС [Honnert et al., 2011; Haupt et al., 2019]. Дальнейшее увеличение разрешения достигается за счёт вложения в область D₃ мезомасштабной модели двух областей D₄ и D₅, в которых выполняется вложение одной области LES в другую, и расчёты проводятся одновременно в обеих областях [Hellsten et al., 2021]. В табл. 1 даны параметры расчётных областей (доменов), пространственное разрешение и диапазон изменения шага численного интегрирования разностных уравнений по времени.

 

Рис. 1. Схема вложения расчётных областей D₁–D₅

 

Таблица 1. Размеры доменов, параметры сетки и величина временного шага WRF и LES

Модель	Домен	Количество ячеек в горизонтальной плоскости Nx × Ny	Количество вертикальных уровней Nz	Длина стороны ячеек в горизонтальной плоскости, км	Высота домена, км	Временной шаг ∆t, c
WRF	D1	200×200	50	9	19.892	45
WRF	D2	120×120	50	3	19.892	15
WRF	D3	90×90	50	1	19.892	5
LES	D4	192×192	64	0.09	3.072	0.284 – 1
LES	D5	192×192	192	0.03	2.304	0.284 – 1

 

В качестве основного набора параметризаций WRF используются рекомендуемые для мелких сеток параметры с уменьшенным временем вызова модулей параметризации радиационного баланса RRTMG (RADT = 10 мин) (табл. 2). Время вызова блока конвекции CUDT = 10 мин. В качестве начальных и граничных условий используются данные глобальной прогнозной модели GFS с разрешением 0.25°: NCEP GFS 0.25 Degree Global Forecast Grids Historical Archive (ds084.1, https://doi.org/10.5065/D65D8PWK). Параметриза-ция поверхности выполнена с использованием Noah Land Surface Model: Unified NCEP/NCAR/AFWA (https://ral.ucar.edu/model/unified-noah-lsm) [He et al., 2023]. Начало расчётов в 18.00 UTC предыдущих суток. Результаты первых 12 ч расчётного периода для анализа не использовались.

 

Таблица 2. Схемы параметризации процессов WRF

Процесс	Параметр WRF	Название параметризации	Ссылка
Облачная микрофизика	MP_PHYSICS = 8	Схема Томпсона	[Thompson et al., 2008]
Длинноволновая радиация	RA_LW_PHYSICS = 4	Схема RRTMG	[Iacono et al., 2008]
Коротковолновая радиация	RA_SW_PHYSICS = 4	Схема RRTMG	[Iacono et al., 2008]
Приземный слой	SF_SFCLAY_PHYSICS = 2	Схема Монина–Обухова с оптимизацией Зилитинкевича	[Janjić, 1994]
Поверхность	SF_SURFACE_PHYSICS = 2	Модель Noah	https://ral.ucar.edu/document-or-file/noah-lsm-users-guide [He et al., 2023]
Планетарный пограничный слой	BL_PBL_PHYSICS = 2	Схема Меллора– Ямады–Яника	[Janjić, 1994]
Конвекция	CU_PHYSICS = 1	Схема Кейна–Фрича	[Kain, 2004]

 

Вихреразрешающая численная модель (LES) турбулентной динамики атмосферного пограничного слоя PALM v. 6.0, с помощью которой в анэластическом приближении находится разностное решение на декартовой сетке фильтрованных уравнений Навье–Стокса, уравнений переноса тепла и водяного пара, используется для получения ансамбля реализаций 3-мерных полей скорости U, потенциальной температуры Θ и удельной влажности Q в доменах D₄ и вложенного в него D₅. Применение процедуры пространственной фильтрации динамических и транспортных уравнений приводит к появлению в них слагаемых, содержащих ковариации компонент турбулентных флуктуаций скорости $R_{ij}=\overline{u_i{u}_j}$ – тензор напряжений Рейнольдса, ковариации скорости и потенциальной температуры $H_i^{SGS}=\overline{u_i\theta }$, а также ковариации скорости и удельной влажности $E_i^{SGS}=\overline{u_{i}q}$, представляющие собой подсеточные турбулентные потоки тепла и водяного пара соответственно. Для моделирования подсеточных переменных используется гипотеза о пропорциональности турбулентных потоков локальным градиентам разрешаемых на сетке компонент скорости U_i, потенциальной температуры Θ и удельной влажности Q, а именно

$r_{ij}=\overline{u_i{u}_j}-\frac{2}{3}e{\delta }_{ij}=-K_m\left(\frac{\partial U_i}{\partial x_j}+\frac{\partial U_j}{\partial x_i}\right)$, (1)

$\overline{u_i\theta }=-K_h\frac{\partial \Theta }{\partial x_i},\overline{u_{i}q}=-K_h\frac{\partial Q}{\partial x_i}$, (2)

где r_ij – девиатор R_ij, e = Tr (R_ij) – подсеточная часть ТКЭ, δ_ij – единичный тензор, K_m и K_h – локальные коэффициенты подсеточной турбулентной диффузии импульса и тепла соответственно, определяемые в случае неустойчивой и нейтральной стратификации АПС следующим образом

$K_m=c_{\ast }{\Delta }_m\sqrt{e},K_h=3K_m,{\Delta }_m=\max \left(\Delta x,\Delta y,\Delta z\right)$, (3)

где вычисление c_* основано на использовании тестового фильтра шириной ∆_T = 2∆_m [Germano et al., 1991]. Подсеточная часть ТКЭ находится из прогностического уравнения

$\frac{\partial e}{\partial t}=-U_j\frac{\partial e}{\partial x_j}-R_{ij}\frac{\partial U_i}{\partial x_j}+$

$+\frac{g}{{\Theta }_v}\overline{u_3\theta }-\frac{\partial }{\partial x_j}\left[\overline{u_j\left(e+\frac{p}{{\rho }_0}\right)}\right]-\epsilon$, (4)

в котором Θ_v – потенциальная виртуальная температура, слагаемое с ковариацией флуктуаций скорости и давления параметризуется аналогично (1) и (2)

$\left[\overline{u_j\left(e+\frac{p}{{\rho }_0}\right)}\right]=-2K_m\frac{\partial e}{\partial x_j}$, (5)

а подсеточная скорость диссипации ТКЭ ε определена выражением

$\epsilon =c_{\epsilon }\frac{e^{3/2}}{\Delta }$, (6)

где c_ε = 0.93, ∆ – среднее геометрическое трёх линейных размеров ячейки сетки.

Сопряжение WRF и PALM выполнялось с помощью мезомасштабного динамического драйвера WRF4PALM v 1.0 [Lin et al., 2021] офлайн методом, то есть моделирование в доменах D₄ и D₅ проводилось после завершения моделирования в доменах D₁ – D₃. Офлайн вложение с динамическим драйвером подразумевает преобразование и интерполяцию географических и сигма-координат WRF к декартовой прямоугольной сетке PALM, интерполяцию на сетку PALM мезомасштабной вертикальной компоненты скорости ветра, постановку граничных условий Дирихле на боковых и верхней границах вложенного домена. Граничное условие Дирихле применяется к трём компонентам скорости, потенциальной температуре и удельной влажности, значения которых из WRF на границе линейно интерполируются по времени на каждом временном шаге для форсинга PALM. Временной ряд значений атмосферного давления у поверхности из WRF используется на самом низком уровне сетки PALM после интерполяции для вычисления вертикальных профилей давления и плотности воздуха. Референсная потенциальная температура, которая участвует в уравнении динамики, рассчитывается, так же как и давление, с помощью горизонтального усреднения по домену.

Незначительное несоответствие между потоками массы в LES-домене между решениями WRF и LES, связанное с использованием в последнем анэластического приближения, устраняется коррекцией вертикальной скорости в каждом узле на верхней границе z_top и на предыдущем слое на величину входящего потока массы воздуха через боковые границы, делённую на площадь верхней границы домена [Hellsten et al., 2021; Kadash et al., 2021]. Высота АПС z_i определяется по положению минимума плотности полного вертикального турбулентного потока тепла. Выше меньшей из высот 1.3z_i и 0.8z_top применяется рэлеевское демпфирование. Значения z_top в доменах D₄ и D₅ установлены 3070 м и 2300 м соответственно.

В качестве нижнего граничного условия между земной поверхностью и первым уровнем модели в соответствии с теорией подобия Монина–Обухова предполагается слой постоянных по высоте потоков импульса, тепла и удельной влажности, определённых локально в условиях неоднородной поверхности. Турбулентные потоки тепла и удельной влажности с поверхности земли определяются на основе встроенной в PALM модели земной поверхности в соответствии с изменяющимся потоком приходящего коротковолнового солнечного излучения, высотой и положениями неоднородностей земной поверхности, включая водные объекты, а также физическими свойствами почвы и растительности [Heldens et al., 2020; Gehrke et al., 2021]. Для определения свойств земной поверхности используются растровые данные из базы данных ESA WorldCover (https://esa-worldcover.org/en). Подготовленные с помощью статического драйвера PALM [Lin et al., 2021] данные ESA WorldCover содержат информацию (рис. 2а), которая используется в модели переноса излучения, моделях растительного покрова и деятельного слоя почвы.

 

Рис. 2. а) Данные ESA WorldCover (https://esa-worldcover.org/en), подготовленные с помощью статического драйвера PALM, б) схема расположения ячеек, в которых применяется метод возбуждения мелкомасштабной турбулентности

 

Топографическая карта земной поверхности в области D₅ задается с использованием интерполированных значений цифровой карты высот с разрешением 30×30 м² (Digital Elevation Model), полученных в ходе спутникового зондирования Shuttle Radar Topography Mission (https://lpdaac.usgs.gov/products/srtmgl1v003, DOI: 10.5067/MEaSUREs/SRTM/SRTMGL1.003, https://srtm.csi.cgiar.org). Сила сопротивления, действующая на течение со стороны возвышений рельефа поверхности, учитывается с помощью дополнительного слагаемого в правой части уравнения динамики фильтрованных горизонтальных компонент скорости [Anderson and Meneveau, 2010]

$f_i\left(x,y,{\Delta }_z/2\right)=-\frac{U_i\left(x,y,{\Delta }_z/2\right)}{{\Delta }_z}{U}_k\frac{\partial h}{\partial x_k}H\left(U_k\frac{\partial h}{\partial x_k}\right),$

$i,k=\mathrm{1,}\mathrm{2,}{x}_1=x,x_2=y$ (7)

где ∆_z – расстояние между вертикальными уровнями, h – высота рельефа, H – ступенчатая функция Хевисайда. В примыкающих изнутри к границам доменов D₄ и D₅ горизонтальных слоёв, состоящих из 6 ячеек нижнего уровня, отсчитываемых вглубь соответствующих доменов (рис. 2б), применяется метод возбуждения турбулентности, основанный на добавлении в уравнение для потенциальной температуры стохастического источника [DeLeon et al., 2019; Senocak and DeLeon, 2023], для параметризации статистики флуктуаций которого использовались временные ряды наблюдений температуры на высоте 2 м от земной поверхности.

4. Результаты и обсуждение

Чтобы исключить влияние облачности на радиационный баланс и турбулентную динамику микрометеорологических переменных в АПС, для моделирования были выбраны два дня 05.08.2022 г. и 17.08.2023 г., в наибольшей степени удовлетворяющие условиям хорошей погоды вследствие отсутствия облаков в выбранный для моделирования временной интервал с 9 до 15 ч (рис. 3). Также в эти даты осуществлялись запуски привязного аэростата с инструментальной платформой для исследования высотного распределения и вариаций атмосферных электрических и сопутствующих аэрофизических величин. Сравнение рассчитанных метеорологических полей и параметров турбулентности в дневном АПС проводилось с результатами наземных и аэростатных наблюдений, акустического и микроволнового зондирования атмосферы. Для сравнения отобраны результаты моделирования WRF и PALM в центральной области домена D₅, которая соответствует местам размещения средств наблюдений и охватываемых ими областей дистанционного зондирования. Отклонения модельного прогноза от результатов наблюдений могут быть оценены с помощью среднеквадратичного отклонения RMSE (root mean square error) и индекса согласованности IOA (index of agreement) [Lin et al., 2021]

$\text{RMSE }=\sqrt{\frac{\text{1}}{N}\sum _{n=1}^N{\left(P_n-O_n\right)}^2}$, (8)

$\text{IOA }=1-\frac{{\displaystyle \sum _{n=1}^N\left|P_n-O_n\right|}}{{\displaystyle \sum _{n=1}^N\left(\left|P_n-\overline{O}\right|+\left|O_n-\overline{O}\right|\right)}}$, (9)

где P_n – рассчитанные в модели значения прогнозируемой величины, O_n – соответствующие значения этой величины, полученные в наблюдениях, со средним Ō. При замене наблюдаемых значений соответствующими данными из другой модели индексы могут быть использованы для оценки отклонений между двумя моделями. Если же рассчитать индексы (8) и (9), используя результаты наблюдений в точках, находящихся на расстоянии L, то можно оценить естественный разброс величин и согласованность их флуктуаций на этом масштабе.

 

Рис. 3. Плотность потока приходящего коротковолнового солнечного излучения

 

На рис. 4 показаны рассчитанные WRF распределения геопотенциальной высоты на уровне 500 гПа и атмосферного давления на уровне земной поверхности в 12 ч по местному времени. Для оценки корректности воспроизведения крупномасштабных метеорологических процессов использовались данные зондирования атмосферы Университета Вайоминга (https://weather.uwyo.upperair) на станциях 27713 – Москва (Долгопрудный), 27038 – Вологда и 27459 – Нижний Новгород в 12 UTC (15 LT). На рис. 5 показаны рассчитанные WRF вертикальные профили потенциальной температуры и скорости горизонтального ветра в сравнении с результатами наблюдений 17.08.2023. Для 05.08.2022 получаются аналогичные результаты: вертикальные профили температуры и скорости ветра, а также положение верхней границы АПС воспроизводятся моделью WRF корректно, показывая соответствие средних градиентов температуры в тропосфере, а также расположение наиболее выраженных минимумов и максимумов в профиле скорости ветра. Максимальное отклонение потенциальной скорости ветра в АПС не превышает 4 м/c.

 

Рис. 4. Рассчитанные WRF распределения геопотенциальной высоты на уровне 500 гПа (верхний ряд) и атмосферного давления на уровне земной поверхности (нижний ряд) в 12 ч по местному времени; точка показывает положение Геофизической обсерватории «Борок» ИФЗ РАН (58° 04ʹ с.ш., 38° 14ʹ в.д.)

 

 

Рис. 5. Наблюдаемые и рассчитанные в модели WRF D₁ высотные профили потенциальной температуры и модуля горизонтальной скорости ветра в 12 UTC 17.08.2023; горизонтальные линии на графиках температуры показывают высоту АПС

 

Вариации температуры, удельной влажности, модуля и направления горизонтальной скорости ветра у земной поверхности с 10-минутным усреднением представлены на рис. 6. В условиях неоднородной поверхности результаты LES зависят от расположения и количества ячеек, по которым проводится усреднение, поэтому на рис. 6 показан диапазон изменения T, q и V_h по всем ячейкам области D₅, в квадрате площадью 1 км² в центре домена, а направление горизонтальной скорости ветра, рассчитанное в LES, усреднено по этой области. Видно, что диапазон рассчитанных в LES величин может не содержать значений, предсказываемых моделью WRF для той же пространственной области, однако отклонения несущественны.

Совокупные количественные оценки качества прогноза моделей WRF и LES температуры и модуля скорости горизонтального ветра, основанные на индексах RMSE (8) и IOA (9), даны в табл. 3 и 4. Для вычисления индексов RMSE и IOA использованы 5-минутные данные микроволнового зондирования, 10-минутные данные акустического зондирования и 1-минутные данные ультразвуковых наземных метеостанций. Кроме этого, выполнена оценка вариабельности наблюдаемых значений модуля скорости горизонтального ветра и температуры на основе индексов RMSE и IOA, рассчитанных по 1-минутным данным, которые были получены с двух цифровых метеостанций, установленных на высоте 2 м от земной поверхности и разнесённых на 2 м и 4 м. Результаты приведены в табл. 5.

 

Рис. 6. Результаты моделирования и наземных наблюдений температуры (T), удельной влажности (q), модуля (V_h) и направления горизонтальной скорости ветра; min и max показывают мгновенные экстремальные значения в квадрате площадью 1 км² в центре D₅, направление ветра LES усреднено по указанному квадрату; высота 42 м отсчитывается от нижней границы D₅, совпадающей с уровнем реки и находится выше всех неоднородностей рельефа в D₅ (перепад высот в D₅ составляет 42 м, медианная высота в D₅ от нижней границы – 10 м)

 

Таблица 3. Отклонения прогноза WRF от результатов наблюдений

Дата	RMSE Vh [м/c]	IOA Vh	N	RMSE T [K]	IOA T	N
05.08.2022	0.65	0.52	798	0.81	0.87	828
17.08.2023	1.11	0.55	798	0.64	0.89	828

 

Таблица 4. Отклонения прогноза LES от результатов наблюдений

Дата	RMSE Vh [м/c]	IOA Vh	N	RMSE T [K]	IOA T	N
05.08.2022	2.36	0.19	1440	0.65	0.89	1656
17.08.2023	1.95	0.34	1440	0.92	0.85	1656

 

Таблица 5. Отклонения результатов разнесённых на расстояние L наземных наблюдений

Дата	L [м]	RMSE Vh [м/c]	IOA Vh	RMSE T [K]	IOA T	N
05.08.2022	2	0.23	0.86	0.50	0.95	1440
17.08.2023	4	0.19	0.84	0.44	0.96	1440

 

Сравнение с индексами, полученными в [Lin et al., 2021] с использованием результатов наблюдений и моделирования с усреднением за 1 ч, показывает, что моделирование АПС с помощью WRF и PALM выполнено корректно, поскольку достигается соответствующий результатам предыдущих работ уровень согласованности с наблюдениями. Так, в [Lin et al., 2021] индексы RMSE для температуры находятся в диапазоне 1.13–2.02 WRF, 0.99–3.61 PALM, для скорости ветра 1.55–2.70 WRF, 1.76–2.45 PALM. Индекс IOA для температуры находится в диапазоне 0.72–0.88 WRF, 0.56–0.81 PALM, для скорости ветра 0.50–0.84 WRF и 0.56–0.69 PALM. Следует также учитывать, что при увеличении временного интервала усреднения значение индекса RMSE уменьшается.

На рис. 7 показаны рассчитанные вертикальные профили температуры в сравнении с результатами микроволнового зондирования атмосферы с помощью профилемера MTP-5. Суперадиабатический градиент температуры, не прогнозируемый моделями, систематически наблюдается в нижней части АПС выше приземного слоя в первой половине дня. Такая форма вертикального профиля температуры может быть связана с испарением микрокапель и высыханием гигроскопичных аэрозольных частиц, приводящих к локальному понижению температуры в слое их повышенной концентрации. В этой области температура в модели WRF приблизительно на 1 K ниже, чем в модели PALM. Выше АПС согласованность модельных прогнозов WRF с наблюдаемыми вертикальными профилями температуры оказывается лучше в отличие от примеров, приведенных на рис. 5. На рис. 8 показаны измеренные и рассчитанные изменения температуры на трёх уровнях. Видно, что максимальное отклонение результатов моделирования от данных наблюдений с 5-минутным усреднением не превышает 4 K и имеет место в области суперадиабатического градиента температуры. Температура и удельная влажность, зарегистрированные с помощью датчиков, размещенных на инструментальной платформе привязного аэростата в сравнении с рассчитанными в моделях показаны на рис. 9. Данные с аэростатной платформы 05.08.2022 и 17.08.20232 получены на высотах 100 м и 400 м соответственно. Максимальное отклонение температуры при скользящем усреднении с окном в 30 мин составляет приблизительно 3 К, а удельной влажности – 4 г/кг, при этом в обоих рассматриваемых случаях модели переоценивают концентрацию водяного пара в воздухе.

 

Рис. 7. Модельные и измеренные MTP-5 высотные профили температуры с усреднением 5 мин

 

Рис. 8. Модельные и измеренные MTP-5 вариации температуры на трёх уровнях высоты с усреднением 5 мин

 

Рис. 9. Вариации температуры и удельной влажности на двух уровнях высоты по данным аэростатных наблюдений с усреднением 40 мин, результатам WRF с усреднением 10 мин и LES с усреднением 1 мин

 

На рис. 10 и 11 показаны рассчитанные в моделях WRF и LES вертикальные профили модуля и направления вектора горизонтальной скорости ветра c 10-минутным усреднением и данные акустического зондирования атмосферы с помощью СОДАРа. В дневное время рассматриваемых дат средняя наблюдаемая скорость ветра до высоты 400 м колеблется около значения 1 м/c, в то время как результаты моделирования находятся в диапазоне 2–4 м/c. Максимальное отклонение результатов наблюдений и моделирования модуля горизонтальной скорости ветра не превышает 4 м/c по данным с временным разрешением в 10 мин. В целом модельный прогноз скорости ветра менее аккуратен, чем прогноз температуры (см. также рис. 5), что находит количественное выражение в различающихся значениях RMSE и IOA (табл. 3, 4).

 

Рис. 10. Высотные профили модуля горизонтальной скорости ветра по данным акустического зондирования (SODAR), результатам WRF и LES с усреднением 10 мин.; серые линии показывают мгновенные профили LES; площадь усреднения результатов LES – 1 км²

 

Рис. 11. Высотные профили направления горизонтальной скорости ветра по данным акустического зондирования, результатам WRF и LES с усреднением 10 мин; площадь усреднения результатов LES – 1 км²

 

Статистика нормированных отклонений рассчитанных на основе моделей WRF и PALM значений температуры, удельной влажности и плотности вертикального турбулентного потока тепла от наблюдаемых в дневное время на высоте 2 м над земной поверхностью по данным одной из двух ультразвуковых метеостанций иллюстрируется диаграммами на рис. 12, на которых также показаны отклонения результатов наблюдений в двух точках (расстояния между метеостанциями указаны в табл. 5). Относительные отклонения между результатами моделирования и наблюдаемыми значениями удельной влажности и плотности вертикального турбулентного потока тепла не сильно отличаются от аналогичных величин, рассчитанных по результатам наблюдений в двух точках, расстояние между которыми соответствует нахождению в инерционном интервале турбулентности в АПС и подфильтровым масштабам LES.

 

Рис. 12. Диаграммы статистики отклонений рассчитанных WRF и LES температуры (T), удельной влажности (q) и плотности вертикального турбулентного потока тепла (H) от измеренных с помощью ультразвуковых метеостанций и отклонений между результатами разнесённых измерений; количество отсчётов и расстояния между метеостанциями приведены в табл. 3–5; периоды усреднения T и q – 1 мин, H – 20 мин

 

Вертикальные профили вторых моментов разрешённых LES явно флуктуаций компонент скорости, разрешённой K и подфильтровой e ТКЭ, усреднённых по домену D₅ и по времени с 11:30 до 12:30 LT для 17.08.2023, показаны на рис. 13. Среднее значение u_* при этом составило 0.29 м/c. В табл. 6 приведены средние за период с 11:30 до 12:30 LT измеренные значения моментов второго порядка турбулентных флуктуаций компонент скорости ветра, ТКЭ и динамическая скорость. На рис. 14 показаны вариации ТКЭ, измеренные ультразвуковой метеостанцией на высоте 2 м от земной поверхности, а также разрешённая и подфильтровая компоненты ТКЭ, рассчитанные с помощью LES на высоте 18 м в одной ячейке сетки D₅, соответствующей месту проведения наблюдений. Средние за период с 11:30 до 12:30 LT значения полной ТКЭ в этой ячейке для 5.08.2022 и 17.08.2023 составили 0.22 м²/c² и 0.28 м²/c² соответственно.

 

Рис. 13. Высотные профили вторых моментов турбулентных флуктуаций компонент скорости ветра, разрешённой (K) и подфильтровой (e) ТКЭ в LES с усреднением по площади D₅ и по времени с 11:30 до 12:30 LT 17.08.2023 

 

Рис. 14. Вариации модельных разрешённой (K), подфильтровой (e) компонент и измеренной полной ТКЭ с усреднением 1 мин

 

Таблица 6. Турбулентные параметры по результатам усреднения наземных наблюдений

Дата	ТКЭ, м2/c2	σu2, м2/c2	σw2, м2/c2	⟨uʹwʹ⟩, м2/c2	u*, м/c
05.08.2022	0.26	0.22	0.08	0.01	0.21
17.08.2023	0.41	0.38	0.09	0.00	0.21

 

С целью сравнения выбранной схемы моделирования дневного АПС, предполагающей сопряжение WRF с вложенным LES-доменом, и схемы с ещё одним вложением WRF-домена с увеличенным горизонтальным разрешением, было выполнено моделирование WRF в четырёх доменах D₁ – D₃ и Dʹ₄, размеры которого 18×18 км², а горизонтальный размер ячеек – 1/3 км. На рис. 15 показаны мгновенные горизонтальные сечения температуры и модуля горизонтальной скорости ветра в D₄ в 11 ч местного времени на основе моделей WRF на сетках D₃ и Dʹ₄ и LES на сетке D₄. Видно, насколько более подробно LES разрешает в полях температуры и скорости характерные для конвективного АПС структуры, позволяя на масштабах ~ 1 км, соответствующих высоте слоя конвективного перемешивания, оценивать средние и экстремальные значения.

 

Рис. 15. Горизонтальные распределения температуры и модуля горизонтальной скорости ветра WRF на сетках D₃ и Dʹ₄ на высоте 160 м и LES на сетке D₄ на высоте 168 м в 11 ч 17.08.2023

 

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

$⟨{\left|f\left(r\right)-f\left(r+\Delta r\right)\right|}^p⟩=⟨{\left|\Delta f\right|}^p⟩$, (10)

где r – модуль радиус-вектора, и усреднение проводится по ансамблю реализаций поля. Если в некотором интервале приращений аргумента ∆r выполняется

$⟨{\left|\Delta f\right|}^p⟩\propto \Delta r^{{\zeta }_f\left(p\right)}$, (11)

то скейлинговую экспоненту ζ_f(p), в общем случае зависящую от порядка структурной функции и масштаба ∆r, можно аппроксимировать, вычисляя логарифмическую производную

${\zeta }_f\left(p\right)=\frac{d\ln ⟨{\left|\Delta f\right|}^p⟩}{d\ln \Delta r}$. (12)

Показатели скейлинга структурной функции второго порядка ζ(2), структурной функции абсолютных приращений первого порядка ζ(1) и показателя β степенной аппроксимации спектральной плотности флуктуаций поля f

$S_f\left(k\right)\propto k^{\beta }$, (13)

где k = 2π /∆r – модуль волнового вектора, связаны соотношением [Lovejoy et al., 2004]

$-\beta =1+\zeta \left(2\right)\approx 1+2\zeta \left(1\right)$. (14)

Вычисление структурных функций (10) первого порядка проводилось для рассчитанных одновременных распределений температуры и модуля горизонтальной скорости ветра внутри каждого из доменов D₁ – D₃, Dʹ₄, D₄, D₅. По данным наблюдений на высоте 2 м с частотой дискретизации 1 Гц были рассчитаны с помощью быстрого преобразования Фурье спектры флуктуаций температуры и модуля горизонтальной скорости ветра в интервале с 09:00 ч до 15:00 ч местного времени. Составные вариограммы, иллюстрирующие скейлинг пространственных разностей, а также частотные спектральные плотности флуктуаций температуры и модуля горизонтальной скорости ветра, показаны на рис. 16. Во всех доменах существуют интервалы самоподобия полей температуры и модуля горизонтальной скорости, при этом скейлиноговые экспоненты ζ(1) структурных функций абсолютных разностей модуля горизонтальной скорости в доменах WRF на масштабах 3–100 км близки к 3/5, что соответствует скейлингу Больджиано–Обухова [Benzi et al., 1998; Lovejoy et al., 2004; Lilley et al., 2008; Глазунов и Дымников, 2013]. Однако поведение структурных функций температуры в модели WRF отличается от скейлинга Больджиано–Обухова ζ(1) = 1/5, рассчитанные структурные функции имеют приблизительно тот же показатель, что и структурные функции модуля горизонтальной скорости. В низкочастотной области спектральной плотности измеренных флуктуаций модуля горизонтальной скорости показатель степенной аппроксимации близок к –11/5, что также соответствует скейлингу Больджиано–Обухова. Однако, как и в случае с показателями структурных функций результатов модели WRF, низкочастотный участок спектра наблюдаемых флуктуаций температуры лучше аппроксимируется степенной функцией с показателем –11/5, чем с соответствующим скейлингу Больджиано–Обухова показателем β = –7/5. Диапазон логарифмических производных структурных функций в области масштабов 300 м – 15 км составил 0.2–0.73 для температуры и 0.17–0.41 для модуля горизонтальной скорости. В интервале масштабов (∆x∆y∆z)^1/3 < ∆r < z_i для обеих переменных существует участок с колмогоровской экспонентой ζ(1) = 1/3, где ∆x, ∆y, ∆z – размеры ячеек LES вдоль соответствующих координатных осей, а z_i – высота слоя конвективного перемешивания [ср. Lovejoy et al., 2004]. Высокочастотная часть спектра с f > 10^–2 Гц измеренных флуктуаций температуры и модуля горизонтальной скорости ветра имеет показатель β = –5/3. В крупномасштабной области (∆r > 5 км), рассчитанного в LES поля скорости ветра заметна тенденция к воспроизведению скейлинга WRF ζ(1) = 3/5. Таким образом, комбинированное моделирование с использованием последовательности вложенных доменов позволяет воспроизводить масштабно-инвариантные свойства атмосферной турбулентности в диапазоне масштабов, охватывающем несколько порядков.

 

Рис. 16. Структурные функции модельных абсолютных разностей (слева), спектральная плотность измеренных флуктуаций температуры и модуля горизонтальной скорости ветра (справа)

 

Полученные недавно с использованием моделей WRF результаты о связи климатических мод с процессами в глобальной электрической цепи [Slyunyaev et al., 2021; Kozlov et al., 2022] и LES о влиянии турбулентности на электродинамические процессы АПС [Anisimov et al., 2020, 2021b] мотивируют использовать методы многомасштабного моделирования динамики АПС для решения задач атмосферного электричества. Так, развитые методы сопряжения мезомасштабной модели WRF и модели PALM могут быть применены для исследования влияния метеорологических процессов на действие электрического конвективного генератора в невозмущённых областях атмосферы и формирование электрической стратификации нижней атмосферы.

5. Выводы

Численное моделирование метеорологических полей и турбулентной динамики дневного АПС выполнено с помощью последовательного вложения двух областей вихреразрешающей модели PALM в три области расчета мезомасштабной модели WRF с использованием нисходящего масштабирования сеток с увеличением горизонтального разрешения от 9×9 км² до 30×30 м². Оценки прогностических возможностей сопряжения моделей WRF и PALM в случае дневного АПС выполнены для двух дней, соответствующих условиям хорошей погоды и отсутствию облачности. Сравнение результатов моделирования с результатами наземных наблюдений, акустического и микроволнового зондирования атмосферы, а также высотных наблюдений с использованием размещенных на платформе привязного аэростата датчиков температуры и влажности показало, что согласованность результатов моделирования и наблюдений соответствует достигнутому ранее уровню. Однако при моделировании не удалось корректно воспроизвести форму вертикальных профилей температуры в средней части АПС и модуля горизонтальной скорости ветра, а также уменьшить отклонения между прогнозируемыми и наблюдаемыми значениями удельной влажности. Отклонения средних за 1 ч динамической скорости и полной ТКЭ, рассчитанных в LES, от результатов наблюдений составляют десятки процентов. Мгновенные поля температуры и модуля горизонтальной скорости ветра, рассчитанные в моделях WRF и PALM, обладают пространственным самоподобием в интервале масштабов, который находится внутри каждого из вложенных доменов. Скейлинговые экспоненты структурных функций первого порядка абсолютных разностей модуля горизонтальной скорости в доменах WRF на масштабах 3–100 км близки к 3/5, что соответствует скейлингу Больджиано–Обухова. Поведение структурных функций первого порядка поля температуры WRF отличается от скейлинга Больджиано–Обухова, рассчитанные структурные функции имеют приблизительно тот же показатель, что и структурные функции модуля горизонтальной скорости. Структурные функции первого порядка абсолютных разностей температуры и модуля горизонтальной скорости LES в области крупных масштабов (больше 300 м) могут быть локально аппроксимированы степенной зависимостью с показателем из диапазона 0.17–0.73, в области меньших масштабов существует участок с колмогоровской экспонентой, равной 1/3.

Учёт изменчивости метеорологических переменных и статистики турбулентных флуктуаций в процессе развития конвекции под воздействием мезомасштабного форсинга позволяет сделать вывод о том, что сопряжённая WRF-LES модель при достаточно подробном задании параметризаций, начальных и граничных условий, а также применении методов учёта неоднородностей рельефа, свойств поверхности, включая растительность, и методов для ускорения адаптации решения к более высокому разрешению корректно воспроизводит динамику АПС и может найти применение для решения задач турбулентного переноса и атмосферного электричества. При этом для верификации статистики турбулентности было бы полезно организовать высотные наблюдения релевантных величин.

Финансирование

Работа выполнена при финансовой поддержке Российского научного фонда (грант № 22-17-00053).

作者简介

S. Anisimov

Schmidt Institute of Physics of the Earth of the RAS

Email: svga@borok.yar.ru

Borok Geophysical Observatory

俄罗斯联邦, Borok, 142, Yaroslavl region, 152742

E. Mareev

Gaponov-Grekhov Institute of Applied Physics of the RAS

Email: svga@borok.yar.ru
俄罗斯联邦, Ulyanova, 46, Nizhny Novgorod, 603950

S. Galichenko

Schmidt Institute of Physics of the Earth of the RAS

编辑信件的主要联系方式.
Email: svga@borok.yar.ru

Borok Geophysical Observatory

俄罗斯联邦, Borok, 142, Yaroslavl region, 152742

N. Ilin

Gaponov-Grekhov Institute of Applied Physics of the RAS

Email: svga@borok.yar.ru
俄罗斯联邦, Ulyanova, 46, Nizhny Novgorod, 603950

A. Prokhorchuk

Schmidt Institute of Physics of the Earth of the RAS

Email: svga@borok.yar.ru

Borok Geophysical Observatory

俄罗斯联邦, Borok, 142, Yaroslavl region, 152742

E. Klimanova

Schmidt Institute of Physics of the Earth of the RAS

Email: svga@borok.yar.ru

Borok Geophysical Observatory

俄罗斯联邦, Borok, 142, Yaroslavl region, 152742

A. Kozmina

Schmidt Institute of Physics of the Earth of the RAS

Email: svga@borok.yar.ru

Borok Geophysical Observatory

俄罗斯联邦, Borok, 142, Yaroslavl region, 152742

K. Aphinogenov

Schmidt Institute of Physics of the Earth of the RAS

Email: svga@borok.yar.ru

Borok Geophysical Observatory

俄罗斯联邦, Borok, 142, Yaroslavl region, 152742

A. Guriev

Schmidt Institute of Physics of the Earth of the RAS

Email: svga@borok.yar.ru

Borok Geophysical Observatory

俄罗斯联邦, Borok, 142, Yaroslavl region, 152742

参考

补充文件