Новый алгоритм расчета инсоляции Земли

Размещено на http://www.allbest.ru//

Размещено на http://www.allbest.ru//

Новый алгоритм расчета инсоляции Земли

И.И. Смульский

Аннотация

Для исследования палеоклимата М. Миланковичем разработана теория инсоляции Земли. Расчеты инсоляции по упрощенным методам используются также в метеорологии, в строительстве и в др. областях. Новый алгоритм расчета инсоляции основан на результатах точного решения задачи 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 будет

. (9)

Зенитный угол z Солнца зависит от угловых координат т. М и Солнца S. В сферическом треугольнике NZS известен угол N = щ и две прилегающие стороны NZ = /2 - и NS =/2- -. Сторона ZS является дугой, которая измеряется центральным углом z (см. рис. 2), т.е. z = ZS. Этот угол и дуга определяется по теореме косинусов (стр. 30 [14])

cos z = cos(ZS) = cos(/2-)·cos(/2-) + sin(/2-)·sin(/2-)·cos . (10)

После подстановки (8) и (10) в (9) получаем закон облучения т. M в зависимости от ее широты и углов и положения Солнца

. (11)

В дальнейшем принимаем, что точки M земной поверхности на всех широтах имеют одинаковое относительное расстояние с до Солнца. В действительности это расстояние меняется в пределах величины радиуса Земли RE. При этом закладывается относительная погрешность в расчет инсоляции порядка RE/a = 4.310-5. Отметим это как упрощение 2. Третье упрощение связано с формой Земли. Зенитный угол Солнца z определяется по отношению к линии зенита ZM (см. рис. 2), которая проведена перпендикулярно поверхности Земли и проходит ее через центр. Для эллипсоидальной поверхности Земли это верно для точек M на полюсах и на экваторе. А на промежуточных широтах зенитный угол для эллипсоида будет немного отличаться от выражения (10). Четвертое упрощение связано с тем, что солнечная постоянная J0 принимается неизменной. Как уже отмечалось, имеются небольшие ее колебания во времени. Возможно, в дальнейшем, будет выявлен закон их изменения, который можно будет заложить в расчет инсоляции. На данном этапе изменение инсоляции, вызванное этими четырьмя обстоятельствами, не учитываются.

5. Суточная инсоляция

Зенитный угол Солнца z согласно (10), а, следовательно, и поток облучения (11) в зависимости от широты точки M и углов и щ положения Солнца может принимать значения больше нуля, равные нулю и отрицательные. Положительные значения инсоляции dW/dt соответствуют дневному времени, нулевые - времени восхода и захода Солнца, а отрицательные - ночному времени.

Как отмечалось ранее, отсчет зенитного угла Солнца z на рис. 2 начинается с полудня. Поэтому наибольшее его отрицательное значение (-90) будет соответствовать восходу, а наибольшее положительное значение (+90) - закату. Из условия cos z = 0 в формуле (10) или условия dW/dt = 0 в выражении (11) получаем зависимость часовых углов восхода и заката от широты точки M и склонения Солнца

щ0 = arcos(-tg tg ). (12)

Итак, дуги SrMd и MdSd на рис. 2 равны 0, а дневное время определяется часовыми углами Солнца -0 0. Поэтому удельное количество тепла, поступающее за световой день, т.е. за сутки на поверхность Земли в т. M, можно определить интегрированием (11) за дневное время. После перехода от времени t к часовому углу щ суточную инсоляцию получаем в виде:

. (13)

Здесь W - количество солнечной радиации, поступающей за сутки на 1 м2 земной поверхности в т. M.

Так как отсчет времени осуществляется по среднему движению Солнца, то времени суток в минутах = 2460 = 1440 мин. соответствует полный оборот = 2 Солнца вокруг Земли, т.е. в среднем dt/dщ = /2. После подстановки dt/dщ и потока солнечной радиации согласно (11) в (13) получаем удельную суточную инсоляцию в следующем виде

. (14)

При записи формулы (14) мы пренебрегли изменением расстояния с точки M до Солнца за время одного оборота Земли вокруг своей оси. За это время также изменяется склонение Солнца (см. рис. 2) за счет его движения по эклиптике ЕЕ. Пренебрежем также этим изменением. Отметим это как пятое упрощение. Его влияние значительно уменьшается, если значения с и относить к середине дня. При этих условиях в подынтегральном выражении (14) только cos щ зависит от часового угла щ. В результате интегрирования (14) получаем удельную суточную инсоляцию в виде:

, (15)

где часовой угол границы дня щ0 определяется выражением (12).

Склонение Солнца можно выразить из прямоугольного сферического треугольника SB (см. рис. 2), в котором угол B прямой, = , а сторона S определяет долготу Солнца , т.е. = S. По теореме синусов: sin /sin = sin / sin /2 находим угол склонения Солнца

= arcsin (sin sin ). (16)

