Геомагнетизм и аэрономия, 2019, T. 59, № 4, стр. 461-473

1. ВВЕДЕНИЕ

Наиболее ионизированной областью верхней атмосферы Земли является слой F2 ионосферы, который существует в любое время суток и имеет высокую изменчивость в зависимости от солнечной и геомагнитной активности. Ключевыми параметрами F2-слоя являются высота h_mF2 и электронная концентрации N_mF2. Электронная концентрации часто выражается через критическую частоту foF2, которая имеет с ней простую связь. Моделирование и прогнозирование временны́х и пространственных изменений параметров F2-слоя необходимо для решения задач распространения радиоволн, радиопеленгации, а также для эксплуатации глобальных навигационных спутниковых систем типа ГЛОНАСС и GPS. Разработка моделей прогноза ионосферных параметров в реальных геомагнитных условиях имеет большое значение. Для этого требуется хорошее знание структуры и динамики области F ионосферы. Важным параметром ионосферы, особенно в задачах распространения радиоволн КВ-диапазона, является критическая частота foF2. Описание даже геомагнитно-спокойного состояния распределения критических частот F2-слоя ионосферы представляет собой трудную задачу, которая решается с помощью эмпирических медианных моделей, таких как, например, Международная справочная модель ионосферы IRI-2016 [Bilitza et al., 2017].

Прогнозирование параметров ионосферы во время ионосферной бури представляет дополнительные трудности. Изменения foF2 во время геомагнитной бури зависят от мирового UT, локального LT времени, продолжительности бури, ее интенсивности, географических и геомагнитных координат, а также сезона и солнечной активности. Реакция ионосферы на большие геомагнитные бури представляет особый интерес, поскольку они приводят к резким возмущениям в ионосферно-термосферной системе в глобальном масштабе. Основными движущими силами таких возмущений в ионосфере являются изменения состава и ветра нейтральной атмосферы, которые возникают в области авроральных широт в результате нагрева термосферы токами Джоуля и высыпающимися частицами во время геомагнитных бурь [Prolss, 2011]. При высокой геомагнитной активности авроральный пик ионизации и вся субавроральная ионосфера, включая главный ионосферный провал (ГИП), оказываются на гораздо более низких широтах, чем при низкой геомагнитной активности. Возможность смещения субавроральной ионосферы на средние и даже в низкие географические широты означает, что для создания глобальной модели критической частоты F2-слоя ионосферы необходимо знание динамики субавроральной и авроральной области ионосферы. Известно, что высокоширотная ионосфера создает значительные проблемы для эмпирического моделирования из-за своей высокодинамичной природы. Международная справочная модель ионосферы (IRI) плохо отражает поведение foF2 в высоких широтах [Bilitza and Reinisch, 2008; Xiong et al., 2013]. Попытки создания модели foF2 для высокоширотной ионосферы предпринимались неоднократно [Chasovitin et al., 1987; Kutiev and Muhtarov, 2003; Bjoland et al., 2016]. В недавно построенной модели E-CHAIM (the Empirical Canadian High Arctic Ionospheric Model) [Themens et al., 2017] был использован региональный подход к моделированию электронной плотности на высотах F2-слоя ионосферы в области геомагнитных широт 50°–90° N. Следует отметить, что в климатологической модели E-CHAIM частично учтены поправки на зависимости пространственной структуры электронной концентрации в арктическом регионе от индексов геомагнитной активности. Однако использование медианной модели в области высоких широт без определения границ этих структур, например, положения минимума ГИП может привести к ошибкам. Как было показано в работе [Деминов и др., 2016], при низкой и средней солнечной активности ГИП не характерен для месячной медианы foF2.

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

2. МОДЕЛЬ ДЛЯ ГЕОМАГНИТНО-СПОКОЙНЫХ УСЛОВИЙ

В качестве фоновой модели foF2 используется глобальная медианная модель ионосферы SDMF2 (Satellite and Digisonde Data Model of the F2 layer) [Шубин, 2017] для спокойных геомагнитных условий. При построении модели SDMF2 были использованы данные 227 наземных ионосферных станций (1950–2014 гг.), а также радиозатменные измерения foF2, полученные с помощью спутников COSMIC (Constellation Observing System for Meteorology, Ionosphere, and Climate), GRACE (Gravity Recovery and Climate Experiment), CHAMP (Challenging Minisatellite Payload) для периодов высокой и низкой солнечной активности (2001–2015 гг.). Входными параметрами модели SDMF2 являются день, месяц, год, мировое время UT, географические координаты. В качестве зависимости foF2 от солнечной активности впервые был применен интегральный индекс F_10.7(day, τ), который представляет собой ряд предшествующих заданному дню суточных значений F_10.7 с экспоненциальным сглаживанием. Этот индекс определяется интегрированием величины F_10.7 на заданном интервале времени T. Для конкретного дня (day) он имеет следующий вид:

