Основные положения и новые результаты астрономической теории изменения климата

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

Российская академия наук

Сибирское отделение

Институт криосферы Земли

УДК 551.583 +521.1

Основные положения и новые результаты астрономической теории изменения климата

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

Тюмень, 2014

Аннотация

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

Введение

В истории Земли известны многократные чередования потеплений и похолоданий [1]-[2]. Десять тыс. лет назад сошел ледник, который покрывал северную и среднюю полосу Европы. С другой стороны, в современных, практически безжизненных полярных районах, некогда была обильная растительность, росли деревья, и обитал обширный животный мир: мамонты, шерстистые носороги, бизоны, лошади и другие животные. Какая причина таких колебаний климата на Земле?

Еще в 19 веке Луи Агассис [3], Ж. Адемар [4], Джеймс Кроль [5] и др. вынашивали идеи о том, что изменение параметров орбиты Земли и оси ее вращения может приводить к изменению количества тепла Солнца, которое поступает на поверхность Земли на разных ее широтах. К концу 19 века достижения в небесной механики стали таковы, что можно было рассчитать изменения орбитальных и вращательных характеристик Земли, и в начале 20ого века Милутин Миланкович [6] завершил создание Астрономической теории ледниковых периодов, которая по существу является астрономической теорией изменения климата. В этой теории на основании трех параметров: эксцентриситета e орбиты Земли, углового положения перигелия цpг и угла наклона е плоскости орбиты к плоскости экватора Земли рассчитывается инсоляция Земли на разных ее широтах.

Процесс облучения Земли Солнцем и величину получаемого ее поверхностью тепла принято называть инсоляцией (in-sol, in - приставка к глаголам, означающим приносить, приводить; solis - солнце). В то же время необходимо помнить, что существуют разные стороны этого процесса, которые обозначены другими терминами: иррадиация, освещенность, облучение и т.д. Инсоляция может быть мгновенная, за сутки, за какой-то сезон или полугодие, за год.

Так как параметры e, цpг и е изменяются и колеблются с периодами десятки тысяч лет, то от эпохи к эпохе может быть рассчитаны инсоляции, например, за летнее полугодие на широте 65°, и по изменению такой инсоляции можно судить об изменении климата. Однако амплитуда колебаний эксцентриситета e и угла наклона е в теории М. Миланковича была небольшой. Например, угол е колебался в диапазоне ± 1°. Эти колебания могли дать колебания температуры также порядка 12°C. Поэтому, как во времена М. Миланковича, так и в последнее время астрономическая теория изменения климата вызывала сомнение ученых [7], занимавшихся изучением палеоклимата.

После М. Миланковича, с периодичностью нескольких десятков лет, его исследования повторяли разные группы ученых [8]-[13]. Они уточняли эволюцию параметров e, цpг, е и первоначальный расчет за 600 тыс. прошедших лет [6] продолжали за более длительное время, например, до 30 млн. лет [10]. Однако главный результат: колебания эксцентриситета e и угла наклона е оставались по-прежнему небольшими.

Во второй половине 20ого века выяснилось, что в морских осадках обнаруживаются колебания изотопа кислорода O18. Были проведены обширные исследования по всему мировому океану. Результаты исследований обобщены в виде стандартных зависимостей относительного содержания изотопа кислорода дO18 по толщине осадков, которая соотносится со временем T [14]. Эти зависимости названы по начальным буквам названия программы CLIMAP [15] или имен авторов LR-4 [16]. Предполагается, что более легкий изотоп O16 в составе воды испаряется и накапливается в ледниках. Поэтому было принято, что избыток изотопа кислорода O18, т.е. превышение дO18 над средним уровнем, пропорционален объему льда, накопленному в ледовых покровах Земли.

Один из больших периодов кривой изменения дO18 равен 100 млн. лет. Он близок к одному из периодов изменения эксцентриситета орбиты Земли e. Поэтому было принято, что эксцентриситет e оказывает существенное влияние на климат Земли. Например, в работе [7] даже утверждается, что М. Миланкович не учитывает прямого влияния эксцентриситета на климат Земли. В действительности теория инсоляции М. Миланковича дает математически строгую зависимость инсоляции от параметров e, цpг и е [17]. Это с одной стороны. А с другой стороны, предстоит еще большая работа по выяснению обстоятельств и причин образования определенных свойств морских осадков, ледовых кернов современных ледников и других палеоклиматических наборов данных. Такая же работа должна быть выполнена по выяснению и надежности выводов астрономической теории климата. Этой работой мы занимаемся в Институте криосферы Земли на протяжении двух десятков лет.

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

Вначале было отмечено, что решения М. Миланковича последовательно повторяли несколько поколений исследователей. Однако все они шли одним и тем же путем, который за столетия сложился в небесной механике. Мы пошли другим путем. Мы не копируем уравнения предшественников, а сами их выводим из первоначальных основ. Во-вторых, при выводе стремимся ввести минимальные упрощения. И, в-третьих, задачи решаем численными методами, стремясь использовать самые высокоточные их варианты или создать новые. В настоящей работе мы излагаем Астрономическую теорию изменения климата в результатах, полученных нами. Вначале речь пойдет об инсоляции Земли, затем об эволюции орбитального движения и об эволюции оси Земли. По первым двум задачам наши независимые исследования подтверждают исследования предшественников. А результаты вращательного движения - другие. Амплитуда колебаний угла наклона е в семь раз больше по сравнению с прежними теориями. Эти колебания дают такие колебания инсоляции, которые могут объяснить имевшие место колебания климата в прошлом. Эта существенная разница по результатам задачи о вращении Земли требует всесторонней проверки. В заключительной части будет уделено внимание проверкам результатов задачи о вращении Земли.

Геометрические характеристики инсоляции

В центре небесной сферы 1 (см. рис. 1) находится наблюдатель M. Его горизонт пересекает небесную сферу по кругу HH. Перпендикуляр к плоскости горизонта пересекает небесную сферу в точке зенита Z. Ось вращения Земли, обозначенная вектором угловой скорости Земли , пересекает небесную сферу в точке северного полюса N. Угол между и плоскостью горизонта является широтой наблюдателя. Напомним, что угол дуги большой окружности сферы равняется центральному углу между радиусами ее концов, например, дуга ц равна HMN.

