Численное моделирование трехмерных задач конвекции в мантии Земли с применением завихренности и векторного потенциала Текст научной статьи по специальности «Физика»

Вычислительные технологии

Том Т, № 1, 2002

ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ТРЕХМЕРНЫХ ЗАДАЧ КОНВЕКЦИИ В МАНТИИ ЗЕМЛИ С ПРИМЕНЕНИЕМ ЗАВИХРЕННОСТИ И ВЕКТОРНОГО ПОТЕНЦИАЛА *

В. В. ЧЕРВОВ

Объединенный институт геологии, геофизики и минералогии СО РАН

Новосибирск, Россия e-mail: [email protected]

The numerical model of convection processes in Earth mantle is constructed on the basis of vorticity and vector potential. The results of numerical experiments are presented.

The comparison with well-known benchmark problem shows good results.

Введение

Численному моделированию трехмерных конвективных течений в мантии Земли посвящено большое количество работ (см., например, [1-7] и приведенную в них библиографию). Особо можно отметить статью Биээе Г. Н. е! а1. [4], в которой содержатся результаты трехмерного моделирования конвективных течений в мантии, полученные разными авторами в применении к одной и той же задаче. При решении трехмерных задач гидродинамики достаточно хорошо зарекомендовали себя переменные “завихренность” — “векторный потенциал” [8-13]. Анализ известных работ показывает, что при численном моделировании конвективных процессов в мантии Земли указанным переменным уделено недостаточное внимание. Так, например, в статье [13] имеется решение только для постоянной вязкости. В настоящей работе предпринята попытка построения численной модели конвективных процессов в мантии Земли, основанной на вышеупомянутых переменных и методе дробных шагов. В качестве первого этапа анализа работоспособности построенной численной модели осуществлено сопоставление с результатами работы [4].

1. Постановка задачи

1.1. Основные уравнения

Для решения задач тепломассопереноса в мантии Земли привлекается хорошо известная система дифференциальных уравнений в декартовых координатах [14, 15]:

ди дь д'ш

дХ + ду + & = ° (1Л)

*Работа выполнена при поддержке Интеграционных проектов СО РАН №27, 30. © В. В. Червов, 2002.

Р _ дТхх + дТху + дТхг х дх ду дг ’ (1.2)

Р _ дТух + дТуу + дТуг у дх ду дг (1.3)

дТгх дТгу дТгг О + о + о + ардг Т' дх ду дг (1.4)

дТ дТ дТ 2 + V— + т— _ XV Т' дх ду дг (1.5)

др дг дТ ~дЬ

где р(х, у, г, £) — плотность; и(х, у, г, £), у(х, у, г, £), т(х, у, г, £) — компоненты вектора скорости и; р(х,у, г,Ь) — давление; ц(х,у, г,Ь) > 0 — коэффициент динамической вязкости; а — коэффициент теплового расширения мантийной жидкости; Т — температура жидкости; дх — г-компонента вектора силы тяжести g = (0,0, —дх), т. е. ускорение свободного падения; х — коэффициент температуропроводности. Компоненты тензора напряжений задаются следующим стандартным образом [16]:

Тх

ди

'дХ'

ху

п

ди дь\ ду + дх )

Тг.

ди дт\

дг + дх ) '

Ту

уу

2 дV 2ПТТ'

ду

дт

>уг

гу

