CC BY
104
УДК 550.311
Быстрое численное моделирование соляной тектоники: возможность оперативного использования в геологической практике
Б.В. Лунев, В.В. Лапковский
Институт нефтегазовой геологии и геофизики им. А.А. Трофимука СО РАН, Новосибирск, 630090, Россия
В статье представлены примеры применения программы расчета развития неустойчивости Рэлея-Тейлора для моделирования соляной тектоники. Благодаря использованию точного фундаментального решения соответствующей краевой задачи удалось достичь высокой эффективности работы программы при расчете различных вариантов эволюции сложных моделей, соответствующих реальным геологическим разрезам. Это делает программу пригодной к широкому использованию при решении практических задач структурной и поисковой геологии в комплексе с другими геофизическими методами.
: соляная тектоника, неустойчивость Рэлея-Тейлора, численное моделирование
Fast numerical simulation of salt tectonics: possibility of the operational application in geological practice
B.V. Lunev and V.V. Lapkovskii
Trofimuk Institute of Petroleum Geology and Geophysics SB RAS, Novosibirsk, 630090, Russia
In the paper the calculation program for Rayleigh-Taylor instability development is used to simulate salt tectonics. Owning to the exact fundamental solution of a given boundary-value problem the program proves to be highly efficient and allows calculating different evolution variants of complex models corresponding to real geological sections. This renders the program suitable for wide use in solving practical tasks of structural and exploration geology in a combination with other geophysical methods. : salt tectonics, Rayleigh-Taylor instability, numerical simulation
Развитие соляных куполов и связанных с ними высокоамплитудных, контрастных структур в осадочных чехлах платформ — явление само по себе любопытное для геологов. Интерес этот сильно подогревается сугубо практическими соображениями, поскольку соляная тектоника развита в большей части нефтегазоносных провинций и зачастую непосредственно контролирует размещение залежей углеводородов как структурообразующий фактор.
В связи с этим, соляные структуры и механизм их образования и развития давно и интенсивно изучались как в полевых условиях, так и методами физического и математического моделирования [1-26]. Однако несмотря на значительные успехи математическое моделирование, также как и физическое, пока остается скорее инструментом изучения общих закономерностей
© Лунев Б.В., Лапковский В.В., 2009
соляной тектоники, а его использование в целях практического исследования конкретных геологических структур ограничено. Дело в том, что моделирование какой-либо конкретной геологической структуры требует расчета, по крайней мере, нескольких альтернативных моделей эволюции с различными начальными условиями, чтобы по реально проверяемым результатам можно было выбрать оптимальный вариант. Поэтому важным требованием к программе моделирования для ее использования в практической геологии является быстродействие. Между тем, медленная сходимость вычислительного процесса при решении дифференциальных уравнений движения (как уравнений Навье-Стокса в случае ньютоновской среды, так и уравнений баланса сил общего вида в комбинации с используемыми реологическими уравнениями) приводит к тому, что даже при
использовании многопроцессорных компьютеров и современных вычислительных методов моделирование ситуаций, близких к реальным геологическим (с большим количеством сложных границ разноплотных слоев и замкнутых тел), требует весьма длительных расчетов для каждой модели.
Благодаря использованию аналитического решения соответствующей задачи для уравнения Навье-Стокса, полученного ранее одним из авторов настоящей работы [27], мы смогли обойти указанные трудности медленной сходимости вычислительного процесса и разработали высокоэффективную компьютерную программу для моделирования развития неустойчивости Рэлея-Тейлора в вязкой жидкости. Программа позволяет даже на обычных персональных компьютерах быстро рассчитывать эволюции сложных моделей с большим количеством слоев и изолированных тел разной плотности и произвольной формы. Это дает возможность оперативного применения численного моделирования соляной тектоники при исследовании реальных геологических структур.
Платой за быстродействие в данном случае являются следующие ограничения: 1) среда рассматривается как несжимаемая вязкая ньютоновская жидкость; 2) вязкость считается константой для всего разреза. Рассмотрим их по отдельности.
Согласно известной теореме Нолла [28], модель ньютоновской вязкой жидкости является (во вполне строгом смысле) первым приближением общего реологического уравнения для необратимо деформирующихся сред. Это приближение адекватно описывает деформирование материала до тех пор, пока оно остается достаточно медленным. Степень «медленности» определяется для данного материала относительно глубины его памяти, характеризуемой его так называемым «естественным временем» . Течение может считаться медленным и корректно описываться уравнениями На-вье-Стокса, если скорость деформации не превосходит некоторой критической величины порядка -1. Для кристаллических горных пород может быть оценено величиной порядка (1-5)-106 лет [29], соответственно критическая скорость деформации оценивается величиной порядка 10-14 с-1. Вероятно, для осадочных пород критическая скорость деформации будет не меньше указанной (а естественное время — не больше). Применительно к процессам соляной тектоники средняя скорость деформации в процессе роста диапира определяется как 10-14-2 -10-15 с-1 [3, 7, 9]. Отсюда следует применимость ньютоновской модели среды для описания процессов соляной тектоники. Даже если в каких-то зонах скорости деформаций могут превышать критические значения, модель, построенная в рамках ньютоновской реологии, будет в целом правильно описывать процесс. Для аномальных зон (если такие имеются) это решение первого приближения может быть уточнено
введением реологического уравнения специального вида. (Конечно, следует помнить, что, помимо прочего, адекватность ньютоновского описания, как и вообще концепции сплошной среды, находится в прямой связи с соответствующим осреднением характеристик по времени и пространству.)
Что касается принятого в настоящей работе предположения о постоянстве вязкости среды, то оно опирается на следующие факты и соображения. При использовании ньютоновской модели как первого приближения для медленных течений фигурирующая в ней вязкость («естественная вязкость») должна пониматься как верхняя асимптота эффективных вискозиметричес-ких вязкостей данного материала [28]. Несмотря на то, что имеющиеся в литературе оценки эффективной вязкости осадочных пород (особенно соли) варьируют в довольно широких пределах (в зависимости от способа оценки и скорости деформации), их верхние асимптоты практически совпадают и имеют порядок 1020 Па-с (за исключением сильно обводненных слоев в зонах аномально высокого пластового давления) [9, 14, 30]. Малые различия (в пределах порядка) в оценке этих асимптотических значений, по-видимому, не следует принимать во внимание, т.к., во-первых, они недостаточно достоверны (в силу понятной проблематичности определения реологических параметров для геодинамических процессов), а во-вторых, согласно численным экспериментам (например [14-16]) такие различия вязкости слоев практически не оказывают влияния на эволюцию. Напротив, значительные различия вязкости всплывающего и перекрывающего слоев (два порядка и более) приводят в экспериментах (как численных, так и физических) к существенным различиям в форме зрелых диапиров: в случае большей вязкости всплывающего слоя формируется «палец», в случае его меньшей вязкости — «пузырь», а в случае близких вязкостей — «гриб» [5-10, 12, 31]. Реальные соляные диапиры на зрелых стадиях там, где они изучены, имеют форму «гриба» [3, 9], что также свидетельствует о малых различиях вязкости внутри системы. Поэтому предположение о постоянной вязкости в данном случае представляется обоснованным. В расчетах будем основываться на оценке ее асимптотического значения порядка 1020 Па-с.
Еще одной важной особенностью исследуемых течений является заведомая малость характеризующего их числа Рейнольдса. (Полагая вязкость 1020 Па - с, характерный размер 104 м, характерную скорость течения 10-9 м/с и плотность 2500 кг/м3, число Рейнольдса Re определится величиной порядка 10-22, т.е. можно считать, что Re 0.) Соответственно уравнение движения принимает стационарную форму. Такие течения, целиком определяемые в каждый момент времени актуальной конфигурацией объемных сил, границ и граничных условий, известны как «ползущие течения». Эволюция ползущего течения представляется как после-
малого промежутка времени, выбранного с учетом полученных значений скорости течения и сложности границ, интегрируется уравнение сохранения массы — определяется новая конфигурация границ. Таким образом получается новое распределение плотности, которое вновь сворачивается с фундаментальным решением для получения нового поля скорости течения и т.д. В результате имеем эволюцию разреза. Поскольку в нашем случае нет нужды каждый раз производить расчет по сетке, покрывающей всю область вычислений, то густота задания точек вычисления на границах тел может быть сколь угодно высокой и изменяться в процессе счета в зависимости от изменения геометрии границы и разности скоростей в соседних точках. Это позволяет добиваться высокой точности расчета движения границ за разумное время. Кроме того, программа позволяет выполнять по заданной сетке расчет полей вектора скорости течения среды, тензора скорости деформации и тензора негидростатического напряжения, а также возмущения свободной границы (рельеф поверхности) и аномалии гравитационного поля. Эти вычисления могут производиться для любого выбранного момента времени в эволюции разреза, для всей исследуемой области или любой ее части.
Использование полученного аналитически фундаментального решения дает важные преимущества при
построении программ расчета ползущих течений, обусловленных действием сил Архимеда, вообще, и моделирования процессов соляной тектоники, в частности:
1. При численном решении стационарного уравнения Навье-Стокса объем вычислений растет пропорционально кубу числа точек, в которых рассчитывается поле течения [20]. В данном же случае он пропорционален произведению первой степени этого числа на число точек, в которых задана аномальная плотность.
2. Кроме того, при численном решении точки, в которых рассчитывается поле течения, и точки, в которых задана плотность, должны быть более-менее равномерно распределены по области расчета. В данном же случае такое распределение необходимо только при плавном («градиентном») распределении плотности, например при расчете тепловой конвекции. В обычных при моделировании соляной тектоники ситуациях, когда плотность в пределах каждого выделенного слоя (или замкнутого тела) постоянна, для расчета эволюции достаточно вычислять скорость только на границах слоев (тел). Аномальная плотность также задается не везде, а лишь там, где она отлична от нуля, т.е. в том случае, если на данной глубине она меняется по латерали (по координатам Хх, Х2). Соответственно при той же точности расчета в данном случае можно обходиться меньшим числом точек либо, сохраняя число точек, можно
Ух,ю-13 м/с
3.5
Уг,1СГ13 м/с
0.0
0.0 4.0
1 МПа
км I
8.0 км -3,5 41. Е
13800
_ 1
1
14.0
0.0
4.0
км I
,0 км -89800
и ■ г
ГГ к
н= к
Т„, МПа
0.0
4,0
138
ш ЛЯ ' я
■ к
0.0 4 .0 1 1 8
|1(х), Ю-2 м
— 0.0
1
4.0
КМ
КМ -30.8
4.39
J 0.0 4.0 КМ 1
0.0 4.0 8.0 км -0.545
Т„, МПа
Л ! 1
■м
■
8.26
0.0 4.0 8.0 км -8.23
д(х), 10"6 м/сг 1.06
0.0 4.0
.0 км -4.36
Рис. 1. Поля горизонтальной ( ) и вертикальной ( ) компонент вектора скорости, скорости вертикальной деформации ( ), горизонтального ( ), вертикального ( ) и касательного ( ) напряжений, возмущения давлений ( ), поверхности ( ) и гравитационного поля ( ), возбуждаемые точечной аномальной массой
и
— ш ■ щцщ ||| ™
0 2 4 в: 6 3
4.19 ■ 10"11 М/С
J
г!
тО.О
НО
I
0.0
4,0
км ____
КМ -4.21 ■ 10~11 М/С
' Окм
' о км
3.22 ■ КГ11 М/С
J
-,0.0
4.0
0.0
4.0
8.0
I
км -4.66 ■ 10-11 м/с
Рис. 2. Классический пример развития неустойчивости Рэлея-Тейлора: 3.5 ( ), 5.8 ( ), 7.3 ( ), 10.1 млн. лет ( ). Разница плотностей всплывающего и перекрывыющего слоев — 0.5 - 103 кг/м3. Поля горизонтальной ( ) и вертикальной ( ) компонент вектора скорости (стадия )
гораздо подробнее описывать границы, что важно с учетом неустойчивости некоторых из них. При этом распределение точек можно оптимизировать, сгущая их там, где форма границы сложнее и выше градиенты скорости течения.
3. Известно [20], что при решении задач рассматриваемого класса на многопроцессорных компьютерах рост эффективности параллельных вычислений падает с увеличением числа процессоров, причем падает быстро. Приведенный в цитированной работе график зависимости эффективности вычислений от числа процессоров имеет логарифмический вид. Это связано с увели-
чением времени на пересылку данных между процессорами. Благодаря тому что в нашем случае используется фундаментальное решение, полученное для точечной единичной силы, его применение может быть особенно эффективным на многопроцессорных машинах, поскольку допускает естественное неограниченное «распараллеливание» задачи, практически без обменов между процессорами. В результате эффективность параллельных вычислений будет расти линейно пропорционально числу процессоров, как в идеальном случае.
Таким образом, использование фундаментального решения позволяет кардинально увеличить производи-
12 16 20 24 28 32 36 40 Рис. 3. Развитие неустойчивости «глубокого слоя»
км
тельность компьютерных программ, моделирующих развитие гравитационной неустойчивости в жидкости. Особенно при расчете больших задач, с большим количеством границ — в этом случае производительность может быть повышена на многие порядки.
Отметим также, что использование решения для сосредоточенной единичной силы позволяет гибко оперировать с любым распределением плотности.
Для обеспечения возможности моделирования на обычных персональных компьютерах была создана программа Диапир-2, которая реализует двумерный вариант решения данной задачи. Двумерные решения получаются интегрированием трехмерных вдоль оси 2. Кроме того, в данном варианте программы в качестве ядер свертки используются не сами двумерные решения, а их определенные интегралы вдоль оси 1, т.е. распределение плотности задается не точками, а горизонтальными отрезками. Плотность каждого выделенного тела или слоя считается постоянной. Программа имеет удобный интерфейс, позволяет задавать исходную модель разреза, просматривать результаты расчетов (эволюцию разреза и поля, характеризующие его напряженно-деформированное состояние), а также редактировать существующую модель на любом этапе ее эволюции. Ниже представлены некоторые результаты работы этой программы.
На рис. 1 показаны поля вектора скорости течения, тензоров деформаций и напряжений, а также рельеф свободной поверхности и аномалия гравитационного поля, возбуждаемые точечной (почти точечной) отрицательной аномальной массой. Подъем вещества над аномальной массой сопровождается его растеканием в стороны и погружением на флангах. С подъемом связана область вертикального сжатия (горизонтального растяжения) и поднятие свободной поверхности. На флангах так называемые компенсационные депрессии, где областям погружения соответствуют зоны вертикального растяжения (горизонтального сжатия), выражены в рельефе «компенсационными прогибами». Эти принципиальные особенности течения видны и в моделях эволюции, получая там соответствующее структурное выражение.
Обратимся к рис. 2, где представлен классический пример, всегда рассчитываемый при моделировании развития неустойчивости Рэлея-Тейлора. Показано развитие неустойчивости километрового слоя плотностью 2 000 кг/м3, перекрытого трехкилометровой пачкой слоев плотностью 2500 кг/м3. Возмущение проходит стадии «подушки», «пальца», превращается в зрелый диа-пир, достигает свободной поверхности и растекается вдоль нее, образуя характерные свесы на краях «шляпки» (обусловленные нисходящими ветвями течения). Для удобства представления во всплывающем слое показана деформация двух систем материальных линий. В скорости роста возмущения наблюдается известная
' ' ■ ■ |__^ ыыы
0 100 200 0 100 2000 '00 200 'км
Рис. 4. Воспроизведение эксперимента из [8, с. 282, рис. 11.93]
стадийность (см., например, [9]): очень медленный рост до формирования «подушки» (рис. 2, ), затем ускоренный (в 2, 3 и до 10 раз) рост в фазе развития «подушки» (рис. 2, ) до начала формирования диапира; небольшое (примерно на 10 %) замедление и установление постоянной скорости роста в фазе «пальца» и на начальных этапах образования «бульбы» (рис. 2, ) и замедление на финальных стадиях при подходе «созревающего» диапира к поверхности (рис. 2, ). Абсолютные значения скорости роста зависят от принятого коэффициента вязкости, контраста плотностей, глубины залегания всплывающего слоя и размеров возмущения.
Те же этапы можно видеть на рис. 3, где представлено развитие неустойчивости глубокого слоя, также обычно рассчитываемое в работах по моделированию диапиризма. Максимальная скорость роста возмущения достигается здесь на начальной стадии формирования «шляпки» диапира (при достижении им глубины 2 км) и составляет при вязкости 1020 Па - с величину 3 м за 10000 лет, что соответствует данным исследования истории реальных соляных диапиров в подобной ситуации [9]. (Можно заметить, что представленные результаты моделирования повторяют для случая постоянной вязкости классические расчеты [5].)
На рис. 4, 5 воспроизведены результаты физических экспериментов, представленных в книге X. Рамберга
во-------------
км
км
'з,20 '^ёб 3.00 Г2~5о' " 2Ю0 270 "
Рис. 5. Воспроизведение эксперимента из [8, с. 336-337, рис. 15.4, 15.5]
Рис. 6. Результат физического эксперимента по книге X. Рамберга [8, с. 282, рис. 11.93]
[8] (рис. 6, 7), с использованием применяемых в книге коэффициентов подобия. Эти вычисления позволяют провести сравнение с экспериментальными данными, хотя и с использованием не вполне корректных (см. [8]) коэффициентов подобия. Можно видеть, что для модели глубокого засасывания субстрата (рис. 6), в которой были использованы однородные материалы (все изготовлено из силиконовой замазки), расчетный результат (рис. 4) совпадет с экспериментальным вплоть до мелких деталей на самых глубоких стадиях эволюции. Для другой модели (рис. 7), в которой были использованы материалы с весьма различными свойствами — силиконовая замазка, грунтовка, модельная глина (различие вязкостей достигает трех порядков [8]), совпадение результатов расчета (рис. 5) с экспериментом не столь точное, но общие черты и основные особенности эволюции модели воспроизводятся. (Этот результат, также как и численные эксперименты других авторов [15, 16], говорит о незначительном влиянии умеренных различий вязкости на основные характеристики развивающегося течения и формирующихся структур.)
Для демонстрации возможностей программы на рис. 8 приведен пример расчета эволюции довольно сложной модели, состоящей из 17 слоев и одного замкнутого тела с различными плотностями. Левая часть модели (до 48 км) содержит 1.5-километровую гравита-
' / / / / / у /././,/ Л (
□
Рис. 7. Результат физического эксперимента из [8, с. 336, 337, рис. 15.4, 15.5]: начальная ( ) и конечная ( ) стадии развития неустойчивости
ционно неустойчивую пачку слоев разной мощности, имитирующих переслаивание легкой соли (22002300 кг/м3) с тяжелыми доломитами и ангидритами (2 700-2 800 кг/м). Пачка залегает на глубине 8.5-10 км и подстилается «карбонатами» с плотностью 2600 кг/м3. Справа «эвапоритовая пачка» выклинивается, упираясь в «барьерный риф», отгораживающий «лагуну» от «океана», где «эвапориты» отсутствуют. Сверху пачка перекрыта сначала «карбонатами» (2 600 кг/м3), а затем несогласно залегающей косослоистой «терригенной толщей», образованной слоями с плотностью 25002300 кг/м3. Самый верхний этаж разреза образует мощный слой (2 км) «слаболитифицированных отложений» со средней плотностью 2300 кг/м3.
Несмотря на наличие тяжелых прослоев «эвапори-товая пачка» вначале деформируется как единый слой (рис. 8, ), что обусловлено значительным контрастом ее средней плотности с перекрывающими «породами». (В ситуациях с малым отличием средней плотности такой пачки от плотности перекрывающего слоя наши расчеты показывают, что сначала развивается мелкая неустойчивость внутри пачки и только затем начинается ее деформация как единого слоя. Здесь мы эти модели не приводим.) Затем, на стадии развитых «подушек» начинает развиваться неустойчивость второго порядка внутри пачки (рис. 8, ). Ее развитие подавляется интенсивным общим течением при переходе растущих возмущений всплывающей пачки в диапировую стадию (рис. 8, , ), где они приобретают характерные грибообразные формы. В дальнейшем (рис. 8, ), после проникновения диапиров в верхнюю часть разреза, где общая плавучесть «эвапоритов» становится отрицательной, снова активизируется течение второго порядка внутри диапиров, приводящее к накоплению легкого вещества («солей»), в их апикальных зонах происходит формирование диапиров внутри диапиров. В результате, общая грибообразная форма больших диапиров искажается — они раздуваются из-за потери общей плавучести (этот эффект наблюдается у реальных диапиров [9]) и приобретают специфическую выпуклость в верхней части «шляпки» (из-за вспучивания внутренними диапирами). Из общих особенностей эволюции можно отметить следующие: развитие возмущений второй генерации из остатков вещества легкого слоя рядом с выросшими диапирами; глубокое засасывание субстрата в диапиры; растаскивание и втягивание внутрь растущего диапира «рифовой постройки» (на рис. 8, в зоне втягивания «рифовой постройки» под крайний диапир формируется идеальная ловушка для углеводородов).
Однако исследование общих закономерностей диа-пиризма не является основным приложением предлагаемой программы. Как нам представляется, основное достоинство разработанного метода — быстродействие в сочетании с корректностью — подразумевает оперативное использование его в целях практической геоло-
гии, а именно подбор такого варианта эволюции какого-либо конкретного разреза, который приводит к формированию структуры, оптимально соответствующей реально наблюдаемым геолого-геофизическим данным. Поясним эту идею следующим примером.
На рис. 9 показана эволюция гипотетического исходного разреза, представляющего собой упрощенный (для наглядности) вариант рассмотренного выше. Разрез включает эвапоритовый слой, залегающий на карбонатных породах, осложненный карбонатными рифовыми постройками; сверху он перекрыт косослоистой пачкой терригенных отложений (плотность последних для наглядности принята одинаковой—2 500 кг/м3). Эвапори-
товый слой имеет переменную мощность около 1 км и выклинивается справа, упираясь в барьерный риф. Кроме того, посередине он разделен другой рифовой постройкой на две части: левая часть — однородная соль плотностью 2200 кг/м3, а правая содержит ангидритовый прослой плотностью 2800 кг/м3. На финальной стадии расчета в левой части разреза образовались зрелые диапиры, один из которых почти достиг поверхности, тогда как в правой части развитие возмущений из-за наличия тяжелого прослоя еще находится в стадии «подушек» (внутри эвапоритового слоя видно развитие неустойчивости второго порядка). В бессолевой фланговой зоне и в зоне первоначального положения внутрен-
ней рифовой постройки сформировались прогибы. Сами рифовые постройки оказались разнесены и частично затянуты вместе с подстилающими карбонатами под растущие диапиры и подушки.
Представим, что имеем дело с реальным разрезом и сосредоточимся на зрелых диапирах в его левой части. (Представленный наклонный диапир с прилегающей к нему слева подушкой похож на реальный соляной купол
Хай-Айленд в округе Галвестон, Техас, США (рис. 10). Впрочем, подобные конфигурации довольно типичны как в расчетных моделях, так и в природе.) Предположим, что в подсолевой части разреза имеются нефтема-теринские породы. Тогда в рассуждении о возможной нефтегазоносности можно обратить внимание на подсо-левые карбонаты (включая рифовые), втянутые под диа-пиры. Действительно, их структурное положение в
Рис. 10. Соляной купол Хай-Айленд, Техас, США [8, с. 262, рис. 11.60]
верхних частях высокоамплитудных антиклиналей под идеальной покрышкой создает весьма благоприятные условия для формирования здесь залежей углеводородов. Кроме того, расчеты показывают (рис. 9, , ), что подкорневым зонам диапиров во все время роста соответствуют зоны пониженного давления. Причем горизонтальный градиент давления в данном примере достигает 0.8 • 104 Па/м, а это составляет примерно треть от литостатического градиента давления, так что разрежение под растущими диапирами весьма значительно и должно служить дополнительным сильным фактором, стимулирующим приток в эти зоны пластового флюида. В этих же зонах определяются максимумы скорости вертикальной деформации, т.е. растяжения вдоль вертикальной оси, что в сочетании с пониженным давлением должно способствовать развитию здесь трещиноватос-ти, улучшая коллекторские свойства пород. Таким образом, под растущими диапирами создаются структурно-динамические ловушки. Притом эти ловушки оказываются высоко поднятыми относительно нормального залегания подсолевых карбонатов, что делает их более доступными.
В надсолевой толще на модельном разрезе также есть высокоамплитудные структуры, связанные с ростом диапиров. Через «окно» между левой и правой частями эвапоритовой пачки углеводородсодержащие пластовые флюиды могут проникать в верхнюю часть разреза и при подходящих литологических условиях могут формировать залежи в висячем крыле опрокинутой складки, образованной наклонным диапиром (имеется в виду зона между 22 и 24 километрами модельного разреза). Заметим, что с этой зоной также сопряжены пониженное давление и вертикальное растяжение.
Если говорить о рассчитываемых компонентах тензора негидростатических напряжений, то на рис. 9, можно видеть, что они достигают довольно значительных величин — порядка 5 МПа для сдвигового напряжения и порядка 10-30 МПа для нормальных напряжений и давления. Помимо упомянутого учета возможного влияния давления и скорости деформаций на миграцию порового флюида, расчет напряжений (особенно сдвиговых) может быть полезен при проектировании буровых скважин. Также, возможно, имеет смысл учитывать их влияние на распространение упругих волн при интерпретации материалов сейсморазведки. (Как уже отмеча-
лось, в связи с квазистационарностью течения абсолютные значения рассчитываемых напряжений, в отличие от кинематических характеристик (скорость, деформации), не зависят от конкретного значения вязкости среды. Они могут рассматриваться и как постоянные упругие напряжения.)
В целом очевидно, что указанные детали структуры и истории деформирования разреза важны в плане поиска месторождений углеводородов. Между тем, если бы разрез, подобный модельному, был исследован с поверхности сейсморазведкой, то мы вряд ли увидели бы на временных разрезах такую картину с опрокинутым залеганием слоев и нависающими карнизами соли. Скорее всего, диапировые структуры предстали бы сглаженными до просто контрастных антиклиналей с крутым падением слоев в крыльях. Тем самым важные с поисковой точки зрения структурные особенности ядер антиклиналей были бы замаскированы. В определении этих особенностей, как мы полагаем, как раз и может помочь предлагаемая программа.
Предположим, что нам неизвестно, из какого начального состояния получен рассматриваемый модельный разрез. Попытаемся воспроизвести его левую часть (со зрелыми диапирами). При этом будем исходить из того, что (как в реальности) нам известны положение на разрезе максимальных и минимальных глубин залегания кровли соляного пласта, а также гравитационная аномалия и осредненный рельеф дневной поверхности. Варьируя из некоторых соображений возможные начальные условия, рассчитаем с помощью нашей программы ряд моделей эволюции с тем, чтобы на основании проверки указанных «известных» характеристик получаемых структур выбрать наиболее близко им соответствующую.
На самом деле, в данном случае мы «подглядывали в решение» и, руководствуясь характерной конфигурацией наклонного диапира с прилегающей подушкой, связанной обычно с выклиниванием соляного пласта или другим резким изменением его мощности, рассматривали различные варианты соответствующего сценария. Тем не менее, нам пришлось просчитать не один десяток вариантов эволюции (для примера несколько результатов приведены на рис. 11), чтобы нащупать наиболее подходящий (рис. 11, ). Невзирая на «подглядывание», можно видеть, что гравитационная аномалия и рельеф поверхности в указанном варианте лучше всего соответствуют таковым у подбираемого модельного разреза.
Конечно, здесь мы имели дело с искусственным и упрощенным для наглядности разрезом. Однако этим примером мы хотели показать, что подобное геодинамическое моделирование может быть полезным при решении реальных задач структурной и поисковой геологии. Понятно, что в силу известных особенностей объекта
. . /
12 1 ( 00 ) 63-74
73
0 I 8 1 2 1 3 2 0 2 28 км
Рис. 11. Варианты подбора эволюции левой части разреза (рис. 9, ) по возмущениям поверхности ( ) и гравитационного поля ( )
гравитационное моделирование структур соляной тектоники давно используется как в дополнение к сейсморазведке, так и самостоятельно. Но ясно и то, что в случае тел сложной формы, какими бывают соляные диапиры (и даже «подушки», учитывая конфигурацию подошвы), прямое гравитационное моделирование в подобной ситуации затруднительно ввиду сложности структуры, категорической недостаточности реперных точек и вытекающего отсюда неограниченного многообразия возможных вариантов подбора. В этой ситуации предложенное геодинамическое моделирование может рассматриваться как естественный регуляризатор, ограничивающий класс возможных решений. В качестве дополнительного средства контроля может использоваться надлежащим образом осредненный рельеф дневной поверхности, в некоторых случаях он может быть более контрастным, чем гравитационная аномалия.
Заметим, что важным условием продуктивности такого использования является возможность просчета большого числа вариантов эволюции, т.е. быстродействие программы геодинамического моделирования, которое обеспечивается в наших программах использованием аналитически полученного фундаментального решения.
Литература
1. . . Соляная тектоника платформенных областей. - М.— Л.: Гостоптехиздат, 1950. — 247 с.
2. s . .Mathematical formulation of salt-dome dynamics // Geophysics. — 1964. — V. 29. — No. 3. — P. 414—424.
3. . Salt-stock families in northwestern Germany // Am. Assoc. Petrol. Geol. Memoirs. — 1968. — V. 8. — P. 261—270.
4. . ., . . Fluid dynamics model for salt-dome evolution // Tectonophysics. — 1978. — V. 47. — No. 1—2. — P. 85—107.
5. .-., . . Finite element models of density instabilities by means of bicubic spline interpolation // Phys. Earth Planet. Inter. — 1980. — V. 21. — No. 2—3. — P. 176—180.
6. . . The propagation of two-dimensional and axisymmetric viscous gravity currents over a rigid horizontal surface // J. Fluid Mech. — 1982. — V. 121. — P. 43—58.
7. . ., . . . Evolution of salt structure, east Texas diapir province. Part 1. Sedimentary record of halokinesis. Part 2. Patterns and rates of halokinesis // AAPG Bulletin. — 1983. — V. 67. — No. 8. — P. 1219—1274.
8. . Сила тяжести и деформации в земной коре. — М.: Недра, 1985. — 399 с.
9. . . ., . . External shapes, strain rates and dynamics of salt structures // Geol. Soc. Am. Bull. — 1986. — V. 97. — No. 3. — P. 305—323.
10. . ., . . . Internal kinematics of salt diapirs // Am. Assoc. Petrol. Geol. Bull. — 1987. — V.71. — No. 9. — P. 1068—1093.
11. . On the relation between initial conditions and late stages of Rayleigh—Taylor instabilities // Tectonophysics. — 1987. — V. 133. — No. 1—2. — P. 65—80.
12. . . ., . . Anatomy of mushroom-shaped diapirs // J. Struct. Geol. - 1989. - V 11. - No. 1-2. - P. 211-230.
13. ., . Numerical simulation of Rayleigh-Taylor instability for single and multiple salt diapirs // Tectonophysics. - 1992. -V. 206. - No. 1-2. - P. 55-69.
14. . ., . ., . ., . . The effective viscosity of rocksalt: Implementation of steady-state creep laws in numerical models of salt diapirism // Tectonophysics. - 1993. -V. 225. - No. 4. - P. 457-476.
15. . . ., ., ., é ., -., . Numerical analysis of how sedimentation and
redistribution of surficial sediments affects salt diapirism // Tectono-physics. - 1993. - V. 226. - No. 1-4. - P. 199-216.
16. ., ., . . . Numerical models of complex diapis // Tectonophysics. - 1993. - V. 228. - No. 3-4. -P. 189-198.
17. . . ., . . Advances in salt tectonics // Continental Deformation. - Oxford: Pergamon Press, 1994. - P. 159-179.
18. é ., . Numerical modelling of salt diapirism: Influence of the tectonic regime // Tectonophysics. - 1994. - V. 240. -No. 1-4. - P. 59-79.
19. . ., - . ., .. Numerical approach to problems of gravitational instability of geostructures with advected material boundaries // Geophys. J. Int. - 1998. - V. 134. - No. 2. -P. 473-483.
20. - . ., . ., ., . Трехмерное моделирование соляного диапиризма: численный подход и алгоритм параллельных вычислений // Вычисл. тейсмология. - 2000. -Вып. 31. - С. 62-76.
21. - . ., .. Гравитационная неустойчивость идеально пластичного слоя, покоящегося на слое вязкой жидкости: следствие для диапиризма // Физика Земли. - 2001. - № 7. - С. 1017.
22. - . ., . ., . . Gravitational and buckling instabilities of a rheologically layered structure: Implications for
salt diapirism // Geophys. J. Int. - 2002. - V. 148. - No. 2. - P. 288-302.
23. . ., . ., - . .Salt structures and hydrocarbons in the Pricaspian Basin // Am. Assoc. Petrol. Geol. Bull. -2003. - V. 87. - No. 2. - P. 313-334.
24. - . ., . ., . ., .. Three-dimensional forward and backward modelling of diapirism: Numerical approach and its applicability to the evolution of salt structures in the Pricaspian basin // Tectonophysics. - 2004. - V. 387. - No. 1-4. -P. 81-103.
25. . ., . . Численное моделирование формирования соляных диапиров в земной коре // Математический журнал. - 2006. - Т. 6. - № 1(19). - С. 67-73.
26. - . ., .. Экструзия и гравитационное течение жидкости: применение к соляной тектонике // Физика Земли. - 2006. - № 12. - С. 34-42.
27. . . Напряжения в океанической литосфере, обусловленные плотностными неоднородностями // Глубинная морская геофизика. - Л.: Недра, 1991. - С. 29-36.
28. ., . Основы гидромеханики неньютоновских жидкостей. - М.: Мир, 1978. - 309 с.
29. . . О природе верхнемантийной аномалии плотности под Срединно-Атлантическим хребтом и ее роли в рифтогенезе и спрединге // Геология и геофизика. - 1996. - Т. 37. - № 9. - С. 87101.
30. ., . Scaling of Newtonian and non-Newtonian fluid dynamics without inertia for quantitative modeling of rock flow due to gravity (including the concept of rheological similarity) // Phys. Earth Planet. Inter. - 1986. - V. 43. - No. 4. - P. 316-330.
31. . . ., . . Dynamics of laboratory diapir and plume models // J. Geophys. Res. - 1975. - V. 80. - No. 5. - P. 705717.
32. .., .. Модели эволюции Земли и планет земной группы // Итоги науки и техники. Физика Земли. - М.: ВИНИТИ, 1980. - 232 с.
Поступила в редакцию 02.09.2008 г.
Лунев Борис Валентинович, ведущий инженер ИНГГ СО РАН, [email protected]
Лапковский Владимир Валентинович, к.г.-м.н., старший научный сотрудник ИНГГ СО РАН, [email protected]