Годовое движение Солнца S прочерчивает на небесной сфере 1 круг эклиптики EE против часовой стрелки. Он пересекает круг экватора AA в точке г. Долгота Солнца отсчитывается от точки г, которая является точкой весеннего равноденствия. Расстояние Солнца от экватора AA определяется склонением .

Вокруг оси MN Земля вращается против часовой стрелки. А небесная сфера совместно с Солнцем вокруг этой оси совершает суточное вращение относительно наблюдателя за часовой стрелкой. Поэтому суточное движение Солнца проходит по кругу SrMdSd, который параллелен кругу экватора. В точке Sr Солнце поднимается над горизонтом, в точку Md оно приходит в полдень, а в т. Sd заходит за горизонт. Часовой угол Солнца отсчитывается от точки полудня Md.

Следует отметить, что на рис. 1 плоскости экватора AA' и эклиптики EE' названы подвижными, т.к. они со временем изменяются в пространстве.



Рис. 1. Основные геометрические характеристики Солнца S при облучении точки M на земной поверхности: 1 - небесная сфера; НН - плоскость горизонта; N - северный полюс; AA - плоскость подвижного экватора; ЕЕ - плоскость подвижной эклиптики, а - угол между плоскостями AA и ЕЕ; Z - зенит точки M, а z = ZMS - зенитный угол Солнца; Дуга HN = - географическая широта точки M; щ = MdNS - часовой угол Солнца, отсчитываемый от полудня Md; = SB - склонение Солнца; = S - долгота Солнца.

Длительность дня пропорциональна длине дуги SrMdSd, а длительность ночи соответствует длине дуги SdMnSr. В показанном положении Солнца S день длиннее ночи. Если Солнце S будет на экваторе в т. г или г, то в течение суток оно будет перемещаться по кругу экватора AA. В этом случае день равен ночи. Если Солнце S будет находиться в южной части небесной сферы, то его путь под горизонтом будет длиннее пути над горизонтом HH, т.е. ночь будет длиннее дня.

На рис. 2а показано положение наблюдателя, находящегося на большей широте . В этом случае суточный круг MdS движения Солнца S не пересекает горизонт и находится над ним. Мы имеем полярный день. При нахождении Солнца S в южной части сферы суточный круг его движения, проходящий параллельно экватору, будет находиться под горизонтом. Для наблюдателя M наступает полярная ночь. При нахождении Солнца на круге ЕЕ в положениях более близких к т. г или г наблюдатель M будет иметь и ночь и день.

В летнее время при нахождении наблюдателя на северном полюсе, ц = 90° (см. рис. 2б), суточное движение Солнца происходит по кругу, параллельному плоскости горизонта HH'. В этом случае круг горизонта на небесной сфере совпадает с кругом экватора AA'. На полюсе нет восходов Солнца (т. Sr) и его заходов (т. Sd). В зимнее время суточное движение Солнца S' на рис. 2б происходит под горизонтом, т.е. на полюсе стоит полярная ночь.

При нахождении наблюдателя на экваторе ц = 0° (см. рис. 2в) северный полюс N лежит в плоскости горизонта HH', а точка зенита Z находится на плоскости экватора AA'. Суточное движение Солнца происходит по кругу SrMdSd, который параллелен кругу экватора AA'. В точках восхода Sr и заката Sd Солнце движется перпендикулярно горизонту, поэтому день начинается и заканчивается практически мгновенно.

Рис. 2. Геометрические характеристики Солнца S при нахождении наблюдателя M на разных широтах земной поверхности: а - в заполярных широтах; б - на Северном полюсе; в - на экваторе: 1 - небесная сфера; соответствующими стрелками показано годичное и суточное направление движения Солнца на небесной сфере; остальные обозначения см. рис. 1.

Инсоляция Земли

Угол падения солнечных лучей z на плоскость горизонта HH (см. рис. 1) отсчитывается от линии зенита ZM. В полдень Md он наименьший zmin, а в точке восхода Sr и заката Sd зенитный угол z = /2. Изменение этого угла в диапазоне от /2 до zmin и от zmin до /2 определяет изменение солнечной радиации в течении светового дня. Тогда количество солнечного тепла на 1 м2 земной поверхности в единицу времени, т.е. мощность инсоляции запишется в виде [19] - [20]:

, (1)

где J0 = 1366.22 Вт/м2 - поток солнечного тепла на расстоянии Солнца от Земли r, равным среднему радиусу земной орбиты a;

= r/a - относительное расстояние до Солнца.

Величина J0 называется солнечной постоянной, а зависимость (1) определяет инсоляцию Земли без учета влияния атмосферы. Если в формуле (1) зенитный угол z выразить через другие углы и солнечный поток dW/dt проинтегрировать за время длительности дня, то суточная инсоляция Земли запишется в виде:

, (2)

где = 243600 - длительность суток в секундах, а часовой угол границы дня щ0, определяемый дугами SrMd = SdMd (см. рис.1), зависит от углов ц и д и рассчитывается так [19] - [20]:

0 = arcos{-tg[arcsin(sinsin)]tg}. (3)

В выражения (2) и (3) для суточной инсоляции входят долгота л Солнца и его расстояние r от Земли. Они зависят от параметров орбиты Земли (см. рис. 3). Земля (E) движется по орбите против часовой стрелки. В точке перигелия PE она находится на ближайшем расстоянии Rp от Солнца, а в точке афелия AphE - на наибольшем Ra. Эксцентриситет орбиты определяется этими расстояниями

e = (Ra - Rp)/(Ra + Rp). (4)

Рис. 3. Схема движения Земли (E) по орбите вокруг Солнца (S): - точка весеннего равноденствия; PE и AphE - перигелий и афелий орбиты Земли, соответственно; о - полярный угол движения Земли по орбите; p - угол перигелия Земли.

