Математическое моделирование эволюции тепломассопереноса в процессе формирования осадочного чехла

Размещено на http://www.allbest.ru/

Автореферат

диссертации на соискание ученой степени доктора физико-математических наук

Математическое моделирование эволюции тепломассопереноса в процессе формирования осадочного чехла

специальность 25.00.10 «Геофизика, геофизические методы поисков полезных ископаемых»

Суетнова Е.И.

Москва - 2008

Работа выполнена в Учреждении Российской Академии Наук

Институт физики Земли им. О.Ю. Шмидта РАН

Научный консультант

Доктор физ.-мат. наук, академик РАН, профессор Глико А.О.

Официальные оппоненты:

Доктор физ.-мат. наук, профессор, член-корр. РАН,

Институт океанологии РАН. Лобковский Л.И.

Доктор физ.-мат. наук, Институт физики Земли РАН. Ребецкий Ю.Л.

Доктор геол.-мин. наук, профессор. Геологический институт РАН.

Хуторской М.Д.

Ведущая организация: Институт геофизики УРО РАН

ОБЩАЯ ХАРАКТЕРИСТИКА РАБОТЫ

Актуальность темы.

Построение математических моделей геофизических процессов, протекающих в недрах Земли, позволяет изучить различные сценарии эволюции геофизических полей и выделить физически обоснованные варианты интерпретации данных наблюдений и сформировать критерии выявления различных глубинных процессов по данным наблюдений. Значительные регионы земной коры фактически представляют собой не сплошную твердую породу, а являются мультикомпонентной средой, содержащей переменное количество флюида (обычно жидкость) рассеянного в трещинах, порах или межзеренном пространстве твердой породы. Наиболее важными примерами таких регионов представляются рифтовые зоны и зоны интенсивного осадконакопления, такие, как осадочные бассейны, континентальные окраины, и глубоководные желоба, в которых поровый флюид участвует в процессах формирования структуры верхней коры. Насыщающий флюид определяет специфику процессов тепломассопереноса в этих зонах и изучение проблем связи движения флюида и деформации породы и роли флюидов в тепломассопереносе, происходящем в процессе эволюции коры этих зон, является очень важным так как взаимодействие механических, тепловых и гидродинамических процессов в течение геологической истории формирования и эволюции верхней коры этих регионов создает условия для образования газогидратных скоплений, рудных и нефтегазовых месторождений. При наращивании мощности осадочного чехла Земли происходит ряд геофизических и геохимических процессов, формирующих структуру осадочных пород, и базовым процессом является уплотнение и связанная с ним фильтрация к поверхности насыщающего флюида. Процессы уплотнения приводят к падению свободной пористости с ростом глубины от поверхности пород вплоть до средней и нижней коры, и прогрессивному обезвоживанию земной коры в течение истории роста мощности осадочного слоя. Однако, при принципиальном сходстве трендов, наблюденные зависимости распределения пористости и порового давления от глубины оказывались различными для различных регионов [Hall, 1993]. В первых работах на эту тему, относящихся к тридцатым годам прошлого века были предложены полуэмпирические законы экспотенциального падения пористости с ростом глубины, причем коэффициенты подбирались путем сопоставления с экспериментальными данными по измерениям пористости в скважинах в конкретном регионе и считалось, что распределение пористости по глубине зависит только от типа осадков [Athy, 1930]. Позднейшие модификации закона Ати уже связывали пористость с эффективным давлением, как это было до этого сформулировано в механике грунтов для гидростатического порового давления [Terzaghi, 1923]. В дальнейшем процесс уплотнения осадков исследовался на основе математического моделирования в предположении упругого или упруго-пластического характера деформирования осадков при потере насыщающего флюида [Parasnis , 1960; Magara, 1978; Audet , McConnell , 1992; Wangen, 1997; Luo, Vasseur, Pouya, Lamoureux-Var, Poliakov, 1998; и другие]. Такой подход позволял удовлетворительно моделировать распределение пористости в приповерхностных горизонтах, но не давал возможности учесть зависящую от времени, необратимую компоненту уплотнения осадков, наблюдаемую как в природе, так и в лабораторных экспериментах [Bjerrum, 1967; Vasseur, 1998; Spiers, Schutjens, 1990]. Зависящая от времени, необратимая компонента уплотнения осадков может численно моделироваться на основании решения задач уплотнения осадков в постановке, которая описывает вязкие деформации среды при уплотнении, аналогично, например, моделированию фильтрации расплава под рифтовыми зонами. В дальнейшем в рамках механизма вязкого уплотнения были исследованы несколько специальных стационарных и квазистационарных случаев уплотнения двухкомпонентной флюидонасыщенной среды коры, которые позволили сделать следующий шаг в моделировании геофизического процесса уплотнения осадков [Mckenzie, 1987; Angevine, Turcotte, 1983; Каракин, Суетнова, 1989; Birchwood, Turcott, 1994; Suetnova, Carbonel, Smithson, 1994; Schneider, Potdevin, Wolf and Faille, 1996; Fowler, Yang, 1999; Podladchikov, Connoly, 2000]. Однако, эти исследования носили во многом фрагментарный характер и проводились в упрощенной постановке и не уделялось внимания сложной реологии двухкомпонентной среды осадков, анализу влияния физических и гидродинамических свойств и роли скорости и времени формирования осадочной толщи. Эти исследования, внеся новый вклад в понимание механизмов процессов уплотнения, не давали возможности описать с единых позиций закономерности и особенности эволюции в течении геологической истории основных характеристик процесса уплотнения, таких как распределение пористости, порового давления, скорости флюида и скорости уплотнения осадков в зависимости от условий осадконакопления. При этом могли неадекватно оцениваться как распределения пористости и давления насыщающего флюида и скорость фильтрации насыщающего флюида, так и их эволюция во времени, что может вносить искажения как в интерпретацию геофизических данных так и в реконструкции процессов, проходящих в осадках в течение их геологической истории. Изучение истории формирования наблюдаемых структур осадков и процессов образования в них полезных ископаемых требует моделирования истории уплотнения, то есть динамики падения пористости и фильтрации насыщающего флюида в зависимости от условий осадконакопления, так как миграция флюидов к поверхности в земной коре может приводить к формированию коллекторов, скоплению флюидов в аккумулирующих резервуарах и оставлять следы в вмещающей среде в виде рудных и гидратных проявлений. Таким образом разработка математических моделей эволюции процессов тепломассопереноса в течении геологической истории формирования осадочных структур на основании современных представлений о вязкоупругой реологии двухкомпонентной среды осадков является актуальной задачей.