В выражение (15) для суточной инсоляции входит относительное расстояние с центра Земли до Солнца. Согласно (2) расстояние r = сa до Солнца определяется углом o, который отсчитывается от перигея.

Угол перигелия Земли p отсчитывается от восходящего узла ее орбиты . Так как перигей Солнца PS (см. рис. 1б) расположен диаметрально противоположно перигелию Земли PE, то угол перигея Солнца будет

= p + . (17)

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

o = - = - p - . (18)

При вычислении полярного угла o необходимо учитывать ряд моментов, связанных с цикличностью изменения углов от 0 до 2. Угол перигелия p от подвижного узла в наших решениях орбитальной и вращательной задач представлен в виде нарастающего по величине ряда чисел. Поэтому угол перигея Солнца необходимо превратить в ряд чисел с циклическим изменением углов от 0 до 2. После выражения долготы в (18) через угол o последний также необходимо привести к диапазону углов от 0 до 2. Эти операции необходимо выполнять при исследовании эволюции инсоляции за большие интервалы времени.

Итак, в расстояние r входит согласно (2) и (18) долгота . Поэтому суточная инсоляция W по выражению (15) зависит от долготы Солнца и широты т. М на поверхности Земли. Долгота является короткопериодическим параметром, цикл изменения которого 1 год. А из долгопериодических параметров, которые входят в суточную инсоляцию W согласно (15), являются эксцентриситет e или параметр траектории б1, угол наклона и угол перигелия p.

6. Изменение долготы Солнца по дням года

Солнце неравномерно движется по орбите, например, в перигее быстрее, а в апогее - медленнее. Поэтому для расчета инсоляции по дням года необходимо определить долготы Солнца, соответствующие этим дням, т.е. необходимо определить зависимость долготы от времени. Время движения tfp от перигея до точки на орбите с углом o определяется выражением (4). Согласно (18) в угол o входит долгота Солнца . Тогда используя зависимость (4) для времени движения tfp по орбите, можно определить ряд долгот, соответствующим дням года. Эта задача решается методом последовательных приближений. Вначале задается равномерный ряд долгот

0i = 0i-1 + 0, i = 1, 2 … 365, (19)

где 00 = 0; 0 = 2 /365.

По формуле (18) определяется соответствующий ряд углов o,i положения Солнца от перигея. Затем по формуле (4) определяется ряд моментов времени tfp,1,i, которые соответствуют долготам (19).

Полученные дискретные величины tfp,1,i/ed, где ed = 243600 - количество секунд в сутках, отличаются от целых чисел дней в году Tdi = 1, 2, … 365. Поэтому рассчитывается уточненный ряд долгот

, (20)

где 10 = 0.

По новым значениям 1i рассчитываются углы o,1,i согласно (18), а затем по формуле (4) - уточненные моменты времени tfp,2,i. Промежутки времени tfp,2,i+1 - tfp,2,i между соседними долготами должны быть равны продолжительности суток ed. Поэтому рассчитывается относительная величина погрешности времени между соседними долготами 1i в виде

, (21)

а также средняя величина погрешности

. (22)

Величина tfp,1,m представляет среднюю относительную погрешность по времени, которую дает уточненный ряд долгот.

Аналогичным образом, начиная с вычисления уточненного ряда долгот по формуле (20), могут быть повторены эти операции для получения новых уточнений 2i, 3i, 4i и т. д. Вычисления показывают, что после первого уточнения величина погрешности tfp,1,m = 7.0910-4, после второго tfp,2,m = 3.4310-5 и после третьего tfp,3,m = 7.1510-7. Величина относительной погрешности порядка 10-5 ч 10-7 вполне достаточна для дальнейших вычислений инсоляции. При необходимости ее уменьшения можно выполнить еще несколько уточнений по вышеизложенному алгоритму.

Далее рассчитывается погрешность определения долготы 365 в последний, 365-й день года. Так как полный оборот Солнца вокруг Земли на угол 2• происходит за сидерический год Psd, то избыток дня сверх 365 дней dsd = Psd - 365. По рассчитанной долготе 365-го дня 365 избыток дня будет

, (23)

где 1d - 1 день, а dsdc вычисляется в днях.

Разница между рассчитанным избытком дня и фактическим dsd представляет накопленную погрешность долготы 365 в днях, т.е. 365 = dsdc - dsd. А относительная погрешность долготы будет 365 = 365/dsd. Получены следующие значения этих величин после каждого из трех уточнений: 365 = -0.2 дня; -161 сек; -11сек; 365 = -0.8; -0.0073; -4.9710-4. Как видно из этих данных, каждое уточнение уменьшает накопленную погрешность более чем на порядок.

7. Изменение суточной инсоляции по дням года

Перепишем выражение (15) для суточной инсоляции, выразив в нем склонение через долготу с помощью выражения (16). Далее под подразумеваем уточненную согласно алгоритму (20) - (22) долготу Солнца. Тогда согласно (15) суточная инсоляция будет

