Математическая модель гляциоизостатических движений 
 
Для входа в суть проблемы математического моделирования некоторого процесса, работу полезно начинать с решения простейшей задачи. В нашем случае такая задача была давно поставлена и решена [Г. Ламб. Гидродинамика. М-Л, 1947]. В ней рассматривается вязкая самогравитирующая сфера. Это значит, что предполагается в отсутствии внешних полей вязкая однородная жидкость, сжатая в шар собственным гравитационным полем. На поверхности этой сферы появляется возмущение некоторой формы. Так если это возмущение можно разложить в ряд по сферическим функциям Лежандра, то даётся решение возвращения этой сферы к невозмущённому состоянию (см. статью на сайте - файл stat74-1.pdf). Именно это решение я и использовал для оценки вязкости материала астеносферы.
 
По литературным источникам были оценены объёмы воды, связанной в ледниковых покровах. По этим сведениям уровень Мирового океана должен был понизиться примерно на 100 м. Вся задача рассматривалась в осесимметричном случае. По современным скоростям поднятия области Фенноскандии (Балтийский ледниковый щит) была оценена вязкость жидкости шара (мантии) в 2*10^21Па*с.
 
 С помощью этой модели была получена зависимость, позволяющая рассчитывать, с учётом гляциоизостатического фактора, изменения объёма материкового оледенения на земном шаре, по известному временному ходу изменений уровня Мирового океана (см. файл stat74-2.pdf) . Из известных кривых по изменениям УМО, в качестве рабочей, для последних 35 тысяч лет, была выбрана кривая Миллимена и Эмери [Millimen J.D. & Emery K.O. Sea levels during the past 35, 000 years. "Science”, Vol. 162, 1968, No 3858, p. 1121-1123]. В вышеупомянутой модели рассматривается влияние гляциоизостатического фактора, вызванного одним ледниковым покровом, причём в осесимметричном случае. В прошлом же на Земле существовала целая группа ледниковых покровов; следовательно, нашу задачу необходимо рассматривать в более общем виде, то есть учитывать на поверхности Земли влияние всего ансамбля ледников. Кроме того, для более или менее реалистичной картины гляциоизостатических движений, необходимо иметь в виду реальные контуры континентов, рельеф их окраин и бассейнов морей. Следовательно, нам нужна была глобальная модель в распределённых параметрах на поверхности сферы. Её построение отражено в файле stat76.pdf.
 
Итак, рассматривается модель Земли в виде твердой сферы, имеющей тонкую трёхслойную оболочку. Нижний слой оболочки - астеносфера – вязкий, пластичный слой постоянной плотности. Где-то внутри астеносферы проходит поверхность компенсации. Второй слой – твёрдая литосфера, включающая в себя надастеносферный субстрат и кору. Третий слой – гидросфера – непрерывен на океанах и включает в себя крупные ледниковые покровы на континентах. Движения в астеносфере вызываются изменениями поверхностной нагрузки в виде изменения уровня моря или толщины континентальных льдов. Известным советским геофизиком-теоретиком Евгением Викторовичем Артюшковым установлено, что литосфера, вследствие её блокового строения, практически не влияет на изостатические движения литосферы. Поэтому он полагает, что все точки её могут двигаться только в вертикальном направлении.
 
Далее работа идёт в стандартном направлении. В сферических координатах рассматриваются три покомпонентных уравнения гидродинамики Навье-Стокса и уравнение неразрывности. Ввиду чрезвычайной медленности астеносферных движений в уравнениях усекаются все «малые» члены и остаются только «весомыеые». Далее, полученная система уравнений решается при известных краевых условиях… и мы получаем одно единственное динамическое уравнение, которое и надо проинтегрировать с учётом изменений ледниковых покровов. А поскольку у нас в модели есть и контуры континентов, и рельеф их окраин, и мощности ледниковых покровов, и их контуры, то проинтегрировать эту систему аналитически невозможно. Остаётся один выход – переходить к численному решению задачи.
 
Численное решение упирается в разностные схемы и в вычислительные способности (мощность, скорость и объём обрабатываемой информации) ЭВМ. Что такое разностная схема? Большинство живых организмов, не говоря уже о людях, постоянно используют разностные схемы. Рассмотрим, например, переход через дорогу в нерегулируемом месте с довольно оживлённым движением. Мы вдалеке видим приближающийся автомобиль, под который нам не хочется попасть – засекаем, где он находится. Спустя некоторый интервал времени, величину которого мы фиксируем в мозгу, бросаем снова взгляд на автомобиль. Автомобиль сместился на некоторое расстояние – его величину мы интуитивно делим на промежуток прошедшего времени и т.о. определяем скорость авто. Далее, деля расстояние до автомобиля на скорость его, мы находим через сколько времени авто может налететь на нас. Оценивая ширину дороги и нашу обычную скорость, мы находим величину времени, необходимую нам для перехода дороги. Если это время значительно меньше времени - авто, то мы начинаем переход. И эти «расчёты» интуитивно производят в уме и очень умные люди и не очень и, при других обстоятельствах, животные. Так же осуществлялся и вычислительный процесс на ЭВМ в старых программах.
 
Если бы мы видели сверху верстовые столбы и имели секундомер, то мы бы точнее могли определять скорость авто, и задача перехода дороги могла бы решаться тоже точнее. Следовательно, нам нужна вычислительная сетка на сфере. Лучше всего если эта сетка будет однообразной – регулярной, тогда и вычислительный процесс можно организовать монотонно однообразным, максимально приспособленным для решения на компьютере. Наиболее простым было бы разделение сферы параллелями и меридианами на боксы (ячейки), одинаковыми в градусной мере но разными по площадям. Но ввиду ограниченных вычислительных способностей компьютеров (ЭВМ) тех времён приходилось экономить на всём, в том числе и на количестве боксов вычислительной сетки. Это делалось укрупнением ячеек регулярной сетки по широте к полюсам, при этом, кстати, площади ячеек были примерно одинаковы. С ростом вычислительных способностей компьютеров, в конце концов, используемая мною сетка стала регулярной 5х5 градусов по широте и долготе (файлы stat79.pdf и stat81.pdf).
 
Как работает вычислительная схема? К каждому рассматриваемому боксу на его границах с четырёх сторон примыкают четыре соседних бокса. Только у боксов, «упирающихся» в северный полюс, отсутствует северный «сосед», а упирающихся в южный полюс – южный). Компьютер поочерёдно опрашивает боксы, соседствующие с рассматриваемым, на следующий предмет, где актуальная меняющаяся поверхностная нагрузка больше: на самом рассматриваемом боксе или на текущем - соседнем. По имеющемуся динамическому закону (при известном коэффициенте вязкости астеносферы) рассчитываемому на определённый интервал времени, для уравновешивания разности этих нагрузок, компьютер «перемещает» часть астеносферного материала из перегруженного бокса в недогруженный. При малом интервале шага по времени, объём перемещаемого материала будет тоже малым. После этого нагрузка на поверхности соседнего бокса пересчитывается за счёт изменения астеносферного слоя. Таким образом обходятся все боксы сферы. Задаётся новый интервал времени и расчёт повторяется.
 
