К вопросу об изменении качества подземных вод на водозаборах Текст научной статьи по специальности «Математика»

ВЕСТНИК САНКТ-ПЕТЕРБУРГСКОГО УНИВЕРСИТЕТА

Сер. 7 2007 Вып. 2

УДК 556.3

В.Г. Румынии, Р.А. Филин

К ВОПРОСУ ОБ ИЗМЕНЕНИИ КАЧЕСТВА ПОДЗЕМНЫХ ВОД НА ВОДОЗАБОРАХ

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

Рис. 1. Откачка из пласта, получающего дополнительное питание посредством: а - перетекания через разделяющий слой из смежного горизонта, б - отжатая из подстилающей толщи. Условные обозначения даны в тексте.

1. Уравнения массового и фильтрационного баланса

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

ЗС 1 о

тЛ + Щ -—г [V! (г)С,)] - С20 = 0, (1)

от г дг

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

© В. Г. Румынии, Р. А. Филин, 2006

горизонте и разделяющих слоях; т, и п1 -мощность и пористость основного пласта; ч^г)- горизонтальная скорость фильтрации в эксплуатируемом пласте (для радиального потока в направлении, противоположном направлению оси г, V, < 0), у0 - вертикальная скорость фильтрации в разделяющих слоях (У0 > 0).

Выполнив дифференцирование по пространственной координате, преобразуем уравнение (1) к виду:

тхпх

дСх dt

+ TWjV,

д£_ dr

1

дг

Q-v0C20= 0.

(2)

На сравнительно длительных этапах развития фильтрационного возмущения в зоне квазистационарного режима вблизи скважины, где понижения напоров £ = Н-к. (Я - начальный напор, - текущие напоры, в частности, в эксплуатируемом горизонте, 5, = Н- - см. рис. 1, а) слабо зависят от времени, уравнение фильтрационного баланса для радиального потока имеет вид:

Г,

г dr

dSx dr

~v„=0,

где 7j = кхтх. Отсюда, учитывая, что Vj = -kxdhx / dr = kxdSx / dr, имеем

mxd r dr

1

тл -v, +

dv j

It

= v0.

Тогда базовое уравнение массового баланса (2) преобразуется к виду:

^L + M^L + 8(C1-C20) = 0,

dt

дг

или

ec+u{r)ec(r)C=Q>

dt дг

(4)

(5)

(6)

где С = Cj - С20, и = u{r) = Vj /п, е = e(r) = v0/mxnx. В общем случае, учитывающем сработку емкостных запасов пласта,

♦ ад

