Новый алгоритм расчета инсоляции Земли
И.И. Смульский
Аннотация
Для исследования палеоклимата М. Миланковичем разработана теория инсоляции Земли. Расчеты инсоляции по упрощенным методам используются также в метеорологии, в строительстве и в др. областях. Новый алгоритм расчета инсоляции основан на результатах точного решения задачи 2-х тел. Он отличается от прежнего метода простотой и большей приспособленностью к компьютерным вычислениям. Приведено обоснование нового алгоритма, и представлена его реализация в среде MathCad. Рассчитаны распределения инсоляции по широте за год и за калорические полугодия в современную эпоху. Вычислено изменение инсоляции в эквивалентных широтах за 200 тыс. лет. Результаты расчетов сопоставлены с результатами других авторов и обоснована достоверность нового метода.
Рассчитана динамика инсоляции в современную эпоху. Показаны перспективы применения ее результатов для выяснения причин изменения природных процессов, которые зависят от инсоляции.
инсоляция земля алгоритм
1. Введение
Количество тепла Солнца, поступающее на Землю, т.е. инсоляция Земли, является важным фактором процессов, происходящих на нашей планете. Одной из наиболее вероятной причиной геологических изменений в истории Земли может являться изменение ее инсоляции. На протяжении почти двух столетий изучается влияние инсоляции на климат Земли прошедших эпох. В этих исследованиях используется методика расчета инсоляции, которую в завершенном виде представил М. Миланкович в своих работах, например, в работе [1]. Она является достаточно сложной. Поэтому ряд операций, как, например, вычисление инсоляции в эквивалентных широтах, обычно не проводится. В последнее время из-за сложности этой теории в литературе наблюдается ее непонимание. Например, в работе [2] утверждается, что “в своих расчетах Миланкович пренебрег прямым вкладом изменений эксцентриситета в инсоляционную кривую…”. В действительности М. Миланкович создал теорию расчета инсоляции Земли с учетом всех факторов, в том числе и эксцентриситета. Поэтому создание более простого алгоритма вычисления инсоляции будет способствовать ее пониманию и адекватному применению.
Изменение инсоляции по поверхности Земли и ее динамика представляет интерес специалистам в различных отраслях деятельности человека и в разных областях науки. Например, в строительстве расчеты инсоляции и ее нормирования играют важную роль для обоснования застройки городов [3]. При исследовании причин изменения различных факторов в биологии, медицине и в метеорологии все чаще привлекается динамика инсоляции для поиска корреляционных связей. В этих работах авторы вынуждены развивать свои подходы по расчету инсоляции земной поверхности и влиянии на нее движения Венеры [4], других планет [5] и Луны [6]. В этих подходах, как правило, из всех факторов, от которых зависит инсоляция, учитывается только изменение расстояния между Землей и Солнцем. Например, в работе [6] получена амплитуда колебания инсоляции 84.9 мВт/м2 за счет месячных изменений этого расстояния при движении Луны. Из приведенных примеров использования инсоляции видно, что создание простого алгоритма ее расчета с учетом всех факторов, от которых она зависит, является своевременным.
В разработанном М. Миланковичем [1] алгоритме расчета инсоляции Земли для определения зависимости долготы Солнца от времени t приближенным аналитическим методом интегрируется дифференциальное уравнение. Оно следует из второго закона Кеплера. На основании решений задачи двух тел мы получили [7] - [8] точную аналитическую зависимость для времени от полярного угла o или долготы . Эта зависимость значительно упрощает выражение для инсоляции. Поэтому появляется большая ясность метода, и упрощаются вычисления.
2. Основные результаты задачи двух тел
Для расчета инсоляции Земли необходимо определять положения Солнца над точкой поверхности Земли. Так как относительно Земли происходит движение Солнца по орбите, идентичной орбите Земли, то параметры орбиты Земли используются для описания движения Солнца в течение года. Аналогично параметры вращательного движения Земли используются для расчета суточного движения Солнца.
Рис. 1. Схемы движения Земли (E) по орбите вокруг Солнца (a) и Солнца (S) относительно Земли (б): - точка весеннего равноденствия; PE - перигелий Земли; PS - перигей Солнца; о - полярный угол Земли; p - угол перигелия Земли; = EPS - угол перигея Солнца; - долгота Солнца; пунктирная линия и точки E, PE и E на ней являются условными изображениями орбиты Земли на схеме орбитального движения Солнца.
На рис. 1 представлены схемы движений Земли относительно Солнца и Солнца относительно Земли. Плоскость орбиты Земли (рис. 1а) пересекает плоскость экватора по линии S. Точку весеннего равноденствия по своей орбите Солнце проходит весной, а Земля - осенью. В т. Солнце находится в плоскости экватора Земли, поэтому по длительности день равняется ночи. Плоскости экватора и орбиты Земли изменяются в пространстве, поэтому точка по орбите Земли (рис. 1а) перемещается. Положение перигелия PE также перемещается по орбите. Угол p измеряется между этими двумя точками и PE.
Будем рассматривать движение планеты с массой m при ньютоновском воздействии на нее со стороны Солнца, масса которого равна М. Дифференциальные уравнения ее движения относительно Солнца в безразмерных величинах имеют вид [7]:
, (1)
где - безразмерный радиус положения планеты относительно Солнца;
= t·vp/Rp - безразмерное время;
- параметр траектории;
G(M+m) - параметр взаимодействия;
G - гравитационная постоянная;
- радиус перигелия;
vp - скорость планеты в перигелии.
В безразмерных переменных , движение, как следует из (1), полностью определяется параметром траектории 1. Следует отметить, что этот параметр идентичен эксцентриситету e, и они связаны следующим выражением: e = - (1 + 1)/1. В дифференциальных уравнениях электромагнитного взаимодействия [7] параметр 1 играет аналогичную роль параметра траекторий, к которым понятие “эксцентриситет” неприменимо. Поэтому предпочтительно использовать параметр 1, а не эксцентриситет e.
В результате интегрирования уравнения (1) получено уравнение траектории и зависимости для времени движения t(r) от расстояния r между телами для различных случаев движения [7]. Уравнение траектории в полярной системе координат (r, о) в размерных переменных имеет следующий вид
, (2)
где о - полярный угол тела на орбите, отсчитываемый от радиуса перигелия r = SPE = (см. рис. 1а).
Уравнение (2) при 1 = -1 представляет окружность, при -1 < 1 < -0.5 - эллипс, при 1 = -0.5 - параболу, при -0.5 < 1 < 0 - гиперболу, а при 1 = 0 - прямую. Следует отметить, что при известной орбите параметр траектории 1 может быть определен через наименьший и наибольший радиусы орбиты следующим образом:
.
Из этой формулы видно, что для окружности = и параметр 1 = -1, а для параболы и 1 = - 0.5.
Скорость в перигелии согласно решениям задачи двух тел [8] рассчитывается по формуле
, (3)
где Psd - период обращения. В задаче двух тел он рассматривается относительно неускоренной системы координат. Для планеты - это период ее обращения вокруг Солнца относительно неподвижных звезд (сидерический период). Для Земли средняя величина периода Psd = 365.25636042 дней. Это значение периода Psd получено в результате обобщения данных наблюдения за всю историю астрономических наблюдений. При решении задачи взаимодействия планет, Солнца и Луны получено [9]-[10], что периоды обращения в разные эпохи колеблются вокруг величины Psd в небольших пределах. Поэтому этим отличием пренебрегаем.
Чтобы найти зависимость времени от полярного угла в решение t(r) задачи двух тел [8] подставляется зависимость r(o) согласно (2). Тогда время движения тела по орбите (см. рис. 1а) от точки перигелия PE до точки E с углом o рассчитывается по формуле
tfp = tfp при o и tfp = 2ta - tfp при < o 2, (4)
где
.
Здесь время движения от перигелия до афелия, как следует из (4), при o = равно
. (5)
При подстановке vp из (3) в (5) получаем ta = 0.5Psd. В решениях полной орбитальной задачи и в наблюдениях время движения от перигелия до афелия ta колеблется в небольших пределах вблизи значения 0.5Psd. В такой постановке мы этими колебаниями пренебрегаем. Отметим пренебрежение колебаниями Psd и ta как первое упрощение.
3. Геометрические характеристики инсоляции
Рассмотрим облучение Солнцем S точки M земной поверхности (см. рис. 2). Используем все основные обозначения величин, принятые в работе М. Миланковича [1]. Плоскость горизонта в т. М на небесной сфере 1 нанесена горизонтальным кругом HH. Перпендикуляр к плоскости HH пересекает плоскость небесной сферы 1 в точке зенита Z. Солнце S совершает вокруг Земли годовое движение по орбите, которая проектируется на небесную сферу в виде круга эклиптики EE. Движение происходит против стрелки часов с началом отсчета долготы Солнца в точке весеннего равноденствия . В этой точке Солнце находится в плоскости экватора, когда из южного полушария переходит в северное.
Суточное движение Солнца происходит следующим образом. Земля совместно с наблюдателем M, кругом горизонта HH и меридианом NZEAH вращается вокруг оси вращения NM, которую называют осью мира. Вращение происходит против часовой стрелки. Поэтому Солнце относительно Земли и, в частности, относительно круга горизонта HH по часовой стрелке перемещается по кругу SrMdSd параллельно экватору AA. В точке Sr оно восходит над горизонтом HH, в точке Md находится в полдень, а в точке Sd заходит за горизонт. Часовой угол Солнца будем отсчитывать от точки полудня Md. Дуги от восхода до полудня SrMd и от полудня до заката MdSd имеют одинаковую длину, которую обозначим как 0. Поэтому часовой угол дня изменяется в пределах -0 0. Длительность дня равна d = 20, а ночи n = 24 - d. Здесь мы использовали часовой угол в часах. Кроме того, далее он будет применяться в угловых единицах.
В представленном на рис. 2 положении наблюдателя M и Солнца S длительность дня d больше длительности ночи n. При нахождении Солнца S в точке E длительность дня будет наибольшая. Это точка летнего солнцестояния. При нахождении Солнца S в точке E наибольшей будет длительность ночи. Это точка зимнего солнцестояния. А при нахождении Солнца S в точках или его суточное перемещение будет происходить по кругу экватора AA. Этот круг пересекает круг горизонта HH по его диаметру, поэтому время нахождения Солнца над горизонтом и под ним одинаково, т.е. длительность дня равна длительности ночи.
Если наблюдатель в точке M на рис. 2 будет находиться на большей широте, т.е. дуга HN будет больше, то окружность SMd не пересечет круг горизонта HH. В этом случае для наблюдателя M наступит полярный день. При нахождении Солнца S в южном полушарии вблизи т. E круг его суточного движения также не пересечет линию горизонта HH. В этом примере широты для наблюдателя в точке M наступит полярная ночь.
Плоскости экватора AA и эклиптики EE изменяются в пространстве, вследствие чего точка весеннего равноденствия по часовой стрелке перемещается в год на 50.25641. Как уже отмечалось, годовое движение Солнца проходит по кругу эклиптики EE против часовой стрелки, что отражается изменением долготы , начиная от точки весеннего равноденствия . Поэтому время прохождения Солнцем двух последовательных точек весеннего равноденствия, т.е. тропический год Ptr = 365.24219879 дней, меньше сидерического года Psd. Для того, чтобы сезоны года не смещались по датам, наш календарь основан на тропическом годе. Поэтому инсоляцию Земли по сезонам года и по полугодиям определяют на основании длительности тропического года.
Рис. 2. Основные геометрические характеристики Солнца S при облучении точки M на земной поверхности: 1 - небесная сфера; НН - плоскость горизонта; N - северный полюс; AA - плоскость подвижного экватора; ЕЕ - плоскость подвижной эклиптики, а - угол между плоскостями AA и ЕЕ; Z - зенит точки M, а z = ZMS - зенитный угол Солнца; Дуга HN = - географическая широта точки M; щ = ZNS - часовой угол Солнца, отсчитываемый от полудня; = SB - склонение Солнца; = S - долгота Солнца.
4. Поток солнечного тепла
В пренебрежении поглощением солнечных лучей в межпланетном пространстве можем принять, что в единицу времени через сферические гелиоцентрические поверхности радиусами r и a протекает одинаковое количество тепла
4r2dWn(r)/dt = 4a2dWn(a)/dt, (6)
где dWn(r)/dt - поток лучистой энергии через единицу поверхности на расстоянии r от Солнца в единицу времени; a - средний радиус земной орбиты, т.е. ее большая полуось a = 0.5 (Rp + Ra), где Ra - радиус афелия.
Поток тепла на расстоянии от Солнца, равном a, называется солнечной постоянной J0
J0 = dWn(a)/dt. (7)
М. Миланкович использовал следующее ее значение J0 = 2 кал/(см2мин). Это значение в других единицах запишется так: J0 = 83.736 кДж/(м2мин) = 1395.6 Вт/м2.
В настоящее время значение солнечной постоянной известно с большей точностью. Радиация Солнца исследуется с 1907 г. в Физико-метеорологической обсерватории Давоса в Швейцарии [11], которая с 1971 г. функционирует под эгидой Всемирной метеорологической организации. В 1996 г. на основе анализа данных различных радиометров было принято в виде мирового радиационного эталона (World Radiometric Reference - WRR) значение J = 1366.784 Вт/м2. С 1978 г. измерение радиации Солнца проводится вне атмосферы - на спутниках. Как в наземных измерениях, так и в спутниковых поток солнечного тепла в каждой серии колеблется в пределах десятых долей процента. А во всех сериях солнечная радиация ступенчато изменяется в диапазоне от 1358 Вт/м2 до 1375 Вт/м2 [12]. В качестве космического абсолютного радиометрического эталона (SARR) принято значение солнечной радиации J = 1366.22 Вт/м2 [13]. С целью сохранения преемственности с расчетами предшествующих авторов ниже во всех расчетах мы используем значение, принятое М. Миланковичем J0 = 83.736 кДж/(м2мин).
Из выражения (6) при обозначении (7) поток солнечной радиации запишется dWn(r)/dt = J0a2/r2. Тогда поток солнечного тепла на расстоянии r от Солнца на площадку земной поверхности, перпендикулярную лучам Солнца, будет:
, (8)
где с = r/a - относительное расстояние до Солнца.
Линия зенита MZ (см. рис. 2) перпендикулярна рассматриваемой площадке. Если угол между Солнцем и зенитом равен z, то поток тепла в т. M будет