где n = 0, 1, 2, …Ts; τ = exp(–1/Ts). Проведенные нами расчеты показали, что период времени Ts = = 27 дней, примерно соответствующий одному обороту Солнца, дает наименьшие ошибки foF2 в глобальной медианной модели SDMF2 [Шубин, 2017]. При этом наибольшие веса относятся к недавнему прошлому, а сравнительно меньшие − к более ранним суточным величинам F_10.7. Здесь индекс F_10.7 − это поток солнечного радиоизлучения на длине волны 10.7 см, который измеряется в единицах 10^–22 Вт/(Гц м²).

Медианные модели пространственно-временнóго распределения foF2 основаны на разложении этих параметров в ряд по системе взаимно-ортогональных функций [Jones and Gallet, 1962; Шубин, 2017]. Обычно при построении эмпирических моделей используется линейная зависимость foF2 от солнечной активности (например, [Zhang et al., 2009, 2014; Altadill et al., 2013]). В работе [Richards, 2001] было показано, что с ростом солнечной активности электронная концентрация в максимуме F2-слоя ионосферы стремится к насыщению. Проведенные нами исследования показали, что логарифмическая зависимость критической частоты от солнечной активности дает лучшую корреляцию foF2 с индексом F_10.7(day, τ), чем линейная зависимость. В общем виде построенную модель SDMF2 критической частоты foF2 можно описать следующей формулой [Шубин, 2017]:

Зависимость foF2 от координат в фиксированный момент мирового времени UT представляется в виде ряда:

где θ = 90°–φ, φ – географическая широта; λ – географическая долгота; $g_{n}^{m}$ и $h_{n}^{m}$ – коэффициенты разложения, а $P_{n}^{m}\left( {\cos \vartheta } \right)$ – присоединенные функции Лежандра. Для интерполяции по времени использовано разложение Фурье с периодом T = 24 ч. В качестве индекса солнечной активности в построенной модели SDMF2 используется интегральный индекс F_10.7(day, 0.96) на заданный день. Для проверки модели SDMF2, кроме интегрального индекса F_10.7(day, 0.96), нами были использованы другие солнечно-ионосферные индексы, которые широко применяются в различных эмпирических моделях ионосферы – R12, IG12 (индексы R, IG, усредненные за 12 месяцев) и F_10.7A (индекс F_10.7, усредненный за 3 оборота Солнца). Рассчитанные по модели SDMF2 величины foF2 с разными входными индексами R12, IG12, F_10.7A и F_10.7(day, 0.96) были сравнены с наблюдаемыми данными ионосферных станций, расположенными в южном и северном полушариях: Scott Base (77.9° S, 166.8° E), Port Stanley (51.7° S, 302.2° E), Townsville (19.7° S, 146.9° E), Boulder (40.0° N, 254.7° E), Moscow (55.5° N, 37.3° E) и Tromsø (69.7° N, 19.0° E) за период с января 2010 г. по декабрь 2016 г. Результаты проверки показали, что наименьшие значения среднеквадратических ошибок RMSE были получены в модели SDMF2 с входным индексом солнечной активности F_10.7(day, 0.96). Для спокойных условий наибольшее улучшение точности модели SDMF2 относительно модели IRI-2016 (URSI) было получено для станций южного полушария: на ~25% для ст. Port Stanley и на ~17% для ст. Townsville [Shubin, 2018].

Преимущество индекса F_10.7(day, 0.96) состоит также в том, что он в отличие от индексов F_10.7A или IG12 и R12, используемых в моделях IRI-2016 [Bilitza et al., 2017], E‑CHAIM [Themens et al., 2017], не требует прогноза на 40 дней или 6 месяцев вперед индексов F_10.7 или IG и R, соответственно. Проведенное тестирование модели SDMF2 [Шубин, 2017] по наземным и спутниковым данным показало, что разработанная модель обеспечивает более высокую точность по сравнению с моделью IRI-2016. Наибольшие различия моделей SDMF2 и IRI-2016 в геомагнитных спокойных условиях наблюдаются при низкой солнечной активности, а также в южном полушарии и над акваториями океанов.

3. МОДЕЛЬ ИОНОСФЕРНОЙ БУРИ