Цель работы.

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

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

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

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

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

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

Для достижения указанных целей были поставлены и решены следующие основные задачи:

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

2. создание алгоритма численного решения сформулированной математической задачи и его программная реализация;

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

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

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

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

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

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

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

Методика исследования.

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

Личный вклад автора.

Постановка задач, выбор объектов исследования, численное моделирование выполнены автором. Интерпретация результатов по моделированию однородного осадконакопления частично выполнена совместно с Г. Васером (Университет Монпелье, Монпелье, Франция). Асимптотическое решение для скоростей уплотнения в краевой задаче для вязкой среды построено совместно с Чернявским В.М. (МГУ, Москва).

На защиту вынесены следующие основные результаты работы:

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

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

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

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

4. Установленные закономерности взаимовлияния процессов уплотнения осадков и процесса осаждения примесей из насыщающего флюида на примере аккумуляции гидратов из растворенного в флюиде газа в зоне РТ стабильности газовых гидратов в субаквальных условиях и сравнительные оценки возможной гидронасыщенности в зависимости от условий формирования осадочного слоя. По анализу результатов моделирования показано, что гидратонакопление определяется значениями безразмерных критериев подобия задачи уплотнения накапливающихся осадков, так, что скорость аккумуляции гидратов оказывается выше при больших значениях критерия подобия V и критерия подобия D.

Научная новизна.

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

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

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

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

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

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

Научное и практическое значение.

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

Апробация работы и публикации.