1 д ( dSA

dt

г dr\ dr

h - упругая водоотдача основного горизонта. Тогда:

т, д , ч . 8S,

г дг

и уравнение (2) преобразуется к виду

dt

dr

mjix dt

(7)

(8)

(9)

2. Основные фильтрационные решения

Для трехслойной системы с перетеканием (рис. 1, а) решение фильтрационной задачи, справедливое для длительных этапов откачки, относительно функций понижения 5, и£2, можно представить в виде [1,3,4, 5]:

5,=

2

4ТС(Г1+Г2)

е

1п

2.25 аг 2

4п(Т1+Тг)

2.25а г

где <2-дебит водоотбора; Тх = кхтх, Т2 =к2тп2 - проводимости слоев (горизонтов); кхик2~ коэффициенты фильтрации в эксплуатационном и смежном горизонтах, соответственно; а = (ТХ +7^)/(ц* + ц2) ~~ эффективная (суммарная) пьезопроводность;

В

™0Т1Т2 К {Тх+Т2у

V =

Г,

(10 а)

(В - параметр перетекания); к0 и т0 - коэффициент фильтрации и мощность разделяющего слоя; К0 (г) - модифицированная функция Бесселя второго рода нулевого порядка.

Вертикальная скорость фильтрации у0 находится из очевидного соотношения:

т0

2пВ{

2

С")

В

где

Подставим У0 из уравнения (11) в уравнение неразрывности фильтрационного потока (4), переписав последнее в виде:

/

(12)

Интегрируя левую и правую части уравнения (12), получаем:

ЯгВ

ГУу = -

2п тхВх

■к,

+ с,

(13)

где К, {г) - модифицированная функции Бесселя второго рода первого порядка.

Константу интегрирования с = —

2

2п тх 1

находим из условия:

й

2п тхг

Г^Л

КВ;

В Г л = — при--> и.

г У В

(14)

Окончательно имеем:

V, =-

0

2п тхг{\+у)

V +—К,

В 1

ГуЛ

(15)

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

Q

4тс Тх

4 « \

1--£ ~ехР

Е1

( 2 У\

Г

4 аЧ

п и=

и=1,3,5... '

Г ПК ^

ц0т0

БШ

Г \

пп г \2то;

(16)

оо -х

•е

с б

где — Е1(— г) = I-йх - интегральная экспоненциальная функция;

1 х

а* =ТХ /(ц* + ц,*) - эффективная (суммарная) пьезопроводность, ц* и Цд - упругая водоотдача хорошо- и слабопроницаемого слоев. Скорость у0 находится из закона Дарси:

= Уо(2 = °) = =0

&

[50 = {) = Н — к0(г, {) ]. В результате имеем:

2 к

тп

Ях £ ехр

я=1,3,5...

Г \г 1 ' пк \ К

ц0т0

(17)

(18)

Ряд в формуле (18) является быстросходящимся, поэтому при больших t, когда справедливо выражение (16) для понижения ^ в хорошо проницаемом пласте, можно записать:

2*о й

т0 4п Тх

г

-Ег

2 \\

4 Л

или

й

ца0

4тс а

' ( ~Е1

ехр

2 У\

2 К

Л

V .4

,-а0'

(19)

(19 а)

где \Х =■

Ио

2к0

■, а0 =—^

Нетрудно заметить, что функция V0(г, *) является

Ц0 + ИГ " №*<>

немонотонной: в любой расчетной точке г у нее имеется временной максимум,

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

Теперь подставим выражения, найденные для функций ^ = 5,1(г,*) и у0 =у0(г,?), в уравнение фильтрационного баланса (8). После несложных преобразований приходим к дифференциальному тождеству:

¿О1) =

871 тхТ^

-^/4«*/ » 2

с1г +■

ка

4тс т0тхТх

К

2

Г

е-2Ло/цЯ ¿гг (20)

//

Интегрируя левую и правую части этого тождества в пределах от г до 0, в предположении, что V, (г —> 0, = —Q / 2к тг, окончательно получаем:

г) = --2-{1 -(1 -¡1X1 -)-ца0t[1-4 Н(Ч)-е-*]}, (21)

2п тхг

(22)

<2

4 а г /

пт1пхг2 — причем у ----, б =

<2 4тг тхпха

Отсюда, в частности, видно, что | ух(г, ?) |^ <2 /2к тхг, т. е. процессы упругой отдачи воды основного горизонта и ее отжатия из водоупорного слоя приводят к уменьшению радиальной скорости потока в сравнении с условиями жесткого режима фильтрации. Такого рода эффект контролируется тремя параметрами - рТ, () и а0.

3. Решение миграционных задач 3.1. Схема с перетеканием

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

Л </г 1 с/С

— --=---; (23)

1 и{г) е(г) С

здесь, с учетом выражений (11) и (12),

и{г) = -

а

1

2п тхпхг 1

е

8 (Г) =

2% тхпхВх

г К0

В г

СВу

В

(23 а) (23 б)

Интегрируя первое равенство (23),

\ <Лг

\dt= . J J,

о

fu{r)

получаем выражение для нахождения времени движения частицы воды от произвольной точки на цилиндрической поверхности Г = г{{) до расчетной точки г = :

t =

КО

dr

I К)!

(25)

Интегральное тождество для определения текущей концентрации С в расчетной точке г = г^ выглядит, очевидно, следующим образом [см. второе уравнение системы (23)]:

:10 Q ^20 r(t) м(г)

(26)

Такая запись предполагает, что на цилиндрической поверхности г — г{{) , определяемой в любой расчетный момент времени I из выражения (25), выполняется условие

С = С M Ио-

Интегрируя левую часть тождества (26), находим окончательное решение задачи

в виде:

— С-С

С = Чо =1-ехр

С -Г

7 8(Г)

i l«W|

dr

(27)

или, в явной безразмерной форме,

С = 1 - ехр

f "(О

V rw

т ыпг *

J v +FKj(F)

(28)

где г = г IВ, ги, = гм11 В. При расчете изменения концентрации в одиночной скважине

можно положить гш = 0.

Графическое представление решения (28) (случай г№ = 0) приводится на рис. 2, где безразмерное время определяется соотношением 1 = 0} I Ъ1тхпхВх (Вх / к0 )•

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

^тах

1

1+V

(29)

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

T'l I IHII| n Hllll|-ГТТТПП]-1 I I IIIIIJ-1 I I lllllj

80 100 0.001 0.01 0.1 1 10 100 t t Рис. 2. Графики изменения относительной концентрации природного индикатора в откачиваемой воде (схема с перетеканием) в линейном (а) и логарифмическом (б) масштабах времени; числа на кривых - соотношение проводимостей (параметр v =ТХ/Т2).

3.2. Схема с отжатием

Перепишем исходное дифференциальное уравнение (9) в форме:

1 дС _, N дС —, ч_ ,, —v 1 -+ «(ПД)—+ е(п,т)С + (1-ц)-С20 = 0, Q от дх\ т

где и(п,х ) = -{l- (1 - ц)(1 - ) - цх е* [l- \ Ei(4) - в"* ]}, _ я aüm\n\r%Q

(30)

(30 а) (30 6)

Л

С = С1-С20, л = " * , Х=а0/, (30в)

Характеристики Т| и х ассоциируются с безразмерной пространственной и временной координатами. У

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

Так, в большинстве случаев, интересных для практики, выполняется условие цЦ » ц* (¡Т »1) (упругая водоотдача слабопроницаемых слоев, представленных обычно глинистыми породами, многократно превышает соответствующий показатель хорошо проницаемых горизонтов, из которых и осуществляется водоотбор). При ц = 1 уравнение (30) преобразуется к виду:

1 дС ч дС ч _ . ■=— + м(г|,т)—+ е(л,т)С = 0, Q от <тп

и(тьт) = -{1-х <Г [14 Ш(Ч)-е"*]} .

Графики функций и(п,х ) и £ (т],т) приведены на рис. 3. а

(31а) (316)

0.01 г- 8

|ш 1з

Рис. 3. Графики функций: (а) - и (Т), 8 (I) при фиксированных Т| и (б) - и (Т|), 8 (Т|) при фиксированных Т. Ы (Т), Ы (Г|) - сплошные линии; 8 (X), 8 (Т)) - пунктирные линии.

Линейному дифференциальному уравнению первого порядка в частных производных (31) соответствует система обыкновенных дифференциальных уравнений:

(х-с)2 ~

Ъх

=/со.

(32)

Решение этой системы с помощью стандартных численных процедур (пакет прикладных математических_программ МАРЬЕ 7) позволило построить графики функции т ад / Я* (рис. 4) и С(х ) = (С, - С10)/(С20 - С10) (рис. 5).

а б

100 -ц

Рис. 4. Графики функции Т|(Х) (а) и Т|(Х) I Qx (6), построенные при различных значениях безразмерного комплекса 0,.

0.2-,

0.16

0.12

Ю

0.08

0.04

0.1

0.05

Точное решение Решение (36) Решение (39)

Моделирование (МТЗОМЭ) 0.01

0 = 0.001

10

Рис. 5. Графики функции С(Т) построенные при различных значениях безразмерного комплекса Q.

Немонотонный характер графиков Г|(т ) / (см. рис. 4) объясняется экстремальным поведением функции скоростей г7(г|,т ) и 8"(г),т) (см. рис. 3, а).

Полученные численные решения для динамического отношения Т| (т ) /т можно аппроксимировать следующей приближенной формулой (см. рис. 4):

Н'-т

ехр

(т-с)2 Ъх

=/ео,

(33)

где коэффициенты а, Ь и с определяются так: а = 3.71£?0'92 - 3.38(2 (для - 1)> Ъ = 28.2, с = 4.2. При такой аппроксимации характеристика Г] является аналитической функцией безразмерного времени:

Л = - ал/Гехр

(т-с)2 Ьх

(33 а)

Второе уравнение системы (32), переписанное в явном виде, можно представить следующим образом:

Г11 I — ПИ— 11 #'1 >1

■<1г\, (34)

[-Ш(-г|/т)]е-

С _ {1 _х [1—СпД)ЕК-пД) - е^ ]}

или, если положить, ¿Л] / ¿/т = т[ = Г|'( т) [производная безразмерной координаты г\ (33 а)],

¿£ =__[-ш(-/(х))]л^-

с

{1 -т [1 - /(т) Е1(-/(Х)) - ]}

с1х ,

(35)

где функция /(т ) определяется зависимостью (33). Интегрирование последнего уравнения позволяет-йолучить решение задачи в интегральной форме:

[ }1 - 9 е-9 [1 - т И(- /(9))- е-/т]

N

(36)

С(т) = (С] (^ - расчетный момент времени). Расчеты по фор-

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

Далее, анализ графишвнарис. 4, б показывает, что в практически важных диапазонах изменения параметра £) (£? < 0.1) безразмерное отношение:

л п т, п,г2 не сильно отличается от 1. Это значит, что можно положить

^ и Г| =£т ; (37 а)

последнее тождество предполагает, что отжатие воды из слабопроницаемого слоя не приводит к отклонению расчетной скорости фильтрации от значений, определяемых по формуле для изолированного пласта: у1 = Q / 2кт1г.

Поэтому другое приближенное решение задачи имеет более простой вид:

С = 1 -ехр| - [-

(38)

_ Если учесть, что знаменатель дроби в подынтегральном выражении (38) при

Q <0.1 мало отличается от 1, то:

С = 1 - ехр [-£ (-Е1 (-£ )) (1- е-)]. (39)

Максимальное значение шкищтращщ. щж х = а0( > 2 -V 4 \

/ С=СЯ„« Ш <0.05 + 0.1). (40)

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

Ъ ЖЛШ Ж, исходить то реальных значении параметров горизонтов и разделяющих их глинистых толщ, воздействие на которые характеризуется низкими значениями безразмерного комплекса () = (214ктхп1а* ~ 0^14тшхп(Гх , то ожидаемые

максимальные относительные концентрации компонентов в откачиваемой воде, первоначально содержащихся в поровом растворе глин, в реальных условиях вряд ли могут превышать 5-10 %. В какой-то степени, данный вывод представляется парадоксальным, поскольку известно [2, 3], что основной объем воды, поступающей в скважины, определяется ее упругими запасами в глинистых толщах. Объяснение этому кажущемуся противоречию можно найти, если вспомнить, что водные запасы пласта, срабатываемые при водопонижении, формируются на периферии депрессионной воронки, причем в силу немонотонного («импульсообразного» - во времени) характера функции v0(r, /) в проводящий пласт в этой зоне из глин поступает ограниченное количество солей. Дальнейшая миграция растворов происходит в области, где запасы вод глинистой толщи уже сработаны, и поэтому рост концентрации солей в транзитном потоке не происходит. Поскольку зона поступления поровых растворов глин в эксплуатируемый горизонт постоянно отодвигается от скважины, то время подтягивания растворов из этой зоны в радиальном потоке все время увеличивает, что также способствует относительно быстрому выполаживанию концентрационной функции. Соответствующее время достижения предельных значений

а#и) > 2 + 4 или i(Cmax) «(1 + 2)ИоЧАо-

Summary

Rumynirt KG., Filin R.A. On ground water quality changes at the wellfield.

Models of subsurface flow and transport describing ground water quality at the wellfields, which withdraw water from stratified systems are developed. These systems initially are characterized by hydrogeochemical heterogeneity. The mathematical set-up of the problem assumes hydrodynamic interactions between the layers which are facilitated by: (a) leakage from the source aquifer towards the discharged aquifer (layer), (b) compression of the semipermeable layer. The obtained solutions describe both steady state and transient stages of the flow and transport processes. These solutions can be used for long-term ground water contamination by environmental component forecasting.

Литература

1. Бочевер Ф. M., Лапшин Н. Н„ Орадовская А. Е. Защита подземных вод от загрязнения. М., 1979. 2. Мироненко В. А., Шестаков В. М. Теория и методы интерпретации опытно-фильтрационных работ. М., 1978. 3. Проектирование, водозаборов подземных вод / Под ред. Ф. М. Бочевера. М., 1976, 4. Синдаловский JI.H. Справочник аналитических решений для интерпретации опытно-фильтрационных опробований. СПб., 2006.5. Huntush М. S. Hydraulics of wells // Chow V. Т. (Ed.). Advances in Hydroscience. 1964. 1. 6. Streltsova T.D. Well testing in heterogeneous formations. New York, 1988.

Данное исследование поддержано грантами CRDF (No.RUG2 2821-МО-Об) и РФФИ (№07-05-

007-96-а) и МНТЦ (№193).

Статья принята к печати 26.12.2006 г.