Для построения динамической модели foF2 нами применен аэрономический подход к описанию ионосферной бури. Он основан на учете физических процессов, происходящих в термосфере. Впервые попытка использовать этот подход для прогноза foF2 была предпринята в работе [Шубин и Аннакулиев, 1995]. В дальнейшем с появлением более совершенных термосферных моделей типа NRLMSISE-00 [Picone et al., 2002] и более надежных лабораторных измерений констант скоростей основных ионно-молекулярных реакции [Hierl et al., 1997] данный подход получил свое развитие [Шубин и Аннакулиев, 1997; Аннакулиев и др., 2005]. Результаты его применения изложены в данном разделе.

Обозначим величины, соответствующие геомагнитно-спокойным условиям, индексом q. Предположим, что изменения критической частоты F2-слоя ионосферы связаны с относительными изменениями нейтрального состава и температуры термосферы следующим соотношением:

где коэффициент C_storm = R(t)/R_q(t), а в качестве значений foF2_q для спокойных условий используется медианная модель ионосферы SDMF2 [Шубин, 2017]. Величины R и R_q являются функциями только параметров термосферы на высоте 300 км. Для их получения была использована аналитическая модель среднеширотной ионосферы [Аннакулиев и др., 2005]: Здесь $n\left( {\text{O}} \right),$ $n\left( {{{{\text{N}}}_{2}}} \right)$ и $n\left( {{{{\text{O}}}_{2}}} \right)$ – концентрации нейтральных составляющих верхней атмосферы; $m\left( {\text{O}} \right),$ $m\left( {{{{\text{N}}}_{2}}} \right),$ $m\left( {{{{\text{O}}}_{2}}} \right)$ – молекулярные веса нейтральных составляющих (в см^–3); T_n –температура нейтральных частиц атмосферы (в K°) на высоте h = 300 км; β – коэффициент линейной рекомбинации ионов O⁺; k₁ и k₂ – коэффициенты диссоциативной рекомбинации O⁺ с N₂ и O₂ [Hierl et al., 1997]. Для определения R используется эмпирическая модель NRLMSISE-00 [Picone et al., 2002], которая является глобальной моделью термосферы, применимой при любом уровне солнечной и геомагнитной активности. В ней косвенно учтено, что основной причиной изменений параметров термосферы в периоды геомагнитных бурь является Джоулев нагрев термосферы высоких широт. В целом термосфера реагирует на такое поступление энергии как слабо пропускающий фильтр, сглаживая отклик параметров термосферы (температуры, плотности, состава), который в каждый момент времени зависит от предыстории поступления энергии в термосферу. Это свойство термосферы учтено в эмпирической модели NRLMSISE-00 с помощью эффективного интегрального индекса геомагнитной активности a_p(t,τ), где τ = exp(–3/T_m) с характерным временем T_m. В модели NRLMSISE-00 характерное время реакции термосферы на изменения геомагнитной активности зависит от широты и различно для температуры экзосферы T_∞ и концентраций каждого из компонентов термосферы. В среднем T_m ≈ 7–12 ч [Hedin, 1987]. Анализ экспериментальных данных показывает, что величина a_p(t,τ) является достаточно точным индикатором геомагнитной активности изменений foF2(t) во время геомагнитной бури [Wrenn and Rodger, 1989]. Использование эмпирической модели термосферы NRLMSISE-00 в качестве составной части полуэмпирической модели бури в ионосфере позволяет отразить важное свойство ионосферы, а именно, зависимость foF2(t) от предыстории изменения геомагнитной активности через соответствующую зависимость параметров термосферы. Статистический анализ показал, что построенная модель воспроизводит зависимость foF2(t) от изменений в составе термосферы на средних широтах во время ионосферной бури любой интенсивности, за исключением начальных периодов геомагнитной бури.

4. МОДЕЛЬ СУБАВРОРАЛЬНОЙ И АВРОРАЛЬНОЙ ИОНОСФЕРЫ

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

Отличительной чертой разработанной модели является использование фонового распределения foF2_q и ряда поправок к этому распределению, связанных с формированием главного ионосферного провала (ГИП), аврорального пика foF2, изменением температуры и состава термосферы. В качестве фонового распределения foF2_q используется медианная модель ионосферы SDMF2 [Шубин, 2017]. Все поправки зависят от изменяющихся гелиогеофизических условий. Модель foF2 субавроральной и авроральной ионосферы включает в себя модели положения минимума главного ионосферного провала Φ_MIT и аврорального пика критической частоты Φ_AVR.