, (24)

где часовой угол восходов и закатов согласно (12) запишется так:

0 = arcos{-tg[arcsin(sin sin )]tg }. (25)

Относительное расстояние до Солнца в выражении (24) с = r/a, где r определяется выражением (2). Так как согласно задаче двух тел [7] - [8] Rp/a = (2б1 + 1)/б1, то с учетом выражения (18) для полярного угла o относительное расстояние запишется в виде

, (26)

где, как уже упоминалось ранее, параметр траектории б1 однозначно связан с эксцентриситетом: б1 = -1/(1 + e).

Итак, выражениями (24) - (26) полностью определяется удельная суточная инсоляция Земли на любой ее широте и в любой день года, соответствующих долготе . В ряде уточненных по дням долгот i, где i = 1…365, долготы начинают отсчет с момента весеннего равноденствия, т.е. с 22 марта и заканчивают 21 марта следующего года. Отметим, что здесь имеется некоторое отличие рассчитанных долгот от тех, которые соответствуют календарным дням. Это обусловлено тем, что момент равноденствия не совпадает с началом дня, календарный год не совпадает с тропическим годом Ptr и в високосный год скачком изменяется на один день. Поэтому далее будет рассмотрен расчет точных дат.

Рис. 3. Изменение суточной инсоляции W в МДж/м2 по дням года на разных широтах от + 90 до - 90. Время Td в днях отсчитывается от момента весеннего равноденствия 22 марта (эпоха 1950 г.). На графиках числами приведены широты в градусах, максимальная Wmax и минимальная Wmin суточные инсоляции на данной широте в МДж/м2: NP - северный полюс; Eq - экватор; SP - южный полюс.

Согласно (24) долгота определяет изменение суточной инсоляции W в течение года, а входящие в выражения (24) - (26) параметры , fp, 1 или e позволяют рассчитать суточную инсоляцию в любую эпоху и исследовать ее изменение. На рис. 3 представлено изменение суточной инсоляции во времени Td в днях на разных широтах от = 90 на северном полюсе до = - 90 на южном полюсе. Счет в днях происходит от момента весеннего равноденствия 22 марта. В средних широтах северного полушария суточная инсоляция W в летнее время достигает максимального значения Wmax и в зимнее время - минимального значения Wmin. Их значения для каждой широты приведены на графиках. В высоких широтах минимальная инсоляция Wmin = 0, т.е. наступает полярная ночь. На экваторе ( = 0) имеется два минимума и два максимума. Максимумы соответствуют моментам прохождения Солнца через экватор, т.е. в период равноденствий, а минимумы - в период солнцестояний.

Наибольшая суточная инсоляция, как видно из графиков, имеет место на полюсах во время полярного дня, например, на Северном полюсе Wmax = 46.4 МДж/м2, а на южном Wmax = 49.6 МДж/м2. На экваторе максимальная суточная инсоляция значительно меньше и равна Wmax = 38.8 МДж/м2. В то же время минимальная инсоляция Wmin = 34.1 МДж/м2 немного отличается от максимальной. Однако, минимальная инсоляция на экваторе, как видно из рис. 3, является наибольшей из минимальных на всех широтах.

Из сравнения Wmax и Wmin видно, что летние инсоляции в южном полушарии больше, а зимние меньше, чем соответствующие инсоляции в северном полушарии. Это обусловлено тем, что в северном полушарии перигелий приходится на зимнее время.

8. Инсоляция за год

Как уже отмечалось, тропический год длится больше чем 365 дней. Определяется избыток дня в тропическом году dtr = Ptr - 365. Тогда полная инсоляция QT за год может быть рассчитана суммированием суточных инсоляций так:

, (27)

где Wi - рассчитанная согласно (24) суточная инсоляция в зависимости от долготы Солнца i, а W365 - суточная инсоляция последнего дня.

Основное отличие нашего расчета инсоляции от метода М. Миланковича заключается в расчете инсоляции за промежутки времени, которые включают ряд дней, например, за сезоны, полугодия и за год. М. Миланкович находит аналитические выражения для инсоляции за необходимые интервалы времени. С этой целью ему приходится выполнять ряд математических преобразований и приближений. Суть их заключается в том, чтобы найти аналитическую зависимость времени от долготы, а затем поток радиации dW/dt, представленный выражением (11), проинтегрировать по необходимому промежутку времени.

В отличие от М. Миланковича мы находим интегрированием потока dW/dt только суточную инсоляцию. Затем для каждого дня на основе решения (4) задачи двух тел мы определяем долготу. Так как отсчет долготы и дней начинается с точки весеннего равноденствия, то любой период времени можно задать по количеству дней. При этом для каждого конца дня суточная инсоляция определяется выражением (24). Поэтому инсоляция за любой отрезок времени определяется суммированием суточных инсоляций. Именно так и была определена в выражении (27) годовая инсоляция.