( ду

+ Ту

Тг

дт

2піт-

дг

1.2. Обезразмеривание

Переменные задачи могут быть обезразмерены с применением масштабных множителей (табл. 1).

Таблица 1

Название и обозначение Масштабный множитель

Длина, ширина и глубина (х, у, г) (х', у', г')б,

Время Ь М /х

Динамическая вязкость ц п'по

Давление р РпоХ/&2

Температура Т Т 'Ид2/рсРх

Компоненты скорости (и,у,,ш) 1з Ъ

Число Рэлея Иа Иа = адг рсРАТ/до х

Здесь Н — тепловыделение радиоактивных изотопов в мантии; ср — удельная теплоемкость при постоянной температуре; й — расстояние от поверхности Земли до ядромантийной границы; АТ = Т0 — разность между температурой на подошве нижней мантии и температурой на поверхности; щ — масштабный множитель для динамической вязкости.

В результате в уравнении (1.4) вместо ардг появится число Рэлея, а в (1.5) “исчезнет” сомножитель х:

др дтгх , дтгу , дтхх

дг

дх

+

дУ

+

дг

— ИаТ,

дТ дТ дТ дТ

~яТ + и------------+ ----+ _ V Т-

ді дх ду дг

(1.4а)

(1.5а)

Из соображений простоты у обезразмеренных переменных оставлены те же обозначения, что и у размерных.

1.2.1. Краевые условия

В рассматриваемой ниже задаче тестирования [4] предполагаются симметричные условия на всех боковых границах:

где Х,У — длина и соответственно ширина вычислительной области (параллепипеда).

На верхней границе температура и компоненты скорости нулевые, на нижней границе температура Т0 = 1 и компоненты скорости нулевые:

2. Описание трехмерных геодинамических течений с применением завихренности и векторного потенциала

Система дифференциальных уравнений, применяемых для расчетов течения мантийного вещества, может быть записана в переменных ш, ф [8, 16, 18].

Завихренность ш определяется из соотношения

дТ ду дт х = 0, 0 < у < У, 0 < г < 1;

дх дх дх , х = X, 0 < у < У, 0 < г < 1;

дТ ди дт у = 0, 0 < х < X, 0 < г < 1;

ду ду ду , у = У, 0 < х < X, 0 < г < 1,

Т = 1, и = у = = 0, г = 0, 0 < х < X, 0 < у < У;

Т = 0, и = у = = 0, г =1, 0 < х < X, 0 < у < У.

1.2.2. Начальные данные

Система уравнений (1.1)-(1.5) устроена таким образом, что при Ь = 0 задаются начальные условия лишь для температуры [17]:

Т(х,у, г, 0) = Т(х,у, г).

В качестве Т в настоящей работе выбиралась функция

Т(х, у, г) = (1 — г) + 0.2(cos(пx/X) + со$,(пу/У)) вт(пг).

(1.6)

ш = Ух И,

которое расписывается покомпонентно следующим образом:

дт ду ди дт ду ди,

ду дг дг дх дх ду

Кроме того, вводится векторный потенциал

ф = гфх + уфу + ,

(2.1)

удовлетворяющий уравнению несжимаемости

У ■ И = 0,

причем

И = Ух ф, (2.2)

т. е.

дф* дфу дфх дф* дфу дфх (2 „)

и ду дг , у дг дх , т дх ду

После подстановки (2.2) в (2.1) получим

Ух (Ух ф) = ш. (2.4)

Если в (2.4) учесть, что векторный потенциал может быть выбран произвольно, например, так, чтобы [8]

Уф = 0,

то из уравнения (2.4) получим векторное уравнение Пуассона

У2ф = ш

или три скалярных:

У2фх = — их, У2фу = — иу, У2ф* = — и*. (2.5)

Применяя операцию “ротор” к уравнениям (1.2)-(1.4) и учитывая, что

их = ту — у*, иу = и* — тх, и* = ух — иу, О = пш,

получим соотношения

У2Сх = Пххих + (пхи% + %иу + П* их) + + ИяТу,

У2Су = Пуу иу + {пх,и% + Пу иу + п* иу) + ^2 — ИаТх, (2.6)

где

_ 2(Пгг ту — Пуу Vz) + 2Щг ^у — Шг) — Пху (иг + Шх) + Пхг ^х — иу); ^2 — 2(пххиг Пгг тх) + 2пхг (шг их) Пуг (vx + иу) + пху (шу vz ); ^3 — 2(пyyvx пххиу) + 2пху (их vy) пхг (шу + vz ) + пуг (иг Шх').

В случае постоянной вязкости запишем

V ш — і ^Иа~д—^ + j ^—Ка~д—^ + к0,

т. е. имеем еще три скалярных уравнения Пуассона (к имеющимся в (2.5)):

дТ дТ

V2иx _ Иа—, V2иy _ — Иа—, VV _0 (2.7)

ду дх

или три уравнения (2.6) в случае переменной вязкости. Причем компоненты вектора завихренности вычисляются простым делением компонент вектора О на п в каждой точке вычислительной области.

Последнее уравнение, замыкающее систему, — уравнение теплопереноса

дТ дТ дТ дТ 2

+ и я---+ v^----+ _ V Т. (2.8)

ді дх ду дг

Таким образом, для решения поставленной трехмерной задачи потребуется интегрирование шести однотипных уравнений Пуассона и вычисление уравнения (2.8) для температуры, которое является общим для любого способа записи основных уравнений.

2.1. Граничные условия для завихренности и векторного потенциала

С учетом краевых условий в п. 1.2.1 граничные условия для функций ш, ф выглядят следующим образом:

на поверхностях 51ед, где х = 0, и 5^^, где х = X (рис. 1), 0 < у < У, 0 < £ < 1:

дфх д„х до

= ФУ = Ф* = 0, = „у = 0, „ = 1Г;

дх дх дх

на поверхностях 5^^, 5ъаск, где у = 0 и у = У (рис. 2), 0 < х < X, 0 < £ < 1:

^ ^ Г = 0, ^ = ЫХ =0; ^ ^

ду ду ду

на поверхностях 5ьо«от, 5;0р, где £ = 0 и £ = 1 (рис. 3), 0 < х < X, 0 < у < У:

^ ^ ^ = Г = =0> „х = _ до „у = ди „ = 0.

дг дг дг дг

(2.9)

(2.10)

(2.11)

Рис. 1. (Зьа, )-плоскости (уг). Рис. 2. (Зьот, Яаск)-плоскости (хг).

р'

/ С!::::::: ой: ——-

X х

Рис. 3. ( б'ьоИот, 6\ор)-ПЛОсКОсТИ (ху).

3. Конечно-разностный алгоритм решения задачи

Для решения поставленной задачи вводилась простейшая равномерная сетка. В прямоугольном параллепипеде П = [0, X] х [0, Y] х [0, Z], который с помощью одномерных сеток

^ : 0 = x < Х2 < ... < Xm = X,

^2) : 0 = Vi < У2 < ... < Vn = Y,

^3) : 0 = zi < Z2 < ... < Zk = Z