Модели положения минимума ГИП и аврорального пика критической частоты F2-слоя основаны на данных зондовых измерений плотности электронов спутником CHAMP с июля 2000 г. по декабрь 2007 г. на высотах от ~450 до ~350 км. С помощью ручной обработки данных каждого анализируемого витка спутника были определены географические координаты минимума ГИП. Затем по международной модели геомагнитного поля IGRF [Maus et al., 2005] и полученным географическим широтам были определены исправленные геомагнитные широты минимума ГИП – Φ_MIT и максимума аврорального пика электронной концентрации – Φ_AVR на высоте 400 км. Модели Φ_MIT и Φ_AVR представлены в виде аналитических зависимостей инвариантной широты минимума провала и аврорального пика критической частоты от местного солнечного времени, геомагнитной активности и географической долготы. В качестве индикатора геомагнитной активности использован средневзвешенный по времени индекс K*(t,τ) [Деминов и Шубин, 2018]:

где коэффициент τ = 0.6 определен из условия минимума отклонения моделей от экспериментальных данных. Из τ = exp(–3/T) следует, что характерное время индекса K*(t,τ) составляет ~6 ч.

Полученная эмпирическая модель положения минимума ГИП имеет вид [Деминов и Шубин, 2018]:

где широты Φ_MIT, Φ(t), Φ(λ) и географическая долгота λ определены в градусах, время t в часах. Величины Φ(λ) = Φ_N(λ) и Φ(λ) = Φ_S(λ) для северного и южного полушарий определяются соотношениями:

Свободный член в Φ(t) подобран так, чтобы в полночь при t = 0 было выполнено условие Φ(t) = 0. Модель не имеет ограничения по местному времени и величине индекса геомагнитной активности. Отметим, что массив минимумов Φ_MIT содержал в основном данные, для которых было выполнено условие K*(t,τ) ≤ 6. Поэтому точность модели для K*(t,τ) > 6 может быть невысокой.

Модель положения аврорального пика критической частоты F2-слоя представим в виде:

Здесь K = 1.2 K*(t,τ) – 1.0. При этом на величину K наложены следующие условия: K = 0, если K < 1 и K = 8, если K > 8.

Для построения модели foF2 в субавроральной и авроральной ионосфере введем поправки к фоновому (медианному) распределению foF2_q. При этом поправку к foF2_q в области ГИП обозначим C_MIT, а поправку для аврорального пика электронной концентрации– C_AVR. Если указанные поправки известны, то для вычисления foF2 в этих областях ионосферы введем следующее уравнение:

Поправка C_AVR вычисляется по известным исправленным геомагнитным широтам заданной точки – Ф, аврорального пика – Φ_AVR и минимума ГИП – Φ_MIT, местному времени SLT для данной даты, данного мирового времени и эффективному индексу геомагнитной активности K:

где

Из уравнения (6) видно, что (C_AVR)_max представляет собой поправку C_AVR на широте Ф = Φ_AVR, где она максимальна. При удалении от этого максимума поправка С_AVR быстро уменьшается c широтой и, например, на широте Ф = Ф_AVR – 2ΔФ₁ = = Ф_MIT эта поправка становится очень малой величиной, которой можно пренебречь – C_AVR = = 0.018(C_AVR)_max. Поправки C_AVR и C_MIT существенны только в области аврорального пика и ГИП, соответственно. Вне этих областей C_AVR$ \ll $ 1 и C_MIT = 1, поэтому foF2 = C_stormfoF2_q для широт, расположенных к экватору от ГИП. Структура уравнения (5) при C_storm = 1 подобна использованной Tascione et al. [1988], но все параметры этого уравнения, включая значения C_MIT и C_AVR, положения аврорального пика foF2 и минимума ГИП, а также используемые индексы солнечной и геомагнитной активности, существенно отличаются от величин, приведенных в работе [Tascione et al., 1988].

На широте Ф = Φ_AVR расположен авроральный пик критической частоты F2-слоя. Из уравнения (5) следует:

В этом соотношении косвенно учтено, что авроральный пик foF2 связан с ионизацией атмосферы потоками авроральных электронов, интенсивность высыпания которых увеличивается с ростом геомагнитной активности и слабо зависит от сезона. Другие факторы, включая ионизацию атмосферы солнечным излучением, зависят от сезона через величину C_stormfoF2_q. В итоге величина foF2 _AVR также зависит от сезона. Например, при низкой геомагнитной активности, когда C_storm ≈ 1, авроральный пик foF2 достаточно отчетливо выражен в ночные часы зимой и практически отсутствует летом, поскольку летом на широте Ф = Φ_AVR во все часы суток выполнено условие – (foF2_q)² > > (C_AVR)_max.

Поправка C_MIT в уравнении (5) вычисляется по известным исправленным геомагнитным широтам заданной точки Ф, аврорального пика Φ_AVR и минимума ГИП Ф_MIT, а также местному времени SLT, зенитному углу Солнца χ, дню года nday и эффективному индексу солнечной активности FF = [F_10.7 + 2F_10.7(day, τ)]/3 в соотношении:

