Прямая и обратная задачи геофлюидодинамики в приложении к прогнозированию зон АВПД в осадочных бассейнах 1. Теоретический аспект Текст научной статьи по специальности «Математика»

Прямая и обратная задачи геофлюидодинамики в приложении к прогнозированию зон АВПД в осадочных бассейнах

1. Теоретический аспект

1 2 А.Г. Мадатов , А.-В.И. Середа

1 НИИ Моргеофизика

2 Судоводительский факультет МГТУ, кафедра высшей математики и программного обеспечения ЭВМ

Аннотация. Разрабатываются постановка и методы решения прямой и обратной задачи флюидодинамики в масштабе геологического времени в приложении к прогнозированию аномально высоких поровых давлений (АВПД). Вводится базис чувствительных к развитию АВПД параметров геологической модели среды. Прямая задача флюидодинамики ставится относительно эволюции уплотнений, УВ-насыщений и поровых давлений при погружении элементарного объема пористой породы в процессе осадконакопления. В качестве входных данных для соответствующей обратной задачи (инверсии) используются измерения давлений, пористостей и УВ-насыщений в скважинах. С целью повышения вычислительной эффективности и обеспечения единственности и устойчивости решения проблемы инверсии размерность задачи сводится к эффективной 1.5D. Обсуждаются частные случаи применительно к однородной и многослойной среде. Приводятся примеры решения прямой и обратной задач флюидодинамики в приложении к прогнозу АВПД.

Abstract. In the paper the numerical solutions of the forward and inverse fluid dynamic problems for geological time scale have been considered in connection to excess pore pressure (EPP) prediction. The set of the most sensitive to the EPP evolution model parameters has been introduced. The forward problem of the fluid dynamics describes the compaction, HC-saturation and overpressure evolution vs. time and depth. The corresponding real well data set is used as an input for the inversion routine. The 1.5D reduction of the 3D inverse problem is developed. Such approach allows to achieve a more unique and stable inverse problem solution in reasonable computing time. The particular forward modelings for a homogeneous and layered medium are discussed. Applications of the approach to the EPP prediction are demonstrated on the practical examples.

1. Введение

Понимание таких явлений геологии осадочных пород, как уплотнение, генерация зон аномальных давлений флюидов и миграция углеводородов (УВ) из генерирующих толщ по разрезу постоянно опережает по фундаментальности описания подходы к их прогнозированию. За последнее десятилетие геологи-нефтяники получили в свои руки мощные инструменты бассейнового моделирования, основанные на численных решениях многофазных 3D задач флюидодинамики (Флюидодинамический фактор..., 1989; Ungerer, 1993; Verweij, 1993), а "воз" прогнозов и ныне стоит в области инженерии и эмпирики.

Так, например, петрофизики давно уже доказали (Авчан и др., 1979; Mann, Mackenzie, 1990), что возникновение и развитие зон аномально высоких поровых давлений (АВПД) в осадочных породах, находящихся на глубинах свыше 2.5-3 км, уже никак не может объясняться исключительно явлением недоуплотнения. Начиная с температур 85-90° (глубины 2.5-3 км), само понятие "уплотнение" для пород, подвергающихся активной цементации и диагенезу, теряет чисто "механистический" смысл, то есть оно является скорее следствием вторичных диагенетических процессов, чем собственно процессов механического уплотнения пород под действием стресса. Тем не менее, методы прогноза давлений до сих пор основаны в основном на анализе отклонений чувствительных к пористости параметров (скорости, плотности и т.п.) от линии "нормального уплотнения".

В чем причина подобного консерватизма? Причин, как нам кажется, несколько.

Первая - традиционный отрыв теоретических исследований от нужд практики.

Вторая - некоторый пессимизм по отношении к возможностям эффективного использования достаточно сложных в математическом отношении и тяжелых в вычислительном аспекте постановок прямых задач для реальных нужд прогнозирования АВПД.

Наконец, третья - четкое разделение прикладных задач по направлениям, гораздо более глубокое, чем требуется в многодисциплинарном подходе к бассейновому моделированию (иначе -отсутствие системного подхода, необходимого для эффективного решения подобных задач). Таким образом, практически анализом данных о первичной миграции УВ занимаются главным образом геохимики-нефтяники, а прогноз АВПД - удел буровиков. Теоретические исследования при этом посвящены в основном проблемам бассейнового моделирования, и практические аспекты прогнозирования АВПД не обеспечиваются. В результате никто никого не понимает.

Между тем общая точка приложения интересов теоретиков и практиков в нефтяной геологии есть, и лежит она в постановке и решении обратной задачи бассейнового моделирования.

На самом деле большинство полуэмпирических способов прогноза уплотнения, генерации зон аномальных давлений флюидов и нефтегазового насыщения коллекторов есть не что иное, как интуитивная инверсия в рамках частных моделей. Например, классические способы прогноза АВПД, такие, как метод эквивалентных глубин или метод Итона (Mouchet, Mitchell, 1989), обращают данные, чувствительные к пористости, на шкалу геофлюидальных давлений. Базируясь на простейшей механической модели уплотнения осадков (Magara, 1978), предложенной Терцаги еще в 1924 г., такая инверсия обеспечивает единственность и устойчивость решения. Но, к сожалению, в силу тех же причин она может приводить к совершенно ошибочным результатам, тем более грубым, чем сложнее история развития и разнообразие литологий изучаемого разреза. Уточнить же решение за счет привлечения других данных невозможно, поскольку простая модель уплотнения не включает более сложные, но не менее значимые феномены миграции порового флюида в процессе осадконакопления.

С другой стороны, в работах (Lurch, 1991; Zhao, Lurch, 1993; Yu et al, 1995) уже сформулированы основные принципы и подходы к инверсии данных о пористости, АВПД и зрелости УВ-генерирующих толщ в область параметров, управляющих динамикой эволюционных моделей уплотнения - УВ-генерации -сверхгидростатических давлений. Однако, будучи пионерскими, данные работы, на наш взгляд, не проработаны достаточно глубоко в теоретическом плане и с технологической точки зрения для реальной интеграции их в практику прогнозов АВПД они не готовы.

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

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

2. Калибровка бассейновой модели на основе рассогласования синтетических и реальных данных. Общая схема

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

В данном разделе мы сформулируем задачу калибровки нелинейных моделей уплотнения-миграции на операторном уровне. Это позволит выдвинуть основные требования к практической реализации подхода.

Как уже было сказано, методы решения ориентированы на минимизацию некоторого функционала J(x), оценивающего степень рассогласования синтетического и реального распределения некоторого целевого фактора в зависимости от значений параметров модели, или иначе - компонент

вектора модельных параметров %. Практически речь всегда идет о нахождении не одного наилучшего вектора модельных параметров, а об определении некоторого множества таких векторов, позволяющего задать область изменения параметров модели, в пределах которой синтетические данные удовлетворительно согласуются с реальными. Критерий рассогласования обычно связывается с качеством исходных данных, т.е., в частности, с состоятельными оценками разброса данных полевых экспериментов. Таким образом, наличие погрешностей в данных выражается в виде некоторой "размазанности" отображения вектора реальных данных на множестве сопоставления его с синтетическим вектором данных, или, иначе, - подмножества возможных векторов реальных данных, ассоциированного с наиболее вероятной реализацией этого вектора.

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

Таким образом, на области сопоставления, выступающей в качестве общего метрического базиса, при подборе вектора модельных параметров % учитываются не точки, но компактные подмножества, во взаимном расположении которых и ищется оптимум ('МаёаО, Ъетейа, 1997; ТгаиЬ, Woznjakovski, 1980).

Сформулируем теперь задачу калибровки более формально, учитывая специфику данных и их модельного описания.

Пусть модель геологического разреза, дающая описание процессов уплотнения-миграции с необходимой полнотой, насчитывает М репрезентативных элементов. Пусть также для каждого репрезентативного элемента определен набор параметров из N компонент для инвариантного по пространству и времени дифференциального оператора прямой задачи. Параметры предполагаются линейно-независимыми и могут быть отнормированы по средним значениям. Тогда в многомерном пространстве параметров моделей определен ортонормированный базис X, в котором вектор параметров подбираемой модели определяется какх = Х21,..., Хм1, Х1 ,..., Хм2, ..., 2а/}Т

По аналогии определим векторное пространство данных Б с элементом й = {й11, й21, ..., dK1,

2 2 3 Ь Т

d1 , ..., dK , d1 , ..., dK } . Размерность этого пространства определяется шагом независимых измерений в пространственной области калибрования К и числом независимо измеряемых физических характеристик Ь. Для простоты будем считать, что пространственная сетка реальных данных одна и та же для различных измерений и, кроме того, согласована с сеткой расчета синтетических данных.

Область сопоставления Ж определим как область векторного пространства, в котором синтетические и реальные данные имеют общие шкалы и выбранную метрику рассогласования. Иными словами, произвольный вектор / = {/^/2, ...,Л}т определяется здесь либо путем трансформации полевых наблюдений, заданных на некоторой пространственной сетке, в шкалу сопоставления с синтетическими данными, либо отображением на эту же сетку вектора модельных параметров, путем применения к нему дифференциального оператора прямой задачи.

Для построения и обеспечения возможности анализа функционала рассогласования введем меру

рассогласования в пространстве сопоставления. С этой целью определим скалярное произведение

Область модельных произвольной пары векторов /¡,/-> е Гв форме: Область <*н,ных параметров