разбит на частичные параллепипеды:

nitj = [xi,xi+i] х [y, Vj+i] х [zfc,zfc+i]

(i = 1, 2,..., M — 1; j = 1, 2,..., N — 1; k = 1, 2,..., K — 1), вводились постоянные hx = xi+i — xi;

hy = Vj+i — Vj, hz = zfc+i — Zk. Кроме того, At = ts+i — ts = Tt = const, где индексом s

обозначен текущий слой по времени: ts = sAt (s = 0,1, 2,...).

На примере задачи с переменной вязкостью алгоритм решения задачи сводится к последовательному интегрированию системы уравнений {(2.3), (2.4), (2.6)} на каждом слое по времени. При этом для решения каждого уравнения Пуассона V2П = Ф используется итерационная схема стабилизирующей поправки [19]:

On+1/3 _ O

— ^xx° + LyyO + Lzz

LxxOra+1/3 + LyyOn + LzzOn _ Ф,

Tf

On+2/3 _ On+1/3

Lyy(On+2/3 _ On), (3.І)

Tf

On+1 _ On+2/3

Lzz (On+1 _ On),

т/

где Ф — правая часть уравнения Пуассона; т/ — итерационный параметр, играющий роль фиктивного времени. Итерационный процесс останавливается при условии

max

On+1 _ On

On+1

< є, (3.2)

где £ > 0 — некоторая малая положительная константа, определяемая в процессе вычислений; верхний индекс в формулах (3.1), (3.2) обозначает номер итерационного шага.

Уравнение переноса тепла (2.6) интегрируется также с применением схемы стабилизирующей поправки:

Т"+1/3 _ Т" - + иХХТ5+1/3 _ £ХХТ5+1/3 = _оЬУТ5 + ЬууТ5 _ ^Т5 + ЬггТ",

Tt

T«+2/3 Ts+1/3

Tt

Ts+1 ts+2/3

+ vLyTs+2/3 _ LyyTs+2/3 = vLyTs _ LyyTs, (3.3)

+ wLzTs+1 _ LzzTs+1 = wLzTs _ LzzTs,

Tt

где Tt — величина шага по времени.

В схемах (3.1) и (3.3) трехточечные разностные операторы Ьхх, Ьуу, и Ьх, Ьу,

аппроксимируют дифференциальные следующим образом:

д2/ = т г = Уг+1 _ 2/ + /i-l

dx2 = xx/ =

hx