Основные результаты работы изложены в научных отчетах, по теме диссертации опубликованы 40 работ, из них 14 в изданиях входящих в рекомендованный ВАК список. Основные результаты работ по теме диссертации докладывались автором на Международном совещании “Новые достижения в области геотермических исследований в скважинах” (Клейн-Корис, Германия, 1993), на Генеральной Ассамблее Американского Геофизического Союза (AGU), 1994, на Генеральных Ассамблеях EGS (Вена, Австрия, 1997; Ницца, Франция, 1998, 2001, 2002), Международных конференциях “Тепловое поле Земли и методы его изучения” (Москва, 2000, 2002), на Третьей международной конференции " Физико-химические и петрофизические исследования в науках о Земле", Москва, 2001, секции “Тепловое поле Земли и методы его изучения” VI, VII, VIII Международных конференциях “Новые идеи в науках о Земле” (Москва, 2003, 2005, 2007), Объединенной Ассамблее EGS-AGU-EUG (Ницца, Франция, 2004), на международной конференции “Structures in the continental crust and geothermal resources” (Сиена, Италия, 2003), ХХ1Х Генеральной Ассамблее Европейской Сейсмологической Комиссии, Потсдам, 2004, V Международной научно- практической конференции “Наука и новейшие технологии при поисках разведке и разработке месторождений полезных ископаемых” ( РГГРУ Москва, 2006г.), на 6х, 7х, 8х, 9х и 10х “Чтениях им. Федынского” (2004, 2005, 2006, 2007, 2008г., Москва), на конференции “Фундаментальный базис новых технологий нефтяной и газовой промышленности” (2007г., Москва), и других.

Структура работы. Диссертация состоит из введения, шести глав, заключения, приложения и списка литературы; содержит 243 страницы текста, 54 рисунка, 2 таблицы. Работа выполнена в лаборатории теоретической геофизики Учреждения Российской Академии Наук Института физики Земли РАН.

Благодарности. Автор считает своим приятным долгом выразить благодарность научному консультанту работы, академику РАН, профессору А.О. Глико за постоянное внимание и поддержку. Работа по теме диссертации поддержана грантами РФФИ № 00-05-65392, № 03-05-65001, № 06-05-652604. Автор выражает искреннюю благодарность соавторам. Автор считает своим долгом вспомнить с благодарностью профессора Е.А. Любимову, под руководством которой автор начинала свою работу в области геофизических исследований и профессора В.П. Мясникова.

КРАТКОЕ СОДЕРЖАНИЕ РАБОТЫ

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

Глава 1. Моделирование процессов тепломассопереноса в флюидонасыщенных областях земной коры (Обзор).

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

Глава 2. Нестационарная задача уплотнения флюидосодержащих осадков вязкоупругой реологии при их накоплении

2.1 Математическая модель задачи.

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

(2.1.1)

(2.1.2)

(2.1.3)

В рамках модифицированного с учетом принципа эффективного напряжения Био-Терцаги для пористой среды, Максвелловского реологического закона, динамическое соотношение для напряжения и скорости деформации записывается как:

(2.1.4.)

Уравнение теплопроводности записывается как

(2.1.5)

где А₁= , А₂= ,

p_f -давление флюида, V_f - скорость порового флюида, V_s -скорость матрицы осадков, t -время, m - пористость, с_f - плотность флюида, C_f - теплоемкость флюида, с_s - плотность осадков, C_s - теплоемкость матрицы осадков, k - проницаемость, м - вязкость флюида, g-ускорение силы тяжести, с - плотность осадoчной среды, f - температура, к- температуропроводность осадков, C -теплоемкость осадков, p_e - эффективное давление p_e = p_tot - p_f , (p_tot-полное давление), з - вязкость осадков при уплотнении [Stevenson, Scott, 1991; Birchwood, Turcotte, 1994; Schneider, et all, 1996], в - эффективная сжимаемость среды. В силу предположения вязкоупругой реологии пористой среды в характеризует сжимаемость пор [Gueguen, Palciauskas, 1994]. В работе далее система уравнений рассматривалась в одномерной постановке, горизонтальный размер слоя предполагается много больше вертикального и накопление осадков на поверхность растущего слоя принимается однородным по поверхности. Рассматриваемая задача не включает рассмотрение в уравнении баланса энергии вклада теплогенерации осадков, и зависимости тепловых свойств от PT условий [Bethke, 1985; Миклашевский, Попов, и др., 2006], что, в рамках сформулированной цели исследования не снижает общности результата. Система (2.1.1-2.1.5) становится полной и замкнутой после формулировки следующих граничных условий (без потери общности рассматривается скомпенсированное осадконакопление): на нижней, непроницаемой границе осадков y =b(t), скорости твердой V_s и жидкой V_f фаз обе равны V₁, - скорости погружения основания. Пористость на верхней границе, являющейся дренажной, постоянна m((y=0),t)=m_b и эффективное давление p_e =0. Температура на поверхности осадков поддерживается постоянной и на нижней границе задается градиент температуры ,

2.2. Физические параметры модели и масштабирование.

Уравнения (2.1.1-2.1.5) включают коэффициенты, которые являются физическими параметрами модели совместно с граничными условиями.

