CC BY
56
В. Г. Румынии, Л. Н. Синдаловский, С. Э. Макашов, А. В. Боронина
НОВЫЕ АНАЛИТИЧЕСКИЕ ЗАВИСИМОСТИ
ДЛЯ ОПИСАНИЯ ПОДТЯГИВАНИЯ ГЛУБИННЫХ РАССОЛОВ
К КОНТУРУ НЕСОВЕРШЕННОЙ СКВАЖИНЫ
Введение
Интерес к исследованию процесса субвертикального подтягивания границы раздела соленых и пресных вод (исходное положение границы — горизонтальное) к несовершенному контуру дренирования (чаще всего — это скважина с коротким фильтром) водоносной толщи не ослабевает у специалистов-гидрогеологов в течение нескольких последних десятилетий. Об этом свидетельствуют регулярно появляющиеся публикации, посвященные проблемам формирования качества воды на водозаборах, функционирующих в условиях неоднородной профильной гидрохимической зональности. Такого рода условия встречаются, например, в прибрежных водоносных горизонтах, при-подошвенная часть которых занята интрузией морских или океанических вод. В артезианских бассейнах эта зональность чаще всего ассоциируется с естественным ростом минерализации подземных вод с глубиной. Еще один пример — линзы пресных вод инфильтрационного происхождения в горизонтах, исходно содержащих соленые воды и рассолы (например, в пределах океанических островов). С похожими проблемами приходится сталкиваться и при анализе формирования качества дренажных вод в платформенных областях на разрабатываемых месторождениях полезных ископаемых [1].
Большинство работ теоретического характера было посвящено поиску аналитических выражений для определения так называемого предельного дебита скважины (ф*), превышение которого приводит к подтягиванию границы раздела разноплотностных вод к водозаборному устройству. В то же время почти не исследованным оказался круг задач, описывающих изменение концентрации соли в воде, отбираемой скважиной (за исключением случаев, когда различиями в плотностях жидкостей можно пренебречь), по времени.
Целью настоящей работы является получение приближенных аналитических выражений, описывающих как равновесный купол соленых вод под скважиной (когда ф < ф*), так и нестационарную фазу процесса (при ф > ф*) —фазу интенсивного засоления водозаборной скважины. Для этого привлекались как аналитические, так и численные методы решения задач фильтрации и миграции подземных вод.
1. Равновесный купол соленых вод под несовершенной скважиной: предельный дебит и критическая величина подъема купола соленых вод
Эксплуатация несовершенной скважины, фильтр которой вскрывает верхнюю часть водоносной толщи (или, в более широком смысле, к несовершенному контуру дренирования водоносной толщи), приводит к формированию «купола» («языка») соленых вод (рис. 1). В реальных условиях язык соленых вод может либо прорваться в скважину, что приводит к нарастанию в откачиваемых водах содержания водорастворенных
© В. Г. Румынин, Л. Н. Синдаловский, С. Э. Макашов, А. В. Боронина, 2010
солей, либо занять относительно стабильное положение, не достигнув ее фильтра. Это зависит от дебита скважины, различий в удельных весах пресных и соленых вод, а также параметров, характеризующих скважину и пласт. Максимально возможный дебит ^*), при котором не происходит прорыва соленой воды в скважину, называется предельным дебитом; этой характеристике процесса отвечает величина Аг* — высота «купола» соленых вод в предельном равновесии (см. рис. 1). Несмотря на повышенный интерес многих исследователей к проблеме куполообразования и большое количество опубликованных работ, точной теории соляного конуса не существует ввиду чрезвычайной сложности задачи.
+ + + + + + + + + + + +
Рис. 1. Концептуальная схема образования «купола» соленых вод.
Так, оценкам предельного дебита скважины, о которых пойдет речь в данном разделе, посвящено большое количество ранних публикаций в нефтяной гидрогеологии [2, 3, 4]. С тех пор выполнены многочисленные исследования в приложении к задачам подтягивания глубинных минерализованных вод к несовершенным контурам дренажа водоносных горизонтов [5-13].
1.1. Постановка проблемы и анализ существующих решений
С математической точки зрения решение задачи о субвертикальном подтягивании профильной границы соленых вод к скважине может основываться на двух подходах (методах).
Первый подход базируется на предпосылке о стационарности (неподвижности) нижней, более тяжелой, жидкости, так что ее контакт с верхней, более легкой, жидкостью является для последней непроницаемой границей. Полагается, что наличие этой границы (в форме конуса) под скважиной не влияет на распределение давления (напора) в области, занятой легкой жидкостью. В соответствии с соотношением Гибена-Герцберга подъем контакта Azr (см. рис. 1) определяется зависимостью [8, 14, 15]:
Azr = Ъ — z(r) = S(r, z = 6), (1)
где Ь — начальная мощность зоны, занятой легкой (например, пресной водой) жидкостью и, г — радиальная координата [Ь], Ар — плотностный градиент (Ар = (ря — Р!)/Р!, Ря и pf —плотность пресной и соленой воды [МЬ-3]), Б — понижение напора [Ь].
Гидростатическая модель (1) определяет пространственные координаты линий тока, вдоль которых частицы пресной воды перемещаются по поверхности соленой воды в предположении, что изменение напора в зоне пресных вод приводит к мгновенному достижению нового равновесного положения границей раздела.
Однако модель (1) адекватно описывает соотношение между Агг и Б (г, г = Ь) только до тех пор, пока не будет превышено некоторое критическое значение Аг* = Агг (г = г«0, гт —радиус скважины. В этом контексте Аг* —высота, при подъеме на которую конус под скважиной становится нестабильным (теряется линейная связь между функцией Агг и величиной понижения Б). Нестабильность контакта приводит к тому, что соленая вода весьма быстро («скачком») подтягивается к водозаборной скважине.
Характеристика Аг* обычно определяется в терминах в (Ь — I), где в — априорно задаваемое отношение между критической величиной подъема и расстоянием между начальным положением контакта (г = Ь) и забоем скважины (г = I). Поэтому уравнение (1) может быть переписано в виде соотношения
Аг* = е(Ъ-1) = -^Б(г1и,г = Ъ), (2)
АР
используемого и для нахождения критической величины ^*. Значение в, согласно различным исследованиям [8, 9, 14-20], может меняться от 0,25 до 0,89.
Критический дебит водоотбора служит одной из интегральных характеристик предельного равновесия поверхности раздела жидкостей. Оценка этой характеристики имеет большое практическое значение: это максимальный дебит водоотбора, происходящего без захвата скважиной соленой воды. Для этого используются разнообразные аналитические решения, полученные для различных условий относительно функции Б(гт, г = Ь) = /(^) —см., например,[8].
Неопределенности расчетов, связанных с выбором параметра в, можно избежать, если при определении пространственных координат поверхности раздела жидкостей в предельном положении допустить необходимым одновременное выполнение двух условий [9, 21]:
первое — давление на контакт разновесомых жидкостей со стороны пресной воды (Р!) при снижении напора в ней на величину Б = Б (г, г) равно давлению на контакт со стороны соленой воды (Ря); использование соотношения Гибена-Герцберга [5] для этого условия приводит к равенству вида
Аг* = Б (г, г = Ъ — Аг*); (3)
АР
второе — вертикальные градиенты давления на поверхности раздела жидкостей в состоянии равновесия также должны быть равны, дР!/дг = дРя/дг, или
аз(г,г = ь-д*;) =
дг
Физически последнее условие означает, что потенциальная скорость гравитационного опускания контакта под влиянием силы тяжести гир = Ар равна скорости вертикаль-
ной фильтрации воды в точке (г, г) гг = — дБ/дг, кг —коэффициент фильтрации в вертикальном направлении.
Второй подход получения соотношений для описания процесса подъема купола соленых вод основывается на решении уравнений Лапласа для функции пьезометрического напора [5, 11, 17] или давления [22] для двух областей, имеющих подвижную границу раздела.
Так, уравнения для описания движения границы раздела Аг(г, £) под влиянием водоотбора из скважины были получены методом линеаризованного пертурбационного анализа. Для случая короткой (в сравнении с мощностью зоны пресной воды) скважины в полуограниченном по мощности пласте положение вершины конуса (г = 0) может быть определено из формулы [17, 19]:
2п Аркг (Ь — I)
1- 1
1 + £
+ АРк* + /К\
‘ = (5)
справедливой, как показано на численных моделях, при выполнении условия Аг/(Ь — I) < 1/4 [23]. Из решения (5) при £ ^ то и Аг* = в(Ь — I) может быть получен критический дебит:
Q* = (в)2пкг А р(Ь — I)2. (5а)
Представленные соотношения для Q* и А г* позволяют описать частное состояние системы — состояние предельного равновесия; некоторые из этих соотношений могут быть модифицированы и преобразованы к виду, представляющему определенный практический интерес.
1.2. Расчет предельного дебита и критической величины подъема купола соленых вод
Как уже отмечено, динамическое равновесие контакта предполагает равенство давлений и их тангенциальных составляющих в каждой точке контакта со стороны двух жидкостей, т. е. одновременное выполнение двух условий (3) и (4), которые имеют следующие безразмерные представления:
Агг = —, (6)
й = “АА (7)
где Б = Б (г, Р*)/Ь, Р = г/Ь, Р = г/Ь, г* = 1 — Аг*, Аг* = Аг^/Ь.
Для определения критической величины Аг* и соответствующего значения Q* достаточно ограничить дальнейший анализ рассмотрением случая г = гт.
Функция безразмерного понижения Б = Б (р ,г) и ее производная дБ/дг для условий стационарной фильтрации в пласте с круговой границей I рода (при г = Д(или, в безразмерном виде, Р = К)) определяются выражениями [24, 25]
5 = 4^<21"7 + л>’ (8>
- = -2-/' (9)
дг А-кТУ31 К ’
где
4 °° 1
/в = /^г, г) = —- — вт (п7г /?) соб(п7г г) Ко(птг г),
п в п
' п=1
4
/' = /з(»’, 8т(шг/?) 8т(шг.г) Ко(пм'г),
в п=1
(10)
(11)
г = -г, Д = 4,
кь ко'
/3 = Т = кгЪ, к = у/кг/к~—коэффициент анизотропии (кг и к2 — коэффициенты фильтрации по напластованию и вкрест напластования породы); здесь и далее Р = Р(р), так что безразмерная величина подъема купола соленых вод под скважиной (см. рис. 1) Ар * = 1 — Р (гw). Заметим, что при малых Р ряды (10) и (11) являются плохо сходящимися.
Подстановка уравнений (8) и (9) для критических условий (когда г = = 1 — Д,г*
и Р = в исходную систему уравнений (6) и (7) преобразует последнюю к виду
Q*
Q*
2п (1 — г )
Т Ъ АП 1д _Д /Лгт>~) 2
Q*
Q*
2п
ть ар
2
(12)
(13)
Кривые на рис. 2 дают представление о решении системы уравнений (12) и (13), полученном численными методами. В частности, выполненные расчеты дали значение в = в диапазоне 0,75 (/3 = 0,1)—0,9 (/? = 0, 7) при Д = 100. Эти значения превыша-
ют часто принимаемые величины в = 0, 3. Кривые на рис. 2 могут интерпретироваться и под углом зрения влияния анизотропии на расчетный предельный дебит. Так, при фиксированных К и I уменьшение проницаемости в вертикальном направлении, ведущее к падению значений безразмерного комплекса К, дает увеличение безразмерного дебита Q*.
Рис. 2. Расчетные значения <3*(а), г* (б) при различных В (ось абсцисс) и в (числа у кривых); = 0, 001.
0.7
0.6
0.3
0.2
0.1 0.05 0.2
__________________0.3
0.5 _______________________0.4
0.5
0.4
1*0 0.6
0.7
0.8
0.1 0.9
0^—.............................. —. —. —.
1Е-005 0.0001 0.001 0.01 0.1
Рис. 3. Изменение относительного предельного дебита 3* в зависимости от при различных в (числа у кривых); В = 100.
Далее, основываясь на решении уравнений (12) — (13), можно показать, что в практически значимой области значений Рw < 0,01 процесс слабо чувствителен к абсолютным значениям собственно параметра Рw, т. е. к радиусу водозаборной скважины (рис. 3). Можно предположить, что такое поведение системы возможно, если функция имеет альтернативную структуру, которая включает логарифмический член ~ 1п(р). Действительно, при в < 0,4—0, 5 и т^/в < 1 функция ]а может быть аппроксимирована конечной аналитической формулой
/8 = 21п + 21п /3 + А(Е), (14)
в
где
0 8 0 18 1
А(Е) = 0,1 + -г- + -V + 4/3 ехр(- - 1) • 1п(1 + 13). (15)
г г 2 г
Это позволяет трансформировать уравнения (12) и (13) к виду
✓=,* 27Г (1 — г)
д =Й1)7Ф ’
«* = -Йт т
2
А'(г) = /' = — [^- -|—• 1п(1 + /3) ехр(- — 1)]. (18)
г 2 г 3 г 2 г
Решение аппроксимирующей системы уравнений (17) и (18) довольно хорошо согласуется с точным решением системы (12) и (13), что подтверждается данными табл. 1. В табл. 1 приведены результаты расчетов по приближенной формуле
(1!1)
где
Таблица 1. Величина относительного предельного дебита (<Э*) при различных параметрах системы для различных методов его расчета (к = 1)
Д = 10 Д = 100
/3 = 0,1 /3 = 0,2 /3 = 0,3 /3 = 0,4 /3 = 0,1 /3 = 0,2 /3 = 0,3 /3 = 0,4
0,01 0,92 0,88 0,81 0,73 0,61 0,57 0,53 0,47
0,001 0,92 0,88 0,81 0,73 0,60 0,57 0,53 0,47
0,0001 0,90 0,87 0,81 0,70 0,60 0,56 0,52 0,47
(17)—(17) 0,94 0,89 0,82 0,75 0,62 0,58 0,55 0,51
(19) 1,12 0,97 0,84 0,72 0,61 0,55 0,50 0,44
точность расчета предельного дебита по формуле (19) в диапазоне 0,2 < в < 0,8 не ниже 10-15%.
2. Изменение концентрации в откачивающей скважине
Описание процесса в упрощенной постановке может отталкиваться от известных решений задачи для потоков постоянной плотности (Др = 0), которые оказываются полезными и для построения приближенных зависимостей, описывающих поступление в скважину растворов повышенной плотности из глубоких частей разреза (Др > 0).
2.1. Аналитическое решение: случай Др = 0
Пусть смещение природного контакта вод, различающихся по качественным показателям и занимающих изначально нижнюю и верхнюю зоны пласта с постоянной мощностью Ь и т — Ь (т — суммарная мощность), происходит под влиянием водоотбора из несовершенной скважины с фильтром конечной длины I, примыкающим к кровле пласта (рис. 4, а).
Рис. 4- Откачивающая скважина в условиях профильно неоднородного исходного концентрационного распределения.
а — две зоны гидрохимической неоднородности, б — аппроксимация переходной зоны кусочно-неоднородной функцией С*.
Первостепенный интерес представляют расчетные зависимости, отражающие: 1) время подтягивания (по кратчайшей линии тока) первых порций вещества к водозаборной скважине от границы раздела вод г = Ь (см. рис. 4, а); 2) изменение содержания вещества (С — безразмерная концентрация некоторого маркирующего компонента) в откачивающей скважине. Влияние гидродисперсии обычно затушевывается разновременным подтягиванием к откачивающей скважине вещества в условиях сильной деформации сетки движения подземных вод, и поэтому имеет смысл рассматривать миграцию в рамках расчетной схемы поршневого вытеснения. В основе математических построений лежат уравнения квазистационарной фильтрации, дополненные кинематическими уравнениями.
Решение для скважины в однородной полуограниченной в разрезе (Ь ^ т) анизотропной толще с исходной (на глубине Ь) горизонтальной границей раздела вод (см. рис. 4, а), различающихся содержанием (С) некоторого маркирующего компонента (Со —в верхнем слое и С1 —в нижнем слое), имеет вид [22]
с(г) = 1 - (20)
где С(£) = [С(£) — Со]/(С1 — Со), а функция /о = / (0, £) находится из трансцендентной формулы
£ = в к2
(1 — /о3) в2
— 3/о (1 — /о)
(21)
/о3
Решение (21а) справедливо при выполнении условия
Е > £о = к2 (1 + 2в3 — —3в2) , (21а)
где Iо — безразмерное время подтягивания границы раздела вод к скважине,
к= \/кг/к2 — коэффициент анизотропии (кг —коэффициент фильтрации в вертикаль-
ном направлении); в = 1/Ь (I- длина фильтра); ^ —объемный расход скважины;
Е = гдо; £о = -£оДо;
. 2ппЬ3 ._.
г° ~ зд ^
— «нормирующая» временная характеристика, отвечающая времени подтягивания границы раздела вод в изотропном пласте (к =1) к скважине с коротким фильтром (в = 0 — точечный источник) [26].
Анализ графических зависимостей С(£) (рис. 5) показывает, что процесс в целом сильно «растянут» во времени: кратковременный период достаточно быстрого нарастания концентрации сменяется длительным периодом весьма медленных изменений функции С(£). Анизотропия фильтрационных свойств проявляется в смещении графиков по оси £" примерно на величину кг/к2.
Полное решение задачи, описывающее динамику смещения контакта вод Е = /(е, Е) (куполообразование), можно представить в следующем трансцендентном виде [22]
£ (г , г ) = 6 в к2
4(1 — Е3) в2 к3 М
“7ТА1- г>
3М3 4к
где М = у/г2 + к2 (г + /З)2 — у^2 + к2(Е — /З)2, г = г/Ь, г = г/Ь.
(
Рис. 5. Графики функций /о (а, в) и С (б, а, б) и анизотропном (к = л/10— в,
г) для солепереноса в изотропном (к =1 — г) пластах. Числа у кривых — параметр в*
В качестве примера на рис. 6 приведены результаты расчетов, иллюстрирующие, в частности, определяющее влияние анизотропии фильтрационных свойств пласта на концентрационные распределения. Хорошо видно, что относительное снижение проницаемости в вертикальном направлении (например, за счет наличия в разрезе небольших глинистых прослоев) существенно снижает скорость подтягивания границы раздела к фильтру эксплуатационной скважины. На больших моментах времени влияние анизотропии сглаживается.
Анализ решений, полученных для вычисления характеристики £о при расположении скважины в ограниченном по мощности пласте [27, 28], показывает, что влиянием подошвы пласта на величину £о, рассчитываемую по (21а), заведомо можно пренебречь, если толщина (Ь) верхнего слоя воды (с концентрацией маркирующего компонента Со) не превышает половины мощности (т) пласта, т. е. когда Ь/т < 1/2.
Наконец, можно показать, что при в = 0 решение (20) преобразуется к виду [29]
С(£)= / (Мо), (23)
г/ ч (£оА1/3 2 I 2пк2пЬ3
/ (#, #о) = 1 - ^Jj , #о = щ—, (23а)
Скважина
Рис. 6. Характер деформации границы раздела вод под скважиной — решение (22).
Сплошные линии—изотропный пласт (к = 1), пунктирные — анизотропный пласт (к = у/5). (3 = 0,3. Числа у кривых — безразмерное время
или, в абсолютных значениях концентрации
С(г) = Со + (С1 - Со) IМо). (24)
При наличии переходной зоны от пресных вод к минерализованным (рис. 4, б) используется суперпозиция решений вида (20) и (24). Так, для скважины с коротким фильтром (в = 0)в соответствии с (24) имеем [29]
I
С (г) = Со + X) (С - С<_1) I (г, г_),
(25)
где I — номер нижнего расчетного слоя (I +1 — общее число слоев с постоянной концентрацией С*, аппроксимирующих изменение содержания компонентов по мощности пласта); I — текущий номер слоя (I = 0 отвечает верхнему слою с фоновой концентраций Со) — см. рис. 4, б; г*_1 —время начала поступления минерализованных вод от верхней границы г-го рас-
четного слоя, залегающей на глубине Ь*_1; оно определяется по второй формуле (23а), в которой Ь необходимо заменить на Ь*_ 1, т. е.
гг— 1 —
2пк2п Ь3_1
зд '
(25а)
Одним из факторов, заметно смещающим представленные здесь оценки, является плотностная конвекция, обусловленная наличием градиента плотности на контакте вод. Влияние этого фактора на миграционный процесс в целом рассматривается в следующем разделе.
2.2. Исследование процесса на численных моделях и аппроксимирующие зависимости: случай Ар > 0
В последние годы математическое моделирование является одним из основных инструментов исследования процессов, связанных с подтягиванием соленых вод к контурам дренажа. В задачах методического характера повышенное внимание уделялось исследованию влияния плотностного градиента и дисперсионных эффектов на концентрационные распределения в водоносных толщах [12], а также оценивалась возможность использования приближенных решений для описания поведения границы раздела жидкостей в «докритическом» состоянии, Q < Q* [13]. Обращение к математическому моделированию позволило показать, что предпосылка о резкой границе раздела вод смещает оценки предельного дебита в сторону его увеличения и, наоборот, недооценивает высоту критического поднятия поверхности раздела при заданном дебите. Наличие переходной зоны требует более низких дебитов водоотбора для того, чтобы качество воды оставалось на приемлемом уровне, — в сравнении с ситуацией резкой границы раздела [6, 12]. Было отмечено, что при наличии переходной зоны обратный
процесс — опускание «купола» соленых вод после остановки скважины — сильно растянут во времени (возвращение границы раздела в исходное положение происходит исключительно медленно, более того, полного совпадения с первоначальным концентрационным профилем никогда не достигается). Эту особенность процесса необходимо иметь в виду, прогнозируя реабилитационные мероприятия на водозаборах: повторное включение эксплуатационной скважины после ее засоления оправдано спустя весьма продолжительные периоды времени. В то же время очевидно, что степень проявления такого рода гистерезиса зависит от исходных различий в плотности вод (Др). Наконец, движение легкой воды по поверхности рассола способствует возникновению циркуляционной ячейки в пределах «купола».
В дополнение к представленному выше анализу в настоящем разделе делается попытка нахождения универсальных (в рамках принятых допущений, о которых речь пойдет далее) аналитических зависимостей, аппроксимирующих концентрационную функцию в случае совместной фильтрации в пласте растворов, имеющих изначально резкую границу раздела и различающихся исходной плотностью (см. рис. 4, а). «Эталонные» решения ищутся численными методами в широком диапазоне параметров, причем рассматривается ситуация, когда дебит несовершенной скважины превышает критический дебит, а гидродисперсионным рассеянием вещества можно пренебречь. Решение для случая Др = 0, рассмотренное в подразделе 2.1, оказывается полезным при обосновании структуры аппроксимирующих формул.
2.2.1. Описание модели, исследование физических закономерностей транспортного процесса и расчетные варианты
В основе исследования поведения концентрационной функции C(t) лежат частные решения соответствующей миграционной задачи, полученные с помощью численного моделирования на программе TOUGH2 [30, 31]—основной объем вычислений и SUTRA 2.1 [32]—контрольные просчеты для частных случаев. Модель имитирует сопряженные фильтрационный и миграционный процессы в однородной по проницаемости полуограниченной по мощности толще (рис. 7). В табл. 2 представлены параметры для базового («В») варианта расчетов. Последующие варианты (1-27) отличались от базового вариациями тех или иных параметров (табл. 3), таких как: 1) расход откачки
Рис. 7. Модельная область (в r — z-сечении).
а — представление одноразмерными блоками (^г и — номера блоков в г- и ^-направлениях); б —ре-
альная сетка (в м): Аr(Nr = 1) = 0,1, Аr(Nr = 30) = 2990,9, Аz(Nz = 90) = 2, Аz(Nz = 1) = 73,4. Штрихпунктирная линия — исходное положение границы раздела пресной и соленой воды.
Коэффициент фильтрации, кг = кг 0,5 м/сут
Пористость (п) 0,05
Плотность рассола при нормальных условиях, р® 1070 кг/м3
Плотность пресных вод при нормальных условиях, р^ 998,8 кг/м3
Плотностной градиент на контакте вод, Ар 0,072244
Расход скважины, <5 1000 м3/сут
Длина фильтра, 1 2,0 м
Мощность пресных вод, Ь 98,97 м
Мощность соленых вод, га3 1700 м
Радиус влияния, К 10000 м
Таблица 3. Характеристика модельных вариантов
Р» кг/м3 До Расход, (), м3/сут (к = 1,0, /3= 0,02, пг = 99 м) Анизотропия, К 03= 0,02, т = 99 м, й = 1000 м3/сут) (к = Параметр /? 1,0, т = 99 м, Є = 1000 м3/сут) Мощность, Ь, м (к= 1,0, 0, = 1000 м3/сут)
о •л СО о о о 1000 1250 2500 о о о сГ о„ о сГ (Ч сч с? сГ о" о *3- о “■> II 2 (Ч О' о" * II д. о ^ о <4 II 2
0 + +
1000 0,002 7 7
1020 0,022 12 8 17 8
1050 0,052 13 9 18 9
1070 0,072 1 2 3 в 4 5 6 14 В 19 22 23 24 25 26 В 27*
1100 0,102 15 10 20 10
1150 0,152 16 11 21 11
* Расход скажины — 5000 м3/сут.
(табл. 4), 2) плотность рассола, 3) коэффициенты фильтрации («прямая» и «обратная» анизотропия фильтрационных свойств, 4) длины фильтра скважины, 5) мощности зоны пресных вод (табл. 4,). Принятые соотношения Ь ^ ш8 (Ь = ко, шя — мощность слоя соленой воды, см. рис. 7, а) позволяют пренебречь влиянием нижней границы модели на характер процесса, что было подтверждено дополнительными расчетами. Поэтому все дальнейшие построения касаются случая полуограниченного по мощности пласта. Соотношение Q/Q* рассчитывалось исходя из численного решения системы уравнений
(12) (13).
Одной из особенностей программы ТОиСИ2 является автоматический учет влияния сжимаемости раствора за счет гидростатического давления на плотность раствора (р = /(г)), что приводило к отклонению значений рf и ря на контакте жидкостей (в точке г = Ь) от нормальных (при атмосферном давлении) величин (р°°и р°). Соответствующий градиент плотности Др, используемый в дальнейшем анализе, рассчитывался по зависимости Др = (ря — рf )/рf (естественно, это значение несколько отличается от величины (р° — р0 )/р0).
Миграционный процесс исследовался в предположении о пренебрежимо малой роли гидродисперсии (6Г =0, 6т = 0). Корректное решение задачи предполагало в этом случае также исключение заметного влияния численной дисперсии на распределение
0(1).
Предварительно расчеты показали, что разбивка области в направлении г и г
№ варианта Я, м3/сут Я/Я* Ра, Кг/м3 (к = 1, О)1 Ра, кг/м3 (к = 10)2 Ра, кг/м3 (к = 0,1)3 Ар 1, м Ь, м /3
Базовый (В) 1000 4,60 1070 0,072 2 98,97 0,020
1 350 1,61 1070 0,072 2 98,97 0,020
2 500 2,30 1070 0,072 2 98,97 0,020
3 750 3,45 1070 0,072 2 98,97 0,020
4 1250 5,74 1070 0,072 2 98,97 0,020
5 2500 11,49 1070 0,072 2 98,97 0,020
6 5000 22,98 1070 0,072 2 98,97 0,020
7 1000 160,6 1000 0,002 2 98,97 0,020
8 1000 14,99 1020 0,022 2 98,97 0,020
9 1000 6,36 1050 0,052 2 98,97 0,020
10 1000 3,25 1100 0,102 2 98,97 0,020
11 1000 2,18 1150 0,152 2 98,97 0,020
11а 1000 6,57 1250 4 0,253 2 98,97 0,020
12 1000 12,5 1020 0,022 2 98,97 0,020
13 1000 5,31 1050 0,052 2 98,97 0,020
14 1000 3,83 1070 0,072 2 98,97 0,020
15 1000 2,71 1100 0,102 2 98,97 0,020
16 1000 1,82 1150 0,152 2 98,97 0,020
17 500 86,2 1020 0,022 2 98,97 0,020
18 500 36,6 1050 0,052 2 98,97 0,020
19 500 26,4 1070 0,072 2 98,97 0,020
20 500 18,7 1100 0,102 2 98,97 0,020
21 500 12,6 1150 0,152 2 98,97 0,020
22 1000 4,69 1070 0,072 10,86 98,97 0,110
23 1000 4,95 1070 0,072 21,28 98,97 0,215
24 1000 6,02 1070 0,072 40,43 98,97 0,409
25 1000 8,72 1070 0,072 60,36 98,97 0,610
26 1000 1070 0,072 2 53,98 0,037
27 5000 1070 0,072 2 211,31 0,009
1 кг = 0,5 м/сут; кг = 0,5 м/сут; 2 кг = 0,5 м/сут; к2 = 0,05 м/сут;
3 кг = 0, 05 м/сут; к2 =0, 5 м/сут; 4 кг = 0,1 м/сут; к2 =0, 1 м/сут.
на блоки, размеры которых увеличиваются от скважины к внешним границам модели по логарифмическому закону, является оптимальной (рис. 7, б). Для выбранных размеров модельной области сетка 30 (по оси г) х 90 (по оси г) позволяет добиться практически полного совпадения модельных результатов с аналитическим решением для схемы поршневого (без дисперсии, Др = 0) вытеснения (рис. 8); некоторое расхождение между решениями за счет численной дисперсии наблюдается только на начальном этапе процесса (когда С < 0,1). На рис. 8 (вариант Др > 0) видно, что более дробное представление сеточной области не улучшает результаты расчетов. Там же отмечается влияние плотности на отклонение модельных кривых от аналитического решения (20), полученного для случая нулевого плотностного градиента
(Др = 0).
Рис. 8. Сравнение численного (сетка 30 X 90) и аналитического (20) решений (случай Ар = 0 — кривые 1 и 2), а также решений для базового варианта (Ар > 0) при различной сеточной разбивке модельной области — кривые 3 (30 X 90), 4 (60 X 60), 5 (60 X 90), 6 (30 X 30).
На заключительном этапе тестирования модели проводилось сравнение распределений C(r, z, t), полученных на программных комплексах TOUGH2 и SUTRA, с поверхностью раздела z(r, t) —см. решение (22).
Как видно на рис. 9, иллюстрирующем результаты такого рода теста, «аналитическая» граница раздела практически точно описывает положение концентрационного фронта С = 0, 5, найденного численным методом. В то же время графики (см. рис. 9) показывают, что моделирование процесса сопровождается формированием переходной (дисперсионной) зоны, имеющей численную природу. Совпадение модельных концентрационных кривых C(t) (концентрация в скважине) с аналитическим решением (20) (см. рис. 8, кривые 1 и 2) говорит о малой значимости отмеченного численного эффекта с точки зрения достижения поставленной цели — нахождение аналитических зависимостей, аппроксимирующих выходную концентрационную функцию. Слабая чувствительность функции С(t) к численной дисперсии объясняется определяющей ролью деформации сетки движения подземных вод. Изменение концентрации происходит в основном за счет разновременного привноса вещества из различных зон по различным линиям тока, а не за счет дисперсионных эффектов.
Говоря о профильном распределении концентрации, необходимо отметить определенный гистерезис процесса (рис. 10), проявляющийся в быстром формировании «купола» на стадии подтягивания воды к скважине при водоотборе (Q > 0) — концентрационные изолинии достаточно быстро достигают положения, близкого к стационарному положению, и довольно медленном опускании «купола» в период остановки откачки из скважины (Q = 0) .
Наконец, анализ графиков (см. рис. 9 и 10) показывает, что фильтры скважины до-
стигают изоповерхности минерализации воды, характеризующейся довольно низкими значениями концентрации солей (C ^ 1). Такой эффект связан с высоким разбавлением соленых вод вблизи скважины потоком пресной воды, занимающей верхнюю зону разреза.
Дальнейший анализ миграцион- о-1
ного процесса ограничивается в основном рассмотрением поведения функции C(t). Рассматриваются модельные ситуации (варианты), в которых дебит водозабора, по крайней мере? в 2 раза превышал предельный дебит Q*(Q/Q* > 2). В противном случае из априорных представлений ясно, что результаты моделирования могут быть в значительной степени искажены численно-дисперсионными эффектами.
Кроме того, при выбранном соотношении Q/Q* (>2) влияние различий в плотностях жидкостей на время достижения скважины первых порций соленых вод (to — см. формулу (21а)) незначительно, что позволяет, как видно из дальнейшего, упростить выбор аппроксимирующих аналитических зависимостей.
Концентрация соли в скважине конечной длины (в > 0), представленной несколькими модельными блоками (длиной Azit), рассчитывалась как средневзвешенная: C =
Y^Ci Qi/J2Q, где Qi —дебит водоотбора из интервала Azit.
На границе модели (r = R) задавалось условие Дирихле: P(t) = pgh = const (h — высота столба жидкости). Для базового варианта влияние границы при R > (ms + b) ~
2000 м (рис. 11) сказывалось на характере функции C(t) лишь на весьма длительных отрезках времени (t = t/t0 > 102 — 103, что в абсо-
Рис. 9. Концентрационные распределения С (сплошные линии — численное решение задачи) и положение границы раздела вод (пунктирная линия — решение (22), в = 0, 02).
а — э = 0, 69, б — э = 1, 99; базовый вариант при Др (см. табл. 2).
О
Рис. 10. Динамика подтягивания границы раздела (С = 0,5) к скважине (а) и опускание (б) сформировавшегося «купола» соленой воды после остановки скважины. Базовый вариант. Числа у кривых—текущее время.
лютных величинах соответствовало десяткам-сотням лет); приближение плановой границы (г = И) к эксплуатационной скважине оказывает заметное выполаживающее влияние на функцию С(£), которое ощутимо на менее коротких временных интервалах. Подобное весьма специфическое поведение функции С(£) в ограниченных по простиранию пластах связывается, по-видимому, с влиянием границы г = И на перерас-
т0
Рис. 11. Влияние расстояния до границы модели (К —числа у концентрационных кривых) на изменение относительной концентрации соли в скважине; А — аналитическое решение (23) для неограниченного пласта.
пределение расходов жидкостей, поступающих в эксплуатационную скважину из двух профильных зон потоков 0 < г < Ь и Ь < г < ш8. Для описания такого рода эффекта в пласте, дренируемого совершенной скважиной, существует аналитическое решение [33].
2.2.2. Аппроксимация численных решений аналитическими функциями
В качестве «ядра» искомого решения предлагается рассматривать функцию (20), преобразованную к виду
, (26)
А = I (ш), (27)
где показатель степени А является функцией безразмерного дебита ш:
“ = Я ** = \кгЪ2Ар. (28)
Как и в случае с миграцией растворов постоянной плотности, значение 1о определяется трансцендентной формулой (21).
На первом этапе модельного анализа рассматривались варианты расчетов (В, 1-21) для изотропного и анизотропного пластов, водоотбор из которых проводился скважиной с коротким фильтром (в = 0,0202). С практической точки зрения — это точечный
источник, для которого выражение (26) трансформируется к виду
с = 1 - (т) ’ (29)
(к = л/кг/кг)- В этом случае удалось установить следующий вид функциональной зависимости:
А=1-— |1 + ^—* (30)
у 4у(1п у)~3/4-1 у ’ У 7
где у = .
Таким образом, показатель степени А(ш)в формулах (26) и (29) контролирует степень выпуклости функции С = С(р). Значение А =1 достигается при ш ^ то, т. е. когда дебит водоотбора Q многократно превышает характеристику Q** (28), зависящую от нескольких параметров (Ь,кг и Др). В этом случае формула (29) трансформируется в решение (23), полученное для случая Др = 0.
Во всем рассмотренном диапазоне Q/Q* > 2 решение (29) при А (ш), определяемом по (30), удовлетворительно (с погрешностью не более 1-2%) описывает всю серию модельных экспериментов, имитирующих миграцию рассолов в изотропных и анизотропных пластах (рис. 12 и 13). Заметные расхождения наблюдаются на больших моментах времени, когда на результатах моделирования начинает сказываться влияние плановой границы пласта (г = Д). Оно хорошо иллюстрируется кривыми на рис. 12, построенными для базового варианта.
Кроме того, из анализа модельных кривых следует, что различия в плотностях контактирующих растворов (вплоть до Др = 0,25, что отвечает минерализации соленой
Рис. 12. Сравнение численного (точки) и аналитического (сплошные линии, формула (29)) решений тестовой задачи. а — варианты 1—6 и В (числа 350—5000 — дебит скважины, м /сут), б — варианты 7—11 и В (числа 1000—1150 — плотность соленой воды, кг/м3. Штриховая линия — решение (29) при
Ар = 0).
Рис. 13. Сравнение численного (точки) и аналитического (сплошные линии, формула (29)) решений тестовой задачи. а — варианты 12—16 («прямая» анизотропия, к = б— варианты 17—21 («обратная»
анизотропия, к = 1/\/10)- Числа 1020—1150 — плотность соленой воды, кг/м3.
0.1 1 10 100 1000 10000 1 10 100 1000 10000
t/tQ */*о
Рис. 14- Сравнение численного (точки) и аналитического (сплошные линии, формулы (а — 26, б — 29)) решений тестовой задачи. а — варианты 22—25 и В—для различных длин фильтра (10—60, м), б — варианты 26—27 и В — при различных исходных мощностях пресной воды (10—211,31, м).
воды С более 300 г/л) практически не влияют на характеристику £о, полученную для случая Др = 0. В то же время видно, что дальнейшие темпы изменения функции С сильно зависят от градиента плотности Др. Так, расчетный вариант 7 (Др = 0,0021, М « 3 г/л) дал примерно 3%-ное занижение концентрации в области С « 0,5 по сравнению со случаем Др = 0, а вариант 8 (Др = 0, 022, М « 30 г/л — морская вода) — около 18%.
Дальнейший анализ показал, что и общее решение (26) оказывается справедливым, если только длина фильтра скважины не превышает половины мощности зоны пресных
вод (в < 0, 5) (рис. 14). Дополнительные расчеты, выполненные для различных b и представленные на рис. 14, подтвердили применимость аппроксимаций (26) и (29).
Литература
1. Ломакин Е.А., Мироненко В. А., Шестаков В. М. Численное моделирование геофильтрации. М., 1988.
2. Muskat M., Wyckoff R. D. An Approximate Theory of Water-Coning in Oil Production, Trans., AIME. 1935.
3. Маскет М. Физические основы технологии добычи нефти. М.; Ижевск, 2004.
4. Щелкачев В. Н., Лапук Б. Б. Подземная гидравлика. М.; Ленинград, 1949.
5. Bear J. Dynamics of fluids in porous media. Dover Publ., INC. New York, 1972.
6. Reilly T. E., Goodman A. S. Quantitative analysis of saltwater-freshwater relationships in groundwater systems — a historical perspective // Journ. of Hydrol. 1985. Vol. 80.
7. Essaid H. I. A multilayered sharp interface model of coupled freshwater and saltwater flow in coastal systems: Model development and application // Water Res. Res. 1990. Vol. 26, N7.
8. Motz L. H. Salt-water upconing in an aquifer overlain by a leaky confining bed // Ground Water. 1992. Vol. 30, N 2.
9. Bower J. W., Motz L. H., Durden D. W. Analytical solution for determining the critical condition of saltwater upconing in a leaky artesian aquifer // Journ. of Hydrol. 1999. Vol. 221.
10. Panday S., Huyakorn P. S., Robertson J. B., McGurk B. A density-dependent flow and transport analysis of the effects of groundwater development in a freshwater lens of limited areal extent: the Geneva area (Florida, U.S.A.) case study // Journal of Contam. Hydrol. 1993. Vol. 12.
11. Lin L.-C., Tsay T. K., Hsu N.-S. Saltwater upconing due to freshwater pumping // Proc. Nat.Sci. Counc. ROC(A). 1999. Vol.23, N2.
12. Zhou Q., Bear J., Bensabat J. Saltwater upconing and decay beneath a well pumping above an interface zone // Transport in Porous Media. 2005. Vol. 16, N 3.
13. Shalabey M. E. E., Kashyap D., Sharma A. Numerical model of saltwater transport toward a pumping well // Journ. of Hydrol. Eng. 2006. Vol. 11, N 4.
14. Bear J. Hydraulics of Groundwater. McGraw-Hill. N. Y., 1979.
15. Motz L. H. Discussion of «A density-dependent flow and transport analysis of the effects of groundwater development in a freshwater lens of limited areal extent: The Geneva area (Florida, U.S.A.) case study», by Pandy et al. (1993) // Contaminant Hydrology. 1995. Vol. 18.
16. Muskat M. The flow of homogeneous fluids through porous media, McGraw-Hill, N. Y., J. W. Edwards, Ann Arbor, 1946.
17. Bear J., Dagan G. Some exact solutions of interface problems by meansof the hodograph method // Journal of Geophysical Research. 1964. 69.
18. Bennett G. D., Mundorf M. J., Hussain S. A. Electric-analog studies of brine coning beneath fresh-water wells in the Punjab Region, West Pakistan // U.S. Geological Survey Water-Supply Paper 1608-J. 1968.
19. Schmorak S., Mercado A. Upconing of fresh water — see water interface below pumping wells, field study // Water Res. Res. 1969. Vol. 5. N 6.
20. Wirojanagud P., Charbeneau R. J. Saltwater upconing in anconfined aquifers // Journ. of
Hydr Eng. 1985. Vol.III, N3.
21. Чарный И. А. Основы подземной гидравлики, М., 1956.
22. Кисель В. А., Абрамов Д. С. Разработка нефтяных залежей с подошвенной водой. М., 1978.
23. Shalabey M. E. E., Kashyap D., Sharma A. Numerical model of saltwater transport toward a pumping well // Journ. of Hydrol. Eng. 2006. Vol. 11, N 4.
24. Huntush M. S. Hydraulics of wells // Advances in Hydroscience, 1 / V.T.Chow (ed.). Aca-
demic Press, New York, 1964.
25. Веригин Н. Н., Васильев С. В., Саркисян В. С., Шержуков Б. С. Гидродинамические и физико-химические свойства горных пород. М.: Недра, 1977.
26. Гольдберг В. М. Гидрогеологические прогнозы качества подземных вод на водозаборах. М., 1976.
27. Карпычев В. А. О перемещении водонефтяного контакта в пластах с подошвенной водой // Инженерный сборник. 1959. Т. 25.
28. Абрамов Ю. С., Кац Р. М. О пространственном движении границы раздела двух несжимаемых жидкостей в пористой среде // Изв. АН СССР. Механика жидкости и газа. 1967. №6.
29. Мироненко В. А., Румынин В. Г. Опытно-миграционные работы в водоносных пластах. М., 1986.
30. Pruess K. EOS7, An equation-of-state module for the TOUGH2 simulator for two-phase flow of saline water and air. Earth Science Division, Lawrence Berkeley Laboratory. Report N LBL-31114, Berkeley, 1991b.
31. Pruess K., Oldenburg C., Moridis G. TOUGH2 — User’s guide. Version 2.0, LBLN-43134, November. Berkeley, 1999.
32. Voss C. I., Provost A. M. SUTRA: A model for saturated-unsaturated, variable-density ground-water flow with solute or energy transport: U.S. Geological Survey Water-Resources Investigations Report 02-4231. Berkeley, 2008.
33. Nordbotten J. M., Celia M. A. Similarity solutions for fluid injection into confined aquifers // J. Fluid Mech., 2006.
Статья поступила в редакцию 25 января 2010 г.