d/ = L / = ji+l _ /i-l dx = x/ = 2hx ■

Для остальных операторов аппроксимации подобны (3.4).

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

1) По закону (1.6) вычисляется начальное распределение температуры.

2) Вычисляется в случае переменной вязкости (в соответствии с работой [4]) п(Т) = ехр[в/(Т + 0) _ 0/(0.5 + 0)], причем в = [225/ 1п(г)] _ 0.251п(г) и, кроме того, г =

п|Т=о/ п|Т =1 = 2°.

3) Вычисляются правые части уравнений (2.6) или (2.7).

4) Из уравнений (2.6) или (2.7) методом установления определяются компоненты вектора завихренности.

5) Из уравнений (2.5) методом установления определяются компоненты векторного потенциала.

6) Вычисляются компоненты вектора скорости по формулам (2.3).

7) Вычисляются там, где это необходимо, граничные условия для завихренности по известной схеме усреднения (Полежаев и др. [20]), которая позволяет повысить устойчивость вычислительного алгоритма:

и

n+1

Здесь у — параметр релаксации, причем у G (0,1), а F(^n+i) — правая часть соотношений, которые аналогичны условиям Тома [21] для двумерного случая: на поверхности (yz), x = const:

2

Ul,j,k = — h2 (^2,j,k), i = 1,

hx

2

WM,j,k = — h2 (^M — i,j,k) , i = M;

на поверхности (xz), y = const:

4z,l,fc = _ —2 (vz,2,fc), j = І,

hy

2

Uz,N,fc = _ —2 (Vz,N -l,fc ), j N;

hy

на поверхности (xy), z = const:

22

^j,! = _ -2 (Vx,j,2), Uy,j,l = _-2(Vy,j,2), k = І,

2

2

и

i,j,K — -2 (Vx,j,K-1) , Uij,K = -2 (Vi,j,K-l), k = K

Остальные условия от векторного потенциала не зависят.

x

8) Проверяется по формуле (3.2) сходимость для векторного потенциала; если условие (3.2) при некотором е > 0 не выполняется, то процедура вычислений повторяется (внутренние итерации).

9) На основе решения уравнения (2.8), определяется распределение температуры внутри расчетной области.

10) Вычисляются вспомогательные параметры задачи, а также величины, необходимые для тестирования (см. п. (1)-(у1) следующего параграфа), после чего начиная с п. 2 вычисления повторяются.

4. Результаты численных расчетов

Поскольку основной целью работы является сопоставление с результатами Биээе е! а1. [4], в задаче вычислялись следующие параметры:

(і) число Нуссельта (Ми) по формуле БІапкепЬаеЬ [22]

Ми = — (ХУ) JJ Тііхііу

^ор

и среднеквадратичная скорость

\

где А — объем параллепипеда со сторонами X, У и Z;

(II) значение вертикальной компоненты скорости ,ш и температуры Т в угловых точках среднего сечения конвективного слоя;

(III) значение теплового потока $ = — дТ/ дг в некоторых точках области;

(1у) интегральный параметр, вычисляемый по формуле

т (х,г) = / о

в точках (0, 0.25), (Ь/2, 0.25), (Ь, 0.25) интегрированием вдоль оси У;

(у) средняя температура Тт, вычисляемая на горизонтальном сечении области на глубине г = 0.75;

(у1) в случае переменной вязкости значение вертикальной компоненты вектора завихренности ^ в точке (0.75, 0.25, 0.75).

Размерные значения (в системе СИ), которые были использованы в [4] и в настоящей работе, принимались следующими:

& = 2 700 000, АТ = 3700, х = 10-6, а = 10-5, р = 3300, дг = 10.

4.1. Постоянная вязкость

Расчеты проводились для следующих значений параметров: длина области X = 1.0079, ширина У = 0.6283, высота Z = 1.000. Масштабный множитель при вязкости п0 = 8.0198 • 1023. Число Рэлея И,а = 30 000.

Система, которая решается при постоянной вязкости, состоит из уравнений (2.5), (2.7) и уравнения переноса тепла (2.8). Граничные условия в точности соответствуют условиям

(2.9) - (2.11).

