Скорости сейсмических волн из пробирки

Bridging the Gap

девиз DARPA

Самая безнадёжная, а потому и самая увлекательная задача при работе с какой-либо моделью (набором уравнений) – это её проверка всеми возможными методами. 

Я занимаюсь, как я писал в прошлых постах (крайний и далее по ссылкам), составлением компьютерной модели плавления вещества в условиях мантии Земли. На каждом этапе необходимо проверять, насколько сильно я вру. Один из неиспользованных мною до сих пор источников данных – это профили скоростей сейсмических волн. Поскольку для Земли есть подобные профили (та же модель PREM), и это единственный прямой источник данных о её структуре, его использовать необходимо. 

Что может дать данное сравнение? Первое – это дополнительная проверка, что вся модель запрограммирована верно. Второе – это выход на какие-то экспериментально измеряемые величины.

Итак, скорости объемных сейсмических волн. Их две - продольные \(V_p\)  и поперечные \(V_s\) (обозначения от слов primary и secondary - в соответствии с порядком, в котором они приходят на сейсмодатчик). Скорости волн определяются как:

\[V_p = \sqrt {\frac {K_s + \frac{4}{3}G}{\rho}} \]

\[V_s = \sqrt {\frac {G}{\rho}} \]

\[V_p > V_s \]

Как видно из формул, они зависят от модуля упругости \(K_S\) и модуля сдвига \(G\) и плотности \(\rho\).

Кстати, про плотность. Тут забавная история. Как нас учили на курсе общей геологии на первом курсе, скорость сейсмических волн зависит от плотности и увеличивается с её ростом. Нас так же продолжили учить и на кратком курсе геофизики на четвёртом курсе, но тут кто-то спросил “а как же тогда скорость сейсмических волн может увеличиваться с глубиной, если с глубиной плотность пород растёт, а она в знаменателе?”. И что я помню совершенно отчётливо – это как зам.декана Степанов П.Ю. был озадачен этим вопросом и отделался фразой “ну модули растут быстрее” без каких-либо пояснений. Искать ни на первом, ни на четвёртом курсе ответ времени у меня не было, но проблема в мозгу осталась. Теперь я могу дать ответ сам: потому что скорость сейсмических волн от плотности, как следует из формул, не зависит!

Применим Истинное Знание (термодинамику) к модулю упругости. Уравнение для него распишется как функция теплоёмкостей (\(C_p\), \(C_v\)), плотности \(\rho\), параметра Грюнайзена \(\gamma\) и коэффициента термического расширения \(\alpha\):

\[ K_S = \gamma \frac{C_p \rho}{\alpha} = \rho \frac{C_p - C_v}{\alpha^2 T} \frac{C_p}{C_v} \]

Остаётся вторая величина - модуль сдвига. Вот она уже не выводится из таких простых формул. Но. Можно сделать фокус. Есть немало работ (напр. [Ledbetter, 1977]), которые показывают, что для изотропных и около-изотропных сред (включая кристаллы кубической сингонии) есть определённое соотношение модуля упругости и модуля сдвига:

\[ G = K_S \left( 9\left(\frac{G}{E}\right) -3\right) = k K_S\]

где \( \left(\frac{G}{E}\right) \), отношение модуля сдвига (shear modulus) к модулю Юнга, является величиной почти постоянной. Теория предсказывает значение 0.375, из экспериментов следует 0.391.  Таки образом,

\[ G = 0.375 \div 0.519 K_S \]
Хорошая это оценка или плохая? Расчёт для мантии Земли по модели PREM даёт значения в пределах этого диапазона. Более того, мы знаем, что мантия Земли состоит из минералов с кубической или около того структурой. Мы можем использовать такие модели, если нас устраивает погрешность в первые проценты - десять процентов. В качестве первого приближения это прекрасный результат.

Итак, обоснование, почему скорость сейсмических волн напрямую не зависит от плотности (на самом деле, немножко сложнее, поскольку \(\alpha\) есть производная логарифма удельного объёма, но это уже высшие сферы):