{/'•■/:), =/'(',''/•■ (1)

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

Введем также норму векторов в соответствии с (1):

\\/ \\р = /,/ )/2. (2)

Пусть также отображение на область сопоставления из пространства данных выполняется оператором трансформации измерений Т(^), а из пространства параметров модели - оператором прямого моделирования М(х). Как подчеркивалось выше, данные операторы порождают в области сопоставления нечеткие подмножества, связанные с неизбежными погрешностями измерений и моделирования (рис. 1). Характеризуя эти ошибки отображения диаметрами (0) соответствующих подмножеств 3/, можно формально записать процедуры моделирования и трансформации измерений в операторной форме:

моделирование - /леЗ/л} = М(%); трансформация измерений - {/*е^/*} = Т(^).

Поскольку мера близости двух векторов в пространстве сопоставления определена, можно в общем виде сформулировать оптимизационную задачу следующим образом:

% = а^шт [./(*)], где 3(х) = \\Т(4) -Ы(Х)\\Р. (3)

Очевидно, что "идеальное" решение, доставляющее нулевое расхождение, недостижимо до тех пор, пока не определены вероятностные критерии при определении области пересечения двух нечетких отображений {/Ле^/Л} и {/*е^/*}, из пространства модельных параметров и пространства измерений, соответственно. Для простоты будем предполагать равномерный закон распределения внутри любого подмножества в пространстве сопоставления. Иными словами, предположим, что любая точка из пространства моделей либо данных отображается в пространство сопоставления в виде некоторого компактного "пятна" с конечными размерами.

При такой постановке критерий допустимых рассогласований отображений, трансформируемых из измерений и синтетических данных, определяется размерами максимальной из двух подобластей неопределенности. Условие, которому должно удовлетворять решение обратной задачи, в терминах проблемы минимизации функционала (3) можно записать в виде:

Х{ е X = а^шт ^(¿)р = \\1Щ -М(*)||р); /ле / п 6/*. (3*)

В этой ситуации возможны три варианта:

<0{$/л} ^ Неоптимальное решение - модель слишком грубая

> 0{$/л} ^ Неоптимальное решение - данные недостаточно точны (4)

0{д/*} = 0{^} ^ Оптимальное соотношение с точки зрения полноты модели и точности данных

Необходимо отметить, что отображения {/Ле^/Л} =М(%) выполняются нелинейным оператором моделирования. Кроме того, единственность решения обратной задачи может быть обеспечена лишь при выполнении дополнительных априорных ограничений в области вариации параметров модели. Область односвязных отображений будем называть областью решения и обозначать Ха с X. Тогда из всех возможных решений в качестве результата будет выбираться вектор %, принадлежащий подмножеству

■ г ¿1>

[Г(</) / М(%) / * I * : / , /

*5Ч Отображение модели ф Отображение данных Рис. 1. Общая схема процедуры подбора

Xa, однозначно отображаемому на область сопоставления (рис. 1). При этом ограничения на множестве Fa с F могут носить "жесткий" характер, когда вектора, не принадлежащие компактному подмножеству Fa, просто отбрасываются, либо "мягкий" - когда приближение изнутри к границам Fa отмечается введением штрафных функций, назначаемых на величину рассогласования. Окончательно стабилизированный по области поиска алгоритм решения задачи подбора может быть схематически записан в форме:

%f е XaсX = argmin(\\f*-fA\\P), f*e / /Ле / n Sf*. (3**)

Стратегия поиска решения (3**) в унимодальной области описывается ниже.

Базовые уравнения флюидодинамики во временной шкале эволюции осадочного бассейна основаны на принципе консервации масс и законе Дарси, удовлетворительно описывающем дивергенцию потока порового флюида в уплотняющейся пористой среде (Bear, Bachmat, 1991) в зависимости от градиента давлений. Здесь мы не будем обсуждать адекватность этой линейной модели реальным геологическим процессам, протекающим во многом катастрофично (Флюидодинамический фактор..., 1989). Постулируя закон Дарси в качестве подтвержденной практикой модели медленно протекающих процессов уплотнения-миграции флюида, мы лишь очерчиваем класс операторов прямой задачи, как дифференциальных операторов, описывающих непрерывно меняющиеся процессы, не имеющие разрывов в пределах области интегрирования. Точно так же кинетические модели УВ-генерации могут быть описаны в классе дифференциальных операторов (Tissot, Welte, 1978). Последние связывают темп и состав фазовых трансформаций органического вещества с абсолютным значением и градиентом температуры по всему ряду генерируемых УВ (от асфальтов до метана).

Решение системы дифференциальных уравнений относительно давлений и насыщений порового флюида углеводородами может быть получено в узлах сетки с точностью, которая должна быть согласована с реальными данными (Madatov, Sereda, 1997). В противном случае калибровка модели уплотнения-миграции, основанная на подборе (3*), окажется неоптимальной в смысле (4). Попутно с основными решениями, в каждом узле сетки рассчитываются пористости и проницаемости, которые формально также должны быть увязаны с соответствующими данными.

Таким образом, в качестве "данных" при калибровке должны выступать пространственные распределения поровых давлений, насыщений УВ, пористостей и проницаемостей, отнесенные к репрезентативным лито-стратиграфическим элементам разреза. Поскольку палеораспределения не доступны измерениям, речь всегда будет идти о настоящем времени - конце временной шкалы интегрирования. За исключением грубой лито-стратиграфической "разбивки" разреза по данным сейсморазведки, ни одно из требуемых измерений не доступно с поверхности. Иными словами, речь всегда идет о сети калибровочных данных, совпадающей с сетью скважин. Вопросам измерения давлений, насыщений, пористостей и проницаемостей посвящена специальная литература (см., например, Авчан и др., 1979; Magara, 1978). Там же можно познакомиться с обработкой комплексов каротажей, промысловых и технологических параметров бурения, используемых для трансформации в соответствующие шкалы, т.е. составить представление об операторах трансформации данных измерений T(d) и их погрешностях (Мухер, Шакиров, 1992).

Геофлюидальные давления доступны для измерений в чрезвычайно ограниченных зонах скважин, связанных с хорошо проницаемыми разностями пород, как правило - целевыми коллекторами. Там же, хотя и значительно реже, производятся измерения проницаемостей. Вся остальная информация о геофлюидальных давлениях (т.е. для плохо проницаемой части пород, которая, как правило, преобладает по разрезу) может поступать из косвенных источников, таких как D-экспонента, ст-каротаж и т.п. (Mouchet, Mitchell, 1989). При аккуратном соблюдении технологии бурении (т.е. при поддержании баланса удельного веса бурового раствора с градиентом формационных давлений), хорошим показателем давлений (по крайней мере, для максимальных значений) служит удельный вес бурового раствора.

Насыщения порового пространства УВ могут быть оценены на основании данных газового каротажа (Мухер, Шакиров, 1992), который, как правило, выполняется по всему забою в процессе бурения. Следует лишь заметить, что необходимые поправки за скорость бурения и промывки, кавернометрию и состав бурового раствора следует дополнить поправками на фазовое состояние смеси на забое. Хроматографический анализ выполняется лишь в зонах отбора проб на нефть и газ. Таким образом, калибровочная информация крайне неравноточна и неравномерна даже в пределах одного ствола. Дополнительно к газовому каротажу, важную информацию о достигнутой стадии зрелости

нефтегазоматеринской породы несут измерения отражательной способности витринита, которые выполняются в лабораторных условиях (Ungerer, 1993).

Открытая пористость также может быть измерена на образцах керна. Однако для массовых расчетов используются пересчеты из акустического, плотностного либо КС каротажей (Котяхов, 1977). Различные формулы, установленные эмпирически для песчанистых, карбонатных и глинистых разностей пород, имеют, как правило, достаточно локальную область применения и требуют поправок за минерализацию поровой жидкости, состав цемента, глинистость и т.д. Таким образом, уместнее говорить не об абсолютных значениях открытой пористости, а скорее о некой трансформанте, коррелируемой с данным параметром.

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

Общий вывод по наличию и качеству доступной информации для целей калибровки моделей уплотнения-миграции состоит в следующем:

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

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

- Точность трансформаций каротажной и технологической информации в требуемые шкалы изменяется от единиц до десятков процентов (относительной погрешности), и зачастую преобразованные данные пригодны лишь для качественных либо ориентировочных оценок предельных значений, как, например, удельный вес промывочной жидкости в применении к данным о геофлюидальных давлениях.

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

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

Очевидно, чем выше размерность модельного базиса, тем, в общем случае, более детально описывается среда. Вместе с тем вся идея прогноза основана на вычленении ключевых факторов, оказывающих доминантное влияние на результат моделирования в области сопоставления. В идеале набор наиболее чувствительных параметров модели среды, которые к тому же являются фундаментальными (т.е. слабо варьируемыми) параметрами разреза, является целевым базисом прогноза, т.е. таким подмножеством Хс с X, которое при минимальной размерности сохраняло бы свойства модельного отображения вектора синтетических данных на область сопоставления. К сожалению, заранее вычленить его, как правило, невозможно. Поэтому в общем случае необходимо начинать анализ чувствительности с наиболее детальной модели, постепенно "загрубляя" ее путем отбрасывания (закрепления) нечувствительных к решению параметров.

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

погружения бассейна. Очевидно, что по мере "старения" осадочной толщи увеличивается неединственность продвижения обратно в геологическом времени, которая становится катастрофической при увеличении числа перерывов в осадконакоплении. Таким образом, трехмерная обратная задача естественным образом не может быть решена единственно.

Однако чисто одномерная постановка не позволяет синтезировать все разнообразие наблюдаемых распределений данных по скважинам. Так, например, спад уровня АВПД, часто наблюдаемый при пересечении границы "УВ-генерирующая покрышка - резервуар", не может быть синтезирован без учета латерального дренажа флюида по проницаемому коллектору.

Нами был предложен и разработан подход к инверсии, основанный на эффективной 1.5D модели уплотнения-миграции геофлюидов, как основа калибрования транспортно-емкостных свойств осадочных пород (Madatov et al., 1996; 1997; 1998). Ниже представлены результаты разработки модели уплотнения-миграции, которая обеспечивает оптимальность подбора в смысле (4) и при этом является относительно простой в вычислительном отношении. Смысл символа "1.5D" состоит в том, что эффекты, вызванные невертикальными компонентами стресса при уплотнении-миграции, учитываются эффективно путем введения среднего за шаг интегрирования темпа оттока-притока геофлюида из одномерной системы в латеральном направлении. При этом практически значимые потоки флюида возможны лишь вдоль хорошо проводящих коллекторов с уровнем проницаемости единицы-сотни миллидарси (Verweij, 1993).

3. Прямая задача флюидодинамики для элементарного объема уплотняющейся пористой породы в масштабе геологического времени

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

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

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

3. Оператор прямой задачи должен непрерывно и однозначно отображать вектор модельных параметров на область сопоставления в шкалах геофлюидальных давлений, УВ-насыщений и пористостей.

4. Пространственная сетка, на которую отображаются синтезированные данные, должна совпадать с пространственной сеткой, на которой эти данные измеряются;

5. Каждый индивидуальный синтетический одномерный (вертикальный) профиль, сопоставляемый с данными по калибровочной скважине, должен эффективно включать 3D эффекты.

Заметим, что последовательное изложение фундаментальных основ теории уплотнения, УВ-генерации и миграции можно найти в ряде известных монографий (Verweij, 1993; Magara, 1978; Bear, Bachmat, 1991; Tissot, Welte, 1978). Поскольку наша задача состоит в описании модели, пригодной для калибровки и удовлетворяющей требованиям "1.5D", мы частично повторим основные отправные пункты и этапы вывода соответствующих уравнений. При этом акцент будет сделан на отличительных особенностях данной модели и ее ограничениях, а конечные формулы будут несколько отличаться от тривиальных.

U( T,hD„)

Время

■ t ■

х Траектория погружения ^ , репрезентативного объема V

Глл'бика |

Стресс,

Температура

Mi-генерация

/Пористость

Рис. 2. Схема уплотнения-УВ-генерации для погружающегося элемента породы U

Для математического описания процесса уплотнения введем две декартовых системы координат: внешнюю и внутреннюю (рис. 2). Внешняя содержит две координатных оси: абсолютного геологического времени и глубины. Ось глубин начинается от уровня моря и измеряет вертикальные движения дна гипотетического бассейна осадконакопления, где положительное направление совпадает с погружением. Временная ось начинается в некоторой стартовой точке истории осадконакопления и

совпадает с возрастом древнейшей формации - условного фундамента.

Во внешней, бассейновой системе координат введем в рассмотрение элементарный репрезентативный объем U - такой квазиоднородный элемент осадочной горной породы, в пределах которого можно пренебречь локальными вариациями свойств минерального скелета и наполняющего его флюида. Как известно, определение таких макросвойств осадочных пород, как пористость, проницаемость и УВ-генерационный потенциал, основано на осреднении по репрезентативному объему (Aziz, Settari, 1983).

Свяжем с элементом U локальную, внутреннюю систему координат, которая вместе с ним следует вдоль траектории погружения (Z, t), подвергаясь воздействию стресса и температуры. Возникающие при этом эффекты уплотнения и фазовых превращений контролируются реологическими (Magara, 1978) и кинетическими (Tissot, Welte, 1978) соотношениями, установленными эмпирически для фиксированного литологического типа. Последние однозначно связывают макрохарактеристики данного литотипа с глубиной (Z), или вертикальным стрессом (о), и/или температурой (T), которых достиг соответствующий репрезентативный объем (рис. 26). Заметим, что эффекты уплотнения, в рамках теории пластичных деформаций горных пород, рассматриваются, как необратимые. Т.е. траектории подъема в бассейновой системе координат соответствует плато на тренде уплотнения. Эффекты вторичной трещиноватости, равно как и цементации порового пространства, описываются, как правило, гораздо более сложными эмпирическими соотношениями, значимыми на локальном уровне (Котяхов, 1977). С нашей точки зрения, они должны включаться в эмпирические соотношения для реологии и кинетики данного литотипа и уточняться в каждом конкретном случае в ходе калибровки моделей уплотнения-миграции. Поскольку ни траектория погружения, ни тренд уплотнения не имеют петель, данная модель имеет тенденцию непрерывного уплотнения репрезентативного объема.

Конверсия части твердой фазы в жидкие и газообразные УВ ведет к повышению пористости и формально может рассматриваться, как "разуплотнение". Однако разница плотностей конвертируемой твердой фазы и нефти незначительна (типичная плотность керогена не превышает 1100 кг/м3, т.е. ненамного выше плотности поровой воды. Таким образом, выведение из матрицы керогена путем первичной миграции УВ, с точки зрения уплотнения материнской породы, эквивалентно выведению поровой воды (Ungerer, 1993)), а в случае генерации газа возникающий при этом дополнительный флюидодинамический потенциал обеспечивает немедленный вывод вновь образующегося флюида из объема, где он генерируется - первичную УВ-миграцию (Жузе, 1986). Поэтому явление УВ-генерации должно учитываться при анализе динамики баланса масс, но не возмущает тренда нормального уплотнения.

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

Прежде, чем записать исходные уравнения, необходимо удовлетворить условию 5, выдвинутому в начале раздела.

Здесь мы вынуждены ввести серьезное ограничение на размерность локальной модели разгрузки порового флюида, чтобы обеспечить возможность ее инверсии в область контрольных физических параметров элемента U. Необходимость ограничения пространственных координат локальной системы только вертикальной осью вытекает еще и из того факта, что законы уплотнения даны нам во внешней системе координат тоже в зависимости лишь от глубины, либо вертикального стресса (Magara, 1978). При этом игнорировать латеральные оттоки из элементарного объема мы не будем, но лишь уйдем от невосстановимых латеральных составляющих горного давления. Для простоты стартовый репрезентативный объем поместим в воображаемый цилиндр с основаниями единичной площади и высотой h. Тогда мера относительного уплотнения AU/U будет совпадать с мерой относительного

линейного сжатия АИ/И по вертикали. Данное ограничение совпадает, по существу, с введением трансверсальной изотропности для пористой среды. При этом, очевидно, линейная плотность вещества численно совпадает с обычной и отличается лишь размерностью. Здесь и далее мы не будем делать различий между линейной (на единицу длины) и обычной плотностью.

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

Заметим, что масса твердой фазы породы потенциально может лишь переходить в эквивалентную массу жидкой или газообразной фазы. В то же время исходный и сгенерированный в погружающемся элементе И поровый флюид может покидать объем в вертикальном и латеральном (радиальном) направлениях в виде потока Дарси. Никаких привносов масс внутрь репрезентативного объема извне не предполагается. Тогда уравнения консервации масс для случая п-компонентного порового флюида могут быть записаны следующим образом:

Здесь dMn - дифференциал массы n-й компоненты флюида, генерируемого из твердой фазы (0-й индекс в уравнениях). Безразмерный параметр v связан с пористостью ^соотношением: v= ^/(1-$. Символы h0, Sn, рп, qn обозначают соответственно: высоту репрезентативного цилиндра, насыщение, плотность (кг/м3) и вектор скорости [m/dt] потока n-й фазовой компоненты, причем h0 = h(1-^ü), а единица времени dt выбирается здесь и далее, исходя из шага квантования шкалы геологического времени.

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

Формально к явлениям конверсии твердой фазы в поровый флюид относятся любые диагенетические преобразования минерального скелета, происходящие в элементе U в процессе его погружения. Однако явлениями дегидратации глинистых минералов (например, конверсии иллита в монтморилонит), которые сопровождаются высвобождением внутрикристаллической воды, мы пренебрегаем. При этом делается допущение, что переход воды из внутрикристаллической в поровую компенсируется уплотнением трансформированного минерала, т.е. уравнение баланса масс не нарушается. Что касается изменения флюидопроводящих свойств матрицы, то они, очевидно, учитываются трендом уплотнения-проницаемости соответствующих глинистых литотипов.

При конверсии же части богатой органикой, нефтегазоматеринской породы в жидкий либо газообразный углеводород в области измеряемых параметров (давления и УВ-насыщения) происходят такие изменения, пренебречь которыми нельзя (Madatov et al, 1998). Более того, УВ-генерация является первопричиной возникновения скоплений нефти и газа и с этой точки зрения является предметом прогнозирования, т.е. соответствующая модель подлежит калибровке.

Число компонент УВ смеси и их плотность, в случае закрытой для привноса извне дополнительных масс порового флюида системы, вообще говоря, определяется фракционным составом первичной УВ-миграции, т.е. типом материнской породы, поскольку тяжелые и легкие фракции УВ взаимно растворяются друг в друге, и анализ плотностей многофазных систем может быть построен на законе "идеальной смеси" (Aziz, Settari, 1983). C точки зрения вторичной миграции удобно различать лишь жидкую и газообразную фазы. Тогда система (5) упрощается до трех компонент: 0 - твердая фаза: скелет осадочной породы, в том числе возможно содержащий конвертируемую в УВ часть - кероген, 1 -жидкая фаза: минерализованная поровая вода, 2 - УВ смесь, содержащая в общем случае газовую и нефтяную компоненту, изменчивая по плотности в связи с разделением легких и тяжелых компонент при падении давления и температуры.

p0dh0/dt = dM0/dt

d( vh0S1p1) /dt = dM1 /dt - (\Н0/ф) div(p1q1)

d(\hcSn-1Pn-1) /dt = dMn-1 /dt- (vh0/ ф) 6iv(pn-fln-i) d( vhßSn Pn) /dt = dMn /dt- (^/ф) div(/>n qn) dM0/dt= dM1/dt+ dM2/dt + ... + dMn/dt S1 + S2 + ... + Sn = 1

(5)

dho/dt = - h0öG,

'0OG0

d(vh0 (1-S)p1)/dt = (Фо/ф)-^^^^^

(5*)

d(vhoSp2)/dt = pohoSGo + (vh^^-divCc^),

где ¿G0 - темп УВ-генерации, являющийся безразмерной (кг сгенерированного УВ на кг материнской породы) функцией температуры и УВ-потенциала материнской породы (подробнее см. раздел 6).

Плотность минералогического скелета р0 принята неизменной, т.е. диагенетические изменения твердой фазы в процессе погружения не допускаются. Поэтому р0 в верхнем уравнении (5) было выведено из производной по времени и сокращено. Таким образом, потеря массы твердой фазы элемента U допускается только за счет УВ-генерации, а его плотность в каждой точке бассейновой системы координат целиком определяется пористостью матрицы, при известных плотностях порового флюида. Плотность жидкой и в особенности газообразной фаз порового флюида меняется под воздействием давлений и температуры, поэтому вместе с удельным насыщением ( vh0S) они введены в производную массы по времени.

С этой точки зрения система дифференциальных уравнений (5*), описывает динамику изменения плотности элемента в связи с механизмами потерь и генерации порового флюида. Для завершения построения математической модели необходимо дополнить конструкцию (5*) уравнениями, связывающими плотность порового флюида и пористость элемента U с его положением во внешней системе координат. Соответствующие тренды уплотнений флюида и матрицы являются эмпирической основой построения прямой задачи, а калибровка соответствующих констант - предметом решения обратной задачи. При этом часть констант может быть априорно известна настолько точно, что в пределах чувствительности обратной задачи не потребует калибровки.

Плотность порового флюида, как функция порового давления (Р) и температуры (7) контролируется в зависимости от фазы соответствующим уравнением состояния (Verweij, 1993; Bear, Bachmat, 1991), которое в линейной аппроксимации выглядит следующим образом:

pj dpj/dt = Yj dP/dt - [ij dT/dt. (6)

Изотермические константы сжимаемости (у) и изобарические константы расширения (J3) достаточно устойчивы и хорошо изучены в лабораторных и пластовых условиях.

Для обеспечения непрерывности решения системы (5*) необходимо, по возможности, исключить особые точки во внешней системе координат, где возможны разрывы 1-го порядка в поведении плотности порового флюида. Разрывные ситуации возникают при фазовых превращениях и отмечаются критическими точками на диаграммах фазового состояния. В частности, они имеют место при крекинге нефти в газ либо при выделении легких УВ фракций из раствора в виде свободного газа (Справочник по геологии..., 1984). Изменения же плотности отдельно жидкой и газообразной компонент УВ на плоскости бассейновой системы координат происходят гладко. Таким образом, для обеспечения непрерывности решения практически необходимо дополнить уравнения (5*-6) диаграммой предельного насыщения поровой жидкости газом (см., например, номограммы Дж. Амикса: Справочник по геологии..., 1984) и контролировать ее пересечение с кривой текущего насыщения газом в любой произвольной точке бассейновой системы координат. Более подробно этот вопрос освещается в следующем разделе.

Пористость элемента U может претерпевать изменения на траектории погружения: при его уплотнении под воздействием стресса а и/или температуры Т пористость сокращается, а при конверсии части минеральной матрицы в геофлюидальную фазу (SG0) - увеличивается. В самом общем виде:

dv/dt = 1/ (1 - ф)2 [(1 - <f)8G0 - dfca,T)/dt]. (7)

Тренды нормального уплотнения матриц, в зависимости от их литотипа, совместно с трендами проницаемостей формируют обобщенные модели уплотнения, которые анализируются в следующем разделе. Темп генерации УВ в аппроксимации модели Тиссо подробнее рассматривается в разделе 6.

4. Модель уплотнения порового флюида

В данной модели поровый флюид будет представляться минерализованной водой, насыщенной в общем случае некоторой массой углеводородов. Растворение последних в воде контролируется кривой предельного насыщения раствора в зависимости от температуры и давления (Справочник по геологии..., 1984). Физические свойства порового раствора, оказывающие влияние на его плотность, в практике определяются по закону "идеальной смеси", т.е. суммой свойств чистых компонентов, взвешенной соответствующими допредельными насыщениями.

Плотность любой "чистой" (без растворенных примесей) жидкой фазы определяется, в зависимости от температуры и давления, через коэффициенты изобарического термального расширения и изотермальной сжимаемости под давлением (см. табл. 1). Для плотности жидкой фазы в "рабочем" диапазоне значений физических характеристик осадочного бассейна (2 = 0-5км; T = 10-200°С; P = 0.150 МПа) достаточно точна линейная модель вида

P1(P,T) = P1(POTG) [1+ЖГ- TO + n(P - Po)],

(8)

Рис. 3. Зависимость плотности поровой воды от температуры при различных аппроксимациях

объемного расширения

1 - параболическая:

Рг = (25.9 - 8.38Т + 0.445Т 2) -10-4

2 - параболическая с поправкой

на минерализацию

3 - линейная: = 4.95

где Р0, Т0 - некоторые нормальные условия, в которых доступны массовые измерения (ниже в расчетах Р0 = 0.1 МПа, Т0 = 15°С).

Для определения плотности "чистой" газовой фазы используется уравнение состояния идеального газа (Ввпв^'И, 1998) с поправкой за "неидеальное" поведение под давлением ^-фактор, по Кацу - см. рис. 5. В формуле (8*) он обозначен через Я):

P2(P,T) = P2(PO,TO)[ (P/PO) (TO/AT ].

(8*)

•10"4T

В природных условиях, однако, "чистые" фазы порового флюида не наблюдаются. Поровая вода с увеличением глубины минерализуется. При этом разуплотнение под действием температуры частично компенсирует этот эффект. В итоге линейная модель (8) оказывается даже точнее, чем параболическая (рис. 3).

Компонентный состав УВ-смеси определяется типом материнской породы (см. раздел 6). При этом жидкая фаза (нефть) содержит тем больше растворенных легких фракций, чем выше пластовая температура и геофлюидальное давление (Aziz, Settari, 1983). Наличие и состав свободного газа определяется при помощи номограмм предельной растворимости (напр., номограммы Дж. Амикса: Справочник по геологии..., 1984). Пример по предельному содержанию метана в "черной" нефти представлен на рис. 4. Таким образом, плотность нефти и газа определяется в зависимости от компонентного состава по закону "идеальной смеси".

Характер поведения плотностей поровой воды, "черной" нефти и газа (метана) с глубиной легко оценить, используя линейное нарастание поровых давлений и температуры с гидростатическим градиентом 10006.2 Па/м и температурным градиентом 0.03°С/м.

,45 Предельная доля газа в нефти .... к

0.3 f' /

>.15 0 ....."

зо

90 ISO

10 10 30 40 50

Геофлюидальное давление [МРа]

Рис. 4. Растворимость метана в тяжелой нефти в зависимости от давления при различной температуре (числа у кривых). Использована парабол. аппроксимация данных по номограмме Дж. Амикса (Справочник по геологии..., 1984)

1.4

; 1.0 -факто]

' "Т* '«QIB - |К:фп ■ Ггкижхдъ -ч, t . 1 1 _____1

фгю¥ № S В-ГТТ1 "

10 20 30 40 50 Геофлюидальное давление [МРэ]

Рис. 5. Зависимость Z-фактора (коэф-та сжимаемости) для метана от геофлюидального давления. Использована параболическая аппроксимация по данным D. Bus, R. Witing, 1962 (см. Denesh, 1998)

ЮОО 3000 ТОО «И) 5000

Глубинам

Рис. 6. Поведение плотностей различных фазовых компонент порового флюида (минерализованная вода, тяжелая нефть, метан) с увеличением глубины при гидростатическом давлении и температурном градиенте 0.03°С/м

Как видно из рис. 6, по мере нарастания давления и температуры плотность нефти приближается к плотности "сухого" газа за счет растворения все большей объемной доли легких фракций. Разуплотнение воды по мере термального расширения и растворения в ней части углеводородов

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

Таблица 1. Физические свойства поровой жидкости (по данным: Справочник по геологии..., 1984)

Свойство Шкала Поровая вода "Черная" нефть

Коэфф. изобарического термального расширения 1/°С 5 • 10-4 7.2 • 10-4

Коэфф. изотермальной сжимаемости под давлением 1/Па 4.78 • 10-10 3.2 • 10-10

Плотность при н.у. кг/м3 1020 850

5. Обобщенная 1.5D модель уплотнения репрезентативного элемента породы

Модели уплотнения осадочных горных пород, как правило, представляются эмпирическими зависимостями пористости выбранного литологического объекта (литотипа) от глубины, стресса и/или температуры (Magara, 1978). Известно, что даже весьма сложные зависимости, комбинирующие эффекты механических (пластическая деформация, растрескивание, и т.д.) и диагенетических (цементация, дегидратация глинистых минералов, растворение под давлением и т.д.) процессов, приводящих к закономерным изменениям плотности пород, заведомо являются упрощением реальных процессов, протекающих на микроуровне (на уровне пор) (Verweij, 1993; Schneider et al., 1994).

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

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

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

Исследования на образцах (Schneider et al., 1994) показывают, что разгрузка вертикального стресса, возникающая вследствие подъема дна бассейна и частичной эрозии вышележащих толщ, не приводит к восстановлению утерянной пористости. Более того, при последующем погружении порода начнет уплотняться вновь лишь после прохождения точки, с которой был начат ее подъем. Сложнее обстоит дело с явлениями вторичных изменений порового пространства. Однако, считая их в основном температурно зависимыми (Котяхов, 1977), можно предположить подобную модель инвариантности к подъему и для диагенетических явлений.

Геологический темп осадконакопления редко превышает 100-150 м вертикальной мощности осадков с начальной пористостью 0.5-0.65 за миллион лет (Verweij, 1993). Это позволяет применить закон Дарси для связи скорости потока флюида из уплотняющегося элемента с действующим на него градиентом давления. Формально горизонтальные составляющие стресса также вызывают уплотнение репрезентативного элемента, т.е. отток флюида. Если распространить закон Дарси и на них, то для дивергенции потока порового флюида можно воспользоваться обобщенным на случай 3D представлением вида (Bear, Bachmat, 1991; Aziz, Settari, 1983):

div(^) = V- ((Kj-ВД fj),

(9)

где q, Kj и jUj - скорость потока, относительная проницаемость и вязкость j-й фазы порового флюида; KL

- проницаемость матрицы, определенная как функция пористости для L-ro литотипа в рабочем диапазоне бассейновой системы координат; fj = {rxVx(0j + pgz), ryVy(0j + pgz), rzVz(0j + pgz)}, где rx, ry, rz

- коэффициенты анизотропии проводимости матрицы, р - плотность поровой воды, g - ускорение свободного падения, Ф}- - флюидодинамический потенциал, имеющий размерность давлений и в зависимости от фазы порового флюида равный (England et al., 1987):

0j = p - gzpj

0j > Capj

(10)

0j = 0 ^ 0J < Capj

Флюидодинамический потенциал определяется, как работа по перемещению единичной массы флюида из бесконечности (практически из точки (Р(0),Г(0)) в данную точку (P(z),T(z)) внешней системы координат. Пренебрегая термальным расширением вытесняемого при уплотнении порового флюида1, можно отнести работу к единичному объему флюида, и тогда сопоставлять флюидодинамический потенциал с геофлюидальным давлением P. В частности, для воды, как следует из (10), динамический потенциал отличен от нуля при геофлюидальном давлении, отличном от гидростатического. Для более легких фракций, движущихся в заполненном водой поровом пространстве, если их рассматривать, как самостоятельную фазу УВ, дополнительно возникает потенциал всплытия. Кроме того, в случае многофазного потока, движение неводной фазы в поровом пространстве может осуществляться лишь с преодолением капиллярного барьера (или "входного давления") по данной фазе - Capj, измеряющегося в шкале давлений (Aziz, Settari, 1983).

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

div(^) = d/dz [ (KjKl/ц) rz (dO/ / dt)] + dQj /dt. (9*)

Здесь член dQj /dt обозначает темп совокупного оттока флюида из репрезентативного объема в латеральном направлении - "сток". Смысл элемента поясняет схема на рис. 7.

Неравномерное погружение дна бассейна и наличие на его краях стабильных зон, не испытывавших погружения, создают предпосылки для появления латеральной составляющей дренажа, направленной от центра к периферии. Эффект латеральной миграции особенно значим для легких УВ фракций порового флюида из-за капиллярных барьеров (покрышек) на пути их вертикальной миграции (Verweij, 1993; Жузе, 1986). Если совместить глобальную зону разгрузки с поверхностью, то, при наличии ненулевого флюидодинамического потенциала в любой точке в пределах области калибровки можно будет

оценить горизонтальную, составляющую градиента потенциала - Ф//Л. Очевидно, что при 8<< А (рис. 7), что, как правило, справедливо даже для глубоких грабенов, ко всей области калибровки можно применить одно и то же среднее удаление А, и таким образом учесть латеральный градиент давлений в составляющей dQj/dt - "стоке" уравнения (9*). Тогда различия в значениях этой составляющей от одной 1.5D модели уплотнения-миграции к другой (от одной калибровочной скважины к другой) будут целиком отражать локальные особенности латеральной проводимости. Эти особенности будут, очевидно, включать в себя также экранирующие свойства разломов, лежащих на пути латерального дренажа. Более подробно вопросы калибровки "стоковых" компонент описаны в (Madatov et al., 1998). В данном разделе мы сосредоточимся на оценке проводимости и параметров вертикального уплотнения матрицы репрезентативного объема. Эти параметры, вместе с эмпирическим законом уплотнения и проницаемости для данного литотипа, определяют важнейшие составляющие базиса пространства параметров изучаемой обобщенной модели уплотнения.

Параметры нормального уплотнения входят в виде констант в соответствующие эмпирические тренды, описывающие уменьшение пористости с ростом глубины, стресса и/или температуры (рис. 9). Проницаемость матрицы KL, в свою очередь, определяется для данного литотипа в зависимости от пористости ф и константы проницаемости AL, характеризующей извилистость поровых каналов, либо

1 Х1аюгощ и Уая^г (1992) показали, что термальное расширение поровой воды существенно влияет на флюидодинамический потенциал только при изоляции репрезентативного объема практически непроницаемыми породами. В этом случае, впрочем, аппроксимация (10) в виде закона Дарси также оказывается недопустимой.

Удаленная точка латеральной разгрузки

Рис. 7. Схема латеральной разгрузки давлений

V» \ / /

0 75 J ' Я , /

0 5 X к, +к2 У

0.25 ^ * X * * К2 * _

О 0.3 04 <¡.6 0 8 1.0

Насыщение порового флюида газом Рис. 8. Зависимость относительной проницаемости матрицы от содержания газа в поровом пространстве (по Викову-Ботсету (Справочник по геологии..., 1984))

удельную площадь внутренней поверхности порового пространства (Котяхов, 1977) (рис. 10). Вязкость ц данной геофлюидальной фазы ] также непрерывно определена в зависимости от температуры эмпирическими соотношениями (Ввпв^'И, 1998) (рис. 11). Наконец, проводимость матрицы по данной геофлюидальной фазе] определяется для данного литотипа Ь в виде ЩКЬ/^ (Увг^ву, 1993; Справочник по геологии..., 1984) и характеризует темп миграции порового флюида из репрезентативного элемента при единичном градиенте геофлюидального давления (рис. 12).

о

0.6 0.5 0.4 0.3 0.2 0.1 0

Пористость

Рис. 9. Эмпирические тренды "нормального уплотнения" в шкале глубины (температурный градиент 0.03 °С/т)

1. Терригенный линотип

ф= фо ехр(-г/Т)

ф0 = 0.6; г/ = 0.0037 - Глинистые матрицы; ф0= 0.5; т]= 0.0025 - Песчанистые матрицы

2. Карбонатный линотип

ф= ф0 — т]2, при Z < Zcr, или же ф= - сг) ехр[-30-^ - Zcr)] Zcr = Тсг/8Г; - глубина коллапса пористости,

где Тсг = 95°С, 5Т = 0.03° С/м фо = 0.45; т] = 0.00006 - Меловые матрицы

10 1 103 104 10 5 10* 10"7 10г* Проницаемость [О] (Юагсу = Ю"12 т2)

Рис. 10. Эмпирические тренды проницаемости в шкале глубины (температурный градиент 0.03°С/т)

1. Терригенный линотип (Глинистые матрицы)

Км = Хм(0.2-ф3^2)) / [1 - ФьШ2),

где Хм = 10-16 м2

2. Терригенный линотип (Песчанистые матрицы)

Кз = ^ ф8з(2), где ^ =10-12 м2

3. Карбонатный линотип (Меловые матрицы)

Кс = 10[10-3^с(2)-76],

где = 3.1510-13 м2

Рис. 11. Эмпирические зависимости вязкости

порового флюида от температуры в шкале глубины (температурный градиент 0.03°С/т)

1. Вода

М2) = 103{5.28 + 3.8А(2 - 0.28[А(2)]3}-1, где А(2 = [Т(2) - 150]/100

2. Нефть

= 1.4186 • 10-10 • ехр[6597/(237 + Т(7))]

3. Газ

Д2 = [10-5 + 1.510-62 - 2.2(Т(2) - 15)] • 10-7,

где Т(2 = Т0 + 2 Т0 = 10°С, ёТ = 0.03°С/м

Удельная проводимость [см/1000лет] Удельная проводимость [см/1000лет]

Рис. 12. Удельная (отнесенная на единичный градиент давлений - 1Pa/m) геофлюидальная проводимость в зависимости от глубины: (а) - для песчанистой матрицы по различному типу порового флюида; (б) - для различного литотипа матрицы по поровой воде. Тренды рассчитаны на основании данных, представленных на рис. 9-11

Ухудшение проницаемости матрицы при наличии газообразной фазы в поровой жидкости учитывается понижающим коэффициентом Kj относительной проницаемости. Он однозначно определяется по каждой из фаз при известном насыщении порового пространства газом (Справочник по геологии..., 1984).

6. Модель генерации УВ флюида

Темп преобразования <G0(t) твердой фазы УВ в жидкую либо газообразную (см. формулу (5) или (5*) в разделе 3) может быть определен, как функция времени, в зависимости от температуры и свойств материнской породы (УВ потенциала), на основе кинетических моделей первичной миграции нефти и газа, впервые предложенных Тиссо (Tissot, Welte, 1978). Классическая кинетическая модель этой реакции описывается следующим дифференциальным уравнением 1-го порядка:

dM(t) /dt = - k(t) M(t), (11)

которое рекуррентно связывает мгновенный темп конверсии фаз с остающейся массой твердого вещества материнской породы, еще способного к конверсии, через кинетическую функцию Аррениуса k(t):

k(t) = A exp[ - E/RT(t)]. (12)

Здесь A - константа Аррениуса, определяющая асимптотический темп реакции (при бесконечно высокой температуре) [1/ед.времени]; Е - энергия активации реакции [ккал/моль]; R - газовая постоянная = 1.986 • 103 ккал/(моль°К); T(t) - температура по Кельвину [°К], данная как функция времени, т.е. в ретроспективе траектории погружения репрезентативного элемента. При линейной связи глубины с температурой через температурный градиент, которая использована в предыдущем разделе, зависимость T(t) однозначно определяется темпом осадконакопления.

УВ-потенциал, как известно (Ungerer, 1993; Tissot, Welte, 1978), определяется долей массы твердой фазы - керогена, потенциально способной к полному преобразованию в УВ флюид в ходе кинетической реакции в единичной массе материнской породы при нормальных условиях (P0, T0), т.е. в условиях практически нулевого темпа реакции. Используя принятые выше обозначения и совмещая начальный момент времени с положением репрезентативного элемента в начале бассейновой системы координат (0;0) О (0;P0,T0), определим это фундаментальное свойство литотипа следующим образом: H0 = M(0)/p0h0. Теперь выразим искомую функцию SG0(t) через параметры кинетической модели и УВ-потенциал данного литотипа. Для этого, в соответствии с определением генерационного члена, перепишем первое уравнение системы (5*) с учетом (11). Получим:

SGo(i) = - k(t)M(t) /poho.

Поскольку решением (11) является M(t) = M(0) exp[-f k(t)dt], то окончательно будем иметь:

SG(t) = - H0 k(t) exp [-Jk(t)dt ]. (13)

Исследования последних лет показали, что для более адекватного описания модель УВ-генерации требует объединения нескольких кинетических реакций (11) в одновременно протекающий параллельный процесс (Akihisa, 1978; Tissot, 1987). Кроме того, необходимо дополнить модель фазовых трансформаций вторичным крекингом нефти в газ, происходящим при относительно более высоких температурах, чем первичный крекинг керогена в УВ (Ungerer, 1993).

Таким образом, вместо выражения (12), для определения темпа k(t) практически используется сумма вида:

k(t) = TrnA, exp [ -E,/RT(t) ], (14)

где i - счетчик фракций УВ, объединенных по общей энергии активации; ш; - весовая доля фракции в керогене на момент начала реакции.

Набор параметров (щ, Ai, Е) характеризует тип керогена и доступен изучению на образцах в лабораторных условиях. В зависимости от распределения весовых долей Шг и соответствующих энергий активации Е;, известны 3 типа керогена, описанные в литературе (Ungerer, 1993; Tissot, Welte, 1978; Жузе, 1986). Это создает предпосылку для локализации зоны поиска, или хороших стартовых условий, при калибровке кинетической модели. Точно так же еще до начала калибровки могут быть высказаны априорные суждения об УВ-потенциале пород. Как правило, нефтегазоматеринскими являются богатые органикой глинистые разности, известные в бассейне заранее. В частности, стандартные характеристики нефтегазоматеринских пород: TOC (Total Organic Content) и HI (Hydrocarbon Index), как правило, колеблются в пределах от 0.1 до 10 % и 100-600 мг конвертируемой фазы/г ТОС, соответственно (Verweij, 1993). Таким образом, параметр Н0 может с высокой степенью вероятности попадать в интервал 10-4 - 610-2 кг УВ/кг тв. фазы породы.

Пример модели УВ-генерации, распараллеленной по шести фракциям, представлен на рис. 13. Здесь темп преобразования <G0(t) твердой фазы в УВ представлен на температурной шкале в нормированном виде (т.е. для единичного УВ-потенциала). Поскольку фазовая диаграмма данной смеси может быть рассчитана заранее, жидкие компоненты объединены в нефтяные фракции, а газообразные - в газовые. Доля вторично генерируемого из нефти газа включена в газовую составляющую генерации и показана отдельно. Стартовый состав керогена представлен на диаграмме энергий активации (рис. 136).

Температура [ С] Энергия активации реакции [ Ккэл/мояь]

Рис. 13. Модель УВ-генерации а - темп преобразования твердой фазы в зависимости от температуры; б - диаграмма энергий активации по исходным фракциям керогена

7. Уравнения для однородной среды. Анализ условий разгрузки геофлюидальных давлений

Пусть во внешней системе координат задана траектория погружения дна бассейна и определен температурный градиент (рис. 2). Условие однородности среды реализуется допущениями о неизменности плотности твердой фазы р0, а также констант уплотнения (фо,г[), проницаемости (Л) и УВ-генерации (Н0, {А, Е, щ}) при любом положении стартовой точки (0,0) истории погружения репрезентативного элемента и при заданных законах уплотнения ф(2) и проницаемости К(ф) матрицы. Таким образом, в любой точке (2, /) на бассейновой шкале может быть определено достигнутое горное

давление L(t,Z) и температура T(t,Z). По определению (Mouchet, Mitchell, 1989), горное (литостатическое) давление на глубине Z равно весу столба породы, приходящегося на единичную поперечную поверхность. Оно может быть также представлено в виде суммы порового давления P(t,Z) и вертикального стресса o(t,Z). Т.е.:

L(t,Z) = {[1 -ф(Щ ро + ф (t)Pl}Zg = P(t,Z) + o(t,Z),

где g - ускорение силы тяжести.

Наша задача в этом разделе будет состоять в получении на основании этих данных уравнений, явно связывающих функции геофлюидальных давлений, УВ-насыщений и пористости матрицы с перечисленными выше фундаментальными константами породы данного литотипа. Исходными будем считать: систему уравнений (5*) баланса масс для репрезентативного элемента U и уравнения (6-8), связывающие гладкое изменение плотностей породы и флюида с траекторией погружения U во внешней системе координат.

Сделанные выше допущения о "репрезентативности" элемента U, однородности среды и практической несжимаемости жидкой фазы порового флюида отразим более формально:

grad (pj) = 0, V j = 0, 1, 2 - päiöägäiöäöeämöü yeaiäiöä U, т.е. постоянство свойства внутри

малого объема;

dp0/dt = 0 - стабильность минерального скелета матрицы в процессе

погружения; (15)

1/р1 dp\/dt = - ß1 dT/dt - несжимаемость и термальное разуплотнение поровой воды;

1/р2 dp2/dt = у2 dP/dt - ß2 dT/dt - сжимаемость и термальное разуплотнение УВ.

Определим темп (производную по времени) насыщения порового пространства УВ-флюидом в процессе уплотнения матрицы, миграции порового флюида и УВ-генерации. Воспользуемся вторым уравнением системы (5*):

dS/dt = 1/vh0 [h0(1-S) dv/dt + v (1-S) dh/dt- vh0(1-S) ß1dT/dt + vh0div q1 /ф\

Учтем первое из уравнений системы (5*) и после упрощений получим:

dS/dt = (1-S) [1/vdv/dt - dG0/dt - ß1 dT/dt + div(q1) /ф ]. (16)

Перепишем последнее уравнение системы (5*), выполняя дифференцирование:

h0Sp2 dv/dt + vSp2 dh0/dt + vh0p2 dS/dt + vh0S dp2/dt = pho dG0/dt - vh0 div q1 /ф.

Выполняя подстановки dh0/dt, dp2/dt, dS/dt из (5*), (15), (16), соответственно, получим после группировок и упрощений:

(1/v) dv/dt - Sy2 dP/dt = [Sß2 + (1-S) ß1] dT/dt + [P0/( vp2) + 1] dG0/dt - (div q1 + div q2) /ф. (17)

Темп относительной потери пористости, с учетом нормального тренда уплотнения и генерации УВ из части твердой фазы, в общем случае определен выражением (7). Для дальнейшего анализа воспользуемся достаточно общим экспоненциальным законом убывания пористости от вертикального стресса:

ф= фо exp[-ro(t,Z)]. (18)

Такая формула оказывается даже более точной для глинистого литотипа (Magara, 1978), а соответствующая константа уплотнения т может быть легко связана с введенной выше константой уплотнения ^ в формуле Anthy (рис. 9):

г= 11/[(1 ~Ф) (Ро-Р1) g]. (19)

Тогда перепишем (7) с учетом (18) и определения литостатического давления (см. начало раздела):

1/vdv/dt = ^dGf/dt- т/(1-^) d(L-P)/dt. (7*)

Подставляя (7*) в (16) и (17) получим после упрощения:

dS/dt = (1-S) /ф [(1-^) dG0/dt - rvd(L-P)/dt - фß1 dT/dt + div q1] (20)

ф [T/(1-<fi)+Sy2] dP/dt = {(1-$ (p/p2 -1) dG/dt + фт/(1-ф) dL/dt + [S^2 + (1-S)^] dT/dt} - (20*)

- [div (q1+q2)].

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

Очевидно, что результат при выбранном литотипе (т.е. зафиксированных законах уплотнения и УВ-генерации) и данной траектории погружения будет однозначно определяться константами уплотнения, проводимости и УВ-генерации, а также параметром "сток" в 1.5D формулировке закона Дарси. Константы вязкости, термального расширения и изотермического сжатия для компонент порового флюида принимаются при этом неизменными (табл. 1).

Оба уравнения не имеют физической размерности и, по сути, описывают динамику баланса нагрузки и разгрузки элемента по содержанию УВ в поровом пространстве и геофлюидальным давлениям. Проанализируем, например, уравнение (20*).

Коэффициент C = ф[т/(1-ф) + S/2] имеет смысл полной (матрицы + флюида) сжимаемости элемента U и измеряется в [1/Pa].

Обозначим первый член правой части (20*) через SU<. Он объединяет в фигурных скобках все источники, вызывающие мгновенный дефицит порового пространства (за единицу времени). К ним относятся, в частности: трансформация части минерального скелета матрицы в менее плотную УВ фазу -1; уплотнение матрицы при погружении - 2 и расширении порового флюида при нагреве - 3. Переходя к конечным разностям (dt ^ St) и принимая шаг по времени St за единицу времени, запишем выражение для источников дефицита порового пространства:

SU< = (1-0 (Р0/Р2-1)Ä?0 + фг/(1-ф)31 + [Sß2 + (1-S) ßi]8T. (21)

Три слагаемых в правой части (21) определяют три важнейших механизма формирования зон аномально высоких поровых и пластовых давлений (АВПД): 1 - газогенерация; 2 - некомпенсированное уплотнение матрицы при вертикальном стрессе за счет быстрого осадконакопления; 3 - температурное расширение заблокированного в поровом пространстве флюида (Mann, Mackenzie, 1990).

Второй член правой части (20 ) обозначим через SU>. Он определяет суммарный (по воде и УВ) вклад относительного оттока порового флюида за единицу времени, происходящего в вертикальном и латеральном направлениях пропорционально градиенту геофлюидальных давлений, достигнутому за тот же промежуток времени.

С учетом перехода к конечным разностям перепишем (20*) в упрощенном виде:

C SP = SU< - SU>.

Очевидно, что SU< >/C имеет размерность геофлюидальных давлений и смысл флюидодинамического потенциала ЗФ накапливаемого (SU</C) и разгружаемого (SU>/C) системой за единицу времени St. Условие полной разгрузки элемента U по геофлюидальным давлениям состоит в том, что дефицит порового пространства, накапливаемый в нем по всем источникам за время интегрирования, не должен превышать суммарного (по фазам и направлениям) оттока из него флюида. Иначе, потенциал нагрузки ¿Ф< не должен превышать потенциала разгрузки ЗФ>. Для численного решения этой проблемы необходимо задать одно начальное и два граничных условия по распределению геофлюидальных давлений на сетке, заданной в бассейновой системе координат. Схема решения этой задачи и результаты для реальных моделей обсуждаются в следующем разделе. Здесь же мы ограничимся приближенной оценкой составляющих дефицита пористости SU< и соотношения SU</SU>, как функций текущей глубины Z элемента U для идеализированной однородной среды.

Прежде всего оценим удельный вклад каждого из трех элементов, формирующих SU< в (21) при следующих допущениях:

- Дефицит порового пространства, вызываемый на любой глубине Z средним за единицу времени темпом осадконакопления Sh(0) [м/MY], целиком определяется трендом нормального уплотнения для данного литотипа (MY - миллион лет. Этот шаг принят за шаг квантования геологического времени).

- Термальное расширение порового флюида рассчитывается для случая однофазного порового флюида (0= 510-4 °С-1 - поровая вода) и градиентной модели нагревания с глубиной T(Z+SZ) = T(Z)+ST-SZ, (ST = 0.03°С/м).

- Дефицит порового пространства связанный с первичной миграцией УВ в уплотненной материнской породе рассчитывается на основании модели УВ-генерации данной в предыдущем разделе (рис. 13).

Как видно из рис. 14, термальный эффект сравнительно невелик и может создавать значительный потенциал нагрузки лишь при низкой (близкой к нулю) проницаемости вышележащей толщи, что вполне

согласуется с (Xiaorong, Vasseur, 1992). Вклад фактора уплотнения в среднем на два порядка значимее. Однако он может уступать значимости фактора УВ-генерации в момент пикового темпа.

На рис. 15 приведены графики соотношений потенциалов нагрузки к потенциалам разгрузки (ЗФ</ЗФ> = ¿и<1§и>) для различных темпов осадконакопления однородной среды при постоянном температурном градиенте 0.03°С/м. Очевидно, что интервал глубин до пересечения кривой 3и<13и> с уровнем 1.0 будет соответствовать области полной разгрузки системы при данных параметрах погружения и константах литотипа. Видно, что при высоком темпе осадконакопления (¿й = 100 м/МУ) система не успевает разгружать элемент и уже с глубины 1000 м. "Включение" с глубины примерно 3100 м газогенерации остается относительно незаметным фактором. Наоборот, при низком темпе (¿й = 1м/МУ) генерация газа в поровом пространстве становится ключевым фактором и может немедленно вызвать АВПД.

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

Естественно, что и прогнозы АВПД, основанные на подобной упрощенной модели среды, могут давать колоссальные ошибки.

8. Численное решение прямой 1.5Б задачи для многослойной среды

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

Для численного решения (20-20*) при заданных начальных и граничных условиях сеточными методами во внешней (бассейновой) системе координат (, /) зададим прямоугольную сетку, на которой может быть реализована любая из известных сеточных схем. В данной работе выбрана четырехточечная неявная сеточная схема, обладающая необходимыми свойствами устойчивости. Численное решение задачи начинается с нулевого временного слоя, соответствующего моменту начала формирования исследуемого глубинного разреза. Для построения решения на следующем временном слое формируется и решается методом прогонки трехдиагональная система линейных уравнений.

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

Пусть флюидодинамический потенциал, возникающий за п-й шаг на временной сетке - "такт накопления", создается мгновенно и определяется, в соответствии с (21), в виде:

¿Ф< = {(1-0 (Ра/р2-1)Я?о + фт/(1-ф)8Ь + [ЗД + (1-5) ¡31]5Г}/ {¿|У(1-0+5й]}. (22)

Тогда уравнение (20*) описывает процесс разгрузки этого потенциала путем миграции двухфазного флюида сквозь п накопленных к этому моменту формаций в вертикальном и латеральном направлениях.

Глубина 1т]

Рис. 14. Дефицит пористости, вызванный различными источниками, для однородной среды при среднем темпе осадконакопления 75 м/МУ

Зона АВПД ¡ !

_I_i_J.-

О 01

0 001

1000 200G 3000 4000 5ÜOO

Рис. 15. Соотношение потенциалов нагрузки и разгрузки элемента по давлениям для однородной модели среды в зависимости от текущей глубины. Числа на кривых - темп осадконакопления, m/MY

Рис. 16. Схема разбивки траектории погружения на такты с переменным темпом осадков

и составом накапливаемых

Тот же прием используем в дифференциальном уравнении (20) для насыщений. Вычленяя дивергенцию потока поровой воды в качестве разгрузки данного флюидодинамического потенциала за такт, определим мгновенное изменение доли УВ в поровом пространстве, произошедшее перед этим тактом вследствие УВ-генерации, уплотнения матрицы и расширения поровой воды в виде:

SS = [(\-ф^0- tvög^ST] (1-ЬУф. (23)

В результате система (20-20*) описывает последовательность релаксационных процессов разгрузки геофлюидальных потенциалов (22) и УВ-насыщений (23), мгновенно генерирующихся на траектории погружения элемента U в каждый такт вычислений.

Схема, представленная на рис. 16, поясняет сказанное. Здесь весь интервал временной шкалы разбит на N тактов нагрузки-разгрузки, в течение каждого из которых происходит мгновенное приращение общей массы литологической колонки единичного поперечного сечения на величину hnpn:

hnPn = honpo n + ho vn[(1-Sn)p1 n + S nP2 n], (24)

где в правой части нижний индекс указывает на фазу вещества, а верхний соответствует номеру такта. При этом каждый из (n-1) уже отложившихся элементов разреза погружается с уплотнением, а также нагревается и мгновенно повышает свой флюидодинамический потенциал (22) и УВ насыщение (23) на величину ЗФ< и SS соответственно.

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

Если по завершении N-ro такта осадконакопления условие полной разгрузки системы выполнено для всех элементов разреза, то, очевидно, давление в вертикальной выводящей трубе будет эквивалентно гидростатическому на любом уровне, а соотношение мощностей твердой и флюидальной фаз по каждому элементу будет соответствовать значениям пористости на трендах нормального уплотнения соответствующих литотипов. Напротив, возникновение флюидодинамических барьеров на пути вертикальной разгрузки приведет к возникновению сверхгидростатического давления, которое ослабит стресс. Формально, если барьер окажется непроницаемым, то разгрузка системы возможна только за счет латерального дренажа из-под покрышки. Если и эта возможность отсутствует, то геофлюидальные давления на уровне барьера и ниже будут возрастать с литостатическим градиентом, а при дополнительной генерации флюида и за счет термального расширения поровой воды существенно быстрее, до тех пор, пока не достигнут предела трещиноватости запирающей дренаж породы (Mouchet, Mitchell, 1989). Последнее условие

выступает в качестве физически обоснованного ограничения роста поровых давлений при реализации разностной схемы.

В ходе численного решения системы (20-20*) необходимо увязать последовательность процессов нагрузки-разгрузки сквозными (т.е. независимыми от номера такта) начальными и граничными условиями по давлениям и насыщениям. Учитывая, что выражения для потенциалов (22-23) легко рассчитываются на произвольной глубине при заданных траектории погружения, температурном градиенте и литотипах осадков, запишем эти условия следующим образом:

Ры'к = 0, ^к = 0 - разгрузка по давлению и насыщению на поверхности для любого такта;

ро.к = р1,к _ 0ТСурСТВИе потока сквозь нижнюю границу, т.е. непроницаемость фундамента;

Р = Р -1,ы-1 + 8Ф - преемственность предшествующего решения по давлению и насыщению для = + любого элемента глубинной сетки;

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

Достигаемые уровни геофлюидального давления и УВ-насыщения в каждой точке глубинной сетки анализируются по окончании каждого такта разгрузки на предмет превышения критических значений. Это необходимо делать, чтобы удовлетворить условиям применимости закона Дарси для двухфазного потока в пористой среде. Ограничения накладываются на значения констант проводимости в течение стольких временных шагов, сколько необходимо для вывода решений из закритических областей. В частности, следующие ограничения вводятся по проводимости УВ-фазы К2к"п:

5к"п < 5ксг ^ К2к"п = 0 - дйдаааёшй ийййшеа аадше оадй а жеаёи Шал офеаа; Рк"n+gzn(pгp2n) < САР2к ^ К2к"п=0 - газодинамический потенциал не превышает капиллярного барьера,

где 5сг - критическое насыщение поровой жидкости газом, САР2 - капиллярный барьер по газу.

172

1132_-92_^_^ 0 | 1-0 ■ ■ ■ 1[5 ■ ■ ■20 ■ ■ ■ ■ ■ ■ 3

.5

3

4

5

УВ Насыщение

Рис. 17. Последовательные этапы численного решения прямой 1.5Б задачи для истории уплотнения - миграции - УВ-генерации, подразделенной на 13 тактов (номера в рамочках). Реальный пример по Центральному грабену Северного моря. а - разгрузка на промежуточном такте 7 (слева схема погружения, справа промежуточное решение); б - разгрузка на промежуточном такте 12 (слева схема погружения, справа промежуточное решение); в - разгрузка на последнем такте (слева схема погружения, справа решения, отвечающие распределению

геофлюидального давления, пористости и насыщения по глубине в настоящее время). Обозначения и индексы совпадают с введенными в тексте для формул (22-24). Дополнительно: ф'мэт - кривая нормального уплотнения комбинированная по различным литотипам;

8П

. - текущий предел растворимости газа в жидком поровом флюиде.

По проводимости матрицы К1 вводится следующая подстановка:

в случае превышения предела трещиноватости породы РЯ

РКп I „ п/ п\ \ Т7пк т^ к"п_ т^ к ^^ т^ кт

+ gz (Р1 - р2 ) > ЕЯ ^ Кр = КРЯ >> Кр , где т - номер такта, на котором произошел "пробой" породы.

1

В общем случае временной шаг для различных тактов нагрузки-разгрузки по оси времени не одинаков, поскольку определяется длительностью отложения данного литотипа. Численное решение (20-20*) в пределах каждого такта осуществляется на прямоугольной сетке (Т, /) с равномерным шагом по времени при соблюдении приведенных выше граничных и начальных (по каждому такту) условий.

Пример расчета прямой задачи для многослойной среды представлен на рис. 17.

9. Обратная задача (задача о калибровке)

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

9.1. Общие положения

Пусть прямая задача может быть решена для некоторого глубинного разреза. Это означает, что имеется информация, достаточная, чтобы априорно задать вектор значений параметров ^ = (х1°, Х2°, ..., Хг°)Т, характеризующих свойства пород и флюида, и позволяющих рассчитать (на основании принятых модельных представлений о геологических, механических и физико-химических процессах - см. предыдущие разделы) модельное (синтетическое) распределение целевого фактора (поровых давлений, пористостей и т.п.) в узлах той же пространственной сетки, в которой представлено "реальное" распределение указанного фактора, получаемое в результате прямых или опосредованных измерений.

В общем случае синтетическое и "реальное" распределения будут различаться. Чем большим будет такое различие, тем больше вероятность, что априорные данные или модельные представления, положенные в основу постановки прямой задачи ошибочны и требуют уточнения (см. Введение).

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

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

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

При рассмотрении формальной постановки и метода решения обратной задачи будем использовать следующие термины и обозначения:

Р = (Р1, Р 2, ..., Р ь)Т - заданное ("реальное") распределение целевого фактора по глубине геологического разреза, Ь - количество точек по глубине разреза, в которых производились измерения значений "полевых" данных. Известны также глубины залегания этих точек. Р еР°, где область Р°сЕсЕЬ (ЕЬ - Ь-мерное евклидово пространство) определяется точностью значений

"реального" распределения целевого фактора, Е- область сопоставления. Х= (Хь Х2, ..., Хг)Т - вектор значений безразмерных варьируемых параметров прямой задачи, г - общее количество таких параметров.

где Х°<^Ег (Ег - г-мерное евклидово пространство) - область априорно допустимых значений варьируемых параметров.

В соответствии с предположением о возможности решения прямой задачи, для любого набора допустимых значений варьируемых параметров, задаваемых вектором % может быть получено (в результате решения прямой задачи) синтетическое распределение Р(^)сЕ - целевого фактора по глубине геологического разреза как функция ^

Р(Х) = С(Х), (25)

где О(х) - сеточный оператор прямой задачи, осуществляющий отображение точек из пространства Х° в точки из области сопоставления Е. Вопросы, связанные с погрешностью, обусловленной конечно-

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

Оператор G(x) является в общем случае нелинейным, и его свойства существенным образом зависят от принятых при решении прямой задачи модельных представлений, используемого сеточного метода и значений компонент вектора х- отображаемой точки в пространстве варьируемых параметров.

Введем в рассмотрение векторную функцию R(x):

R(X) = n{P(x) - P*)/ ||P*||, (26)

где: Q - диагональная LxL матрица заданных весовых коэффициентов для компонент вектора (P(Z) ~ P ); ||P || _ евклидова норма вектора P . Элементы матрицы Q определяются в зависимости от степени "доверия" значениям соответствующих компонент вектора (P(x) — P ) и могут, вообще говоря, меняться в зависимости от х.

Тогда обратную задачу можно сформулировать следующим образом:

Пусть по глубине геологического разреза задано "реальное" распределение P eP0. Требуется определить в рамках заданных модельных представлений такие значения варьируемых параметров прямой задачи — компонент вектора % которые доставляют минимум 9(х) — функции рассогласования синтетического и "реального" распределений целевого фактора по глубине данного геологического разреза:

min = |ВД||2 (27)

^

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

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

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

Тем не менее подбор в результате решения обратной задачи "эффективных" значений варьируемых параметров, позволяющих получить в результате моделирования (решения прямой задачи) синтетическое распределение P(x)eP°, несомненно имеет важное практическое значение для прогнозирования этих распределений в пределах региона, в котором осуществлялась калибровка модели.

9.2. Методы решения

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

]£+1 = / + к = 0,1,2,..... (28)

То есть, начиная с некоторого начального (стартового) вектора £°eX°, они будут строить в области X0 по рекуррентной формуле (28) бесконечную в общем случае последовательность векторов ^, х1, ..., ^ ..., сходящуюся к некоторому х - точке локального минимума функции рассогласования ¥(£). Вектор Sk в формуле (28) задает направление движения из точки в точку 1, а число qk определяет длину шага в выбранном направлении.

Свойства решений, получаемых в результате работы методов типа (28), будут существенно зависеть, среди прочего, от адекватности модельных представлений, положенных в основу прямой задачи, от точности полевых исследований (область P0) и априорных представлений о возможных значениях варьируемых параметров задачи (область X0). Если полученное решение по тем или иным причинам неудовлетворительно, или требуется найти и другие локальные минимумы функции рассогласования, то другое решение может быть, в принципе, найдено если начать с другого начального

приближения j^&X». В целом следует отметить, что, в связи с недостаточным знанием в настоящее время свойств функции рассогласования реализация любой стратегии движения в пространстве

параметров не должна исключать участие эксперта.

Главное отличие одного метода вида (28) от другого будет заключаться в способе выбора вектора Sk и числа qk на каждом шаге итерационного процесса (28). Важной характеристикой используемого для решения задачи (27) метода является также критерий окончания его работы. Будем считать приемлемым с практической точки зрения критерием попадание модельного распределения P(%) в область P0.

Одним из простейших подходов к численному решению обратной задачи (27) является градиентный метод наискорейшего спуска. При реализации этого метода в качестве вектора Sk на каждом шаге итерационного процесса (28) выбирается вектор антиградиента функции рассогласования вычисленный в точке а число qk определяется посредством решения одномерной задачи минимизации:

¥/(^k+qkSk) = min ¥(yk+qSk) (29)

q>0

В процессе реализации метода наискорейшего спуска значения компонент вектора антиградиента - частных производных функции рассогласования ¥(£) - определяются в конечно-разностном виде. Таким образом, на каждом шаге процесса потребуется дополнительное решение r прямых задач для определения значений компонент вектора Sk. Решение одномерной задачи оптимизации (29) для определения оптимальной длины шага также потребует решения некоторого конечного числа прямых задач при любом выборе оптимизационного алгоритма.

Недостатки метода наискорейшего спуска общеизвестны. Тем не менее, он может оказаться полезным для построения достаточно близкой к решению точки ^eX» - стартовой для более эффективного алгоритма, используемого для уточнения решения. С этой целью будем использовать метод, который может быть отнесен (Dennis, Scnabel, 1983) к методам типа метода Гаусса - Ньютона для решения нелинейной задачи о наименьших квадратах (27).

Переформулируем обратную задачу (27) как задачу поиска решения системы нелинейных уравнений:

Пусть для некоторого глубинного разреза задано "реальное" распределение PeP0. Требуется определить в рамках заданных модельных представлений такие значения варьируемых параметров прямой задачи - компонент вектора x'^X0, что выполнены условия:

Rk(*) = 0, k=1, 2, L, (30)

где Rk(x), k = 1, 2, ..., L — компоненты вектора R(%), определяемого согласно (26).

Очевидно, если существует точное решение нелинейной системы уравнений (30), то оно отвечает нулевому значения функции рассогласования ¥(£).

Начинаем с некоторого начального вектора априорных данных, = (rf>1, ..., Z»r)T, X?GX°. Строим последовательность точек ^, ^,..., ^ ..., принадлежащих X0, по следующему правилу:

J(/) Sk = -R(/), (31)

= /+ qkSk, k = 0,1,2,..., (32)

где: Sk = (Sk1, S2, ..., Sr)T - по-прежнему вектор поправок к вектору а число q задает длину шага в выбранном направлении; J(jf) - матрица вычисленных в точке $ частных производных компонент векторной функции R(jf) по компонентам вектора варьируемых параметров задачи (30). Как и в методе наискорейшего спуска, частные производные находятся в конечно-разностном виде, и для их определения на каждом шаге процесса требуется дополнительное решение r прямых задач.

Реализация метода (31-32) для решения задачи (30) предполагает практическое разрешение ряда достаточно известных и в тоже время непростых проблем как вычислительного, так и теоретического характера, связанных с решением системы уравнений (31) и определением значения qk на каждом шаге процесса. Кроме того, необходимо иметь в виду, что направление движения из точки ^, задаваемое вектором Sk, в общем случае может не быть направлением спуска для функции рассогласования ¥(#), т.е. не будет приводить в этом случае к уменьшению значения функции рассогласования в точке по сравнению с ее значением в точке ^ Детальное рассмотрение этих проблем не является целью настоящей работы, поэтому ограничимся ниже лишь кратким комментарием по поводу основных из них.

В общем случае матрица J($) прямоугольная Lxr, с числом строк, большим, чем число столбцов (L>r). То есть для определения направления Sk= (Sk1, Sk2, ..., Scr)T на каждом шаге метода (31-32) решается, в общем случае, переопределенная система линейных уравнений (31). Для таких систем, как правило, ищется обобщенное решение, которое отвечает решению в смысле метода наименьших квадратов. При этом, если указанное решение не единственно, обеспечивается выбор решения с минимальной нормой - min ||Sk||. Такой выбор целесообразен, так как при линеаризации нелинейной задачи, что имеет место на каждой итерации метода (31-32), естественно стремление к малым приращениям значений компонент вектора решения. В реализованной авторами вычислительной процедуре обобщенное решение системы (31) ищется с использованием сингулярного разложения матрицы J(jf).

Если J(jf) - матрица системы (31) - плохо обусловлена, то при решении этой системы возникают принципиальные проблемы. В этом случае предлагается проводить тщательный сингулярный анализ с целью выявления неустойчивой части решения. Техника проведения сингулярного анализа достаточно хорошо изучена (Lawson, Hanson, 1974). Его практическая реализация в настоящей работе в целом основана на общеизвестных рекомендациях.

Использование сингулярного разложения матрицы J(jf) для нахождения обобщенного решения переопределенной системы (31) оказывается удобной при решении обратной задачи (30) еще и потому, что сингулярное разложение и его анализ позволяет получать на каждом шаге метода (31-32) полезную информацию о системе (31) в целом. Эта информация может быть продуктивно использована для изучения и улучшения свойств метода решения обратной задачи, в том числе, для анализа проблем, связанных с оценками погрешностей и организации различных физически обоснованных процедур регуляризации. Так, например, наличие близких к нулю сингулярных чисел свидетельствует о почти линейной зависимости между столбцами матрицы J(^k), а следовательно, и о зависимости между соответствующими компонентами решения системы в (31). Коэффициенты этой почти линейной зависимости могут быть достаточно просто получены и проанализированы при известном сингулярном разложении. Эта информация может быть чрезвычайно полезна эксперту при принятии решения о корректировке априорных данных. В частности, в этом случае известны базисные векторы нуль-многообразия матрицы J(£:), что в принципе позволяет проводить анализ этого множества решений системы (31), а также организовывать при необходимости корректировку решений системы (31) без ухудшения качества этих решений с точки зрения суммы квадратов невязок.

Кроме того, для улучшения свойств плохо обусловленной системы линейных уравнений (31) еще до ее решения целесообразно осуществление ряда предварительных операций. К их числу относятся:

- масштабирование варьируемых параметров задачи;

- взвешивание уравнений системы (31);

- присвоение отдельным варьируемым параметрам фиксированных значений.

Главной целью таких операций является стремление к повышению устойчивости решений и снижению вычислительной погрешности до уровня соизмеримого (или меньшего), чем неопределенности в решении, вызванные погрешностями в исходных данных. Обоснование и способы организации указанных операций также хорошо известны (Lawson, Hanson, 1974). Схема реализации зависит от особенностей конкретной задачи и может варьироваться на каждом шаге метода (31-32), исходя из анализа свойств системы (31) на данном шаге.

Наконец, нахождение на каждом шаге метода (31-32) значения qk, определяющего длину шага в выбранном направлении Sk, следует осуществлять, исходя из требований, предъявляемых к длине шага в методах ньютоновского типа, учитывая при этом наличие ограничений на изменение значений варьируемых параметров задачи. Последнее обстоятельство является принципиальным, поскольку рассматриваемые методы решения обратной задачи в классическом варианте не ориентированы на наличие ограничений на изменение значений варьируемых параметров и относятся к так называемым безусловным методам. Область X0, задающая интервалы возможного изменения этих параметров, как правило, существенно ограничена в пространстве, что отражает естественное стремление к обоснованному заданию интервалов как можно меньшей длины. Это соответствует большей априорной информированности о свойствах изучаемого геологического разреза и объективно должно способствовать получению более качественных решений обратной задачи, но осложняет процесс нахождения этих решений.

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

случаях, когда точка ^gX° находится достаточно близко от границы области X0, указанное условие может оказаться нарушенным, и потребуется выполнение корректной процедуры корректировки

направления Sk и значения qk. В результате фактически речь идет о необходимости построения и реализации некоторых модификаций стратегий безусловного поиска с учетом ограниченной допустимой области изменения значений варьируемых параметров задачи.

9.3. Особенности программной реализации

На основе приведенных выше положений была разработана и программно реализована процедура формирования и решения обратной задачи. Процедура оформлена в виде программного модуля и позволяет осуществлять решение обратной задачи как методом наискорейшего спуска, так и методом типа Гаусса - Ньютона с использованием сингулярного разложения при решении системы (31). Можно отметить следующие характерные особенности реализованных процедур решения обратной задачи:

- Оба указанных режима работы программного модуля предполагают решение одномерной оптимизационной задачи (29) методом последовательного дробления шага с кубической аппроксимацией минимизируемой функции в окрестности текущего значения длины шага. При этом блокируется возможность выбора чрезмерно малой или чрезмерно большой длины шага (Dennis, Scnabel, 1983). При реализации метода типа Гаусса - Ньютона обеспечивается, по возможности, полный ньютоновский шаг.

- В процессе работы программного модуля исключается возможность построения очередного приближения yf+'íX0. В случае необходимости осуществляется корректировка выбранного на

k-м шаге процесса направления Sk. Полученное в результате корректировки направление гарантированно остается направлением спуска для функции рассогласования.

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

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

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

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

20-25 итераций метода (31-32) для достижения достаточно близких к нулю значений функции рассогласования. При этом качество калибровки параметров прямой задачи оказывалось достаточным для надежного прогноза поровых давлений в сходных условиях.

10. Заключение

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

калибровка модели ^ интерполяция параметров в точку прогноза ^ моделирование,

представляется более надежным, чем простая интерполяция данных либо прогноз на основе эмпирической калибровки.

Технология, развитая на основе описанного подхода реализована в виде пакета сервисных программ PANDA® в среде WINDOWS. Пакет PANDA® успешно опробован и получил широкое

распространение при прогнозе плотности бурового раствора при бурении глубоких скважин в шельфовых зонах Северного и Печорского морей.

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

Литература

Aziz K. and Settari A. Petroleum reservoir simulation. Applied Science Publisher, New York, p.475, 1983. Bear J. and Bachmat Y. Introduction to modeling of transport phenomena in porous media. Kluwer Academic

Publishers, London, p.553, 1991. Denesh A. PVT and phase behaviour of petroleum reservoir fluids. Elsever, p.400, 1998. Dennis J.E. and Scnabel R.B. Numerical metods for unconstrained optimization and nonlinear equations.

Prentice-Hall, Inc., Englewood Cliffs, New Jersey, 440 p., 1983. England W.A., Mackenzie A.S., Mann D.M., and Quigley T.M. The movement and entrapment of petroleum

fluids in subsurface. Journal of the Geological Society, v.144, p.327-347, 1987. Lawson Charles L. and Hanson Richard J. Solving least squares problems. Prentice-Hall, Inc., Englewood

Cliffs, New Jersey, 230 p., 1974. Lurch I. Inversion of dynamical indicators in quantitative basin analysis models. 1. Theoretical considerations.

Mathematical Geology, v.23, 1 6, p.817-832, 1991. Magara K. Compaction and fluid migration. Elsevier Scientific Publishing Company, p.319, 1978. Madatov A.G. and Sereda A.-V.I. Computing aspects of the solution of some inverse problems in geophysical investigations. Presented at the International Conference "Northern Universities". (October 16-18, 1997). MSTU, Murmansk, p.57-58, 1997. Madatov A.G., Sereda V.-A.I., Doyle E.F. and Helle H.B. The "1.5-D" Inversion approach to the pore pressure evaluation. Concept and application. Presented on Workshop "Compaction and Overpressure Current Research", December 9-10, 1996, Institute Francais du Petrole, Paris, p.6, 1996. Madatov A.G., Sereda V.-A.I. and Doyle E.F. Pore pressure prediction by using inversion before and during drilling. Presented on the workshop "New methods and technologies in petroleum geology, drilling and reservoir engineering", June 19-20, 1997, Krakow, Poland, p.12, 1997. Madatov A.G., Sereda V.-A.I. and Doyle E.F. Integration and inversion of well data into the basin evolution model: Way to new generation of pressure prediction technologies. Presented on American Association of Drilling Engineers (AADE) Forum "Pressure regimes in sedimentary basins and their prediction", September 2-4, 1998, Houston, Texas, USA, p.7, 1998. Mann D.M. and Mackenzie A.S. Prediction of pore fluid pressure in sedimentary basins. Marine and

Petroleum Geology, v.7, № 2, p.55-65, 1990. Mouchet J.P. and Mitchell A. Abnormal pressures while drilling. Manuels techniques, Elf Aquitaine, Boussens, p.286, 1989.

Schneider F., Potdevin J.L., Wolf S. and Faille L. Modele de compaction elastoplastique et viscoplastique

pour simulateur de bassins sedimentaires. IPFRevue, v.49, 1 2, p.141-148, 1994. Tissot B.P. and Welte D.H. Petroleum formation and occurrence. Springer-Verlag, New York, p.538, 1978. Traub J.F. and Woznjakovski H. Theory of optimal algorithms. Academic Press, New York, p.220, 1980. Ungerer P. Modelling of petroleum generation and expulsion - an update to recent reviews, in basin modelling:

Advances and applications. NPFSpecialPubl., v.3, p.219-232, 1993. Verweij J.M. Hydrocarbon migration systems analysis. Development in Petroleum Science, v.35, p.276, 1993. Waples D.W. and Kamata H. Modelling porosity reduction as a series of chemical and physical processes.

NPF Special Publication 3, Elsevier, Amsterdam, p.303-320, 1993. Xiaorong L. and Vasseur G. Contributions of compaction and aquathermal pressuring to geopressure and the

influence of environmental conditions. AAPG bull, v.76, 1 10, p.1550-1559, 1992. Yu Z., Lurch I. and Bour Q. Inversion of dynamical indicators in quantitative basin analysis models. 3. Multiwell information and two-dimensional case history. Mathematical Geology, v.27, 1 1, p.41-68, 1995.

Zhao K. and Lurch I. Inversion of dynamical indicators in quantitative basin analysis models. 2. Synthetic tests and a case history using dynamical indicator tomography. Mathematical Geology, v.25, 1 2, p.107-123, 1993.

Авчан Г.М., Матвеенко A.A., Стефанкевич З.Б. Петрофизика осадочных горных пород в глубинных условиях. М., Недра, с.222, 1979.

Жузе Т.П. Миграция углеводородов в осадочных породах. М., Недра, с.188, 1986. Котяхов Ф.И. Физика нефтяных и газовых коллекторов. М., Недра, c.287, 1977.

Мухер A.A., Шакиров А.Ф. Геофизические и прямые методы исследования скважин. М., Недра, с.335, 1992.

Справочник по геологии нефти и газа. Под ред. Еременко H.A. М., Недра, 480 с., 1984. Флюидодинамический фактор в тектонике и нефтегазоносности осадочных бассейнов. Сб. ст., АН СССР, Научный совет по проблемам геологии и геохимии нефти и газа. М., Наука, 320 е., 1989.