На рис. 1 плоскость орбиты Земли расположена в плоскости эклиптики EE', а линия Sг орбиты совпадает с линией Mг. По кругу эклиптики EE' перемещается Солнце S относительно Земли, находящейся в т. М. Поэтому на рис. 1 изображения Земли будут зеркальными Солнцу. В точке г весной находится Солнце, а Земля в этот момент находится в т. г', если в т. М расположено Солнце S. Плоскости экватора и орбиты Земли изменяются в пространстве, поэтому точка по орбите Земли (рис. 3) перемещается. Положение перигелия PE независимо от точки также перемещается по орбите. Угол p измеряется между этими двумя точками и PE.

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



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

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



Рис. 5. Параметры орбиты и оси Земли в неподвижных экваториальной x,yz и эклиптической x1y1z1 системах координат.

1 - небесная сфера; неподвижные плоскости в эпоху Т0: 2 - экватора Земли, 3 -орбиты Земли (плоскость эклиптики); подвижные плоскости в эпоху Т: 4 -экватора Земли, 5 - орбиты Земли; единичные вектора: - оси Земли, - оси орбиты Земли; - вектор момента количества движения Солнечной системы; - точка весеннего равноденствия в эпоху Т0; B - положение перигелия на небесной сфере; = 02 - угловое расстояние восходящего узла орбиты; р = 2B - угловое расстояние перигелия; i - угол наклона плоскости орбиты к плоскости неподвижного экватора.

Угловые параметры и расстояние , входящие в формулу суточной инсоляции (2) зависят от параметров орбиты Земли и оси ее вращения, и они, как уже отмечалось, изменяются со временем. На рис. 5 кругами на небесной сфере 1 показаны плоскости экватора Земли 2 и ее орбиты 3, застабилизированные на определенную эпоху, например 2000 года. Угол между этими плоскостями равен е0. В другую эпоху плоскость экватора переместится в положение 4, а плоскость орбиты - в положение 5, и угол между ними будет е. Для определения законов изменения параметров, входящих в инсоляцию (2), необходимо решить две задачи об изменении: 1) орбитального и 2) вращательного движений.

Эволюция орбитального движения

Согласно закону всемирного тяготения тело с номером k притягивает тело с номером i, и это притяжение выражается силой:

, (5)

где G - гравитационная постоянная;

- радиус-вектор от тела с массой mk до тела с массой mi.

На i-тое тело все остальные k-тела оказывают воздействие силой (5). Тогда просуммировав силы всех воздействий, и поделив на массу i-того тела, получаем его ускорение в виде:

i = 1,2,…,n, (6)

где радиус-вектор тела относительно некоторого центра в неускоренной системе координат (в данном случае относительно центра масс Солнечной системы).

Соотношение (6) представляет систему 3n нелинейных дифференциальных уравнений, где п = 11 (девять планет, Солнце и Луна). Для их решения разработан метод высокой точности, на основании которого создана система Galactica [23], которая имеется в свободном доступе [24]. Эффективность метода показана в результате решения ряда задач современной небесной и космической динамики [25]-[29].

Эволюция орбиты рассматривается в неподвижной системе координат xyz, связанной с неподвижной плоскостью экватора 2 (см. рис. 5). Исследуется эволюция эксцентриситета орбиты e, угла ее наклона i к плоскости экватора 2, угла ц = г0г2, восходящего узла орбиты 2 и углового положения перигелия р = 2B, где B - проекция перигелия на небесную сферу 1.

На рис. 6 точками представлена динамика изменения элементов орбиты Земли: эксцентриситета e и углов , i и р. На промежутке 7 тысячелетий (от -3.4 тыс. лет до 3.6 тыс. лет) эксцентриситет е и наклон орбиты i уменьшаются, а угол перигелия р увеличивается. В то время как угол восходящего узла вначале этого интервала уменьшается, а затем возрастает. При этом минимум наступает за 1000 лет до 30.12.49 г. Угол перигелия р в отличие от других параметров изменяется неравномерно. Линиями 2 и 3 на рис. 6 приведены средние изменения этих параметров, которые получены С. Ньюкомбом [30] и Дж. Л. Симоном и др. [31] по наблюдениям. Из графиков видно, что результаты расчетов подтверждаются наблюдениями.



Рис. 6. Вековые изменения орбиты Земли 1 и сопоставления с аппроксимациями 2 и 3 данных наблюдения С. Ньюкомба [30] и Дж.Л. Симона и др. [31], соответственно: e - эксцентриситет; i -наклон плоскости орбиты к плоскости экватора 2000.0 г.; - угловое положение восходящего узла орбиты от оси x на эпоху 2000.0 г.; p - угловое положение перигелия в плоскости орбиты от восходящего узла; углы даны в радианах, а время Т - в столетиях от 30.12.1949 г.; интервал между точками - 200 лет; 1 cyr - столетие.

Эволюция параметров орбиты Земли за 3 прошедших млн. лет показана на рис. 7. Эксцентриситет испытывает короткопериодические изменения со средним периодом Те1 = 94.5 тыс. лет вокруг его среднего значения еm = 0.028. Кроме того, наблюдаются более длительные колебания с периодом Te2 = 413 kyr, которые приводят к достижению крайних значений эксцентриситета е = 0.0022 и е = 0.062. Долгота восходящего узла цЩ меняется со средним периодом ТЩ = 68.7 тыс. лет вокруг среднего значения цЩт = 0.068 радиан. На основной период накладываются колебания с большей длительностью.

Угол наклона плоскости орбиты і испытывает колебания с таким же периодом Ті = 68.7 тыс. лет вокруг среднего значения іm = 0.402 радиан. Колебания угла i происходят в пределах от 0.36 до 0.45 радиан. Диапазон колебаний составляет 5є.

Угол положения перигелия цр увеличивается со временем. Перигелий перемещается в направлении обращения Земли вокруг Солнца, совершая в среднем за 3 млн. лет один оборот за 147 тыс. лет. Вращение перигелия происходит неравномерно. Наряду со средним вращением против часовой стрелки существуют его возвратные движения по часовой стрелке.



Рис. 7. Эволюция орбиты Земли за 3 млн. лет в прошлое. Обозначения те же, что и на рис. 6. Величины с индексом m - средние значения параметров за 50 млн. лет; Тe, Тi, ТЩ - основные периоды колебаний соответствующих параметров в тысячах и миллионах лет, а Тp - средний за 3 млн. лет период обращения перигелия; интервал между точками - 10 тыс. лет; kyr - тысячелетие; Myr - 1 миллион лет.

