Учет капиллярных эффектов при проникновении жидкого углеводородного загрязнителя в увлажненную почву

УЧЕТ КАПИЛЛЯРНЫХ ЭФФЕКТОВ ПРИ ПРОНИКНОВЕНИИ ЖИДКОГО УГЛЕВОДОРОДНОГО ЗАГРЯЗНИТЕЛЯ В УВЛАЖНЕННУЮ ПОЧВУ

К.А. ПОТАШЕВ*, Н.Д. ЯКИМОВ**

*Казанский государственный университет **Казанский государственный энергетический университет

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

Введение

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

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

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

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

© К. А. Поташев, И.Д. Якимов

Проблемы энергетики, 2007, № 7-8

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

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

1. Постановка задачи

Будем рассматривать задачу о проникновении НВЖ в почву с ее поверхности (при заданном объемном расходе). До попадания загрязнителя почва находится в состоянии неполного равновесного насыщения водой. Требуется определить распределение НВЖ по глубине почвы в течение процесса загрязнения.

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

При моделировании будем исследовать слой почвы от дневной поверхности до глубины расположения зеркала грунтовых вод Ь . Этот слой в каждый момент времени можно разделить на две подобласти - первую (I), незагрязненную, и вторую (II), в которую проник загрязнитель. Первая область будет насыщена двумя фазами - водой и воздухом, а вторая будет находиться в условиях трехфазного насыщения водой, воздухом и НВЖ.

Подвижной границей, разделяющей верхнюю (II) и нижнюю (I) области, будет являться фронт проникновения загрязнителя, находящийся в момент времени * в точке г = | (*):

^ = 1, вП = 0, при г ,Ь],

= 1, > 0, при г е[0,1],

где - фазовые насыщенности (нижний индекс (г = , п, g) характеризует

соответствующую фазу, верхний индекс (N = I, II) - соответствующую область).

2. Описание капиллярных взаимодействий

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

wg

В качестве зависимости капиллярного давления вода-воздух Pc = Pg — Pw

(при заполнении порового пространства водой и воздухом) от водонасыщенности примем уравнение Ван-Генухтена [6]:

1

PcwgI (0w )= Pc (0w ) = "

a

0 w

1

M

\ — n

— i

где a, n, M - параметры, определяемые из условия наилучшей аппроксимации уравнением (1) экспериментальных данных; Рi - давление в i - фазе;

Si - Sir

0 i =-------- - эффективная насыщенность пористой среды i - фазой (Sir -

1 - Sir

величина остаточной насыщенности, при которой данная фаза в пористой среде теряет связность). Обычно предполагается, что M = 1 — 1/n [7]. Параметр a является обратным по отношению к величине давления проникновения несмачивающей жидкости Pd - значения капиллярного давления, при котором начинается дренаж пор (“bubbling pressure”) [6]. Часто предполагается, что параметр a можно оценить, зная величину поверхностного натяжения на границе вода-воздух у wg и контактный угол смачивания в wg, с помощью зависимости

1/a = о wg-Jm/k . Здесь m, k - пористость и проницаемость пористой среды, и введено обозначение о wg = у wg ^sP wg.

Параметр n связан с индексом распределения пор по размерам n = X +1 и является мерой в отношении размера каналов, связывающих поры [6]. Типичное значение n для природных грунтов лежит в диапазоне 2...5.

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

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

Важную роль во взаимном распределении НВЖ и воды в пористой среде играет явление растекания органической жидкости по поверхности воды [8]. Рассматривая наиболее распространенные типы жидких углеводородных загрязнителей, можно предполагать, что вся поверхность «вода-воздух» будет покрыта слоем НВЖ.

В загрязненной области (II) при описанном расположении фаз в порах необходимо задавать две функции капиллярного давления PCm = Pn — Pw и

Паркером и др. в работе [5], согласно которому предполагается, что функция Р™п

Здесь Pc (0) - функция капиллярного давления Ван-Генухтена,

определяемая уравнением (1). Разность давлений в воде и воздухе в загрязненной

Этот подход согласуется с результатами экспериментов в условиях 3-фазного насыщения пористых сред. Основным его достоинством является выражение капиллярных свойств сложных систем через капиллярные свойства канонической системы вода-воздух.

Для твердой матрицы почвы, сильно смачиваемой водой и НВЖ, в условиях термодинамического равновесия можно получить следующее соотношение [8]: ° ^ = ° ^ + ° ^. Тогда из трех величин поверхностных натяжений

независимыми остаются лишь две.

3. Описание транспорта воды и НВЖ

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