При достаточно малом шаге по времени и при большой точности вычислений этот процесс очень устойчив. Так в своё время, но вручную, рассчитывал планетные орбиты покойный Исаак Исаакович Ньютон.
 
Перейдём к рассмотрению самих задаваемых поверхностных нагрузок, то есть к динамике последних ледниковых покровов Плейстоцена. Данные по максимальным мощностям ледниковых покровов по боксам брались с карты проекта CLIMAP [ McIntyre A., Moore T.C., Andersen B. et al. The surface of the ice-age Earth. – Science, 1976, v.191, p. 1131-1137]. Они увеличивались в полтора раза, для получения максимального понижения уровня Мирового океана в ходе оледенения примерно на 135м и таким образом для согласования с кривой Миллимена и Эмери. Также учитывалось положение о том, что в течение последнего оледенения объёмы Антарктического и Гренландского ледниковых покровов увеличивались на 20% от современного [Суетова И.А. Основные морфометрические характеристики Антарктиды. – Гляциология . 168, №19].
 
Сами боксы сетки подразделялись на три типа: морские, континентальные и полуморские (с наличием в них и суши и моря). Нагрузка на морских боксах изменялась только за счёт изменения уровня Мирового океана (УМО). Нагрузка на континентальных – только за счёт ледниковых покровов, если они там были в прошлом. Нагрузка на полуморских складывалась за счёт изменений: УМО; ледниковых покровов; самой высоты, например, острова или края континента над изменяющимся УМО. В последнем имеется ввиду, то сколько воды добавится/убавится в боксе за счёт залива/отлива прибрежной/шельфовой зоны. Эта величина аппроксимировалась линейными функциями для каждого полуморского бокса отдельно.
 
Численная реализация первоначально проводилась на ЭВМ БЭСМ-6 а далее на машинах серии ЕС: 1020, 1030, 1033, 1060. Верификация модели производилась следующим образом. Проводились численные эксперименты с различными коэффициентами вязкости астеносферы во временных пределах от 35000 лет назад до нашего времени. В результате были получены данные по вертикальным движениям земной коры, в том числе и по территориям, нагруженным в прошлом крупными ледниковыми покровами: центральная часть Северо-Американского покрова и  Балтийского ледникового щита. Для этих районов по литературным источникам известны современные скорости поднятия литосферы. Путём сравнения результатов моделирования с этими скоростями, а также с их изобазами для максимума оледенения, был определён коэффициент эффективной вязкости материала астеносферы – 4*10^18Па*с [см. файл stat84.pdf]. Таким образом, вязкость материала астеносферы в экспериментах была «подогнана» под реальные современные скорости движений коры.
 
Эта модель должна была войти одной из компонент в глобальную климатическую модель «ледники-океан-атмосфера» [Сергин В.Я., Сергин С.Я. Системный анализ проблемы больших колебаний климата и оледенения Земли. Л., Гидрометеоиздат, 1978, 280 с.]. К сожалению, эта модель так и не была построена.
 
Немного раньше чем у меня, примерно аналогичная (по ожидаемым результатам) модель, была построена группой американских учёных [Clark J.A., Farrell W.E., Peltier W.R. Global changes in postglacial sea-level: A Numerical Calculations. – Quaternary Res., 1978, v.9, p. 265-287.] . Их модель была полуаналитической – методом функций Грина, получисленной. Их было трое вначале, потом компания по крайней мере удвоилась за счёт итальянцев. Мне же всё приходилось делать одному, и это на грани выживания. Трудно быть учёным в СССР а потом и в России. Результаты у них, в принципе – те же, только в США они были как-то востребованы. Никаких важных глобальных выводов американцы не сделали. Представляю, какой бы они подняли шум, если бы разглядели в результатах своих экспериментов и Атлантиду, и некоторое другое. А у нас упомянуть в трудах АН СССР Атлантиду, даже если ты в этом твёрдо уверен – мракобесие. Однако, диалектический материализм. Именно этому «мракобесию» я и посвящу следующие страницы сайта.