При определении астрономических сезонов и полугодий для северного полушария так, как определял их Миланкович [1], номера дней и их количество рассчитываются по долготе : весна - = 0 90; лето - = 90 180; осень - = 180 270; зима - = 270 360; летнее полугодие - = 0 180; зимнее полугодие - = 180 360. Суммированием суточных инсоляций (24) по рассчитанным номерам дней и их количествам определяются инсоляции за сезоны и полугодия. Следует отметить, что вышеупомянутые астрономические периоды определяются значением долготы , поэтому количества дней будут не целые числа, т.е. будут включать в себя части дня.

Этот путь оказался возможным, потому что мы воспользовались решениями для времени задачи двух тел (см. п. 5.2 [7]). М. Миланкович связь долготы Солнца со временем t устанавливает на основании дифференциального уравнения (см. (34) [1]):

. (28)

Уравнение (28) М. Миланкович интегрирует приближенными аналитическими методами. Как уже отмечалось, в работе [7] приведены точные аналитические решения задачи двух тел для четырех возможных случаев движения: по эллипсу, параболе, гиперболе и для случая прямолинейного движения. Формулой (4) представлена зависимость времени от полярного угла o или долготы для эллиптического движения. Таким образом, формулами (4) и (28) выражена основная математическая разница нашего метода от метода М. Миланковича.

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

9. Инсоляция за калорическое полугодие

В астрономической теории палеоклимата важную роль играет инсоляция за летнее и зимнее полугодия. Однако в связи с движением перигелия Земли относительно восходящего узла длина астрономических летнего и зимнего полугодий изменяется. Поэтому были введены калорические полугодия одинаковой длины. Согласно М. Миланковичу [1], летнее калорическое полугодие определяется так, чтобы инсоляция W любого его дня была больше инсоляции любого дня зимнего полугодия.

Как уже отмечалось, год Ptr определяется количеством дней между двумя прохождениями Солнца точки весеннего равноденствия . Величина Ptr почти на четверть дня превышает число 365. Для нахождения калорических полугодий рассматривается инсоляция за два равных промежутка по 182 дня (см. рис. 4). Их сумма меньше Ptr, но это отличие неважно для нахождения полугодий.

Мы разработали два способа расчета калорических полугодий. При первом способе определяется наибольшая суточная инсоляция Wmax в году (см. рис. 4) и соответствующий номер дня Kmax. Задается начальный номер для летней инсоляции

Kstp = Kmax - 91 - Kstp, (29)

где Kstp = 2 при широте || 45; Kstp = 5 при 25 < || < 45; Kstp = 70 при || 25.

Величина смещения Kstp начального номера дня была установлена в результате расчетов летней инсоляции на разных широтах.

Рис. 4. К определению калорических полугодий: 1 - летнее калорическое полугодие; 2 - зимнее калорическое полугодие; W - суточная инсоляция, Td - время в днях.

Чтобы упростить вычисление летней инсоляции, создается массив инсоляций Vi2, i2 = 1, 2… 730, который включает два массива Wi инсоляций за год:

Vi = Wi; Vi+365 = Wi. (30)

Расширенный массив суточных инсоляций Vi2 необходим, чтобы выбранный интервал для полугодия (см. рис. 4) не разрывался на части.

Рассматривается несколько вариантов летних полугодий с начальными индексами, изменяющимися в диапазоне jS = Kstp…Kstp + 2Kstp. Рассчитываются суммы инсоляций за эти полугодия с длительностью 182 дня:

, (31)

где j3 = 0, 1…181; здесь js = js и j3 = j3, т.е. по-разному записанные индексы - идентичны.

Максимальная величина ряда величин VSjS с индексом jp0 принимается как основное слагаемое инсоляции за летнее калорическое полугодие. А полная инсоляция определяется выражением:

QS = max(VSjS) + 0.5(1 + dtr)Vjp0+182. (32)

В формуле (32) дополнительное слагаемое 0.5(1 + dtr)Vjp0+182 обусловлено прибавлением инсоляции за половину разности тропического года и рассматриваемого периода 2182 = 364 дня.

Итак, при этом способе летнее калорическое полугодие определяется максимальной инсоляцией за 182 дня.

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

При известной инсоляции летнего полугодия рассчитывается инсоляция за зимнее калорическое полугодие

QW = QT - QS. (33)

Рис. 5. Распределение по широте Земли удельного тепла в ГДж/м2 в современную эпоху (1950 г.): QS - за летнее калорическое полугодие; QW - за зимнее калорическое полугодие; QT - за весь год: на графике величина QT уменьшена в два раза; > 0 - северное полушарие; < 0 - южное полушарие; Mil - расчеты М. Миланковича для эпохи 1800 г.