где На величину ΔФ₂, вычисленную из уравнения (9), наложено ограничение: На величину C_MIT наложено дополнительное ограничение, связанное с величиной поправки C_storm над данным пунктом для данной даты и данного мирового времени: Дополнительное ограничение, наложено на foF2 в области ГИП: т.е. величина foF2 в этой области не может быть меньше 1.3 МГц.

Из приведенных соотношений (8–12) можно видеть, что поправка C_MIT определяет глубину и форму главного ионосферного провала. На широте минимума ГИП поправка C_MIT достигает минимального значения C_MIT = 1 – (C_MIT)_max, и это значение является глубиной провала. При удалении от этого минимума величина С_MIT увеличивается вплоть до С_MIT = 1, где провал отсутствует, а параметр ΔФ₂ является характерным масштабом, определяющим ширину экваториальной и полярной стенок провала. В соответствии с данными наблюдений провал обычно наблюдается в неосвещенное время суток, что отражено в уравнении (9). Глубина провала максимальна зимой в послеполуночные часы при низкой солнечной активности (см. явный вид коэффициентов T_d, T_hr и T_F), ширина экваториальной стенки провала обычно больше ширины его полярной стенки (см. явный вид ΔФ₂), но ширина экваториальной стенки редко превышает 5° (см. уравнение 10).

Полярная стенка ГИП формируются не только процессами, приводящими к образованию минимума провала, но и процессами, связанными с высыпаниями авроральных электронов. В модели это учтено согласованием характерных масштабов ΔФ₁ для аврорального пика foF2 и ΔФ₂ для ГИП: ΔФ₁ = ΔФ₂ = (Ф_AVR – Ф_MIT)/2 в интервале Ф_MIT < Ф < Ф_AVR (см. явный вид ΔФ₁ и ΔФ₂ в уравнениях (7) и (9)). Это обеспечивает также непрерывность и плавность изменения foF2 при переходе от минимума ГИП к авроральному пику foF2 для всего возможного диапазона изменений уровней солнечной и геомагнитной активности.

Величины Ф_AVR (4) и (C_AVR)_max (7), зависят от геомагнитной активности через средневзвешенный индекс K*(t,τ) c характерным временем ~6 ч. Величина C_storm также зависит от геомагнитной активности, но с характерным временем ~11 ч. Тем самым косвенно учтено, что реакция аврорального пика foF2 на изменение геомагнитной активности происходит с двумя характерными временами. Вначале изменяется положение и величина этого пика, связанная с изменениями потока высыпающихся авроральных электронов и области локализации этих высыпаний, затем происходит дополнительное изменение величины аврорального пика foF2, связанное с изменениями температуры и состава термосферы.

Вклад относительного изменения температуры и состава термосферы в формирование ГИП увеличивается с ростом геомагнитной активности, что косвенно учтено через поправку C_storm, которая уменьшается не только с ростом геомагнитной активности, но и с широтой. Вклад поправки C_storm в формирование ГИП становится преобладающим при C_storm < 0.75 (11). Тем самым косвенно учтено, что провал редко бывает очень глубоким, т.е. чаще среди разных процессов, связанных с образованием провала, один из процессов становится преобладающим, определяя формирование провала. Это, по-видимому, связано с разными характерными временами протекания процессов, влияющих на образование провала.

5. РЕЗУЛЬТАТЫ ПРОВЕРКИ МОДЕЛИ И ОБСУЖДЕНИЯ

Для наглядности, на рис. 1 приведен пример последовательного вычисления широтного распределения foF2 по уравнению (5) (с шагом по широте 0.1°) на долготе λ = 45° E в 21 UT 14.12.2006 г. Эффективные индексы солнечной и геомагнитной активности для этого мирового времени и данной даты следующие: F_10.7(day, 0.96) ≈ ≈ 90 и a_p(t,τ) ≈ 80 нТл, что соответствует низкой солнечной и повышенной геомагнитной активности. Решение задачи разбивается на четыре этапа.

1. По модели SDMF2 [Шубин, 2017] вычисляется фоновое значение foF2_q.

2. Из уравнения (1) определяется распределение foF2, связанное с относительными изменениями нейтрального состава и температуры термосферы. Величина R в уравнении (2) вычисляется по эмпирической модели термосферы NRLMSISE-00 [Picone et al., 2002].

3. Из уравнения (6) вычисляется поправка C_AVR, и по уравнению (5) определяется величина foF2 в субавроральной ионосфере с учетом аврорального пика электронной концентрации.