Следует заметить, что для задачи в представленном виде разноименные компоненты векторов ш и ф не связаны между собой ни через соотношения (2.5), (2.7), (2.3), ни через граничные условия. Благодаря этому обстоятельству можно разбить процедуру вычислений на три независимых потока, что значительно сократит время расчетов на многопроцессорной ЭВМ.

Результаты. В работе Busse [4] наиболее полно представлены результаты Кристенсена (Chr), которые можно принять за эталонные. Сопоставление результатов расчетов Кристенсена и автора (Che) представлено в табл. 2. Соответствие достаточно хорошее.

Т а б л и ц а 2

Векторный потенциал — вектор завихренности. Постоянная вязкость

№ Наименование параметра Результаты Кристенсена на сетке 32 x 32 x 64 Результаты автора на сетке 32 x 32 x 32 Относительная ошибка, %

1 Nu 3.Б44 3.Б77 0.9312

2 ^rms 41.00 41.41 1.0000

3 w(0,0,0.5) 11б.б 121.2 3.94Б1

4 w(0,S,0.5) 40.Б0 40.00 1.234б

Б T(0,0,0.5) 0.В013 0.В099 1.0733

б T(0,S,0.5) 0.б188 0.б174 0.22б2

7 т(0,0) -0.3647 -0.3648 0.0274

В т (L/2,0) -0.1292 -0.1196 7.4303

