Разработка программных средств моделирования процессов региональной миграции веществ в подземных водах

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

Современное состояние моделирования массопереноса.

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

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

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

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

Для регионального моделирования геофильтрации используются достаточно "грубые" конечно-разностные сетки. Например, при региональной переоценке эксплуатационных запасов и ресурсов подземных вод Московского артезианского бассейна используется равномерная квадратная сетка с шагом в плане 1000м. Количество расчетных блоков сетки в плане достигает 10⁵. Число расчетных водоносных горизонтов достигает 8-10. Таким образом, общее количество расчетных блоков геофильтрационной модели составляет порядка 10⁶.

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

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

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

В данной разработке (Сonvective Mass transport Simulation in Multilayer Systems, CMSMS) предлагается использовать новый подход, свободный от описанных вычислительных эффектов, позволяющий решать задачу массопереноса в точности на той же сетке, что и геофильтрационную.

1. Основные идеи, положенные в основу моделирования конвективного переноса и используемые допущения.

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

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

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

Пространственное построение сечений трубок тока проводится на основе общих предположений о структуре фильтрационного потока:

В силу выполнения предпосылки Дюпюи сечение трубки тока на грани блока можно аппроксимировать областью с вертикальными и горизонтальными границами (рис.2).

Рис.2 Распределение трубок тока на грани расчетного блока.
S_i-сечение трубки тока с типом воды T_i

Наличие наклонных границ между трубками тока свидетельствовало бы о существовании наклонных нейтральных линий тока, что противоречит предположению о постоянстве напора по вертикали (предпосылка Дюпюи). Таким образом, пространственное распределение воды разных типов по граням блоков можно описать матрицей составленной из сечений трубок тока S_i={q_i; T_i}, где q_i-расход, текущий через данное сечение, T_i-тип воды текущей по этой трубке тока (рис.2).

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

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

Рис.3 Распределение расходов и типов воды на боковых гранях блока модели,
а - блок геофильтрационной модели,
б - распределение типов воды по входным граням блока.

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

где q_w, q_s, q_n, q_e - абсолютные значения расходов потоков через западную, северную, восточную и южную грани соответственно.

Условие баланса в таком блоке имеет вид q_w+q_s=q_n+q_e.

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

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

В силу выполнения предпосылки Дюпюи о постоянстве напора по вертикали, граница раздела между водой типов T_w и T_s будет вертикальной. Положение границы раздела на выходных гранях можно найти следующим образом. Составим из двух выходных граней одну обобщенную грань, с суммарным расходом Q_c=q_w+q_s=q_n+q_e. Спроектировав на такую выходную грань втекающие в блок потоки, получим распределение, показанное на рис. 5а.

Рис. 5 Определение распределения типов воды на выходных гранях блока модели,
а - распределение входных потоков на обобщенной грани, q_w=q_s,
б - выходные потоки, "вырезаемые" из потока через обобщенную грань,
q_n>q_e, Q_c=q_w+q_s=q_n+q_e.

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

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

Рассмотрим случай, когда внутри блока находится совершенный источник.

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

Рис. 7 Определение распределения типов воды на выходных гранях блока модели при наличии в нем совершенного источника.
а - входное распределение расходов, спроектированное на обобщенную грань с суммарным расходом потока,
б - выходные потоки, "вырезаемые" из потока через обобщенную грань,
в - распределение сечений трубок тока с водой разных типов на выходных гранях.

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

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

В рамках модели CMSMS не рассматриваются источники/стоки несовершенные в пределах водоносного пласта, поскольку модель строится с учетом выполнения предпосылки Дюпюи. При таких допущениях поток от несовершенного источника лежащего внутри блока, на гранях блока будет эквивалентен потоку от совершенного источника, расположенного в той же точке и имеющего тот же расход. При моделировании массопереноса в многослойных системах учитывается несовершенство источника относительно всей водоносной толщи.

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

Рассмотрим перетекание в блок из верхнего слоя потока воды одного типа T_u. При этом в рассматриваемый блок есть приток воды в плане типов T_w и T_s (рис.8).

Для простоты будем считать, что в блоке есть только одна выходная грань с расходом q_e=q_w+q_s+q_u, где q_w, q_s - расходы потоков втекающих через западную и южную грани, q_u - расход поступающий в блок из вышележащего горизонта. Обобщенная выходная грань в этом случае совпадает с единственной выходной гранью блока (E).

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

Рис. 9 Распределение сечений трубок тока на выходной грани при наличии перетекания в блок из вышележащего горизонта.

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

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

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

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

2. Возможности, реализуемые в данной разработке.

Основные возможности.

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

Задачи могут решаться как в плановой, так и в многослойной постановке.

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

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

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

Схема процесса моделирования.

Процесс моделирования миграции предполагается организовать по следующей схеме:

3. Результаты некоторых тестов и сравнение с другими методами.

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

Рассмотрим однослойную плановую задачу с однородной изотропной проводимостью T=10. В модели заданы две питающие и две откачивающие скважины с одинаковым расходом Q=100, расположенных, как показано на рис. 10б. Поскольку расходы скважин равны, модель разбивается на две области, между которыми не должно существовать обмена водой. Пометим воду, закачиваемую одной из скважин. Будем считать, что вода поступающая из этой скважины имеет концентрацию пассивного трассера C₁=1. Вода закачиваемая другой скважиной имеет концентрацию трассера C₂=0. Стационарное поле концентраций полученное из очевидного аналитического решения показано на рис. 10б. На рис. 10а представлен результат моделирования методом конечных разностей (для моделирования использовался пакет MT3DMS [1]). В данном тесте ярко виден результат проявления эффекта численной диффузии, который существенно искажает решение.

Рис. 10 Тестовый пример. Однородная область, два источника и два стока.
а - метод конечных разностей, б - CMSMS

Результат моделирования CMSMS показан на рис. 10б. Поле концентраций полностью совпадает с аналитическим решением.

Сравнение CMSMS с другим популярным методом - методом характеристик можно провести на следующем тесте рис11а. Рассматривается изотропная плановая задача с однородной проводимостью T=1. По левой и правой границе сетки заданы граничные условия с постоянным напором. На левой границе H_l=800, на правой H_r=0. Симметрично относительно верхней и нижней границы сетки задан дублет из откачивающей и закачивающей скважин. Расход откачивающей скважины равен расходу закачивающей скважины Q=800.

Рис. 11 Тестовый пример. Дуплет скважин в бытовом потоке подземных вод.
а - метод характеристик, б - CMSMS

Скважины довольно слабо возмущают поле напоров, и содержащие их блоки являются проточными, т.е. в них есть приток или отток от границ первого рода. Пометим воду, поступающую от закачивающей скважины, задав концентрацию трассера C_src=1, а для воды поступающей от левой границы первого рода зададим концентрацию C_ch=0. Поскольку расход поглощающей скважины равен расходу питающей и задача симметрична, вся вода из источника должна поглощаться стоком и в любом блоке на правой границе концентрация должна быть равной C_r=0.

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

На рис. 11б показано поле концентраций, полученное с использованием CMSMS. В данном случае эффект численной диффузии полностью отсутствует.

Заключение.

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

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