Проницаемость, являющаяся по экспериментальным наблюдениям, нелинейной функцией пористости, обычно в моделях рассматривается как степенная функция пористости [Николаевский, Басниев, Горбунов, Зотов, 1970] и показатели степени варьируют от 2 или 3 [Birchwood, Turcott, 1994] до 8 [Fowler, Yang, 1999], и вопрос вида зависимости проницаемость-пористость продолжает быть предметом исследований. В данной работе принимается зависимость проницаемости от пористости k=k₀ m^l , l=2;3. Величина k₀ зависит от типа осадков. Экспериментальные значения проницаемости для осадков могут варьировать от 10^-12 до 10^-21м² [Файф, Прайс, Томсон, 1981; Gueguen and Palciauskas, 1994]. Используя данные имеющиеся в литературе, в данной работе был принят для k₀ диапазон от 10^-15 до 10^-17м². Вязкость осадочной среды принимается не зависящей от температуры и давления на масштабах исследуемых глубин и варьируется в диапазоне 10¹⁹Па с -10²²Па с [Birchwood, Turcott 1994; Schneider, Potdevin, Wolf, Faille, 1996]. Параметр оценивается как 10^-10 Па^-1 - 10^-8 Па ^-1 [Hart and Wang, 1995; Schneider et al, 1996]. Скорость погружения основания осадков (при скомпенсированном осадконакоплении она определяет скорость роста мощности слоя) V которая, в общем случае, непостоянна в ходе формирования осадочной толщи, согласно литературным данным оценивается как 10^-10мс^-1-10^-12мс^-1 [Audet, 1996; Fowler, Yang, 1999]. Остальные принятые физические параметры модели соответственно равны:? с=2.6?10^-4Па с, с_f=1.0?10³кг/м^-3, с_s=2.65?10³кг/м^-3.

Используя приведенные выше величины входящих параметров и следуя классической теории вязкоупругости [Nadai, 1963, Мельхиор, 1976], мы можем оценить время Максвелловской вязкоупругой релаксации , как 10⁹ - 10¹² с , и следовательно, можно ожидать заметного проявления вязких эффектов уплотнения накапливающихся осадков на временах порядка от сотен до тысяч лет в зависимости от свойств осадочной среды. Предполагая, что механические параметры задачи не зависят от температуры, мы можем исследовать механику и гидродинамику уплотнения. Для построения и исследования решений был произведен переход в систему координат, движущуюся со скоростью погружения основания осадков, и произведена замена переменных и необходимая процедура масштабирования и приведения к безразмерному виду переменных и коэффициентов для выявления параметров подобия [Седов, 1977] так что, масштабом длины является ; масштабом времени является . Масштабом пористости является m_b и масштабы давления и скорости, записываются как P=?сgL, V=L/T.

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

(2.2.1)

(2.2.2)

(2.2.3)

m, V_s и p_e - безразмерные неизвестные и a = 1 / m_b .

В системе уравнений (2.2.1-2.2.3) все коэффициенты и переменные представляют из себя безразмерные величины: z=z/L, t=t/T, p_e=p_e/P,. V_s= V_s /(L/T), a = 1 / m_b , з=з/з_o, D= в_p=в_p Pm_b= з₀ в_p/T. Граничные условия записываются как: V_s(z=0,t)=0 ; p (z=h( t),t)=0; m(z=h(t),t)=1; m(z,t=0)=m_initial(z). Скорость движения верхней границы области V₀= V₀ /(L/T). Безразмерные параметры V₀ = V и D являются нетривиальными критериями подобия сформулированной системы уравнений.

2.3. Численное решение

Для решения полученной нелинейной системы уравнений в частных производных в области с движущейся границей была написана программа для персонального компьютера (приводится в приложении). Система уравнений (2.2.1-2.2.3) решалась численно модифицированным методом конечных разностей. Алгоритм решения задачи в наращиваемой области был построен по полунеявной схеме с контролем переменного шага по времени [Флетчер, 1991; Press et al, 1992 ], обновлением сетки в соответствии с ее уплотнением на каждом шаге по времени, и наращиванием элементов сетки в соответствии с дискретной аппроксимацией скорости роста области и контролем устойчивости [Флетчер, 1991; Press et al, 1992].