В результате анализа полученных решений было установлено, что ось орбиты Земли (см. рис. 5) вращается за стрелкой часов, т.е. против орбитального движения Земли с периодом TS = -68.7 тыс. лет. Напомним, что ось орбиты - это перпендикуляр к ее плоскости 5. Благодаря этому вращению происходят колебания углов i и цЩ, показанных на рис. 7. Вращение оси орбиты происходит вокруг вектора (рис. 5), который является моментом количества движения всей Солнечной системы.

Дифференциальные уравнения орбитального движения (6) были проинтегрированы с помощью программы Galactica за 100 млн. лет. На этом интервале представить результаты в прежних параметрах невозможно, т.к. их колебания сливаются в сплошную полосу. Поэтому на рис. 8 показано эволюция параметров орбиты Земли за 50 млн. лет в несколько ином виде. Скользящие средние эксцентриситета es осреднялись на интервале, равном 2Te2 ? 850 kyr. Как видно из графика здесь наблюдается третий период изменения эксцентриситета Te3 = 2.31 Myr. Вместо углов орбиты i и цЩ на рис. 8 показана эволюция углов Ss и S оси орбиты относительно вектора (см. рис. 5). Угол наклона Ss, также представленный в виде скользящего среднего, является углом между векторами и .



Рис. 8. Эволюция параметров орбиты Земли за 50 млн. лет впрошлое: es - скользящие средние эксцентриситета; цp - угол перигелия; s - угол прецессии и Ss - скользящие средние угла наклона оси орбиты . Углы р и s даны в радианах, ?Ss - в градусах.

Максимальное отклонение оси от момента составляет ?Smax = 2.94є. Периоды ее колебаний следующие: Т97.4 kyr, Т2 1.16 Myr и Т2.32 Myr. Колебания с этими периодами видны на рис. 8. Угол прецессии S идентичен углу цЩ, но только измеряется не в плоскости экватора 2, а в плоскости, перпендикулярной моменту количества движения . Как и угол перигелия цр, угол прецессии S изменяется неравномерно, и существуют его возвратные движения. Однако в таком масштабе, как показано на рис. 8, этих колебаний цр и S не видно. Вследствие неравномерности вращения средние периоды перигелия и оси орбиты на интервале -50 млн. лет немного другие и равны: Tp = 150 тыс. лет, и TS = -72.9 тыс. лет, соответственно.

На интервале от -50 млн. лет до -100 млн. лет изменение параметров не отличается от рис. 8 [32]. Аналогичные исследования были выполнены для остальных планет и Луны [33] - [34]. Параметры их орбит так же, как и Земли совершают установившиеся колебательные движения. Поэтому можно сделать вывод, что, кроме установленных периодов и амплитуд колебания параметров орбит, других не существует. Из полученных результатов также следует вывод об устойчивости Солнечной системы. Необходимо отметить, что в работах других исследователей [35]-[37], которые решали задачу другими методами, были получены изменяющиеся колебания. Поэтому они пришли к выводу о неустойчивости орбит и Солнечной системы в целом [37].

Эволюция вращательного движения Земли

Вращающаяся Земля под действием центробежных сил растягивается в плоскости экватора (см. рис. 9). Рассмотрим воздействие тела B на две половины Земли: на ближнюю - силой и дальнюю - силой . Если бы Земля была центрально-симметричным шаром, то равнодействующая сил и прошла бы через центр Земли O. Для сплюснутой Земли центр масс ближней части Земли приблизится к телу B, а дальней - удалится. Поэтому сила увеличится, а сила уменьшится, в результате чего возникнет момент сил mO, направленный по часовой стрелке. Этот процесс воздействия на Землю определяется теоремой изменения момента количества движения:

, (7)

где - момент количества движения Земли относительно центра О в невращающейся системе координат x1y1z1, - сумма моментов сил, действующих на Землю тел.

Рис. 9. Прецессия оси Земли при воздействии на нее тела B: 1 и 2 - плоскости экватора Земли и ее орбиты, соответственно; 3 - плоскость орбиты воздействующего на Землю тела B.

Задача о вращении Земли рассматривается (см. рис. 5) в невращающейся системе координат x1y1z1, связанной с неподвижной плоскостью орбиты Земли 3. Подвижная плоскость экватора Земли 4 определяется углом наклона к плоскости 3 и углом прецессии = г0г1. Кроме того рассматривается скорость собственного вращения Земли относительно ее подвижной оси . На основании теоремы (7) мы вывели [38]-[39] дифференциальные уравнения вращательного движения Земли в следующем виде:

(8)

(9)

, (10)

астрономический теория изменение климат

где Jx, Jy и Jz - моменты инерции вращающейся Земли на оси системы координат, связанные с вращающейся Землей;

- динамическая эллиптичность Земли;

= const - проекция абсолютной скорости вращения Земли на ее ось (см. рис. 5)

n - количество действующих на Землю тел, а Mi - их масса, а x1i, y1i, z1i - их координаты.

Уравнениями (8) и (9) определяются углы положения и подвижной плоскости экватора 4 относительно плоскости неподвижной орбиты 3 (см. рис. 5). В результате интегрирования уравнений (8)-(9) определяются углы наклона и прецессии подвижного экватора 4 относительно неподвижной плоскости орбиты 3. А в результате интегрирования уравнений орбитального движения (6) рассчитываются углы наклона i и восходящего узла цЩ подвижной плоскости орбиты 5 относительно неподвижной плоскости экватора 2, а также угол перигелия цр. Эти параметры вращательного и орбитального движения позволяют определить углы наклона е и перигелия цpг подвижной плоскости орбиты 5 относительно подвижного экватора 4. Именно эволюция этих углов е, цpг и эксцентриситета e определяют эволюцию инсоляции.

Рассмотрим основные результаты задачи (8)-(9) о вращении Земли. На рис. 10 показано изменения угла наклона ? на пяти разных интервалах времени In. На графиках указаны основные периоды Tni и амплитуды (ai и) колебаний угла наклона: полумесячные Tn2, полугодовые Tn3 и Tn4 = 18.6 лет. Эти колебания называются нутационными. Угол прецессии имеет аналогичные периоды колебаний, а амплитуды - в 2 - 3 раза большие.