Будем считать, что сечение такого канала можно характеризовать эффективным радиусом гс по величине капиллярного давления Р^, необходимого для удержания мениска в данном месте (гс = о /Ре^ ), причем гс не зависит от природы жидкостей. Так что при определенном давлении Pc имеется значение гс (Рс ) такое, что сечения канала с радиусом г < гс заняты водой, ас г > гс - воздухом. Тогда такую геометрическую характеристику структуры пористой среды, как относительный объём пор меньшего гс размера

ng

Рс = Pg — Рп. Для их определения воспользуемся подходом, предложенным

зависит только от эффективной водонасыщенности 0 w, а функция Рп - только от эффективной насыщенности пористой среды газом 0 g :

рwn = р — р = о wn р

рс = рп рw = рс V0 w ),

о wg

wg

области можно определить из (2) и (3): Р^П = Pg — Pw = Рп + Р™.

уг = Уг

9

можно отождествить с насыщенностью воды Уг = Sw при двухфазном

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

сечения. То есть в местах каналов с г <гС находится вода, при г™ <г <гр -

углеводород, и при г < - воздух, причем «граничные» значения г™ , гр

определяются по насыщенностям воды и углеводорода в соответствии с (4).

При > 8СГ имеются каналы для транспорта воды и воздуха. Но для

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

гС < г < Гп. Ясно, что идеально ровные каналы в реальной пористой среде невозможны и не может найтись связных каналов со сколь угодно близкими между собой наименьшим и наибольшим сечениями (здесь г™ , гр). Поэтому и

модель среды должна отвечать следующему свойству. Если рассмотреть множество всевозможных каналов с наиболее узким сечением, равным некоторому г, то во всех них наиболее широкое сечение будет не меньше некоторого г. Каналов «более ровных», в которых по всей длине Г < г < г - 8, г > 0, не существует.

Следовательно, если насыщенности воды в данном месте среды

соответствует по

значение гр (так что Sw + Бп = Уг больше г .

Таким образом, для описания транспорта в системе «вода-углеводород-воздух», в отличие от случая «вода-воздух», требуется дополнительная характеристика структуры порового пространства, описывающая степень «неровности» каналов, в виде, например, связи

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

Вместо (5) можно использовать зависимость вида

насыщении и определить, например, по зависимости Sw = Sw \РС ) (см. (1)), вместе с указанной связью гс \РС ).

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

насыщенностью воды . Здесь = Уг (г), Бт1П = Уг (г )-.

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

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

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

Тогда давление Р% в водной фазе перед фронтом и за фронтом одинаково, а

насыщенность углеводородом соответствует минимально возможной Бт1П по (6) для его транспорта.

4. Основные уравнения

Уравнения баланса масс НВЖ и воды в двух областях имеют вид [9]:

дтБ% Т

-------+ йlvqw = 0,

д *

дшБ„ 1Т

----— + Шу?" = 0

д *

дшБЩ ТТ ’

п + ¿тП1 = о

д *

где Чщп - скорость фильтрации воды, НВЖ.

Считаем, что фильтрация флюидов подчиняется закону Дарси:

1^1 1

ч™ =-нw,

д ™

=-----н%

д

%

чП =-------/п нП

д п

где д%,п - вязкость воды, НВЖ; /щп - относительные фазовые проницаемости почвы для воды, НВЖ; Нщп - гидронапор в водной и НВЖ фазах.

Согласно уравнениям (2), (3) и представлению функции напора, получим следующие выражения для напоров воды и НВЖ в областях (I) и (II):

нЧ (г) = -Рс (в% )-Р,

Н% (г) = -Рс (1 - в8 )-—Рс (в% )-Р,

нп (г) = -Рс (1 - вg )- Рngz,

- wg

где р I - плотность I -фазы; g - гравитационное ускорение.

Функции относительных фазовых проницаемостей воды и НВЖ в загрязненной области выберем в виде, предложенном Ван-Генухтеном [10]:

1

\2