На рисунках 2.1 (а,b) показаны результаты расчетов для следующих значений физических параметров задачи: V₀= 10^-11 м с^-1, з_o =510²⁰ Па с, в =10^-9 Па ^-1, k₀ =10^-15м, и конечная толщина бассейна принималась равной 4 км. Такой набор параметров модели приводит к следующим значениям масштабных параметров T и L : T = 710¹³ с, L =4.310³м и масштабированной скорости роста границы V= V₀ = 0.016.



Рисунок 2.1а ,b . Результаты расчетов распределения пористости и порового давления по глубине в различные моменты времени формирования бассейна. Все величины на рисунках безразмерные.

На рисунке 2.2а, б. показаны результаты расчетов для различных моментов времени для следующих значений физических параметров задачи: V₀= 10^-11 м с^-1, з_o =510²⁰ Па с , в =10^-9 Па ^-1 , k₀ =10^-16м , и конечная мощность осадков принималась равной 4 км. Такой набор параметров модели приводит к следующим значениям масштабных параметров T и L : T = 2.210¹⁴ с, L =1.410³м и безразмерной скорости роста границы V₀ = 1.6.

Рисунок 2.2. а ,b . Результаты расчетов распределения пористости и порового давления по глубине в различные моменты времени формирования бассейна. Все величины на рисунках безразмерные.

Численное решение (рис. 2.1, 2.2) показывает, как процесс уплотнения накапливающихся осадков развивается во времени. В результате в верхнем, приповерхностном слое, пористость монотонно уменьшается, а поровое давление растет по глубине, в формирующемся нижнем слое, примыкающем к основанию осадков пористость практически не меняется, оставаясь равной минимальному значению и градиент порового давления близок к литостатическому. В промежуточной области происходит резкое убывание градиента пористости и возрастание градиента давления. В целом, наличие такой переходной зоны по давлению на глубине около 3 км. часто обнаруживается в геофизических исследованиях осадочных бассейнов, например в прибрежном бассейне Луизианы [Mello, 1994], но не может быть объяснено в рамках моделей упругого уплотнения [Hubbert, Rubey, 1959; Mello, 1994].

Рисунок 2.2.3. Обобщенное типичное распределение порового давления по глубине осадков в прибрежном бассейне Луизианы [Mello, 1994]

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

В любой заданной точке осадков в приповерхностной зоне пористость уменьшается с ростом значения параметра безразмерное время. Так как значение параметра подобия ґбезразмерное времяґ определяется не только реальным временем процесса, но и физическими свойствами осадков, то осадки меньшей вязкости или (и) большей проницаемости будут демонстрировать большее падение пористости в течение заданного реального времени, чем более вязкие и (или) менее проницаемые. Из анализа результатов численных решений видно, что градиенты падения пористости и роста давления флюида по глубине зависят от значений безразмерного характеристического критерия V, являющегося нелинейной комбинацией физических параметров осадков и скорости их аккумуляции. А именно: градиент давления насыщающего флюида в приповерхностном слое оказывается выше при большем значении V. Давление насыщающего флюида повышается с увеличением расстояния от поверхности осадков и градиент развивающегося в накапливающихся и уплотняющихся осадках давления флюида тем больше, чем больше скорость осадконакопления и тем меньше, чем больше пористость и проницаемость поступающих осадков. При формировании осадочной толщи в течении заданного времени, градиент давления флюида будет большим для слоя, на поверхность которого поступали менее проницаемые осадки. Безразмерный параметр подобия D в уравнении (2.2.3) характеризует вклад упругого эффекта в процесс уплотнения. Этот параметр также влияет на распределение порового давления и пористости по глубине [Суетнова, 2008 ], так, что градиент падения пористости по глубине осадков тем больше, чем больше значение D. Для наглядного выделения роли критерия D на рис.2.2.4 приведены в безразмерном виде результаты расчетов распределения эффективного давления по глубине накапливающихся осадков для одинаковых значений параметра подобия V (V=0.06), и различных значений параметра подобия D.



Рис.2.2.4.Распределение эффективного давления по глубине накапливающихся осадков для одинаковых значений параметра подобия V (V~0.06 ) и различных значений параметра подобия D (D~ 0.002-сплошная линия, D ~2.- пунктирная линия).