4. Из уравнения (8) c учетом условий (9–12) вычисляется поправка C_MIT и определяется величина foF2 (уравнение 5) в области ГИП. Исправленные геомагнитные широты в уравнениях поправок C_MIT и C_AVR находятся из моделей IGRF, которые входят в IRI-2016 [Bilitza et al., 2017].

Из рисунка 1 видно, что зависимость foF2 от геомагнитной активности в модели GDMF2 обусловлена в основном двумя факторами – уменьшением foF2 с ростом геомагнитной активности на средних, субавроральных и авроральных широтах из-за изменения температуры и состава термосферы через поправку C_storm и увеличением foF2 в авроральной области через поправку C_AVR.

В качестве примеров на рис. 2а–2в приведены результаты моделирования изменения критической частоты слоя F2 с широтой во время геомагнитных возмущений разной интенсивности для трех сезонов (равноденствия, лета, зимы) и их сравнения с измеренными данными foF2. На этих рисунках сплошной толстой кривой обозначены результаты расчетов foF2 по модели GDMF2, тонкой кривой – по модели SDMF2 для спокойных геомагнитных условий. Результаты расчета foF2 по модели IRI-2016 с опцией STORM = ON представлены штриховой кривой. Крестиками обозначены значения foF2 наземных ионосферных станций, черными точками обозначены величины foF2, полученные на спутнике. На рис. 2а, 2б представлены расчеты foF2 на λ = 15° E, а на рис. 2в на λ = 20° E. Выбор долготы и мирового времени в приведенных примерах определялся наличием как можно большего числа спутниковых и наземных данных foF2 вблизи этих λ и UT. Экспериментальные данные foF2, представленные на рис. 2 были получены в диапазоне ±7.5° от указанной долготы λ и в интервале ±0.3 ч от указанного времени UT. Отметим, что данные 2004 г. не были использованы для построения фоновой медианной модели SDMF2, так как в этом году величина сглаженного индекса солнечной активности F_10.7A находилась в интервале 80 < F_10.7A < 120 (см., Шубин [2017]). Таким образом, фактически осуществляется проверка прогностической способности модели, а не только ее соответствия набору использованных при построении данных.

На рисунке 2а приведен пример широтной зависимости foF2 в равноденствие. Как видно из этого рисунка, уменьшение foF2 относительно фоновой модели достигает географических широт ~±30° в северном и в южном полушарии. При значении K_p≈5.3 ст. Juliusruh (φ = 54.6°, λ = 13.4°, Ф = 50.7º) оказывается в области ГИП. В южном полушарии области ГИП соответствуют данные foF2, полученные методом RO с помощью спутника CHAMP. На рис. 2б приведен пример широтный зависимости foF2 для лета (северное полушарие) во время очень сильной и продолжительной геомагнитной бури. Она началась 09 UT 22 июля 2004 г. и продолжалась ~7 дней (максимальный значение индекса D_st ~ 200 нТл в 13 UT 27 июля). В 23 LT положение минимума ГИП на долготе ст. Juliusruh достигло Φ_MIT = 48.6°. Отметим, что максимальное положение минимума ГИП на этой долготе по модели (3) составило Φ_MIT = 45.6° 27 июля в 03 UT. В данном примере в летнем северном полушарии наблюдаются очень большие различия между моделью для спокойных условий (тонкая линия) и моделью с учетом изменений состава и температуры термосферы (толстая линия). Это связано с высокой геомагнитной активностью: в 23 LT величина индекса K_p ≈ 7.3. На рис. 2в приведен пример широтный зависимости foF2 для зимы (северное полушарие). В данном примере K_p ≈ 4.7, поэтому отличия между фоновой SDMF2 и динамической моделью GDMF2 не настолько большие, как на рис. 2б. Как видно из приведенных на рис. 2 сравнений с экспериментальными данными, модель GDMF2, по крайней мере, качественно ведет себя адекватно.