На интервале In = 0.1 года видны полумесячные колебания и просматриваются суточные; на интервале In = 1 год появляются полугодовые колебания; на интервале In = 10 лет появился тренд колебания с периодом 18,6 года, а на интервале In = 100 лет колебания с этим периодом преобладают.

На интервале In = 100 лет видно, что рассчитанный 1 угол наклона колеблется вокруг среднего 2 угла наклона по С. Ньюкомбу [30] и Дж. Симону и др. [31]. Амплитуда колебаний a49.2Ѕ периода Tn4 лет также совпадает с наблюдениями. Эта амплитуда в астрономии называется постоянной нутации. Рассчитанный угол прецессии также колеблются относительно осредненного угла прецессии по наблюдениям, а средняя динамика также совпадает с наблюдениями.

Рис. 10. Динамика угла наклона и (в радианах) плоскости экватора Земли к плоскости ее орбиты на пяти интервалах In времени: yr - год; - 0; 0 - угол наклона в начальную эпоху 30.12.1949 г. Tn2, Tn3, Tn4 и a2, a3, a4 - периоды и амплитуды колебаний угла наклона ; 1 - по нашим результатам численного интегрирования; 2 - аппроксимация данных наблюдений по С. Ньюкомбу [30] и Дж. Симону и др. [31]; 3 - по результатам интегрирования Дж. Ляскара и др. [12]; 4 - по результатам интегрирования Ш. Г. Шараф и Н. А. Будниковой [10].

Как видно из рис. 10, на интервале 10 тыс. лет такое совпадения с аппроксимацией данных наблюдения и результатами решений других авторов [10] и [12] имеется до 2000 лет. Далее рассчитанный нами угол наклона начинает отличаться от результатов решений [10] и [12].

На рис. 11 показано, что со временем отличия нарастают, и дальнейшая эволюция рассчитанного нами угла наклона ? значительно отличается от эволюции, полученной другими авторами, которые решали упрощенную задачу о вращении Земли. Как видно из графиков, колебания угла ? по нашим решениям происходят в диапазоне от 16.7є до 31є, в то время как по прежним решениям от 22.26є до 24.32є, т. е. диапазон колебаний больше в 7 раз. Мы выполнили расчеты инсоляции Земли с новой эволюцией угла. Ее колебания также в 7 раз превышают колебания инсоляции в прежних астрономических теориях палеоклимата.

Второй угол в уравнениях (8) - (10), угол прецессии , изменяется несколько по другому, по сравнению с углом . Он непрерывно уменьшается, но это уменьшение совершается с такими же колебаниями, как и угла , только амплитуда их в 2-3 раза больше. Такое изменение угла свидетельствует, что ось Земли (см. рис. 5) вращается за стрелкой часов, и средний период ее вращения PN = -25740 лет.

Выполненные исследования по этим двум задачам позволяют представить следующую картину движений в Солнечной системе. В плоскости орбиты происходит неравномерное вращение перигелия против стрелки часов со средним периодом Tp = 147 kyr и деформация орбиты с колебаниями эксцентриситета с периодами 94.5 kyr, 413 kyr, 2.31 Myr. Аналогичные изменения орбиты со своими периодами испытывают Луна и другие планеты, за исключением Плутона. Перигелий Плутона вращается по часовой стрелке.



Рис. 11. Эволюция угла наклона (в радианах) плоскости экватора Земли к плоскости ее орбиты на интервале 200 тыс. лет: 1 - по нашим результатам численного интегрирования; 2 - по результатам интегрирования Дж. Ляскара и др. [12]; 3 - по результатам интегрирования Ш. Г. Шараф и Н. А. Будниковой [10]. В градусах приведены максимальные и минимальные значения угла .

Ось орбиты Земли (см. рис. 5) вращается по часовой стрелке вокруг вектора момента со средним за 3 млн. лет периодом TS = 68.7 kyr. Ось также совершает колебания с разными периодами 97.4 kyr, ?1.16 Myr и 2.32 Myr. Ось орбиты Луны (на рис. 5 не показана) вращается по часовой стрелке вокруг подвижной оси орбиты Земли. Период этого вращения равен 18.6 лет. Кроме того, ось орбиты Луны совершает колебательное движение с периодом 0.4745 года.

Ось вращения Земли (см. рис. 5), как уже отмечалось, вращается по часовой стрелке, как и ось орбиты Луны, вокруг подвижной оси орбиты Земли со средним периодом 25740 лет. При этом ось совершает колебания с периодами: полумесячным, полугодовым и 18.6 лет. Кроме того, существуют колебания оси большой амплитуды (см. рис. 11) и с периодами десятки и сотни тысяч лет. В процессе этих колебаний угол наклона оси Земли изменяется от 16.7° до 31°.

Эволюция инсоляции

На рис. 10 и рис. 11 было представлено изменение угла наклона оси Земли в будущем. На рис. 12 линией 1 показана эволюция угла наклона за прошедшие 200 тыс. лет. По сравнению с расчетами по прежним теориям (линия 2) наши решения так же, как при расчетах в будущее, дают в 7 раз большие амплитуды колебаний. В начальный период несколько тысяч лет в прошлое, начиная с T = 0, угол наклона 1, как и угол наклона 2, увеличивается. Затем он начинает уменьшаться и достигает минимума, в то время как угол наклона 2 по прежним теориям находится в максимуме. На остальном интервале времени также максимумы и минимумы колебаний углов наклона 1 и 2 не совпадают. Однако более существенным являются величины этих экстремумов. По прежним теориям на этом интервале времени угол наклона изменяется от 22.21 до 24.43. А по нашим решениям наклон плоскости экватора Земли к плоскости ее орбиты изменяется от 14.8 до 32.1. Такой же диапазон изменений угла е получен при решении задачи за будущие 200 тыс. лет (см. рис. 11).