\[V_p = \sqrt {\frac {1}{\rho} \left( K_S + \frac{4}{3} G \right)} = \sqrt {\frac {1}{\rho} \left( K_S + \frac{4}{3} kK_S \right)} = \sqrt {\frac {1}{\rho} \gamma \frac{C_p \rho}{\alpha} \left( 1 + \frac{4}{3} k \right)} \]

Или

\[V_p = \sqrt {\gamma \frac{C_p }{\alpha} \left( 1 + \frac{4}{3} k \right)}= \sqrt {\frac{C_p - C_v}{\alpha^2 T} \frac{C_p}{C_v}\left( 1 + \frac{4}{3} k \right)} \]

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

Стоит отметить, что модуль объёмного сжатия убывает с ростом количества пор (в т.ч. заполненных жидкостью), например, [Russell, Smith, 2007]. Модуль сдвига падает по степенной функции до нуля примерно при 40% жидкости (напр. [Kováčik, 2008]). 
Модуль объёмного сжатия зависит от текстуры и структуры агрегата зёрен. По этому вопросу есть замечательная статья [Man, Huang, 2011], называющаяся "простая явная формула ...", от которой у меня начинала болеть голова примерно на третьей странице с суммами тензоров.

Теперь же о том, как всё это внедрялось. На этом графике рыжая линия - это то, что должно быть, синяя - это то, что получилось по старой методике расчёта, зелёная - по новой методике расчёта.

В чём разница?

В "голубой" модели я использовал простую старую простую параметризацию:

\[ V = (a_1T+a_2)P^2 + (b_1T + b_2)P + (c_1T+c_2) \]

\[ E = (d_1T+d_2)V^2 + (e_1T + e_2)V + (f_1T+f_2) \]

Как я уже писал в прошлых постах, эта схема неплохо описывала многие параметры для мантии типа её плотности и вполне вписывалась в допустимые диапазоны по теплоёмкости и коэффициенту термического расширения.

Почему она "сломалась" и стала рисовать эту странную голубую параболу?

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

Собственно, дифференцирование полинома дало параболу, зависящую от давления. Расчёт модуля упругости "усилил" эту зависимость, а окончательно она нарисовалась при расчёте скоростей сейсмических волн.

Когда я это понял, я поплакал в ванной, после чего потратил день и заменил параболы для объёма и энергии на уравнения Birch-Murnaghan (о них я немного писал, подробнее есть в той же википедии). И мы сразу получили "зелёный" вариант, который я считаю невероятно крутым результатом.  Особенно если учесть, что исходные коэффициенты мы берём из экспериментов над миллиметровыми количествами вещества или вообще в численных опытах на компьютере, а в итоге можем нарисовать реалистичную картину тех недр планет, которых мы никогда в жизни не видели и в руках не держали. Погрешность менее 5% при переходе через 9 порядков масштаба - это ничто.

Почему я сразу не сделал хорошо?

Интересный вопрос. Во-первых, в корректном решении новых уравнений есть много своих подводных камней, начиная с корректной итерационной схемы и кончая подходом Mie-Gruneisen (см. прошлый пост). Когда я начинал работать над моделью, я про них не знал. Это сейчас я потратил день. Тогда я бы потратил бы неделю-другую-третью, и не факт, что успешно. И это было правильное решение - начать с простого, перейти к сложному. Сам себя не похвалишь - весь день как оплёванный.

Мораль сей басен такова:

1. Всё ложь. Если в уравнении [сейсмических волн] нарисована буковка [плотность], это ещё не значит, что она имеет значение.
2. Всё ложь. Из того, что значения полинома описывают набор значений некоторого параметра, не следует то, что производная этого полинома будет описывать производную этого параметра.
   К слову, частично эта проблема решается сплайнами (см. напр. тут).
3. Мы смело можем шагать из пробирки в планеты через 9 порядков масштаба, получая ошибки менее 5%. И это можно считать действительно потрясающим результатом.