На рисунке 3 приведено моделирование foF2 на ст. Moscow (φ = 55.5°, λ = 37.3°, Ф = 51.0°) во время очень сильной геомагнитной бури, которая произошла 13–17 марта 1989 г. На верхней панели рис. 3 приведены изменения со временем геомагнитных индексов K_p (толстая) и D_st (тонкая) кривая. На нижней панели представлены изменения значений foF2 (отмечены крестиками), рассчитанные по модели GDMF2 – сплошной и IRI-2016 – штриховой линиями. Индексы геомагнитной активности во время этой бури 14 марта в 03 UT достигли своих максимальных величин K_p = 9, D_st = –589 нТл. Солнечная активность была также очень высокая F_10.7 ≈ 266. При этом модель GDMF2 довольно точно отражает резкое падение foF2 до ~2 МГц, что характерно для отрицательной фазы ионосферной бури. Минимум ГИП, равный Φ_MIT = 41.0°, в 03 UT по модели (3) оказался южнее ст. Moscow почти на 10°. Поэтому уменьшение foF2 могло быть обусловлено как изменениями в составе (foF2 ~ O/N₂) и температуре термосферы, так и нахождением ионосферной станции в области ГИП. В течение 00–12 LT 14.03.1989 г. ст. Moscow находилась вблизи минимума ГИП, вычисленного по уравнению (3). В эти часы модель GDMF2 в отличие от модели IRI-2016 более точно воспроизводит foF2, так как в модели IRI-2016 нет поправки на foF2, связанной с влиянием ГИП.

Для проверки модели в области высоких широт были выбраны ст. Dikson (φ = 73.5°, λ = 80.4°, Ф = 67.9°), которая обычно находится на полюсной границе аврорального овала и ст. College (φ = = 64.9°, λ = 212.2°, Ф = 65.1°), которая обычно находится вблизи экваториальной границы аврорального овала. На рисунках 4 и 5 представлены зависимости foF2 во время геомагнитного возмущения, во время которого в ионосфере наблюдалось последовательно две бури. Первая началась 13.05.1991 г., и практически сразу за ней последовала вторая буря 17.05.1991 г. Как видно из этих рисунков, построенная модель достаточно хорошо воспроизводит foF2 в разных долготных секторах. Значительное увеличение foF2 12.05.1991 г. с 16–23 LT и 14.05.1991 г. с 14–19 LT на ст. Dikson (рис. 4), вероятно, связано с ее нахождением в области аврорального пика электронной концентрации. При этом, хотя модельное значение Ф_AVR из соотношения (4) практически совпадает с геомагнитной широтой Ф ст. Dikson, коэффициент (C_AVR)_max в (7) не обеспечивает достаточный рост foF2. Но в остальные часы всюду модель GDMF2, в отличие от модели IRI-2016, достаточно хорошо воспроизводит изменения foF2. На ст. College 15.05.1991 г. модель GDMF2 плохо воспроизводит дневные значения foF2 (рис. 5). Это может быть связано как с недостатками фоновой модели SDMF2, так и с более сложной зависимостью границ субавроральной ионосферы от геомагнитной активности в уравнениях для Ф_MIT (3) и Ф_AVR (4).

Из приведенных рисунков можно видеть, что динамическая модель GDMF2 количественно правильно воспроизводит экстремумы отрицательной фазы ионосферных бурь (минимальное значение foF2) и время наступления этих экстремумов. Видны и недостатки модели. В дневные часы после максимума главной фазы ионосферной бури модель дает завышенные значения foF2 (14.03.1989 г. на рис. 3) относительно экспериментальных данных. Этот эффект требует дальнейшего детального исследования. Кроме того, модель не воспроизводит положительную фазу ионосферной бури (12.03.1989 г. на рис. 3). Приведенные предварительные результаты в целом показывают перспективность данной модели, поскольку ее точность не сильно ухудшается при переходе от спокойных условий к периодам очень интенсивных геомагнитных бурь. Главным достоинством представленной модели является ее возможность воспроизводить отрицательную фазу ионосферной бури во время интенсивной геомагнитной бури. Точность данной модели определяется главным образом точностью используемой медианной модели foF2 (SDMF2) в качестве фона для геомагнитных спокойных условий и модели термосферы (NRLMSISE-00) для возмущенных условий.

Соотношения (2) были получены из решения уравнения непрерывности для электронной концентрации дневной среднеширотной F-области ионосферы. Однако, как показали результаты статистической обработки данных foF2 среднеширотных станций Moscow и Иркутск [Шубин и Аннакулиев, 1995], эти соотношения приближенно можно использовать для всех часов суток во время геомагнитных бурь. В созданной модели зависимость foF2 от геомагнитной активности обусловлена двумя факторами: уменьшением foF2 с ростом геомагнитной активности на средних и более высоких широтах из-за изменения температуры и состава термосферы через поправку C_storm, а также увеличением foF2 в авроральной области через поправку C_AVR (см. уравнения (2) и (5)).