В астрономических теориях палеоклимата рассматривается инсоляция не за астрономические, а за равные калорические полугодия. Начало и конец летнего калорического полугодия определяется так, чтобы инсоляция за любой его день была больше инсоляции за любой день зимнего полугодия. Далее будем рассматривать инсоляцию на широте 65° северного полушария, которое обозначается буквой N. Мы рассчитали изменение инсоляции QS65N за 200 тыс. лет в прошлое как по определенным нами параметрам: e, е и цpг (линия 1 на рис. 12), так и по тем же параметрам, рассчитанных Дж. Ляскаром и др. [12] (линия 2). Как видно из графиков, инсоляция QS65N за летнее калорическое полугодие на северной широте 65 по нашим решениям имеет амплитуду колебаний также в 7 раз большую, нежели по прежним теориям. Моменты потепления и похолодания по нашим расчетам 1 и по прежним теориям 2 также не совпадают. Вначале, начиная с T = 0, как видно из QS65N на рис. 12, летняя инсоляция растет в течение 4-х - 5-и тысяч лет, затем начинает падать и достигает минимума 16 тысяч лет назад. После этого минимума идет потепление, которое завершается большим максимумом инсоляции 31 тыс. лет назад.

Средний период колебаний угла наклона по прежним теориям (см. 2 на рис. 12) равен 41.1 тыс. лет. Этот же период имеется в колебании инсоляции QS65N (линия 2). Период колебаний 41.1 тыс. лет является наиболее упоминаемый в астрономической теории палеоклимата. Как видно из новой зависимости для угла наклона е (линия 1) характерный период колебаний имеет в 1.5 - 2 раза меньшую длительность.



Рис. 12. Эволюция угла наклона е, а также летних инсоляций QS65N и I за 200 тыс. лет в прошлое. Сравнение новых результатов 1 с результатами прежних теорий 2 на примере работы Дж. Ляскара и др. [12]. е - угол наклона в радианах экватора Земли к плоскости ее орбиты; QS65N - инсоляция в ГДж/м2 за летнее калорическое полугодие на северной широте 65; I - инсоляция в эквивалентных широтах за летнее калорическое полугодие на северной широте 65. В градусах приведены максимальные и минимальные значения угла е.

По нашим решениям получены в 7 раз большие колебания инсоляции. Насколько они существенны? На этот вопрос может дать ответ представление инсоляции в эквивалентных широтах I, которая рассчитывается следующим образом. Если в эпоху T летняя инсоляция на широте ц была такая, как в современную эпоху на широте ц0, то инсоляция в эквивалентных широтах будет I = ц0. На рис. 12 представлена инсоляция I в эквивалентных широтах, рассчитанная нами для широты 65° северного полушария как по нашим данным (линия 1) так и по данным Дж. Ляскара и др. [12] (линия 2). Начиная с T = 0, инсоляция I по нашим данным 1 с широты 65° уменьшается на несколько градусов, т.е. на широте 65° становится теплее. Затем после максимума величина I начинает уменьшаться до широт 80° и 90°. При T = -15 тыс. лет летняя инсоляция на широте 65° меньше, чем современная летняя инсоляция на полюсе, поэтому на графике I проведена горизонтальная площадка. Таким образом, горизонтальная площадка примерно 15 тыс. лет назад означает, что инсоляция на широте 65° была меньше, чем сейчас на полюсе. Такое малое количество тепла могло привести к оледенению территорий на широте 65°.

При дальнейшем изменении времени до T = -30 kyr инсоляция I в эквивалентных широтах достигает широт 50°, 40°, 30°, т.е. на широте 65° становится значительно больше солнечного тепла. Горизонтальная площадка в эпоху T = -30 kyr означает, что на широте 65° поступает больше тепла, чем сейчас на экваторе.

Линией 2 приведена инсоляция I в эквивалентных широтах по прежним теориям. Как видно, на широте 65 летняя инсоляция I на рассмотренном отрезке времени в 30 - 50 тыс. лет изменяется в пределах от 60 до 70. Маловероятно, что изменения количества тепла на широте 65 до таких величин, которые сейчас имеются на широтах 60 и 70, могут привести к существенным потеплениям климата или его похолоданиям. Такие небольшие изменения инсоляции всегда вызывали сомнения [7].



Рис. 13. Распределение по широте ц° Земли удельного количества тепла в ГДж/м2 в самую жаркую в высоких широтах эпоху T = 31.28 тыс. лет в прошлое: Qs - за летнее калорическое полугодие; Qw - за зимнее калорическое полугодие; QT - за весь год: на графике величина QT уменьшена в два раза; > 0 - северное полушарие; < 0 - южное полушарие.

Выше мы исследовали изменение инсоляции во времени на широте 65° северного полушария. Теперь рассмотрим изменение инсоляции по широте в отдельные моменты времени. На рис. 4 была представлено распределение инсоляции по всему земному шару в современную эпоху. На рис. 13 показано изменение по широте ц° летней Qs, зимней Qw и уменьшенной в два раза годовой инсоляции QT в эпоху T = -31.28 тыс. лет. В это время на широте 65° северного полушария была наибольшая за 200 тыс. лет летняя инсоляция (см. рис. 12). Как и в современную эпоху (см. рис. 4), инсоляция за год QT - наибольшая на экваторе и убывает к полюсам. Зимняя инсоляция Qw имеет максимальное значение вблизи экватора, а на полюсах стремится к нулю. Летняя инсоляция Qs имеет максимумы на тропиках, которые находятся на больших широтах по сравнению с современной эпохой (см. рис. 4). Кроме того, летняя инсоляция Qs на высоких широтах близка к максимальной и существенно превышает летнюю инсоляцию на экваторе.

Инсоляция за год QT в эпоху T = -31.28 тыс. лет на высоких широтах больше (рис. 13), чем в современную (рис. 4). В то же время годовая инсоляция на экваторе меньше, чем в современную эпоху. Но это уменьшение не такое значительное, как увеличение в высоких широтах. Зимняя инсоляция Qw в эпоху T = -31.28 на всех широтах меньше по сравнению с современной эпохой (рис. 4).

Рис. 14. Распределение по широте ц° Земли удельного количества тепла в ГДж/м2 в самую холодную в высоких широтах эпоху T = 46.44 тыс. лет в прошлое. Обозн. см рис. 13.