9 ,0) (L, Т( -0.1104 -0.1052 4.7101

10 Tm(0.75) 0.Б21Б 0.Б197 0.34Б2

Относительная ошибка вычислялась по формуле

Err

Che - Chr

Chr

100%.

(4.1)

Близкие результаты были получены и на сетке (24 x 24 x 24).

4.2. Переменная вязкость

Расчеты проводились для следующих значений параметров: длина области X = 1.00, ширина У = 1.00, высота X =1.00. Масштабный множитель при вязкости: По = 1.20165 • 1024. Число Рэлея И,а = 20 000, г = п1т=0/п1т=1 = 20-

Система, которая решается при переменной вязкости, состоит из уравнений (2.5), (2.6) и уравнения переноса тепла (2.7). Граничные условия в точности соответствуют условиям

(2.9)-(2.11).

Кроме этих шести уравнений Пуассона в расчетах используются соотношения (2.3) для вычисления вектора скорости.

Результаты. Как и в предыдущем случае с постоянной вязкостью, полученные автором результаты расчетов с вязкостью, зависящей от температуры, сопоставлялись с результатами расчетов Кристенсена (табл. 3). Установлено их хорошее соответствие.

Таблица 3

Векторный потенциал — вектор завихренности. Переменная вязкость

№ Наименование параметра Результаты Кристенсена на сетке 32 х 32 х 64 Результаты автора на сетке 32 х 32 х 32 Относительная ошибка, %

1 3.03927 3.040 0.0240

2 35.132 35.18 0.1366

3 ■ш(0,0,0.5) 165.91 167.2 0.7775

4 ад(0,Я,0.5) -26.72 -26.43 1.0853

5 ад(Ь,Я,0.5) -58.23 -58.23 0.0000

6 Т(0,0,0.5) 0.90529 0.9067 0.1558

7 Т (0Д0.5) 0.49565 0.4965 0.1715

8 Т (ЬД0.5) 0.23925 0.2379 0.5643

9 0(0,0) 5.83390 5.8121 0.3737

10 0(0,8) 1.71360 1.7322 1.0854

11 0(Ь,Я) 0.7684 0.7731 0.6117

12 т(0,0) -0.5059 -0.5141 1.6209

13 5? (0, т( -0.1921 -0.1891 1.5617

14 т (Ь,я) -0.1388 -0.1356 2.3055

15 Тт(0.75) 0.56593 0.5636 0.4117

16 Тт(0.50) 0.58158 0.5799 0.2889

17 0(0.75,0.25,0.75) -11.125 -11.25 1.1236

Результаты расчетов сравнивались с данными Кристенсена как наиболее полными из имеющихся в статье [4]. Ошибка вычислялась по (4.1). Аналогичные результаты получены на сетке (48 х 48 х 48).

Выводы

Полученные результаты позволяют сделать оптимистический вывод о жизнеспособности подхода с использованием переменных “завихренность” — “векторный потенциал” для задач тепломассопереноса в мантии Земли.

Автор благодарит Г. Г. Черных и С. А. Тычкова за постоянное внимание к работе.

Список литературы

[1] Рыков В. В., Трувицин В. П. Трехмерная модель мантийной конвекции с движущимися континентами // Вычисл. сейсмология. 1994. Вып. 27. С. 21-41.

[2] Рыков В. В., Трувицин В. П. Численное моделирование трехмерной мантийной конвекции и тектоника литосферных плит // Там же. Вып. 26. С. 94-102.

[3] GABLE W., O’Connell J. Convection in three dimensions with surface plates: generation of toroidal flow // J. Geophys. Research. 1991. Vol. 96, No. B5. P. 8391-8405.

[4] BussE F.H., Christensen U., Clever R. et al. 3D Convection at infinite Prandtl number in cartesian geometry — a benchmark comparison. Geophiys. Astrophys // Fluid Dynamics. 1993. Vol. 75. P. 39-59.

[5] GLATZMAIER G. A., Schubert G. Three-dimensional spherical model of layered and whole mantle convection // J. Geophys. Research. 1993. Vol. 98, No. B12. P. 21969-21976.

[6] MACHETEL P., Thoravan C., Brunet D. Spectral and geophysical of 3-D spherical mantle convection with an endothermic phase change at the 670 km discontinuity // Phis. Earth and Planetary Interiors. 1995. Vol. 88. P. 43-51.

[7] Zhong S., ZUBER M. Role of temperature-dependent viscosity and surface plates in spherical shell models of mantle convection // J. Geophys. Research. 2000. Vol. 105, No. B5. P. 11063-11082.

[8] РОУЧ Х. Вычислительная гидродинамика. М.: Мир, 1980.

[9] АНДЕРСОН Д., ТАННЕХИЛЛ Дж., ПЛЕТЧЕР Р. Вычислительная гидромеханика и теплообмен. М.: Мир, 1990.

[10] ФЛЕТЧЕР К. Вычислительные методы в динамике жидкостей. Т. 1, 2. М.: Мир, 1991.

[11] БЕССОНОВ О. А., БРАйЛОВСКАЯ В. А., ПОЛЕЖАЕВ В. И. Пространственные эффекты конвекции в расплавах: концентрационные неоднородности, возникновение несимметрии и колебания // Механика жидкости и газа. 1997. №3.

[12] Aziz K., Hellums J. D. Numerical solution of the three-dimensional equations of motion for laminar natural convection // Phys. Fluids. 1967. Vol. 10. P. 314-324.

[13] Sparks D.W., Parmentier E. M., Morgan J.P. Three-dimensional mantle convection beneath spreading centers spreading centers Implications for along-axis variations in crustal thickness and gravity // J. Geophys. Research. 1993. Vol. 98, No. B12. P. 21977-21995.

[14] ТЫЧКОВ С. А. Конвекция в мантии и динамика платформенных областей. Новосибирск: Наука. Сиб. отд-ние, 1984.

[15] ДОБРЕЦОВ Н. Л., КИРДЯШКИН А. Г. Глубинная геодинамика. Новосибирск: Изд-во Сиб. отд-ния РАН, НИЦ ОИГГМ СО РАН, 1994.

[16] ЛОйЦЯНСКИй Л. Г. Механика жидкости и газа. М.: Наука, 1978.

[17] ФЕДОРЮК М. В. Характеристики течений несжимаемой жидкости в гравитационном поле // Мат. сб. 1988. №137(179), 4(12).

[18] АРФКЕН Г. Математические методы в физике. М.: Атомиздат, 1970.

[19] ЯНЕНКО Н. Н. Метод дробных шагов решения многомерных задач математической физики. Новосибирск: Наука. Сиб. отд-ние, 1967.

[20] ПОЛЕЖАЕВ В. И., ГРЯЗНОВ В. Л. Метод расчета граничных условий для уравнений Навье — Стокса в переменных вихрь, функция тока // Докл. АН СССР. 1974. Т. 219, №2.

[21] ТОМ А., ЭйПЛТ К. Числовые расчеты полей в технике и физике. М.: Энергия, 1964.

[22] BLANKENBACH B., BussE F. ET AL. A benchmark comparison for mantle convection codes // Geophys. J. Int. 1989. Vol. 98. P. 23-38.

Поступила в редакцию 2 октября 2001 г., в переработанном виде — 24 октября 2001 г.