На рис. 5 приведены изменения по широте летней QS, зимней QW и уменьшенной в 2 раза годовой QT инсоляции в современную эпоху. Инсоляция за год QT от экватора монотонно убывает к полюсам. При этом на полюсах она в 2.4 раза меньше, чем на экваторе. Зимняя инсоляция QW вблизи экватора имеет максимальное значение, а на полюсах стремится к нулю. Летняя инсоляция QS имеет максимальные значения вблизи тропиков ( = = 23.4) и наименьшие значения на полюсах. При этом в современную эпоху летняя инсоляция QS на Северном полюсе, как видно из графика на рис. 5, меньше, чем на Южном. Также летняя инсоляция на северном тропике в 1.04 раза меньше, чем на южном.

10. Инсоляция в эквивалентных широтах

Состояние климата в данном географическом месте, кроме инсоляции, зависит от других факторов. Чтобы иметь возможность по величине инсоляции определять изменения климата в разные эпохи, М. Миланкович [1] предложил выражать инсоляцию Земли в эквивалентных широтах. Если в эпоху T летняя инсоляция QS на широте была такая, как в современную эпоху T0 летом на широте 0, то инсоляция в эквивалентных широтах будет I = 0. По величине такой инсоляции I в пункте с географической широтой можно сказать, насколько произошло потепление или похолодание климата.

Представленная на рис. 5 летняя инсоляция QS в эпоху 1950 г. была принята за стандартную. Обозначим ее значения на разных широтах как QS,n,i3 и n,i3, где n - значок стандартной широты, i3 - индекс конкретного значения широты. Итак, имеется функциональная зависимость для стандартной летней инсоляции QS,n,i3(n,i3). Для определения стандартной широты, на которой инсоляция QS в эпоху T будет равна QS,n в эпоху T0, необходимо использовать обратную зависимость n,i3(QS,n,i3). Однако, т.к. инсоляция QS,n, как видно из графика на рис. 5, имеет несколько максимумов, то зависимость n,i3(QS,n,i3) неоднозначна. Поэтому необходимо ее разбить на монотонные участки. Например, для северного полушария рассматривается участок от минимального значения QS,n,min при = 87.5 и с индексом imin до максимального значения QS,n,max при = и с индексом imax.

Так как зависимость n,i3 (QS,n,i3) для стандартной инсоляции дискретна, то инсоляция в эквивалентных широтах при любом значении удельной инсоляции QS находится с помощью параболической интерполяции

I = AiQ-1Q2S + BiQ-1QS + CiQ-1; (34)

где iQ - индекс начала параболического участка аппроксимации.

Входящие в формулу (34) коэффициенты Ai3, Bi3 и Ci3 для любого индекса i3 рассчитываются так:

; (35)

; (36)

, (37)

где i3 = imx … imin-2.

Индекс iQ в выражении (34) подбирается так, чтобы рассматриваемое в эпоху T значение инсоляции QS находилось между значениями QS,n,iQ-1 и QS,n,iQ.

Итак, по известной удельной летней инсоляции QS в кДж/м2 в эпоху T летняя инсоляция I в эквивалентных широтах определяется по выражениям (34) - (37) в зависимости от стандартной летней инсоляции QS,n(n) в эпоху 1950 г. Ее вид представлен графиком QS на рис. 5. Для летней инсоляции в северном полушарии на монотонном участке 2587.5 в табл. 1 представлены распределения стандартной инсоляции по широте и коэффициенты A, B, C.

Таблица 1. Распределение стандартной летней инсоляции QS,n в кДж/м2 по широте северного полушария (эпоха 1950 г.) и коэффициенты A, B, C интерполяционной параболы при J0 = 83.736 кДж/(м2мин).

i3	n	QS,n	A	B	C
1	25	7198511	-1.586825E-08	0.22787858	-818091.62
2	27.5	7193482.5	-1.305063E-09	0.018607205	-66291.062
3	30	7176406.7	-3.5779392E-10	0.0050386608	-17702.802
4	32.5	7147447	-1.4359084E-10	0.0019853727	-6822.3648
5	35	7106727.9	-7.3370248E-11	0.00099094642	-3301.7795
6	37.5	7054735.1	-4.192297E-11	0.00054922519	-1750.6619
7	40	6991670.4	-2.6076644E-11	0.00032880899	-984.20762
8	42.5	6917939.3	-1.7191908E-11	0.00020662631	-564.15994
9	45	6834036.5	-1.1777458E-11	0.00013312729	-314.74172
10	47.5	6740568.2	-8.226025E-12	8.5613197E-5	-155.83001
11	50	6638283.5	-5.7080254E-12	5.2460194E-5	-46.711181
12	52.5	6528122.1	-3.7293019E-12	2.6856673E-5	36.105695
13	55	6411290.4	-1.9029696E-12	3.6610174E-6	109.74904
14	57.5	6289386.7	2.8443607E-13	-2.3580964E-5	194.55854
15	60	6164627.4	4.2047288E-12	-7.1427882E-5	340.53552
16	62.5	6040307.7	2.0009759E-11	-0.00026049237	905.89164
17	65	5921990.6	6.0889024E-11	-0.00074065184	2315.7568
18	67.5	5823804.4	8.5878878E-11	-0.0010299741	3153.1396
19	70	5753783.9	1.3723453E-10	-0.0016181826	4837.383
20	72.5	5699845.5	2.3455868E-10	-0.0027235597	7975.9714
21	75	5657840	4.9616553E-10	-0.0056753239	16302.244
22	77.5	5625366.5	1.1992349E-09	-0.013568618	38456.459
23	80	5601538.9	4.3162616E-09	-0.048438005	135975
24	82.5	5585207.6	4.5065099E-08	-0.50323412	1404964.3
25	85	5575751.7	4.5065099E-08	-0.50323412	1404964.3
26	87.5	5572727.9	4.5065099E-08	-0.50323412	1404964.3