На рис. 14 представлено изменение этих же составляющих инсоляции в наиболее холодную за 200 тыс. лет эпоху T = -46.44 тыс. лет. В высоких широтах, например, на широте ц° = 65°N летняя инсоляция Qs = 4.72 ГДж/м2 существенно уменьшилась по сравнению с Qs = 7.43 ГДж/м2 в эпоху T = -31.28 тыс. лет и Qs = 5.92 ГДж/м в современную эпоху T = 0. Инсоляция за год, например QT = 3.58 ГДж/м2 на Северном полюсе, уменьшилась с QT = 7.43 ГДж/м2 в эпоху T = -31.28 тыс. лет и с QT = 5.58 ГДж/м2 в современную эпоху T = 0. Зимняя инсоляция Qw в эту холодную эпоху немного возросла на всех широтах.

Итак, от теплой эпохи к холодной уменьшение летней инсоляции Qs на широте 65° происходит в 1.57 раза. А на полюсе годовая инсоляция QT уменьшается еще сильнее, а именно - в 2 раза. В то же время в экваториальной зоне происходят обратные по знаку изменения, но величина их в 4 раза меньшая. При этом на широте ц = 45° годовая инсоляция QT практически не изменяется, т.е. изменения летней инсоляции Qs компенсируются изменениями зимней инсоляции Qw.

Полученные в наших расчетах колебания инсоляции могут вызвать наблюдавшиеся изменения климата. Уменьшение летней инсоляции I 19 - 12 тыс. лет назад до значений меньших, чем на полюсе (линия 1 на рис. 12), могло привести к оледенению. Это, по-видимому, подтверждается существованием в это время покровного ледника в северных частях Евразии и Америки. А увеличение летней инсоляции I 35 - 28 тыс. лет назад до значений больших, чем в экваториальной области, могло привести к благоприятным условиям существования мамонтовой фауны. По-видимому, это подтверждается преимущественной датировкой ископаемых представителей этой фауны.

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

Проверка достоверности решения задачи о вращении Земли

Как уже отмечалось во Введении, все три задачи: инсоляция Земли, ее орбитальное движение и вращательное движение Земли мы решали по-другому, начиная с вывода уравнений и заканчивая обработкой и анализом результатов. Например, дифференциальные уравнения (6) для орбитальной задачи ранее непосредственно не решали. Начиная с И. Ньютона, эту задачу рассматривали в рамках так называемой теории возмущений. Исследовалось движение одного тела под воздействием материнского тела: Солнца - для планет, и Земли - для Луны, а воздействия остальных тел рассматривались как возмущающиеся факторы. В итоге, в результате работы поколений таких выдающихся ученых как Эйлер, Лагранж, Лаплас и мн. др. были получены уравнения возмущенного движения для шести элементов орбиты относительно неподвижной эклиптики (см. 3 на рис. 5). Эти уравнения приближенные. А затем они также решались приближенными аналитическими методами разложением в ряды. В настоящее время эти методы переведены на численный счет, благодаря чему количество членов в рядах достигает несколько сотен. Тем не менее, приближенный характер уравнений при решении задач за сотни тысяч и миллионы лет всегда оставляет сомнение в справедливости результатов. Поэтому полученные нами решения неупрощенной орбитальной задачи (6) позволили подтвердить многие результаты, полученные прежними методами, а некоторые результаты, как неустойчивость орбит и Солнечной системы в целом, были опровергнуты.

Задача о вращении Земли (8)-(10), при ее решения за большие промежутки времени в работах прежних авторов упрощалась более существенно, по сравнению с орбитальной. В уравнениях (8)-(9) отбрасывались вторые производные и , а также отбрасывались произведения первых производных и . Поэтому такое решение задачи не позволяло получить короткопериодических колебаний, которые получены нами и показаны на рис. 10. А полученные колебания угла е с периодом 41 тыс. лет (см. линия 2 на рис. 11 и на рис. 12) не имели других подтверждений.

Задачу (8)-(10) о вращении Земли за сотни тысяч лет мы решили впервые. В рамках принятой технологии ее решения все необходимые проверки были выполнены. Например, задача решалась последовательно, при воздействии одного из тел [38]. Полученные периоды колебаний оси Земли подтвердились общетеоретическими выводами, а также результатами других авторов [39]. При воздействии всех тел задача решалась за разные интервалы времени, и полученные результаты подтвердились наблюдениями. Интегрирование уравнений (8)-(10) за 200 тыс. лет проводилось с разными начальными условиями и с разными шагами интегрирования. Вид результатов от этого не изменился. Однако полученные амплитуды колебаний оси Земли в 7 раз превышали амплитуды по прежним теориям. Поэтому требовалось подтверждение полученных результатов. По-видимому, нужно решать эту задачу каким-то другим способом. На протяжении ряда лет мы занимались поиском таких способов и их реализацией.

Задача о вращении Земли является одной из самых сложных задач механики. Это также подтверждается видом уравнений (8)-(10). Их вывод связан с рядом переходов из одной системы координат в другую, с принятием определенных упрощений и приближений. Поэтому кардинальной проверкой полученных результатов являлось бы их получение без решения дифференциальных уравнений (8)-(9).

При исследовании орбит выяснилось, что эволюция орбиты Луны подобна эволюции оси вращения Земли. Этот результат привел нас к составной модели Земли (см. рис. 15), в которой часть массы Земли равномерно распределена между периферийными телами, обращающимися вокруг центрального тела по круговой орбите. Под воздействием Луны, Солнца и планет орбиты периферийных тел начинают изменяться. Эволюция оси орбиты одного из этих тел моделирует эволюцию оси вращения Земли. В первоначальной серии исследований [26] было исследовано три модели, и была подтверждена возможность моделирования эволюции Земной оси. В этих моделях периоды прецессии осей орбит моделей были 170 лет и 2604 года, в то время как средний период прецессии оси Земли TprE = 25740 лет. В последующем было разработано еще 11 моделей, пока не был достигнут необходимый период прецессии.

Рис. 15. Составная модель вращения Земли. Масса Земли распределена между центральным телом и периферийными телами: a - радиус орбиты периферийных тел.

