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

Этот путь оказался возможным, потому что мы воспользовались решениями для времени задачи двух тел (см. п. 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мин).

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

Рассмотренный алгоритм расчета суточной 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], а также графическим характером результатов последних авторов.