Аналогичным образом в эквивалентных широтах могут быть рассчитаны годовая и зимняя инсоляции.

Рассмотренный алгоритм расчета суточной W, годовой QT, летней QS, зимней QW инсоляций, а также инсоляции в эквивалентных широтах I был реализован в среде MathCad в виде программы Insl2bd.mcd (см. Приложение). Для расчета эволюции инсоляции дополнительно задаются данные об эволюции орбитального и вращательного движения Земли в виде ряда величин: T, e, p и . Они считываются из файла с именем, например (см. Приложение), INSO_LA2004.txt или OrAl1c_8.prn. Для расчета летней инсоляции I в эквивалентных широтах в файле InsCvSNJ.prn задается ряд величин: i, QS,n, n, A, B и C из табл. 1. В программе используется также ряд не упомянутых здесь условий, которые обеспечивают работоспособность алгоритма при возможных сочетаниях используемых параметров.

11. Проверка достоверности алгоритма

М. Миланкович в табл. 14 [1] приводит распределение инсоляции по широте за калорические полугодия в эпоху 1800 г. Эти результаты приведены в канонических единицах, введенных М. Миланковичем. С помощью коэффициента KKn = 10-5PtredJ0/60 = 440 их можно пересчитать в кДж/м2. По вышерассмотренному алгоритму для эпохи 1800 г. были рассчитаны летние и зимние инсоляции для северного и южного полушарий. Они практически совпадают с расчетами М. Миланковича. Отличия не превышают величины 0.1%, и, как правило, обусловлены единицами последнего разряда чисел, которые приведены в табл. 14 [1].

На графиках рис. 5, рассчитанных для эпохи 1950 г. по новому алгоритму летней QS и зимней инсоляции QW, ромбиками нанесены результаты расчета М. Миланковича. Несмотря на то, что эти результаты относятся к разным эпохам 1950 г. и 1800 г., они практически не различаются на графиках. Относительное отличие этих результатов для разных эпох находится на уровне сотых долей процента и достигает наибольшего значения 0.6% для широты = 80.

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

Были выполнены расчеты по проверке алгоритма расчета инсоляции в эквивалентных широтах. В работе Шараф и Будниковой [15] инсоляция в эквивалентных широтах приведена в графическом виде. Однако в этой работе не приведены числовые данные об эволюции орбитального и вращательного движения Земли. Эти данные приведены в работе Ляскара и др. [16]. В этом случае в программе (см. Приложение) использовался файл INSO_LA2004.txt с параметрами об эволюции орбитального и вращательного движения Земли за 21 млн. лет. Эти данные доступны на сайте [17]. Основываясь на данных этого файла, мы рассчитали инсоляцию I в эквивалентных широтах за 200 тыс. лет в будущее. На рис. 6 наши расчеты сопоставлены с расчетами Шараф и Будниковой. Как видно, они практически совпадают. Небольшие отличия в экстремальных точках могут быть обусловлены двумя обстоятельствами: различием исходных данных Ляскара и др. [16] и Шараф и Будниковой [15], а также графическим характером результатов последних авторов.



Рис. 6. Сравнение инсоляций в эквивалентных широтах, рассчитанных по двум методикам: новой (1) и Миланковича (2). Инсоляция 1 рассчитана по параметрам эволюции орбитального и вращательного движения Ляскара и др. [16]- [17], 2 - графические результаты расчетов Шараф и Будниковой [15].

Итак, все выполненные проверки подтвердили идентичность результатов нового алгоритма расчета инсоляции результатам алгоритма М. Миланковича.

12. Динамика инсоляции в современную эпоху

Как уже отмечалось вначале, расчеты инсоляции Земли в современную эпоху используются в разных областях. По изложенному алгоритму с помощью представленной в Приложении программы Ins12bd.mcd инсоляцию можно рассчитать, задавшись параметрами орбиты Земли и плоскости экватора для эпохи 2000.0 г. [18]

e0 = 0.01670863; 0 = 0.4090926; p0 = 1.7965956; JDp0 = 2451548. 541, (38)