Такое моделирование вращательного движения Земли включает несколько стадий решения орбитальной задачи (6) с помощью программы Galactica. Выяснилось, что получить модель с необходимым периодом прецессии возможно в случае, если ослабить или усилить гравитационное взаимодействие между телами модели. Поэтому была создана программа Galactica с измененным взаимодействием между определенными телами.

На рис. 16 точками представлены результаты 13-ой составной модели Земли на интервале 300 лет. Видно, что угол наклона ? и угол прецессии колеблются с периодом 18.6 лет. Решения на меньших интервалах времени дали колебания с полумесячными и полугодовыми периодами, т.е. эти результаты совпали с результатами интегрирования уравнений (8)-(10). Амплитуды этих колебаний также совпали. Прямыми линиями на рис. 16 нанесены средние изменения углов ? и . Они также совпали с результатами прямой задачи о вращении Земли, а также с данными наблюдений. Идентичность результатов видна также при сопоставлении рис. 16 с графиком для угла ? на интервале 100 лет на рис. 10.

Такое совпадение результатов модельной задачи с результатами прямой происходит до 3 тыс. лет. В процессе интегрирования уравнений (6) программой Galactica происходит накопление погрешностей для тел модели, вследствие чего ее размеры начинают изменяться. Изменяется также период прецессии орбиты. Например, к концу интервала интегрирования 13.763 тыс. лет период прецессии с 25740 лет уменьшился до 14840 лет. Поэтому модель перестает представлять эволюцию Земной оси. Этот процесс обусловлен весьма напряженными динамическими параметрами модели. Например, радиус орбиты периферийных тел (см. рис. 15) a равен радиусу Земли, их период обращения 0.142 часа, а взаимодействие между телами модели усилено в 9.6 раза по сравнению с гравитационным. Таким образом, тела модели в 170 раз вращаются быстрее Земли. Поэтому шаг интегрирования задачи (6) приходится уменьшить в 1000 раз по сравнению с шагом при решении орбитальной задачи. Это приводит к большому времени счета. Например, решение задачи за 13.763 тыс. лет заняло 2.13 месяца. Чтобы модель не изменялась на интервале 200 тыс. лет, нужно уменьшить шаг интегрирования до таких величин, при которых требуется невыполнимое время счета.



Рис. 16. Эволюция углов наклона и прецессии оси вращения составной модели Земли № 13 за 300 лет. Точками приведены результаты интегрирования уравнений (6) с помощью программы Galactica. Расстояние между точками 3 года. Линиями показаны средние изменения углов угла е и ш со скоростями = - 2.28·10-4 1/cyr = - 0.470 Ѕ/год. = - 2.42·10-2 1/cyr = - 49.9 Ѕ/год.

Итак, на данном этапе составная модель Земли на интервале 3000 лет подтвердила полученные результаты интегрирования дифференциальных уравнений (8)-(10) вращения Земли. Это свидетельствует, что принятые допущения и упрощения при выводе уравнений (8)-(10), их вывод, метод решения и преобразование результатов интегрирования в конечный вид также подтвердился.

Вторая независимая проверка заключалась в использовании более точного метода при интегрировании уравнений (8)-(10), а именно метода Рунге-Кутты 8-ого порядка в реализации Дорманда и Принса [40]. До этого использовался метод интегрирования Рунге-Кутты 4-ого порядка в реализации [41]. Мы использовали его в течение нескольких десятков лет для решения разнообразных задач, и всегда получали удовлетворительные результаты. При интегрировании уравнений вращательного движения (8)-(10) за интервалы времени порядка 200 тыс. лет, мы столкнулись с неожиданным недостатком метода. В решениях этих уравнений существуют суточные колебания производных и . Амплитуда этих колебаний к концу вышеупомянутого интервала интегрирования увеличивалась на несколько порядков. Несмотря на то, что были проведены проверки по влиянию суточных амплитуд и на конечные результаты и приняты определенные меры по их устранению, оставалась опасность непредусмотренного их влияния.

Интегрирование уравнений (8)-(10) было выполнено новым методом в реализации Дорманда и Принса на разных интервалах времени, включая 200 тыс. лет. Все полученные ранее результаты подтвердились. При этом амплитуды суточных колебаний производных и не увеличиваются и остаются на одном уровне. Итак, метод интегрирования уравнений не влияет на полученные результаты, и применение более точного метода подтверждает их.

Третья независимая проверка заключалась в смене технологии решения задачи. В дифференциальные уравнения вращательного движения (8)-(9) входят привязанные к плоскости орбиты координаты x1i, y1i, z1i воздействующих на Землю тел. При решении орбитальной задачи (6) с помощью программы Galactica мы получаем координаты тел xi, yi, zi, привязанные к плоскости неподвижного экватора. Они пересчитывались в координаты x1i, y1i, z1i. Однако при интегрировании задачи (8)-(9) за большие промежутки времени массив этих данных будет занимать нереализуемый объем памяти. Поэтому мы создали математическую модель Солнечной системы [21], которая выдает в необходимый момент времени координаты тел: x1i, y1i, z1i на основании результатов задачи двух тел: материнского тела и его спутника. При этом параметры орбиты тела: e, i, цЩ, цр, Rp. и др. в каждый момент определяются по данным, предварительно рассчитанных с помощью программы Galactica. В процессе решения этой задачи математическая модель Солнечной системы была всесторонне проверена. Тем не менее, оставалась вероятность того, что на больших интервалах времени несущественные отличия результатов математической модели Солнечной системы от значений координат, полученных с помощью программы Galactica, могут повлиять на эволюцию параметров вращательного движения е и .

В течение работы над задачей о вращении Земли мы неоднократно возвращались к вопросу о смене вышеописанной технологии. Однако все попытки были безуспешными, пока не появилась определенная идея эффективного объединения этих двух задач в одну. В итоге была разработана новая программа совместного решения орбитальной задачи и задачи о вращении Земли. В ней за один шаг по времени происходит решение методом программы Galactica орбитальной задачи (6), затем методом Дорманда и Принса за этот шаг по времени решается задача о вращении Земли (8)-(10). С помощью новой программы были выполнены решения этих двух задач за разные интервалы времени, включая интервал 200 тыс. лет. Все ранее полученные результаты подтвердились. Этой проверкой была также подтверждена математическая модель Солнечной системы на большом интервале времени.