Такие значения параметров подобия V и D соответствуют следующим значениям реальных параметров исследуемого процесса: 1) V₀= 10^-11 м с^-1, з_o =110²⁰ Па с , в =10^-10 Па ^-1 , k₀ =10^-15м , m₀=0.3 и финальной мощности осадочного слоя равной 8км; 2) V₀= 10^-11 м с^-1, з_o =110²¹ Па с , в =10^-7 Па ^-1 , k₀ =10^-15м , m₀=0.3 и той же финальной мощности осадочного слоя равной 8км. На рисунке видно, что распределение эффективного давления по глубине зависит от значения параметра подобия D. Чем больше значение D, тем эффективное давление оказывается меньше, и, следовательно, давление насыщающего флюида - больше, для одинаковых значений гидродинамического параметра подобия V и реальной мощности накопленных осадков. Кроме того, глубина перехода давления насыщающего флюида к около литостатическому давлению оказывается большей для большего значения D. Это означает, что число Деборы (De=D/t) оказывается меньше для меньшего значения D и, следовательно, вязкие эффекты в этом случае проявляются больше, что и демонстрирует рисунок 2.4.7. На следующих рисунках 2.2.5а,б представлены графики распределения пористости и порового давления по осадочной колонке для различных значений V и D при одинаковой финальной мощности осадков, которые наглядно показывают различие в градиентах падения пористости и роста порового давления по глубине в зависимости от значений этих критериев подобия. Значения критериев подобия определяются следующими значениями физических и гидродинамических параметров осадков: V₀= 10^-10 м/с, з= 5 ?10²¹ Па· с, k=10^-15 м², в=1/K_p= 10^-9 Па ^-1 , m₀ =0.3, V= 0.576, D =0.0649, t=0.766 (обозначение ромбы); V₀=10^-10 м/с, з=5? 10²¹ Па· с, k=10^-15 м², в=1/K_p= 10^-8 Па ^-1 , m₀ =0.3, V= 0.576, D =0.649 , t=0.78 (обозначение квадраты); V₀= 10^-10 м/с, з= 5? 10²⁰ Па· с, в=1/K_p= 10^-9 Па ^-1, k=10^-14м² , m₀ =0.3, V= 0.0576, D =0.0649, t=7.77 (обозначение треугольники); V₀= 10^-10 м/с, з= 5 ?10²⁰ Па· с , в=1/K_p= 10^-8 Па ^-1 , k=10^-14 м² , m₀ =0.3, V= 0.0576 , D =0.649, t=7.80 (обозначена сплошной линией).

Рисунок 2.2.5а,б. а) Распределение пористости по глубине для скорости осадконакопления 10^-10 м/с. и различных значениях физических и гидродинамических свойств осадков. Обозначения в тексте. б). Распределение порового давления по глубине осадков для скорости осадконакопления 10^-10 м/с. и различных значениях физических и гидродинамических свойств осадков. Обозначения: крест-литостатическое давление флюида, косой крест- гидростатическое давление, остальные обозначения как на рисунке 2.2.5а.

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

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

Глава 3. Асимптотическое решение задачи уплотнения слоя осадков в процессе роста его мощности.

Методы построения решения в виде асимптотических разложений являются эффективным инструментом построения и анализа решений дифференциальных уравнений, однако, существуют классы задач, в которых в уравнениях или в граничных условиях присутствует малый параметр, что приводит к непригодности на больших временах решений, полученных с помощью прямых асимптотических разложений по времени [Найфе, 1976]. Аналитические методы построения асимптотических решений, равномерно пригодных на всем диапазоне исследуемых времен задачи, чрезвычайно важны при решении многопараметрических нестационарных задач, в которых нахождение характерного масштаба времени уравнений или системы уравнений и сопоставление его с характерным временем изменения граничных условий или роста области может выявить принципиальные характерные свойства и особенности поведения решений. В этой главе приводится построение асимптотического решения задачи о течении свободного флюида в вязкодеформируемой пористой среде с подвижной границей полученное с использованием идей построения равномерно пригодных асимптотических решений для дифференциальных уравнений с малым параметром [Коул, 1972; Найфе, 1976; Бэтчелор, 1973].

3.1 Математическая постановка задачи и решение.

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

(3.1)

(3.2)

(3.3)

(3.4)

p= -(с_f + Дс)gy-p_f (3.5)

Здесь используются следующие обозначения для неизвестных величин: v - скорость осадков, v_f - скорость насыщающего флюида, и для независимых переменных: времени - t и вертикальной координаты с положительным направлением вверх - y, остальные обозначения как в главе 2. Уравнения (3.1-3.5) должны удовлетворяться в области, заключенной между нижней границей y=0 и верхней границей y=U?t, движущейся вверх с постоянной скоростью U. Граничные условия: неподвижная нижняя граница является непроницаемой и на ней скорости флюида и матрицы равны нулю, и на верхней границе пористость и давление равны заданным постоянным значениям.