где JDp0 - определенный по нашим данным юлианский день прохождения среднего Солнца через перигелий; ему соответствует календарная дата 05.01.2000 г. и время 1:30:00 по Гринвичу.

Параметры e0, 0 и p0 являются осредненными элементами орбиты Земли и положениями ее оси, которыми в астрономии характеризуется движение “среднего” Солнца и “средней” оси вращения Земли. Эпохе 2000.0 г. соответствует юлианский день JD0 = 2451545 и календарная дата 01.01.2000 и время 12:00:00 по Гринвичу.

По алгоритму в п. 6 программы, описанному формулой (15), рассчитывается распределение суточной инсоляции W по широте и дням года. Счет дней года начинается с точки весеннего равноденствия, которой соответствует номер дня j = 0. А номер дня jp, соответствующий моменту прохождения Солнцем перигелия, определяется долготой перигелия p. Юлианская дата каждого дня JDd определяется, если в этом есть необходимость, количеством дней kd, прошедших после юлианской даты перигелия JDp,

JDd = JDp + kd. (39)

В случае, если день предшествует дате перигелия, kd < 0. Календарный день года определяется по переводным таблицам соответствия юлианских дней JDd определенной дате (см., например, [19]). Промежуток времени движения Солнца от перигелия до момента равноденствия определяется по формуле (4) при долготе = 0. В п. 4 программы количество времени от перигелия до дня равноденствия выражено параметром tpd, т.е. kd = tpd = 76.84. Тогда определенный согласно (39) юлианский день равноденствия в эпоху JD0 = 2451625.6. Ему соответствует календарная дата 22 марта 2000 г. и время 02:24 по Гринвичу.

По рассмотренному в п. 5 программы ряду долгот i рассчитываются дни, которые соответствуют сезонам года, границы которых уже отмечались ранее: весна 0 < < 90, лето 90 < < 180, осень 180 < < 270 и зима 270 < < 360 или соответствуют астрономическому лету 0 < < 180 или астрономической зиме 180 < < 360. Суммированием суточных инсоляций определяются инсоляции за эти периоды времени.

По п.8 программы в Приложении рассчитывается инсоляция за год и за калорические полугодия. Как уже отмечалось ранее, распределения инсоляции по широте за 1950 г. и 1800 г. отличаются преимущественно в сотых долях процента. Поэтому, когда не важны короткопериодические колебания инсоляции, то рассчитанную с неизменными параметрами e0, 0, p0 и JDp0 инсоляцию можно использовать на интервале порядка 100 лет.

Для исследования причин короткопериодических изменений в метеорологии, медицине, биологии и др. областях, которые обусловлены инсоляцией, необходимы более точные ее расчеты. На рис. 7 представлена динамика инсоляции за 100 лет. Она основана на результатах решения орбитальной задачи с помощью системы Galactica [20] и задачи о вращении Земли [21] при одновременном воздействии на нее планет, Солнца и Луны. Динамика орбитального движения рассматривается относительно неподвижного экватора Земли, а вращательное движение Земли - относительно неподвижной орбиты Земли. Затем полученные параметры пересчитываются в относительные величины и p, определяющие положение подвижной орбиты относительно подвижного экватора. Поэтому все колебания орбиты Земли и плоскости экватора включаются в колебания углов , p и эксцентриситета e.

Таблица 2. Динамика параметров орбиты Земли и ее оси вращения за 100 лет от 30.12.1949 г. с JD49 = 2433280.5; T - в сидерических столетиях; JD - юлианский день эпохи; e - эксцентриситет орбиты; - угол наклона между подвижными орбитой и экватором, p - угол между перигелием орбиты и ее подвижным восходящим узлом на плоскости подвижного экватора. JDP - юлианский день прохождения Землей перигелия.

T, столетия	JD	e		p	JDP
0	2433280.500	0.01674220136	0.4091616307	1.787141716	2433284.737
0.05	2435106.782	0.01672566559	0.4091879678	1.781436815	2435109.814
0.1	2436933.064	0.01676227585	0.4092257260	1.788220450	2436936.972
0.15	2438759.345	0.01673867665	0.4091642935	1.788629228	2438764.131
0.2	2440585.627	0.01673105105	0.4091133198	1.800851582	2440590.997
0.25	2442411.909	0.01671231153	0.4091644133	1.797431716	2442415.087
0.3	2444238.191	0.01672990037	0.4091759892	1.782235510	2444242.099
0.35	2446064.473	0.01676031787	0.4091001848	1.793227470	2446069.258
0.4	2447890.754	0.01669882216	0.4090794136	1.794799473	2447894.371
0.45	2449717.036	0.01673052444	0.4091332328	1.787641468	2449721.529
0.5	2451543.318	0.01672283742	0.4091164493	1.808234553	2451548.542
0.55	2453369.600	0.01671511876	0.4090402810	1.805005623	2453372.485
0.6	2455195.882	0.01669494496	0.4090512378	1.803246385	2455199.498
0.65	2457022.163	0.01670168259	0.4091002365	1.801029825	2457026.802
0.7	2458848.445	0.01675111524	0.4090517953	1.807771818	2458853.815
0.75	2460674.727	0.01669781681	0.4089876396	1.791918327	2460678.928
0.8	2462501.009	0.01671860947	0.4090259849	1.808659189	2462506.086
0.85	2464327.291	0.01669587424	0.4090541814	1.822946059	2464332.806
0.9	2466153.572	0.01669390754	0.4089876672	1.799729397	2466157.042
0.95	2467979.854	0.01671288342	0.4089473460	1.804742513	2467984.055
1	2469806.136	0.01667219608	0.4089995992	1.808594771	2469809.460