/ЧТ=в 2 (1 -(1 - в ч(п-1 ’Г1 у п

1 /(п-1 )Л2 (п-1 )/ п

/ЧТ = вЧ ( 1 -(1 - вп)п

В незагрязненной области (I) в связи с незначительным сопротивлением фильтрации воды со стороны воздуха будем использовать линейную зависимость фазовой проницаемости воды от ее насыщенности [11]:

/Ч = в ч .

5. Начальные и граничные условия

До попадания загрязнителя почва находится в состоянии равновесного насыщения водой в ч (г, * = 0) = в Ч (г). Эта зависимость определяется из условия гидростатического распределения давления в воде и условия на нижней границе -на зеркале грунтовых вод, являющейся свободной поверхностью (давление воды

здесь совпадает с атмосферным давлением Рg, которое примем нулевым):

Н% (Ь) = -Рс (1)-РwgL = -РwgL, Рс (в% )= Рwg(L - г), (при *=0).

Поровый объем почвы на глубине г = Ь всегда полностью занят водой.

Величина объемного расхода загрязнителя на поверхности почвы задана:

чЩ (г = 0,*) = ?п.

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

НЧ

Уравнение (6), связывающее насыщенности воды и углеводорода на фронте, было задано для примера в форме

0П

2

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

1^/

где [Яч ] = (^)-ЯЧ И) и [Чч ] = чЧ(^)-чЧ1 ).

6. Численный алгоритм решения задачи

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

1. Расчет идет по насыщенностям воды и НВЖ в узлах: [г ], Яо [г ] - на

нижнем, 1[г ], Яо 1[г ] - на верхнем слое по времени. Кроме того, требуется

подсчитать приращение насыщенностей , йЯо в каждом узле. Если узел г и соседние с ним не являются узлом фронта или прилегающими к нему узлами, то приращения рассчитываются как =(чш- - (г-1 ))й*/йг,

йЯ о =(ч о - ч о (г-1 ))ж/йг, где потоки , ч ог рассчитываются по

насыщенностям в соответствующих узлах с расчетом по ним давлений:

/ Р™г - Рч (г-1) 4

-----+ Р

(— ) Рш Рп(і-1) (— ) Роі Ро(і-1) ^

Чм>і = к (ЯП ) + р Ч оі _ к ( о ./I + р о 3

dz

dz

где = (і + (і-1 ))/2, Я о =( Я оі + я о (¿-і))/ 2 - средние между узлами

насыщенности. Формулы, используемые для подсчета давления, зависят от того, в какой зоне идет расчет: на фронте, до или после него. Фронт находится в элементе с номером іі/. Его узловые значения используются для хранения значений насыщенностей на верхней стороне фронта. Насыщенность на нижней стороне фронта хранится в отдельной переменной.

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

2.1. Точная координата фронта учитывается в виде: номер элемента сетки іі/, в котором находится фронт, плюс относительное положение фронта в

элементе \ (обычно 0 < \ < 1 ), так что = dz(іі/ + \ - 0,5).

2.2. Соотношение для насыщенностей на фронте предполагает его движение таким, что насыщенность Я о на фронте - предельная (при имеющейся

насыщенности воды Я%), так что Я% = Я% (я% ) (см. (7)) и Я0 = Я0 (я% ) (см. (8)).

2.3. При определении насыщенностей воды и НВЖ в элементе й/ считаем, что НВЖ с насыщенностью Я о располагается выше линии фронта, а вода - выше

фронта с насыщенностью Я% , ниже - Я%.

2.4. Далее требуется найти перемещение фронта и значения

насыщенностей после шага < по времени. Рассматриваются только случаи с движением фронта вниз. При этом в элемент поступает НВЖ чо > 0 и из элемента

должна уходить вода 4% - < 0.

Это вызывает изменение насыщенностей и движение фронта. Насыщенности связаны (7), (8) и как бы сами перераспределяются в элементе в соответствии с этими соотношениями. При этом будет наблюдаться баланс по загрязнителю Я0 (| )<7 - Я0 %<7 = Ч0 & и по воде

Я% (+)<11+Я% (1 -1 - )<11 - Я % - Я% (1 -1 )<7 = (ч% - <.

Здесь символом V обозначены величины на начало, Л - на конец шага по времени. Таким образом, для сдвига фронта имеем соотношения:

Ч 0 & Я 0 - Я 0

= -°-----1-0------ по НВЖ,

Я 0 <7 Я 0

=

по воде.

Следовательно, можно получить условие совпадения (тождественности) движения фронта по НВЖ и воде (аналог соотношения (9)):

<7

— -1

' <1 = Ч 0

V <7 /

Таким образом, получаем уравнения (8), (10) для определения и Я% .

2.5. Далее следует расчет насыщенностей на фронте.

2.6. При расчете значений потоков от элементов, соседних с фронтовым, градиенты давления считаем по расстояниям от узлов до фронта (например, («/ -1)- 7/), а не до узла в центре фронтового элемента. Это справедливо и для

баланса на фронте при решении (11), (10) при расчете самих соседних элементов.

3. Остается предусмотреть переход фронта через границу фронтового элемента, когда должно быть изменено значение и/. Когда | станет больше единицы, нужно перевести фронт в следующий элемент, то есть сделать элемент с © Проблемы энергетики, 2007, № 7-8

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

7. Результаты расчетов

Результаты численных расчетов представлены на (рис.1-4). В качестве основного набора параметров при расчетах использовались следующие значения: m=0,5; n=3; ¿=1,5 м; k=10-11 м2; д w = 0,001 Па • с; дn = 0,002 Па • с;

р w = 1000 кг / м 3; р n = 800 кг / м 3; чП = 1 • ¿¿j, где T = (l2 д w )До wg*Jmk ).

Сравнение распределения насыщенностей на (рис. 1) и (рис. 2) показывает, что величина скачка водонасыщенности при переходе через фронт загрязнения растет с понижением значения о wn при фиксированном о wg.

Сопоставление рис. 1 и рис. 3 демонстрирует влияние гравитационных сил на скорость проникновения загрязнителя в почву. На рис. 3 приведены результаты расчета процесса загрязнения почвы тяжелыми НВЖ (ТНВЖ). При фиксированных значениях поверхностных натяжений заданная при расчетах плотность модельного загрязнителя превышает плотность легких НВЖ (ЛНВЖ) из первого расчета (рис. 1) более чем в два раза. В результате скорость фронта ТНВЖ также превышает скорость фронта ЛНВЖ, однако в меньшей степени, чем соотношение их плотностей.

Рис. 1. Распределение фаз в различные моменты времени_______ю2 = 0,5, 1,1,5, 2, 2,5, 3

при о wg = 0,07, о ng = 0,025, о wn = 0,045 Па • м

ч T ' wn

(

Рис. 2. Распределение фаз в различные моменты времени 102 = 0,5, 1,1,5, 2, 2,5, 3

при о ^ = 0,07, о ^ = 0,035, о ^п = 0,035 Па • м

к Т

’ wn

1 2

— 10^ = 0,5, 1,1,5, 2, 2,5,3

IТ

3

Рис. 3. Распределение фаз в различные моменты времени

при о ^ = 0,07, о ^ = 0,025, о п,„ = 0,045 Па • м, р п = 2000 кг / м

Рис. 4. Распределение фаз в различные моменты времени

IT

102 = 0,5, 1,1,5,2,2,5,3

при о ^ = 0,07, о ^ = 0,025, о ^ = 0,045 Па • м, L = 3 м

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

Выводы

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

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

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

t

Работа выполнена при поддержке гранта МНТЦ #2419 Summary

The one-dimensional mathematical model of non-aqueous phase liquids (NAPL) vertical infiltration in the humidified ground is constructed. As an example of NAPL hydrocarbons as the most widespread contaminant are considered. The front of NAPL © Проблемы энергетики, 2007, № 7-8

penetration is mobile border between areas of two and three-phase saturation. Due to presence of NAPL in the polluted area changes water retention ability of ground that results to gallop of a water-saturation at the front pollution. The numerical algorithm of the decision of a problem is developed.

Литература

I. Современные проблемы гидромеханики и гидрогеологии. Сборник докладов конференции. - СПб. 2002. - 572 с.

2 Pinder G.F., Abriola L.M. On the simulation of nonaqueous phase organic compounds in the subsurface // Water Res. Res. 1986. - V. 22. - No. 9. - P. 109-111.

3. Костерин А.В., Бреус И.П., Неклюдов С.А., Бреус В.А. Защита грунтовых вод от углеводородного загрязнения путем усиления барьерных функций почвы // Фундаментальные проблемы воды и водных ресурсов на рубеже 3-го тысячелетия: Материалы Всеросс. науч. конф.. - Томск: Изд. НТЛ. - 2000. -С. 240-243.

4. Поташев К.А., Якимов Н.Д., Бреус И.П. Транспортные процессы в

пористых средах: капиллярные эффекты (обзор современной зарубежной

литературы) // Вопросы атомной науки и техники, Математическое

моделирование физических процессов. - 2005. - Вып. 3. - С. 72-79.

5. Parker J.C, Lenhard R.J., Kuppusami T. A parametric model for constitutive properties governing multiphase flow in porous media // Wat. Resour. Res. 1987. -No. 23(4).- P. 618-624.

6. Van Genuchten M.Th. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils // Soil Sci. Soc. Am. J. 1980. - No. 44. - P. 892-898.

7. Chen J., Hopmans J.W., Grismer M.E. Parameter estimation of two-fluid capillary pressure-saturation and permeability functions // Adv. in Water Resour. 1999.

- No. 22(5). - P. 479-493.

8. Zhou D., Blunt M. Effect of spreading coefficient on the distribution of light non-aqueous phase liquid in the subsurface // J. Contam. Hydrol. 1997. - No. 25. - P. 1-19.

9. Баренблатт Г.И., Ентов В.М., Рыжик В.М. Движение жидкостей и газов в природных пластах. - М.: Недра, 1984. - 211 с.

10. Helmig R. Multiphase Flow and Transport Processes in the Subsurface - A Contribution to the Modeling of Hydrosystems. Springer-Verlag. - 1997.

II. Басниев К.С., Кочина И.Н., Максимов В.М. Подземная гидромеханика: Учебник для вузов. - М.: Недра, 1993. - 416 с.

Поступила 29.03.2007