На рисунке 6 представлены среднеквадратические ошибки моделей GDMF2 и IRI-2016, полученные при сравнении с радиозатменными данными спутников GRACE в 2016 г. Отметим, что эти данные не использовались при построении медианной фоновой модели SDMF2. Линии проведены через RMSE (МГц), полученные в каждом часе местного времени. Модель GDMF2 – толстые линии (точки) IRI-2016(URSI) – тонкие линии (крестики). Буквы обозначают широтные зоны: а – экваториальные (30° S–30° N), б – средние (30° S–60° S, 30° N–60° N) и в – высокие широты (60° S–90° S, 60° N–90°N). Как видно из этого рисунка, наибольшие ошибки модели достигаются на экваториальных широтах. Это связано как с небольшим числом наземных ионосферных станций в этой области, использовавшихся при построении фоновой модели SDMF2, ошибками в определении foF2 радиозатменным методом, так и со сложностью динамических процессов на этих широтах. Но, тем не менее, как видно из приведенных на рис. 6 среднеквадратических ошибок, построенная глобальная динамическая модель foF2 демонстрирует заметное улучшение по сравнению с IRI-2016 в среднем на ~16% в высоких и на средних и на ~7% на экваториальных широтах.

Созданная динамическая модель GDMF2 предназначена для вычисления foF2 как для спокойных условий, так и периодов геомагнитных бурь, главным образом, на средних, субавроральных и авроральных широтах. Более высокая точность разработанной модели по сравнению с IRI-2016 обусловлена, в первую очередь, более точной фоновой моделью SDMF2 для спокойных условий [Шубин, 2017], а также введением в разработанную модель GDMF2 поправок foF2, связанных с образованием ГИП и аврорального пика foF2, которые отсутствуют в IRI-2016. Значения foF2, вычисленные по этой модели на экваториальных широтах и в области полярной шапки, являются приближенными и определяются в основном фоновой моделью.

Отметим прогностические возможности предлагаемой модели. Зависимость от солнечной активности в модели определяется интегральным индексом F_10.7(day, 0.96). Для его определения необходимы данные индексов F_10.7, предшествующие заданному дню, и индекс F_10.7 в данный день. В качестве характеристики геомагнитной активности в термосферный параметр R(t) уравнения (1) входит индекс a_p(t,τ), который также учитывает предысторию изменения уровня геомагнитной активности. Для среднего значения характерного времени τ ∼ 10 ч текущее значение a_p(t) входит в a_p(t,τ) с относительно малым весом. Поэтому точного знания планетарных индексов a_p(t) или K_p(t) к моменту времени t для определения foF2(t) не требуется. Следовательно, простейшим способом прогноза foF2(t) является прогноз на 3 ч вперед при предположении, что измеренные за предыдущие 3 ч значения индексов K_p(t) или a_p(t) соответствуют последующим трем часам.

6. ЗАКЛЮЧЕНИЕ

Построена глобальная динамическая модель, предназначенная для вычисления foF2 как в спокойных, так и геомагнитно-возмущенных условиях. Термин “динамическая модель” подразумевает зависимость изменения foF2 на средних, субавроральных и авроральных широтах от текущего уровня геомагнитной активности с учетом предыстории ее развития. Отличительной чертой разработанной модели является использование медианного (для спокойных геомагнитных условий) распределения foF2 и ряда поправок к этому распределению, связанных с формированием главного ионосферного провала, аврорального пика электронной концентрации, изменением температуры и состава термосферы. При этом все поправки зависят от изменяющихся гелиогеофизических условий. Проведенное тестирование динамической модели показало, что разработанная модель обеспечивает более высокую точность по сравнению с моделью IRI-2016.

Использование в качестве входных параметров модели GDMF2 интегральных индексов солнечной F_10.7(t,τ) и геомагнитной a_p(t,τ) активности, которые учитывают предысторию их изменения, позволят осуществлять краткосрочный прогноз электронной концентрации максимума F2-слоя ионосферы. При этом продолжительность прогноза foF2 будет определяться прогнозом солнечных F_10.7 и геомагнитных (K_p или a_p) индексов.

7. БЛАГОДАРНОСТИ

Авторы благодарны спонсорам и операторам миссии CHAMP, GRACE (GFZ) и COSMIC (UCAR). Радиозатменные данные foF2 были получены из профилей электронной концентрации, размещенных на сайте COSMIC Data Analysis and Archival Center (CDAAC, http://cdaac-www.cosmic. ucar.edu/cdaac/tar/rest.html), величины foF2 наземных ионосферных станций были взяты с сайта Space Physics Interactive Data Resource (SPIDR, http://spidr.ngdc.noaa.gov/). Данные индексов геомагнитной активности были взяты с сайта World Data Center for Geomagnetism, Kyoto (http://wdc.kugi.kyoto-u.ac.jp/).

8. ФИНАНСИРОВАНИЕ РАБОТЫ

Работа частично поддержана Российским фондом фундаментальных исследований (грант № 17-05-00427) и Программой 28 Президиума РАН.