Параметры орбит планет определялись с интервалом 1 год, поэтому в эти расчеты вошли их короткопериодические изменения с периодами несколько лет. Решения для вращательного движения содержат суточные, полумесячные, полугодовые колебания, а также колебания с периодом 18.6 лет. Чем больше период, тем больше амплитуда этих колебаний.

По полученным параметрам , p и e с помощью программы в Приложении была рассчитана динамика инсоляции. В табл. 2 эти параметры приведены с интервалом 5 лет. Кроме того в табл. 2 приведены юлианские дни эпохи JD и моментов прохождения Землей перигелия JDp. В программе Insl2bd.mcd все эти параметры считываются из файла OrAl1c_8.prn. Следует отметить, что файлы данных OrAl1c_8.prn и InsCvSNJ.prn, а также файл программы Insl2bd.mcd доступны по адресу [22].

На графике годовой инсоляции QT для широты 80 северного полушария (см. рис. 7) хорошо видны колебания с периодом 18.6 лет. Они обусловлены воздействием прецессионного движения Луны на вращательное движение Земли: ось Земли колеблется с этим периодом. Амплитуда колебания инсоляции в этом случае составляет 532 кДж/м2. В годовой инсоляции QT имеются еще полугодовые колебания с амплитудой 31 кДж/м2, которых на графике из-за их малой величины не видно.

В инсоляции QS за летнее калорическое полугодие на широте 80 наряду с колебаниями с периодом 18.6 лет видны колебания более короткого периода. В большей мере они просматриваются в зимней инсоляции QW. Эти колебания инсоляции обусловлены колебаниями угла перигелия p орбиты Земли с периодами 2.75 и 11.86 лет, а также ее эксцентриситета e с периодами 3.98 года и 11.86 лет. Если для летней инсоляции колебания перигелия незначительны, то для зимней инсоляции они даже превосходят колебания с периодом 18.6 лет.

На широте 80 все три инсоляции: QT, QS и QW имеют общий ниспадающий тренд. Он обусловлен долгопериодическим уменьшением угла наклона между плоскостями экватора и орбиты.

Как видно из рис. 7 с изменением широты = 65, 45, 25 и 0 динамика инсоляции качественно изменяется. На широте 65 в летней инсоляции QS сильнее проявляются колебания перигелия и эксцентриситета. На широте 45 эти колебания становятся главными для летней QS и зимней QW инсоляций. При этом летняя и зимняя инсоляции изменяются в противофазе друг к другу. Кроме того в летней и зимней инсоляции отсутствует нисходящий тренд, который был на предыдущих широтах.

На широтах 25 и 0 появляется восходящий тренд для годовой инсоляции QT, а также летней - на широте 25 и зимней - на широте 0. На этих широтах летняя и зимняя инсоляции также колеблются в противофазе. В их колебаниях проявляются периоды изменения перигелия 2.75 года и эксцентриситета 3.98 года, а период 18.6 лет колебания угла наклона отражается только в изменении годовой инсоляции QT.

Рис. 7. Динамика инсоляции за 100 лет, начиная с 30.12.1949 г., на разных широтах северного полушария. Обозначения см. рис. 5.

Представленная на рис. 7 динамика инсоляции представляет широкий спектр ее изменения во времени и по поверхности Земли. Эти изменения получены в рамках строго детерминированного рассмотрения взаимодействия тел Солнечной системы и геометрических характеристик освещения Земли Солнцем. Детерминированный подход дает разные частоты, разные тренды, которые по-разному проявляются для соответствующих инсоляций QT, QS и QW. Если мы обратимся к какому-либо процессу, происходящему в природе, например, к изменению погоды или климата, то также обнаруживаем широкий спектр изменения. Так как одни причины изменения неизвестны, а механизм действия других причин не установлен, то исследователь вынужден использовать статистические методы для установления вероятностных связей. Поэтому установленное знание об этом процессе также носит вероятностный характер, который допускает появление непредусмотренных ситуаций. Тем не менее, со временем ряд причин удается установить, и казавшиеся ранее случайные процессы становятся полностью объяснимы.