y=0: v _f =0 , v=0 (3.6)

y=U?t: m=M, p=0 (3.7)

Соотношения (3.1) - (3.7) представляют полную формулировку задачи. В соответствии с современными представлениями принимается, что проницаемость является степенной функцией пористости. Система уравнений (3.1) - (3.5) заметно упрощается при использовании обычного предположения теории пористых сред: m<<1. Из уравнений (3.1), (3.2) и условий на границе (3.6), (3.7) следует интеграл

mv_f +(1-m)v=0 (3.8)

Подстановка выражения для давления флюида p_f из (3.5) в уравнение (3.2) с учетом соотношения . (3.8), приводит к системе уравнений

(3.9)

(3.10)

(3.11)

с краевыми условиями

y=0: v _f =0 , v=0 (3.12)

y=U?t: m=M, p=0 (3.13)

Решив задачу (3.9-3.13) можно определить скорость и давление флюида из соотношений (3.1) и (3.5). Исходные параметры задачи могут варьироваться в широких пределах, некоторые величины в интервале до двух порядков (см. главу 2). Для построения решения выбираем функцию k(m)=k₀?m^l, при l=2, и характерные значения параметров в единицах СИ: 0.1<M<0.2, 10^-12<U<10^-10, 10²⁰< з₀<10²¹, м=2.5?10^-4, 10^-15<k₀<10^-14, с _f=10⁺³, с=2.6?10³, g=10. Время, в течение которого развивается процесс, определяется глубиной бассейна H или задается непосредственно и составляет 1-10 млн лет. Характерные масштабы вводятся с помощью величин: скорости движения границы U, ускорения свободного падения g, разности плотностей матрицы и флюида с-с_f, а также пористости на верхней границе M. Масштабы длины, давления и времени, записываются соответственно, в виде L=(з₀U/g(с-с_f))^1/2, P=(с-с_f)gL, T=L/U. Величины, измеряемые в этих масштабах, связаны с исходными величинами соотношениями

y=L?y`, t=T?t`, m=M?m`, v=M?U?v`, p=P?p`, е=мU/(k₀(с-с_f)gM)

В новых переменных, у которых для краткости записи опущены штрихи, задача при t>>1 записывается в виде

(3.14)

(3.15)

(3.16)

с краевыми условиями

y=0 :v=0 ; (3.17)

y=t: m=1, p=0 (3.18)

Величина е, оцениваемая по значениям исходных данных, оказывается заключенной в интервале 0.8?10^-4<е<1.6?10^-2 и во всяком случае можно считать, что е<<1, то есть е является малым параметром сформулированной задачи.

Представим решение s=m,v,p в виде асимптотического ряда по степеням е

s=s₀+е?s₁+... (3.19)

Подставляя разложение (3.19) в задачу (3.14-3.18 ), и ограничиваясь нулевым приближением, получаем аппроксимацию решения в виде

p₀ = t-y

m₀ = exp(-0.5(y-t)²) (3.20)

v₀ = -m₀ + exp(-0.5t²)

Решение (3.20) при е=0 удовлетворяет краевым условиям (3.17-3.18)), но не является равномерной аппроксимацией при е не равном 0 и t>>1. Для того, чтобы убедиться в этом, рассмотрим уравнение (3.15)) в следующем приближении - порядка е, в котором оно записывается как

?p₁/?y = -v₀/m₀² (3.21)

Если правая часть (3.21) оказывается величиной порядка 1/е, то величина е*p₁ не является асимптотически малой по сравнению с p₀ и разложение (3.19 ) неравномерно. При t>>1 и /y-t/~t значения величины m₀~е, и v₀~е, и правая часть (3.21) оценивается как ~ 1/е. Таким образом, решение (3.20 ) является нулевым приближением внешнего разложения по е задачи (3.14-3.18), которое аппроксимирует точное решение вблизи верхней границы, пока пористость остается величиной порядка единицы и следовательно, необходимо построить соответствующее решение вблизи нижней границы и произвести процедуру сращивания (сшивания) решения. Исходя из предыдущих оценок, следуя теории [Коул, 1972; Найфе, 1976], для нахождения решения введем новые масштабированные переменные

...