CC BY
129
УДЕ 624.131
. РАЗРАБОТКА ПРОГРАММНОГО ОБЕСПЕЧЕНИЯ ДЛЯ РЕШЕНИЯ
ГИДРОГЕОЛОГИЧЕСКИХ ЗАДАЧ
Букаты М.Б.
Рассмотрены тенденции и современное состояние в развитии гидрогеологического программного обеспечения» некоторые методы автоматизированной обработки геогидродинамической и гидрогеохимической информации. Сформулированы рекомендации по разработке моделирующих программных комплексов и АРМ, Описаны теоретические положения, методы и алгоритмы, реализованные в авторском программном комплексе 1ТО32. Обоснованы принципы построения сеточных моделей геомиграции, учитывающих гидрогеохимические процессы.
Последние десятилетия прошедшего века ознаменованы началом массовой компьютеризации фундаментальных и прикладных наук, связанной главным образом с развитием принципиально новых методов хранения и обработки информации, а также способов доступа к ней. Уже сейчас достижения в этой области, применительно к геологии вообще и гидрогеологии в частности, весьма впечатляющи: локальные, региональные и глобальные базы и банки данных, поисковые и геоинформационные системы, почти не ограниченная сложностью математическая, в том числе статистическая, обработка информации, 8 В и динамическая её визуализация, АРМ и моделирующие программные средства [1, 2]. В то же время, состояние наработок в этой сфере всё ещё находится больше в зародышевом, чем пригодном для полноценного применения, виде. Перевод накапливавшейся десятилетиями гидрогеологической информации на машинные носители, её рациональное с точки зрения доступа представление и организация стандартизованными средствами и в форматах СУБД и ГИС, использование локальных и глобальных сетей, современное состояние теории и средств численного моделирования, интерфейс имеющихся крайне разнородных программных продуктов - эти и многие другие вопросы сегодня всё ещё далеки от разрешения»
Тем не менее, с постепенным преодолением множества возникающих на этом пути проблем, информационные технологии постепенно становятся тем фундаментом, без которого полноценная профессиональная деятельность гидрогеолога становится уже невозможной. Важно подчеркнуть, что их разработка и применение не просто по-другому организуют обработку и представление информации, упрощая и облегчая рутинную повседневную работу специалиста (впрочем, "упрощение и облегчение" при этом часто сомнительно), но переводят её на качественно иной, более высокий уровень, позволяя получить новые знания на основе применения всё более трудоемких расчетов и численного моделирования, которые без соответствующей компьютерной базы часто в принципе неосуществимы.
Характерной чертой компьютеризации гидрогеологических исследований, на фоне уже достаточно большого числа имеющихся программных разработок, стало сохранение дефицита многофункциональных программных средств, пригодных для массового практического применения. В этой связи предметом настоящей работы стал краткий обзор применяемых в гидрогеологии прикладных компьютерных программ и освещение многолетнего опыта 348
авторской разработки программного комплекса HydroGeo (ПК HG) на кафедре гидрогеологии и инженерной геологии ТХГУ, который может рассматриваться в качестве одной из версий АРМ, как традиционного специалиста-гидрогеолога, так и гидрогеолога-нефтянника.
ПРОГРАММНЫЕ СРЕДСТВА ГИДРОГЕОЛОГИИ
Если в отношении программных средств общего назначения Россия далеко отстаёт от развитых стран и вряд ли догонит их в близком будущем, то в профессионально ориентированном программном обеспечении (ПО), во всяком случае гидрогеологическом, отечественные программы и программные комплексы, как правило, не только не уступают зарубежным аналогам, но часто и превосходят их по своим возможностям. Достаточно обширные обзоры компьютерных программных средств, в той или иной мере применяемых в гидрогеологии, опубликованы в работах [3-6]. Это дает возможность остановиться лишь на тех из них, которые могут быть отнесены к наиболее известным или используемым в России. Все они условно подразделяются на изначально гидродинамические и геохимические, которые, сближаясь, формируют класс наиболее сложных программ моделирования массо- и теплопереноса, обладающих возможностями более или менее полноценного учета процессов геохимического взаимодействия в системе раствор - порода.
Среди отечественных преимущественно гидродинамических пакетов программ наибольшее распространение получили программные разработки ЗАО Теолинк Консалтинг", ВСЕГИНГЕО, МГУ (г. Москва) и С ПбГИ(У) - ВИМС (г. Санкт-Петербург).
Основными продуктами компании 'Теолинк"® являются "Географическая информационная система GeoLink 2.13", "Программная система гидрогеологического моделирования" (включает программы: "Моделирование геофильтрации GWFS" и DOS-версии программ "Моделирование массопереноса MTS", "Моделирование фильтрации водонефтяной смеси": OWFS, TFDD, TFDD-rt и NAPL, "Расчет систем гидрогеологических скважин WSC", "Моделирование подземного и поверхностного стока на водосборе BASIN"), систему программ-модулей "Мониторинг геологической среды" (Мониторинг 6-11, включая VB Teis 12.00 P.C.ХПтенгелова - МГУ), "Мобильный программный комплекс AquaPalm" (для Palm-совместимых карманных компьютеров). Перечисленные программные продукты, охватывающие, пожалуй, наиболее широкий круг производственных и научных задач в области традиционной гидрогеологии, построены по модульному принципу и могут быть объединены в произвольно компонуемые пакеты программ, известные ранее как АРМ ГЕОо Наиболее сильными сторонами продукции, поставляемой "Геолинк", являются, насколько известно автору, блоки интерпретации результатов режимных наблюдений, обработки данных опытно-фильтрационных и опытно-миграционных работ, а также сеточного моделирования геофильтрации.
Во ВСЕГИНГЕО, наряду с рядом других гидрогеологических программных продуктов, развиваются 2 системы сеточного моделирования задач геофильтрации и массопереноса (пока практически без учета взаимодействий в системе вода-порода.) в среде подземных вод, близкие по назначению и своим возможностям: TOPAS А.Плетнёва и "Система специального программного обеспечения - автоматизированных сеточных шоделей гидрогеологических
• www.geolink-ltd.com; руководители проектов Л.В.Боревский, И. С. Пашковекий, A.A.Рошаль и др.
объектов" - ССПО Модель Е.А. Лолшкова. Данные системы, являющиеся примером наиболее перспективных отечественных гидродинамических программ гидрогеологического назначения, с большей или меньшей полнотой, охватывают основные используемые на практике расчетные схемы фильтрации вод, позволяя учитывать их неоднородность по плотности и вязкости, 3-х мерную изменчивость ^ фильтрационно-ёмкостных свойств многопластовых и/или .многослойных водоносных систем, напорный и безнапорный режимы фильтрации, разнообразные стационарные и динамически меняющиеся внешние и внутренние граничные условия, переменное разрешение (разномасштабность) и различную топологическую дискретизацию расчетной сетки моделей, решение краевых задач геомиграции инертных или изменяющихся по явно заданному закону компонентов состава вод.
■К этому же типу программного; обеспечения относится система Stimul Е.А.Ломакина и др., -разрабатывавшаяся в Санкт-Петербурге, сначала в ВИМС - С-ПбГИ(У), а затем в фирме "Водные ресурсы" [7], и, хотя и мало пока известная, но весьма перспективная, система моделирования,геофильтрации и массо- и теплопереноса В.И.Гунина и А.М.Плюснина (ЧИПР, г. Чита) [8]. Ряд интересных гидродинамических программных продуктов разработан в МГУ: TEIS и REGIM (обработка данных откачек и режимных наблюдений; Р.С.Штенгелов, М.И.Казаков), MCG (сеточная модель геофильтрации; С.О.Гриневский), TRANSFER (сеточная геофильтрация и геомиграция; А.В.Лехов и др.) [5].
Особый вид гидродинамического ПО представляют программные разработки А.В.Кирюхина (С-ПбГИ - ДВО СО РАН) [9] и О.П.Полянского, В.В.Ре-вердатто (ОИГГМ СО РАН, г. Новосибирск) [10], предназначенные для моделирования тепло- и паропереноса в пределах гидротермальных систем, а также обособленно пока развивающееся программное обеспечение геокриологии.
Наиболее известными зарубежными аналогами перечисленных программных продуктов* является моделирующая система GMS (Groundwater Modeling System - система, объединяющая модули MODFLOW, MODPATH, MT3D, RT3D, FEMWATER, SEAM3D, SEEP2D, PEST, UTCHEM и UCODE) лаборатории Brigham Young University США, пожалуй, одна из наиболее развитых и сервисных в настоящее время, MODFLOW геологической службы США (USGS) и VisualMODFLOW фирмы WHI Software. Как самостоятельный продукт, поставляется также мощная система визуализации результатов моделирования Visual Groundwater. Для обработки опытных данных применяются программы AquiferTest, AquiferWin, Adept и др. В целом, можно говорить о большей сервисности и, особенно, тестированности зарубежных программ по сравнению с отечественными, хотя с содержательной точки зрения российские продукты, как правило, обладают более широкими наборами используемых методов и, соответственно, возможностями корректного решения наиболее распространенных задач геогидродинамики.
Подобная ситуация наблюдается и по гидрогеохимическому ПО. Практически не имеют близких по возможностям зарубежных аналогов геохимические моделирующие системы, основанные на методе минимизации свободных энергий Гиббса, GIBBS■ Ю. В.Шварова.,(МГУ) [11, 12] и Селектор-С И.К.Карпо-
• http://www.scisoftware.com (Groundwater Geoenvironmental Software. Scientific Software Group),
http://water.usgs.gov/software, h ttp://www.chemcad.fr/en/index.htmi и др.
ва, К.В.Чудненко и др. (Институт геохимии СО РАН, г. Иркутск) [13], а также разработки Н.Н.Акинфиева (МГГА, г. Москва) [14]. Последние версии этих программ позволяют в упрощенной постановке моделировать геохимические процессы в потоке подземных вод [15-17].
Не уступают лучшим зарубежным системам, а во многих случаях и превосходят их, и программные средства, базирующиеся на методе "констант равновесия". Наиболее интересны среди них генератор гидрогеохимических моделей В12 и полученная с его помощью серия узкоприкладных программ-имитаторов SOXXXX фирмы СофДек В.Н. и С.В.Озябкиных (г. Санкт-Петер-бург), в т.ч. учитывающих геомиграцию [18], и программная система MIF Г.А.Соломина (ВСЕГИНГЕО) [11, 19].
К этому ряду примыкает и разработанный в ТПУ программный комплекс HydroGeo [20], объединяющий одновременно гидродинамические и гидрогеохимические модули, что позволяет отнести его к одному из вариантов АРМ-гидрогеолога. Последние его версии (HG32 для Windows), наряду с аналитическими гидродинамическими расчетами, включают полноценные одно- и двумерные сеточные модели геомиграции [21]. В отличие от других гидрогеологических программных продуктов, данный ПК ориентируется как на традиционные гидрогеологические задачи, так и на специфику глубокозалегающих подземных вод и методов нефтегазовой гидрогеологии, и может широко использоваться в этой области.
Среди зарубежных программ гидрогеохимического направления наиболее интересны HYDROGEOCHEM, также, судя по рекламным описаниям, позволяющая решать задачи геомиграции с учетом протекания геохимических процессов, PHREEQE и PHRQPITZ. Для первичного анализа гидрохимических данных широко применяется программа AquaChem, которая может использовать PHREEQC.
Специализированное ПО, применяемое для изучения подземных вод, использует в принципе близкие методы обработки информации, поэтому ниже более подробно рассматриваются методы, реализованные в программном комплексе HydroGeo (ПК HG32).
МЕТОДЫ ОБРАБОТКИ ГИДРОГЕОЛОГИЧЕСКОЙ ИНФОРМАЦИИ
Реализованные в ПК HG32 гидрогеологические расчеты условно объединены в три группы:
1) вспомогательные, предназначенные для оценки фильтрационно-ёмко-стных свойств пород, пересчетов данных по химическому составу воды, газа, минеральной фазы и др.,
2) прямого расчета и оптимизации объектов добычи и подземного захоронения вод на основе аналитических методов и
3) имитационного моделирования фильтрационного движения и физико-химического моделирования геохимических процессов, изменяющих вещественный состав подземных вод и вмещающих их пород, на численных 0-, 1- и 2-мерных моделях. •
Определение фильтрационно-ёмкостных свойств водоносных горизонтов
В качестве основных расчетных зависимостей для определения фильтрационно-ёмкостных свойств по результатам гидрогеологических откачек и испытания глубоких скважин, используются широко известные формулы геогидродинамики, записанные применительно к системе СИ и специфике вод глу-
боких горизонтов (с учетом их плотности и вязкости):
1) логарифмического приближения уравнения Тейса, описывающего плоско-радиальную фильтрацию флюида к скважине:
О и ахс
, (1)
47СтКп г"
где Р0 - пластовое давление (Па), Р3 - текущее забойное давление (Па), <5 - текущий дебит (м3/с), ш - вязкость флюида (Па с), ¡1 - эффектная мощность пласта (м), Кц - коэффициент проницаемости (м2), а - коэффициент пьезопроводиости (м2/с), х - время замера Р3 (с), г - расстояние от оси скважины до внутренней границы изучаемой зоны дренирования: в случае одиночной скважины - ее радиус в зоне залегания пласта, а для наблюдательной - расстояние до нее (м), с=2.24584 - константа, вычисленная как 4/еу (V - постоянная Эйлера);
2) обработки кривых восстановления давления по методу Хорнера:
(2)
- 47ГшКГ[ т ' 1 ;
где Qcp - средний дебит за время притока (м3/с), Т - продолжительность притока, х - текущее время КВД (с).
При этом давление и напор (в метрах) связаны соотношением Р^рН, а ^ коэффициенты фильтрации (в м/с) и проницаемости зависимостью К^КП где g - ускорение свободного падения (в среднем §—9.80665 м/с2) и р - плотность воды (кг/м3).
Алгоритм определения фильтрационно-ёмкостных параметров использует стандартные графоаналитические методы обработки кривых понижения и восстановления давления или уровня Тейса-Джейкоба и Хорнера-Сейза [22], полученных при откачках/закачках флюида. Последний применяется также для расчета параметров по кривым восстановления давления (КВД) в глубоких скважинах, опробованных при помощи испытателя пласта (ИП).
При использовании этих алгоритмов в процессе визуализации графиков, для более корректной их интерпретации, выполняются проверки времени наступления квазистационарного режима хк>2.5г2/а и времени хе>20-ТСг28/©,
начиная с которого погрешность расчетов за счет влияния ёмкости ствола скважины не превышает 5% .
Ввиду необходимости учета одновременного изменения дебита и забойного давления в ходе опробования глубоких скважин с помощью ИП или испытания в колонне, данные по получаемым при этом кривым притока (КП) обрабатываются аналитическим методом, базирующимся на численном интегрировании КП и КВД на основе принципа суперпозиций (наложения течений). Используемые при этом расчетные зависимости вытекают из уравнения (1), в соответствии с которым разность пластового и текущего забойного давления для каждой точки N кривой притока, начиная с 3-й, составляет
Р0 - К'
Х<эп -VI Инч 1пс1Л
п=0 Хп и * у
(3)
В случае, когда Р0 заранее не известно, используются последовательные разности выражений, записанных в соответствии с уравнением (3) для каждой пары расчетных точек:
Рц - Рнч = ^N-1 """ Т. ХТн V, ).....,........^ ^ ¡п (^N-1 ~~ ^N-2 Х^М ыч) _
Коэффициент пьезопроводности во всех случаях, кроме кустовых откачек (или гидропрослушивания), рассчитывается, методом последовательных приближений, используя равенство, примененное при программной реализации метода суперпозиций:
Кп
" " ^оРф+Э)5 (5)
где п0 - открытая пористость коллектора в д.е., - коэффициент сжимаемости флюида (при отсутствии лабораторных определений для воды может быть задан равным 3 10~10 Па"1), (3 - коэффициент сжимаемости пористой среды (породы), при необходимости определяемый для гранулярных коллекторов по уравнению регрессии, описывающему график Холла (точность ±6%; Па-1):
Р=4.15 1012п0+1.9910"11. (6)
В первом приближении задается а=10 м2/сут и по уравнениям (3) или (4) рассчитывается значение Кп в расчетных точках, затем они подставляются в уравнение (5), по которому определяются уточненные значения а, используя которые уточняются величины Кп и т.д., пока а и Кп в ходе итераций не перестанут изменяться в пределах заданной точности их оценки. Возможные вариации действительных значений п0, и Р по сравнению с принятыми в расчете мало влияют на достоверность оценки Кп, но могут изменить величину а в пределах одного-двух порядков, что вполне допустимо для одиночных опытов.
Ориентировочная величина скин-эффекта 8, суммарно характеризующего несовершенство скважины по характеру и степени вскрытия пласта, для одиночной скважины вычисляется по зависимости [23]:
'К-
где КпД и Кп?0 - проницаемости в удаленной и призабойной зонах пласта, соответственно (например, полученные по последним и первой сериям расчетных точек), г3 - радиус зоны призабойного изменения фильтрационных свойств
пласта при проходке скважины (обычно около 5 см).
В случае если давление насыщения пластовой воды газом (Рг) попадает в интервал между начальным и конечным давлениями КП, то во время притока происходит скачкообразное увеличение коэффициента продуктивности (д)
при смене двухфазной фильтрации на однофазную, когда Р3 = Рг. В отсутствие выраженного эффекта смыкания трещин это может позволить оценить Рг по характеру зависимости а = ДР3).
Приведенные выше формулы записаны для условий напорного режима фильтрации. Для безнапорных водоносных горизонтов в них производится переход от давлений к напорам и понижениям уровня и замена 28 т на 8{2Н-8), где Н - мощность безнапорного горизонта. Соответственно, а в этом случае заменяется на ау - коэффициент уровнепроводности.
Расчет и оптимизация водозаборов или объектов захоронения вод Аналитические расчеты выполняются для конкретных схем размещения водозаборных или нагнетательных скважин (или иных горных выработок) с учетом граничных условий в плане и сводятся к оценке величины понижения уровня в каждой из них на конец срока эксплуатации. Кроме того, в зависимости от полученных результатов предусмотрена возможность автоматизированной оптимизации дебитов и размещения скважин (из подобных программных разработок известна программа расчета и оптимизации неупорядоченного размещения скважин на основе ЛП-поиска, которая использует смещение координат по равномерно распределенной последовательности точек, разработанная М.С.Шутовым в ЗапСибНИГНИ). Все расчеты проводятся для условий однородного (водопроводимость и пьезо- или гидропроводность не изменяются) горизонта, изолированного от выше- и нижележащих водоносных толщ. Переход к безнапорному режиму фильтрации, описанный выше, является относительно несложным, поэтому описание расчетов для безнапорных условий далее приводится в сокращенном виде.
Понижения уровня в скважинах и любых других заданных расчетных точках (формально они отличаются от скважин лишь нулевым дебитом) вычисляются по принципу сложения фильтрационных течений, в соответствии с которым они определяются как сумма срезок уровня, обусловленных независимой работой каждой из эксплуатируемых скважин. Основная расчетная зависимость, базирующаяся на формуле логарифмического приближения уравнения Тейса (1), с учетом различного времени пуска скважин и их дебита, в этом случае принимает вид
1 п ( ~ 8' " 47ГГ ^^
1п Ь
(8)
где 8| - суммарное понижение уровня в расчетной точке 1, вызванное работой
водозабора, м; Т - водопроводимость горизонта (Т-Кт), м2/сут; п - число расчетных эксплуатационных скважин (включая зеркально отраженные при задании границ пласта); С^ - дебит расчетной скважины ], м3/сут; т- продолжительность откачки из скважины ] до расчетного момента времени где Тр - расчетное время, а ~ время запуска скважины), сут; - расстояние между расчетной точкой 1 и скважиной м; Ь - для однородного горизонта является постоянной величиной, зависящей от пьезопроводности.
При квадрат расстояния между расчетной точкой и расчетной скважиной вычисляется по формуле: ¿у^УгУ^+^-х^), где х, у - координаты, а когда 0 и 1 совпадают, обозначая одну и ту же скважину) - /у2^2, где г - ра-
диус скважины. Параметр Ь определяется как Ь=1п(а-с), где а - пьезопровод-ность, а с=2.24584.
В случае безнапорного водоносного горизонта, имеющего мощность Н0, зависимость (8) принимает вид
н2 У ^ 0
/ ' \ + Ь
Ч У
Для слабо неоднородных по фильтрационным свойствам горизонтов параметры К, т и а определяются как средневзвешенные по площади развития воронки депрессии, а при существенной неоднородности вместо их обобщенных величин можно использовать конкретные значения в каждой расчетной точке. Например, для напорных условий:
8 = 1.....
' \ Т,1
\
1п_ц 1 + \пс
'5
где Ту и а^ вычисляются как среднее арифметическое по расчетным точкам, попадающим в эллипс с осями /у, 0.5-/^, вписанный между точками 1 и }.
В расчетах учитываются два основных вида границ пласта в плане: постоянного напора (например, гидравлически связанный с эксплуатационным водоносным горизонтом водоём или водоток) и непроницаемые (выклинивание горизонта, изолирующее разрывное нарушение и т.п.), для чего используется метод зеркальных отображений. При этом границы рассматриваются как бесконечные прямо л инейные, а число расчетных скважин в зависимости (8) определяется выражением:
п = пс(г + 1)9 (9)
где пс - число реально существующих скважин и ж - число учитываемых границ. .
В случае непроницаемой границы расчетные дебиты отраженных скважин приравниваются к дебитам соответствующих им реальных скважин, а в случае границы с постоянным напором - равны дебитам реальных скважин с обратным знаком. -
Если задать каждую из границ координатами двух лежащих на их линии точек х1?У! и х2,У2 (рисунок), то координаты отраженных скважин определяют уравнения:
-х] ~2(Ах1 + В-у^-зта-соза 2
Уко = У* - Т" (Ах, + В - у,) • 81Г1 а • соэ а, (10)
у у
где А = -1 --, В - у, - Ах,, а = агс!ц А .
х, — х2
Очевидно, что при у]/=у2, х[9 а при х1=х2? У|?0=У| (границы параллельны осям координат). В этих случаях У^о^Уг ЖУтУх) и х1>0^х1-2(х1-х1), соответственно.
Расчетные значения в точках, совпадающих с эксплуатационными скважинами, сравниваются с величиной максимального допустимого пони-
.•355
У
о
Ма \
В
X
Рис. Размещение реальной и отраженной скважин относительно границы
с,
х
А8СК = ^ 2 ИТ
жения вм. Последняя при существенном несовершенстве скважин должна предварительно уменьшаться на величину дополнительной срезки уровня (А8СК) за счет скин-эффекта (з):
СМ
(И)
Размещение скважин водозабора может считаться оптимальным, если понижения уровня во всех скважинах на конец расчетного срока эксплуатации будут близки к величине допустимого понижения, поскольку при 81<8М запасы избыточны, а при 81>3М - недостаточны, В первом случае оптимизация размещения скважин заключается в уменьшении расстояния между ними, а во втором - в его увеличении, либо в удалении водозабора от непроницаемых границ или же приближении к границам постоянного напора, для чего используется процедура обработки изменения координат каждой из скважин, не удовлетворяющих условию: : •
<8М±5, (12)
где 8 - заданная точность, при их перемещении в поле графического окна с помощью мыши или клавиатуры.
Необходимые для проведения оптимизации расстояния скважины от границы определяются по зависимости
' 1ьг=(Ах1+В~у|)со8 а. " (13)
При изменении координат каждой скважины все расчеты повторяются до выполнения условия (12) по всем скважинам,, Критериями возможности перемещения каждой конкретной скважины также являются: её сближение с другими скважинами или с границей постоянного напора на минимально допустимое расстояние, а также удаление от соседней скважины или от непроницаемой границы до исчезновения их взаимодействия« При невозможности выполнения условия (12) за счет изменения размещения эксплуатационных скважин, требуется принципиальный пересмотр исходной схемы водозабора (его суммарной производительности, числа, размещения, времени запуска/остановки и деби-тов скважин) и, возможно, конструкции скважин (диаметров, степени гидравлического совершенства, типа водоподъёмного оборудования и т.п.).
Аналогично, кроме/вместо оптимизации размещения, выполняется оптимизация дебитов скважин, для чего при заданном (или полученном в результате оптимизации) их размещении выполняется пошаговый подбор величины. оптимального дебита по каждой из скважин водозабора с заданной точностью. Дебиты в этом случае ограничиваются сверху величиной максималь-
.ной водозахватной (водоприёмной) способности, задаваемой дополнительно, исходя из конкретной конструкции фильтров эксплуатационных скважин.
Численное моделирование совместно протекающих гидродинамических и
гидрогеохимических процессов Основой совместного численного моделирования гидродинамических и гидрогеохимических процессов является расщепление вычислений на условно независимые гидродинамическую и геохимическую составляющие. При этом осуществляется разбиение области фильтрации на относительно однородные по фильтрационно-ёмкостным свойствам и минеральному и химическому составу пород и флюидов блоки, образующие пространственную сетку, между ячейками которой с использованием методов конечных разностей, конечных элементов или граничных интегральных уравнений моделируется фильтрация водных растворов и связанные с ней процессы гидродинамической и диффузионной дисперсии вещества вод и смешение растворов. На каждом шаге или через заданное число временных шагов геофильтрации для всех расчетных блоков выполняется моделирование действия внешних и внутренних источников-стоков вещества. Под последними обычно подразумеваются собственно гидрогеохимические процессы: растворение-осаждение, сорбция-десорбция и радиоактивный распад, причем каждый из таких процессов также рассматривается как условно независимый.
По сравнению с традиционно применяемыми аналитическими методами, практически всегда требующими более или менее значительной схематизации природных условий, такой подход в большинстве случаев позволяет лучше учесть пространственную неоднородность фильтрационных свойств и состава пород и флюидов, сложный характер и разнообразие граничных условий, инфильтрацию, перетекание, изменение Р,Т-условий и другие факторы, реально присущие сложным гидрогеологическим объектам. Это обусловливает возможность получения с их помощью более точных и достоверных оценок поведения водоносной системы, хотя, в то же время, современный уровень как теоретического, так и практического развития методов, программ и технического обеспечения численного моделирования геофильтрации и геохимических процессов во многих случаях пока еще явно недостаточен для их применения в весьма сложных гидрогеологических условиях без существенного их упрощения и схематизации. Значительный интерес представляет также использование численного моделирования для создания и эксплуатации постоянно действующих моделей реальных гидрогеологических объектов. Такие модели могут постепенно уточняться и/или изменяться по мере накопления гидрогеологической информации.
Моделирование геофильтрации Общая конечно-разностная форма дифференциального уравнения нестационарной плановой фильтрации, например в изотропном напорном пласте, имеет вид
а
х Эх
(14)
■ Эх , , .
. Эх
\ /
где Т - водонроводимость, Н - напор, Ш - расход источников/стоков, {1 -удельная емкость пород (водоотдача), а "..." здесь и далее заменяет аналогичные предшествующему члены по пространственным координатам у и г.
Переход от дифференциальной записи к разностной может быть выпол-
нк =
нен после разбивки области фильтрации на сеть элементарных блоков/узлов и записи для каждого блока уравнения баланса воды. В ПК Н032 для этого использован метод конечных разностей в форме, предусматривающей итерационный расчет напоров /давлений с использованием методов релаксации и переменных направлений по зависимости (при итерационном "прогоне" матрицы блоков в направлении роста индекса 1)
ш^з^Н^ -ьр^Н^1 +...-нЕ/АтхНк~1 + \¥к~|) + (1 ^ ^ ^ (15)
где р - межблочная проводимость среды между текущим и смежными блоками по координате с индексом блоков 1, Е - текущая емкость блока, со - ускоряющий коэффициент релаксации (используется значение 1.7).
Массоперенос
В общем случае массоперенос вещества в движущемся потоке вод включает несколько различных процессов, основными среди которых являются [5]:
1) Конвективный перенос фильтрующимся потоком (вынужденная конвекция) под действием гидравлического градиента, протекающий со скоростью и=у/п, где V - скорость фильтрации Дарси, а п - эффективная пористость. Удельный поток вещества (масса вещества, перемещающаяся в единице объема породы за единицу времени) за счет конвективного переноса 1к=иС, где С - содержание компонента в растворе (или его концентрация).
2) Плотностная конвекция, связанная с наличием вертикального градиента плотности растворов. Ее скорость Ур=К2Др. В этом выражении ур - вертикальная составляющая скорости фильтрации за счет градиента плотности, К2 - коэффициент фильтрации в вертикальном направлении, а Др=(рв-рн)/рн, где индексы "в" и "н" отвечают выше- и нижезалегающим водам. Из-за эффекта начального градиента фильтрации движение в этом случае может возникать лишь тогда, когда минерализация более плотного вышезалегающего раствора на 2-5 г/л выше, чем у менее плотного нижезалегающего. Удельный поток
1р==УрС' аг
3) Диффузионный перенос, подчиняющийся закону Фика 1т = Ом ,
а/
где Вм - коэффициент молекулярной диффузии в водонасыщенной среде, а / -расстояние в направлении вектора градиента концентрации.
4) Гидродинамическая дисперсия, контролируемая коэффициентом дисперсии В¥=8уо Используемый здесь коэффициент дисперсивности 8 определяется размерами и структурой пустотности породы, а удельный поток вещест-
дС
ва за счет дисперсии 1Эр = Ву — . Поскольку диффузия и дисперсия являются
физическими аналогами, разделить действие которых по данным натурных
экспериментов крайне сложно, может быть введено понятие суммарного ко-
Э С
эффициента дисперсии Б=ВМ+ВУ. В этом случае \т + 10р = В -- определяются
: Э/
совместно, например, по данным геомиграционных опытов.
; Общее уравнение баланса вещества записывается как
•■«>
где 1 - суммарный удельный поток вещества в направлении координаты, указанной нижним индексом, - удельмшя объемная интенсивность источников-стоков. С учетом перечисленных процессов миграции, суммарный удельный поток по каждой из рассматриваемых пространственных координат 1=1к+1р+ГО£+Шр.
В настоящее время в ПК НШ2 для расчета дисперсионной составляющей массопереноса использован подход, принципиально отличающийся от описанного выше, который заключается во введении для учета гидродинамической дисперсии модели двойной пористости [18]. В этом случае общая пустот-ность породы (открытая пористость) условно разделяется на проточную (эффективная пористость) и непроточную (разность открытой и эффективной пористости) составляющие, а расчет сводится к замещению части проточной пу-стотности расчетного блока втекающим раствором, конвективному смешению содержащихся в проточной части блока добавившегося и остающегося растворов и полному, либо частичному, если задаются соответствующие кинетические ограничения, диффузионному смешению растворов проточной и непроточной зон. Ввиду крайне медленного протекания в большинстве реальных ситуаций молекулярной диффузии, её самостоятельный вклад в суммарную дисперсию пока не учитывается.
Величина ^^ включает в себя действие как внешних по отношению к системе источников/стоков вещества (например, привнос вместе с закачиваемым в скважину раствором или инфильтрующимися атмосферными осадками), так и внутренних, главным среди которых является физико-химическое взаимодействие раствора с вмещающими породами, ск л ад ывающееся, как уже отмечалось, из процессов растворения/осаждения минералов породы и ионного обмена между раствором и поглощенным комплексом породы, то есть собственно гидрогеохимические процессы. В ряде случаев, кроме того, в число внутренних источников/стоков необходимо также включать изменение состава раствора и породы вследствие радиоактивного распада.
Таким образом, алгоритм моделирования геомиграции включает несколько этапов:
1) Геофильтрационный расчет, в результате которого определяется поле напоров/давлений в каждом из блоков сеточной модели на момент времени тк.
2) Независимый расчет объемов перемещения растворов и содержащихся в них компонентов между смежными блоками - расчет конвективной и дисперсионной составляющей массопереноса на тк за Дт.
3) Самостоятельный расчет действия источников/стоков вещества, включая расчет взаимодействий в системе вода-порода, на том же или кратном нескольким шагам Дт временном шаге геохимического моделирования Дтгх.
В результате расчетов в соответствии с пп. 1-3 определяется распределение в расчетных блоках/узлах текущих напоров/давлений, а также составов и других расчетных параметров жидких и твердых фаз системы, после чего повторяется выполнение всех расчетов для последующего временного шага и т.д., до достижения заданного конечного момента времени.
Следует отметить, что в принципе, при сохранении заданной цикличности, данные расчеты могут сменять друг друга в произвольном порядке, а при наличии многопроцессорных вычислительных систем или возможности рас-
пределенных сетевых вычислений выполняться параллельно,.Также параллельно , одновременно для нескольких расчетных блоков по числу, доступных процессоров вычислительной среды, могут выполняться наиболее ресурсоёмкие расчеты действия внутренних источников/стоков вещества.
Источники/стоки вещества В основу моделирования внутренних источников/стоков вещества в соответствии с применяемым в ПК НС32 "методом констант равновесия" положено понятие элементарных реакций
г-ЬД+^ + ^В^В^..^, . <17>
где г - элементарная реакция; В^ Ъг - 1-я частица раствора и её стехиометричес-кий коэффициент в реакции г; В1Ь1...ВпЬ^ образующийся в результате г минерал, либо ионный ассоциат, совокупность которых исчерпывающе описывает анализируемые природные процессы, а также методы равновесной термодинамики и химической кинетики. Для учета неидеальности раствора применен метод активности Льюиса [24]. Расчет коэффициентов активности компонентов раствора и активности растворителя-воды основан на методике Питцера [25].
При этом в качестве параметров элементарных процессов рассматриваются:
оТЛг=о?+(298-т)8;+о;вГ+о;>г,
ЬКТЛг=-ЦЗ т,рДЕ1)Л (18)
-i
щ '
Т . Т £0 ■ • р' • • • .
где параметры G^ = -j C°r -dT-T и о;?г = jvT°r vdP' вычисляются,
298 298 ^ P¿ . .
по формулам интегрирования после подстановки функций €P°=f(T) и VT°=f(P') в виде полиномов четырех видов: Майера-Ке л л и, Карпова или Тангера-Хелге-сона, либо принимая их постоянными; Р - давление в термодинамических Дж; &т,р,г> G°r, S°r, С°Р)Г, V°r - мольные изменения термодинамических параметров в ходе реакции г при заданных и стандартных Т-Р условиях; R - универсальная газовая постоянная; КТ>Р?Г - термодинамическая константа равновесия, Ра г -произведение активностей а компонентов раствора и минерала, участвующих в реакции г, получаемое в соответствии с законом действия масс.
Моделирование растворения-осаждения проводится по зависимостям
kv=k?^, lr = 1пктг -1праг,
vr =0.3(T-298Xnr/mr)drLrBDr, r = rw=max,
(19)
v = yd, ittbn и kv = рт,п к + xdr г'.
i i
Здесь aI r и KT?r - текущие условное время и константа равновесия; Lr и vr -параметр насыщенности раствора и условная скорость осаждения минерала Dr в процессе г; dr и br - относительная скорость диффузии и число составляющих минерал Dr независимых ("базовых") компонентов; шв - молярность частицы
раствора В1? входящей в минерал Вг; Ру г - произведение коэффициентов активности для реакции г в соответствии с законом действия масс; хВг- недостаток насыщения раствора по отношению к минералу Вг или шаг реагирования (число молей породы, переводимых в единицу объема раствора на каждой итерации по минералам породы; определяется итерационным путем по методу дихотомии).
Расчет модели комплексообразования, необходимый для изучения форм миграции и определения действительных (а не валовых, получаемых при химическом анализе) концентраций компонентов в растворе из Ег ионных ассо-циатов и комплексных соединений, включённых в систему моделирования, проводится по формулам
mDr =кт,р,л:'птв,> *v =хкт,рл>лтв,>
(20)
где i = l,n, г = l,Rr, jei, icr, и
t \
тв, = mBj,i 1
£кТ>РдР?>^"'EK', | l = i\j. (21)
При этом текущие содержания иона водорода, гидроксил-иона и активность электронов, используемая при моделировании окислительно-восстановительных взаимодействий, определяются на основе условия электронейтральности тн\онг 1>.т. ^ и зависимости = j, где z - заряд, F - постоянная Фарадея и Eli - окислительно-восстановительный потенциал раствора.
Сходимость расчета комплексообразования обеспечивается путем итераций по методам среднего геометрического, а затем Ньютона раздельно сначала по катионам, потом по анионам и нейтральным частицам, после чего оценивается электронейтральность и вычисляются величины рН и Eh раствора.
Аналогичный подход применяется для расчета водно-газовых равновесий, но вместо приведенных термодинамических зависимостей используется система уравнений регрессии, полученная на основе эмпирических (экспериментальных) данных по растворимости газов в водных растворах [26]. Примененная методика позволяет обрабатывать данные по составу водорастворен-ных и свободных тазов и выполнять расчеты фазовых равновесий в системе вода - природный газ, в том числе при совместном моделировании геохимических процессов в системе вода-газ-порода.
Особый случай физико-химического взаимодействия водного раствора с породой представляет также ионообменная сорбция, для моделирования которой в настоящее время использован метод, основанный на применении констант фазового распределения, которые определяются по данным лабораторных экспериментов или полевых (натурных) геомиграционных опытов,
С В тв
К-р, =-;■--"-........, CB. TR = Кр^г -СВ(1Ж , (22)
В; ,Ж
где индексы "тв" и "ж" обозначают твердую и жидкую фазы.
Очевидно, что в данном случае КР)Г=Кт^р>гРу>г и при наличии термодина-
мических параметров и способа оценки активности частиц раствора сорбированных на поверхности минералов для описания сорбции может использоваться обычный физико-химический подход, описанный выше»
В модели ионообменной сорбции принято, что ёмкость обмена, сосредоточенная на поверхности водонасыщенной твёрдой фазы, является электронейтральной за счёт компенсации разности заряда сорбированных из раствора катионов и/или анионов ионами ОН или Н+, являющимися продуктом диссоциации растворителя - воды. Величина ёмкости обмена на данном этапе считается постоянной и всегда полностью заполненной.
ЗАКЛЮЧЕНИЕ
Приведенное рассмотрение свидетельствует об активном развитии рынка профессионально ориентированных гидрогеологических программных средств. Их разработка потребовала модификации известных и создания новых методик и алгоритмов обработки информации, многие из которых отличает своеобразие и принципиальная новизна теоретических положений, положенных в их основу. Развитие гидрогеологического ПО проводится в двух основных направлениях: гидродинамическом и гидрогеохимическом, по каждому из которых за два последних десятилетия накоплен большой объём весьма разнообразного и разнопланового программного обеспечения, в том числе и первые программные разработки, объединяющие гидродинамические и гидрогеохимические модули в виде единых геомиграционных моделей, либо АРМ, в рамках которых они могут функционировать и как единое целое и самостоятельно.
Применительно к ПК Н032 гидродинамическое направление включило стандартные процедуры обработки результатов гидродинамических исследований скважин с целью расчетной оценки фильтрационно-ёмкостных параметров (проницаемости, пьезопроводности, пластового давления, скин-эффекта, коэффициента продуктивности, давления насыщения жидкой фазы газом, плотности флюида, сжимаемости пористой среды), алгоритмы подсчета эксплуатационных запасов подземных вод и оптимизации размещения проектных скважин водозаборов или производительности систем нагнетания, а также процедур, обеспечивающих достаточно эффективный расчет многомерных сеточных гидродинамических моделей.
Второе, гидрогеохимическое направление, занимающее в ПК основную часть его объёма, реализует алгоритмы наиболее характерных типов гидрогеохимических расчетов (пересчеты единиц измерения концентраций, определение активностей, внутренняя балансировка анализов с учетом рН и ЕЬ, расчеты комплексообразования, смешения, испарения и т.п.), связанных с ионно-солевым и газовым составом подземных вод и гидрогеохимическими процессами в системе вода - газ - порода. Особенностью последних является отсутствие принципиальных ограничений применяемых методов в зависимости от минерализации или ионной силы реальных природных растворов, что достигается использованием для определения активности компонентов ионно-солевого состава авторской модификации метода Питцера, а для водорастворенных газов - эмпирических уравнений, полученных на основе обработки экспериментальных данных. Это обусловливает возможность расширения области применения гидрогеохимического блока ПК на глубокозалегающие воды, характеризующиеся весьма широким разбросом возможной минерализации и состава.
Другой особенностью данного блока, является наличие вспомогательных
процедур, с помощью которых осуществляется выбор минимально необходимого в каждом конкретном случае списка участвующих в расчете компонентов раствора и минералов из имеющихся в постоянно пополняемой базе термодинамических данных. Наличие такой "настройки" позволяет сохранить универсальность гидрогеохимического блока ПК в отношении различных геохимических систем и оперативно управлять быстродействием алгоритмов и расходом оперативной памяти, являющихся наиболее критичными для сложных гидрогеохимических моделей. В связи с недостаточной согласованностью значительной части справочных термодинамических параметров, в системе гид-рогеохимических. расчетов предусмотрена возможность их автоматической корректировки в результате обработки экспериментальных данных по растворимости минералов (решается обратная задача). При этом невязка расчетных и экспериментальных данных устраняется за счет изменения величины стандартной свободной энергии образования минерала и, при наличии данных по его растворимости при разных температурах, - за счет изменения стандартной энтропии, что повышает достоверность результатов моделирования процессов растворения - осаждения для изученного интервала температур.
Наличие в ПК большого числа относительно самостоятельных гидродинамических и гидрогеохимических модулей обусловило относительно малую трудоёмкость разработки и включения в его состав алгоритма моделирования геомиграции, учитывающей наиболее важные из протекающих при этом гидрогеохимических процессов. Следует отметить, что включение таких алгоритмов неизбежно обусловливает возрастание объёмов затрачиваемой памяти и числа необходимых вычислений ЭВМ, тесно взаимосвязанных друг с другом, сразу на несколько порядков, требуя привлечения даже для относительно простых моделей наиболее мощных персональных компьютеров, а в идеале -высокопроизводительных многопроцессорных рабочих станций или сетей для распределенных вычислений, либо суперкомпьютеров.
Рассмотренный подход к одновременному моделированию геофильтрации и геомиграции, безусловно, не охватывает всех сторон этой достаточно сложной проблемы, но уже в данной постановке и реализации позволяет решать достаточно сложные реальные гидрогеологические задачи. Накопленный при работе над ПК НО опыт позволяет сформулировать несколько общих принципов, которые могут стать полезными при разработке аналогичных программных средств:
1. При проектировании ПО необходим тщательный отбор решаемых с его помощью задач. Предпочтительными среди них являются, с одной стороны, наиболее распространенные (даже если они вполне тривиальны), а с другой -те, для которых применение ЭВМ дает существенный выигрыш во времени и трудоемкости, или в принципе не решаемые другими методами.
2. Программные средства должны быть выполнены в виде самостоятельных модулей или внешних процедур, объединяемых в специальные библиотеки, но в то же время позволять использование общего интерфейса пользователя, иметь при необходимости пересекающееся пространство констант и переменных, и быть взаимоувязанными в подчиненную общей логике единую структуру.
3. Обязательными требованиями к программной реализации должны стать достаточно гибкие возможности текущей настройки моделей и сценариев моделирования и оперативное динамическое распределение и перераспределение памяти, отводимой каждой из расчетных моделей.
4. Сервис ПО желательно создавать исходя из требований минимальной необходимости и достаточности по мере тестирования и возникновения в нем действительной потребности пользователей.
Работа выполнена при поддержке Российского фонда фундаментальных исследований (гранты № 02-05-64623 и 02-05-81014).
ЛИТЕРАТУРА
1. Цветков В.Я. Геоинформационные системы и технологии. - М.: Финансы и статистика, 1998, - 288 с.
2. Коновалова Н.В., Капралов Е.Г. Введение в ГИС. Учебное пособие. -М.: ГИС-Ассоциация, 1997. - 160 с.
3. Мироненко В* А*, Румынин ВТ. Проблемы гидрогеоэкологии. Т. 1. Теоретическое изучение и моделирование геомиграционных процессов. - М.: Изд-во МГУ, 1998. ~ 611 с.
4. Крайнов С.Р. Обзор термодинамических компьютерных программ, используемых в США при геохимическом изучении подземных вод. Система компьютеризации научных лабораторий США. // Геохимия. - 1993, - № 5. -с с 685-695.
5. Гидрогеодинамические расчеты на ЭВМ. / Под ред. Р.С.Штенгелова. -М.: Изд-во МГУ, 1994 - 335 с.
6. Langmuir D. Aqueous Environmental Geochemistry. - New Jersey: Prentice-Hall Inc. Upper Saddle River, 1997. - 600 p.
7. Ломакин E.A., Мироненко B.A., Шестаков B.M. Численное моделирование геофильтрации. ~ МНедра, 1988. - 228 с.
8. Плюснин A.M., Гунин В.И. Природные гидрогеологические системы, формирование химического состава и реакция на техногенное воздействие (на примере Забайкалья). - Улан-Удэ: Изд-во БНЦ СО РАН, 2001. - 137 с.
9. Кирюхин A.B. Теплогидродинамическая модель: гидротермальная система - неглубоко залегающий магматический очаг. // Вулканология и сейсмология. - 1984. - № 3. - с. 25-34.
10. Полянский О.П., Ревердатто В.В. Конвекция флюида в коллекторах осадочного бассейна при термическом воздействии даек и си л лов. // Геология и геофизика. - 2002. - Т. 43. -№ 1. - с. 27-41.
11. Методы геохимического моделирования и прогнозирования в гидрогеологии. / Под ред. С.Р.Крайнова. - М.: Недра, 1988. - 254 с.
12. Борисов М.В., Шваров Ю.В. Термодинамика геохимических процессов. - М.: Изд-во МГУ, 1992. - 256 с.
13. Минимизация свободных энергий при расчете гетерогенных равновесий. / И.К.Карпов, К.В.Чудненко, В.А.Вычинский и др. // Геология и геофизика. - 1995. - Т. 36. - № 4. - с. 5-21.
14. Акинфиев H.H., Умрихин В.А. Алгоритм и программа SOL для расчета равновесий в природных и техногенных водах. // Геоинформатика. - 1994. - № 2. - с. 60-66.
15. Шваров Ю.В. Алгоритмизация численного равновесного моделирования динамических геохимических процессов. // Геохимия. - 1999. - № 6. - с. 646-652.
16. Кулик Д.А., Чудненко К.В., Карпов И.К. Алгоритм физико-химического моделирования эволюции локально-равновесных резервуаров, связан-
ных потоками подвижных групп фаз. // Геохимия. - 1992. - № 6. - с. 858.
17. Романова М.А., Акинфиев Н.Н. Термодинамическое моделирование как метод прогнозирования гидрогеохимических процессов в системе "подземные воды-порода-С02" при закачке техногенной углекислоты. // Известия ВУЗ. Геология и разведка. - 1999. - № 1. - с. 79-88.
18. Озябкин В.И., Озябкин С.В. Программные имитаторы для моделирования геохимической миграции неорганических загрязнений. // Геоэкология. - 1996. - № 1. - с. 104-120.
19. Соломин Г.А., Крайнов С.Р. Геохимические условия приложения расширенного уравнения Дэвиса для расчета коэффициентов активности ионов в рассолах. // Геохимия. - 2000. - № 5. - с. 510-515.
20. Букаты М.Б. Разработка программного обеспечения в области нефтегазовой гидрогеологии. // Разведка и охрана недр. - 1997. - № 2. - с. 37-39.
21. Букаты М.Б. Численное моделирование гидрогеохимических процессов в фильтрационном потоке. / Гидрогеология и инженерная геология. Геоэкология и мониторинг геологической среды. - Томск: Изд-во ТПУ, 2001. - с. 10-13.
22. Мироненко В.А., Шестаков В.М. Теория и методы интерпретации опытно-фильтрационных работ. - М.: Недра, 1978. - 325 с.
23. Карнаухов M.JL, Рязанцев Н.Ф. Справочник по испытанию скважин. - М.: Недра, 1984. - 269 с.
24. Льюис Г., Рендалл М, Химическая термодинамика. - М.: ОНТИ, Химтеорет, 1936. - 390 с.
25. Термодинамическое моделирование в геологии: минералы, флюиды и расплавы. / Р.К.Ньютон, А.Навротеки, Б.Дж.Вуд и др. - М.: Мир, 1992. - 534 с.
26. Букаты М.Б. Методика моделирования водно-газовых равновесий в связи с прогнозом нефтегазоносности // Геология нефти и газа. - 1992« - № 1. -с. 7-9. ~
SOFTWARE ENGINEERING FOR THE SOLUTION OF HYDROGEOLOGICAL
PROBLEMS
M.B.Bukaty
The tendencies and state of the art of development of the hydrogeological software, some methods of the automated processing geohydrodynamic and hydrogeochemical of the information are reviewed. The guidelines on mining modelling programmatic complexes and hydrogeological automated workstation are formulated. The idealized positions, methods and algorithms realised in an authoring programmatic complex HG32 are described. The principles of construction of net models of geomigration which is taking into account hydrogeochemical processes are justified.
