Геоэкология. Инженерная геология, гидрогеология, геокриология, 2021, № 5, стр. 72-86

ВВЕДЕНИЕ

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

Для решения прогнозных геомиграционных задач широко используются методы численного моделирования с помощью программных средств, позволяющих рассчитывать распространение контаминантов в ЗА и насыщенной зоне (НЗ), как в рамках единой модели [10, 18], так и отдельно для каждой зоны [9, 19]. Расчет распространения загрязнения в рамках единой модели целесообразен при сопоставимых масштабах переноса в ЗА и НЗ. Однако данный подход сопряжен с объективными сложностями, среди которых можно назвать сильную нелинейность связи параметров, определяющих влагоперенос в ЗА и входящих в уравнение Ричардса [5], что вызывает проблемы со сходимостью решения при расчетах влагопереноса. Кроме того, в перечисленных программах отсутствуют средства создания моделей со сложным геологическим строением, а также возможности автоматической калибровки геофильтрационной модели по уровням. Разработка двух моделей миграции загрязнения отдельно для каждой зоны, с одной стороны, позволяет преодолеть описанные выше сложности. С другой стороны, возникает необходимость согласования результатов, полученных на каждой модели, для их совместного использования.

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

Применение методики рассмотрено на примере района, где выявлено хлоридное загрязнение подземных вод. Загрязнение в течение 12 лет поступает от вторичного источника в ЗА после ликвидации первичного источника.

МЕТОДИКА СОЗДАНИЯ КОМПЛЕКСА МОДЕЛЕЙ ЗОНЫ АЭРАЦИИ И НАСЫЩЕННОЙ ЗОНЫ

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

Перед началом разработки комплекса моделей необходимо решить вопрос о целесообразности рассмотрения миграции загрязнения в ЗА для конкретного района исследований. Это можно сделать на основе качественной оценки защитных свойств отложений ЗА [2]. Важно подчеркнуть, что все характеристики ЗА (мощность, фильтрационные, водно-физические и сорбционные свойства) нужно рассматривать совместно, поскольку ЗА, имея мощность даже первые метры, может обладать значительными задерживающими свойствами по отношению к определенным видам загрязнения. Также возможна количественная экспресс-оценка времени миграции загрязнения в ЗА, например, с использованием известной формулы Н.Н. Биндемана [7].

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

В противном случае, за основу предпочтительнее брать модель НЗ, в которой массовый поток контаминанта задается в границах проекции источника на уровень грунтовых вод (УГВ) и калибруется по данным мониторинга (рис. 1, 2). Результатом модельных расчетов НЗ является ореол контаминанта в потоке грунтовых вод в реальных концентрациях. Затем разрабатывается модель ЗА, необходимая для оценки ее задерживающих свойств по отношению к загрязнению, поступающему от приповерхностного источника. Мощность ЗА принимается в соответствии с положением УГВ, полученном на модели НЗ. При этом миграционные модели обеих зон не связываются друг с другом граничными условиями.

Рис. 1.

Геологическая карта территории (а) и карта-врезка (б). 1–4 – аллювиальные отложения надпойменных террас р. Енисей; 5–6 – породы: 5 – мезо-кайнозойской коры выветривания, 6 – архейского возраста; 7 – границы области моделирования (по контуру АБ – непроницаемая граница, по контуру БВ – граница с постоянным напором, по контурам ВГ и АГ – границы по ручью 2 и р. Енисей); 8 – граница площадки цеха; 9–10 – границы зон с повышенными фильтрационными свойствами: 9 – гравийно-галечников в русловой фации аллювия, 10 – песков в коре выветривания; 11 – ручей 2; 12 – золоотвал; 13 – скважины и их номера: а – режимные, б – прочие; 14 – линия гидрогеологического разреза; 15 – границы карты-врезки; 16 – солехранилище: а – ликвидированное, б – действующее.

Рис. 2.

Схематический гидрогеологический разрез по линии 1–1'. 1 – техногенные отложения; 2 – аллювиальные отложения V надпойменной террасы р. Енисей; 3 – породы мезо-кайнозойской коры выветривания; 4 – гравийно-галечники с песчаным заполнителем; 5 – пески, местами глинистые; 6 – супеси; 7 – суглинки; 8 – выветрелые гнейсы (щебень, дресва, суглинки, супеси); 9 – УГВ, м; 10 – границы: а) литологические, б) стратиграфические; 11 – направление распространения ореола хлорид-иона; 12 – скважина с указанием положения фильтра, вверху – номер, слева – абсолютная отметка УГВ, м, справа – осредненные концентрации хлорид-иона, мг/л; 13 – скважина, снесeнная на разрез.

Расчеты миграции загрязнения в ЗА ведутся в относительных концентрациях. При такой постановке их основной результат – время, на которое отложения ЗА, с одной стороны, задерживают начало поступления контаминанта на УГВ (при активации приповерхностного источника), а с другой, – начало деградации ореола в грунтовых водах (при его ликвидации). В обоих случаях выделяется по два временных интервала (рис. 3а, 3б).

Рис. 3.

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

При активации источника интервал $t_{1}^{A}$ соответствует времени, когда источник уже функционирует, но поступающие на УГВ концентрации контаминанта не превышают фоновых значений (С_{фон з.а.}), а интервал $t_{2}^{A}$ – времени, когда происходит увеличение поступающих концентраций от фоновых до максимальных значений (С_{max з.а.}). После ликвидации источника в течение интервала $t_{3}^{A}$ сохраняются максимальные поступающие на УГВ концентрации, а в течение интервала $t_{4}^{A}$ они снижаются до фоновых значений.

Для модели НЗ принципиально выделяются те же интервалы времени, что и для модели ЗА (рис. 3в, 3г). При этом важно подчеркнуть, что в данном случае влияние ЗА не учитывается. Интервалы $t_{1}^{Н}$ и $t_{2}^{Н}$ связаны с распространением ореола от проекции источника на УГВ до условной наблюдательной скважины, расположенной ниже по потоку, а интервалы $t_{3}^{Н}$ и $t_{4}^{Н}$ – с его деградацией при ликвидации этого источника. Аналогичный вид будут иметь графики изменения массового потока контаминанта, поступающего в зону разгрузки, на которых также выделяются соответствующие временные интервалы.

Таким образом, зная временные интервалы $t_{1}^{A}$–$t_{4}^{A}$, полученные при разработке модели ЗА, и суммируя их с соответствующими интервалами модели НЗ (например, $t_{1}^{Н}$–$~t_{4}^{Н}$ – для конкретной наблюдательной скважины), можно провести совместный анализ результатов моделирования обеих зон (рис. 3д, 3е). На представленном далее примере это позволило оценить время начала снижения концентрации контаминанта в наблюдательных скважинах после ликвидации источника с учетом задерживающих свойств отложений ЗА.

В данной работе разработка комплекса моделей начата с НЗ. На основе принятой геофильтрационной схематизации последовательно созданы геологическая, геофильтрационная и геомиграционная модели НЗ в программном комплексе GMS [14]. Комплекс включает в себя инструменты для подготовки исходных данных, разработки геологической модели, а также обработки полученных результатов. Основные модельные расчеты проводятся в GMS с помощью программных кодов MODFLOW-2005 (геофильтрация) и MT3DMS (геомиграция) [9]. Для моделирования влагопереноса и миграции загрязнения в ЗА применена программа HYDRUS-1D [19].

Важно отметить, что разработка большинства моделей имеет итерационный характер, что продемонстрировано в работе [4]. Данный подход чаще применяется для моделей НЗ и позволяет постепенно детализировать отдельные модельные аспекты, в частности связанные со схематизацией фильтрационных неоднородностей [1], получая на каждом этапе более полное их соответствие фактическим данным, но избегая необоснованного усложнения.

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

КРАТКАЯ ХАРАКТЕРИСТИКА РАССМАТРИВАЕМОГО ОБЪЕКТА И ДАННЫХ МОНИТОРИНГА

В качестве объекта, на примере которого рассматривается представленная методика, выбрано солехранилище (СХ), которое располагалось на территории предприятия ФГУП “ГХК” и эксплуатировалось с 1967 по 2007 г. (см. рис. 1). Данный объект представлял собой заглубленное на 6 м сооружение из бетонных плит, облицованных изнутри металлом. Приготовление раствора NaCl происходило непосредственно в СХ. В ходе ликвидации оно было промыто и засыпано грунтом до дневной поверхности.

В 2016 г. на площадке цеха, в 300 м к северо-западу от ликвидированного СХ, была организована наблюдательная сеть (рис. 4). Мониторинг состояния грунтовых вод в ряде наблюдательных скважин выявил устойчивые превышения предельно допустимой концентрации (ПДК) для питьевых вод [6], а также фоновых концентраций по хлорид-иону (350 и 50 мг/л). По данным разовых замеров в скв. 91, не входящей в режимную сеть, концентрации хлорид-иона составляли более 10 ПДК. По схематической карте гидроизогипс, построенной по результатам мониторинга (см. рис. 4), и анализу потенциальных источников хлорид-иона было установлено, что загрязнение поступает от ликвидированного СХ.

Рис. 4.

Схематическая карта гидроизогипс грунтового водоносного комплекса (осредненные данные за период 2016–2019 гг.). 1 – гидроизогипсы с шагом 1 м; 2 – граница площадки цеха; 3 – скважины, использованные при построении карты гидроизогипс, и их номера: а – режимные, б – с разовыми замерами УГВ; 4 – СХ: а – ликвидированное, б – действующее; 5 – направление геофильтрационного потока от ликвидированного СХ; 6 – ручей 2; 7 – золоотвал.

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

ГЕОЛОГИЧЕСКОЕ СТРОЕНИЕ И ГИДРОГЕОЛОГИЧЕСКИЕ УСЛОВИЯ ИЗУЧАЕМОЙ ТЕРРИТОРИИ

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

На площадке цеха распространены четвертичные аллювиальные отложения, в которых выделяется нижняя (русловая) и верхняя (пойменная) фации (см. рис. 2). Русловая фация аллювия представлена гравийно-галечниками и песками. Пойменная фация сложена супесями и суглинками. На участке СХ русловая и пойменная фации не выделяются, а разрез аллювиальных отложений представлен чередованием песков, супесей и суглинков.

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

ОЦЕНКА ИНФИЛЬТРАЦИОННОГО ПИТАНИЯ И ЭВАПОТРАНСПИРАЦИИ

Перед разработкой комплекса математических моделей проведена оценка ИП и ЭТ. По фондовым материалам среднемноголетнее ИП, оцененное по модулю подземного стока ручья 2 (см. рис. 1) при допущении полного дренирования грунтового водоносного комплекса, составляет 105–115 мм/год.

Для оценки ИП применен балансовый метод. ЭТ рассчитывалась с помощью зависимости М.И. Будыко [8]. Испаряемость определялась по методу Шуттлевота-Валласа [17] с помощью программы SurfBal [3]. Поверхностный сток, рассчитанный по методу нумерованных кривых стока [8], составил около 55 мм/год. В результате при ЭТ 345–360 мм/год и среднемноголетнем количестве осадков 490 мм/год за период 1967–2018 гг. (по данным метеостанции Опытное поле, г. Красноярск) ИП составило 75–90 мм/год.

В дальнейшем при моделировании рассмотрено два сценария расчетов: со среднемноголетним ИП равным 75 мм/год (2 × 10^–4 м/сут) и 110 мм/год (3×10^–4 м/сут).

РАЗРАБОТКА ГЕОФИЛЬТРАЦИОННОЙ И ГЕОМИГРАЦИОННОЙ МОДЕЛЕЙ НАСЫЩЕННОЙ ЗОНЫ

Модели НЗ охватывают юго-западную часть территории предприятия, на которой расположено ликвидированное СХ (см. рис. 1а). При создании геологической модели геологическое строение района схематизировано четырьмя модельными слоями. Первые два слоя охватывают четвертичные отложения пойменной и русловой фаций (см. рис. 2). К третьему слою относятся породы коры выветривания, к четвертому – породы архея. Мощность последнего слоя, который на V надпойменной террасе не вскрыт, задана 10 м в соответствии с предполагаемой глубиной трещиноватости архейских гнейсов.

На основе геологической модели разработана геофильтрационная модель, которая, в свою очередь, в процессе вычислительной схематизации разбита на 20 расчетных слоев с мощностью ячеек от 0.1 до 2 м. В плане разделение на модельные ячейки проведено с равномерным пространственным шагом 10×10 м. Режим геофильтрационного потока на модели принят стационарным. Границы модели показаны на рис. 1а.

Калибровка геофильтрационной модели проводилась с целью достижения максимальной сходимости модельных и натурных УГВ. Для этого использованы осредненные УГВ 15 наблюдательных скважин за период 2016–2019 гг., а также разовые замеры уровней в 60 скважинах. Калибровка параметров осуществлялась в автоматическом режиме с помощью модуля PEST [11]. В качестве калибруемых параметров выступали К_Ф отложений четырех модельных слоев, а также коэффициенты, характеризующие фильтрационное сопротивление ложа ручья 2 и р. Енисей. К_Ф отложений каждого модельного слоя задавались одинаковыми в пределах области, где фиксируется хлоридное загрязнение.

Процессы испарения и транспирации моделировались с помощью задания максимальной ЭТ со свободной поверхности УГВ при нулевой глубине их залегания. На модели ЭТ задана линейно уменьшающейся до глубины 3 м (критическая глубина, ЭТ = 0). Максимальная ЭТ, как нечувствительный параметр, зафиксирована на уровне 300 мм/год.

Полученный график сопоставления модельных и натурных уровней показан на рис. 5а для сценария с ИП 75 мм/год. Стандартное отклонение модельных уровней от натурных составило 1.57 м.

Рис. 5.

Схематическая карта модельных гидроизогипс грунтового водоносного комплекса территории и графики сопоставления модельных и натурных УГВ для моделей до (а) и после (б) детализации геологического строения. 1 – модельные гидроизогипсы: а) – основные, с шагом 20 м, б) – дополнительные, с шагом 5 м; 2 – модельное направление геофильтрационного потока от ликвидированного СХ; 3 – граница площадки цеха; 4 – скважины, использованные при построении карты гидроизогипс, и их номера: а) режимные, б) с разовыми замерами УГВ; 5 – модельный стабильный ореол хлорид-иона; 6–7 – границы зон с повышенными фильтрационными свойствами: 6 – гравийно-галечников в русловой фации аллювия, 7 – песков в коре выветривания.

На основе откалиброванной геофильтрационной модели проведены геомиграционные расчеты. Хлорид-ион рассматривался как нейтральный (несорбируемый, стабильный) контаминант. Учитывались его конвективный перенос с фильтрационным потоком и диффузионно-дисперсионные явления. Расчеты проводились без учета плотностной конвекции, влияние которой, по-видимому, имело место только в непосредственной близости от СХ. Массовый поток хлорид-иона задавался в границах проекции источника на УГВ (см. рис. 2) с начала эксплуатации СХ (1967 г.) и калибровался по данным мониторинга.

В процессе геомиграционных расчетов проведена верификация результатов калибровки геофильтрационной модели по данным мониторинга концентраций хлорид-иона в грунтовых водах. Она показала, что направление распространения модельного ореола контаминанта в потоке грунтовых вод не согласуется с данными мониторинга (см. рис. 4, 5а), несмотря на хорошую сходимость модельных и натурных УГВ. В связи с этим возникла необходимость возврата к этапу разработки геологической модели с целью выявления эффективной фильтрационной неоднородности геологической среды [1], определяющей направление распространения ореола загрязнения.

Переинтерпретация данных геологического строения. Дополнительный анализ карты гидроизогипс (см. рис. 4) позволил выявить область с недостаточной детализацией геологического строения, признаком которой была смена направления потока с северо-западного на участке СХ на преимущественно северное на площадке цеха, а также падение его градиента. Такие изменения в структуре потока оказались связаны с наличием двух зон повышенной проницаемости в северо-восточной части площадки цеха (см. рис. 1б, 2), которые были оконтурены по данным бурения. Одна зона, сложенная гравийно-галечниками, выделена в русловой фации аллювия (во 2-м слое геологической модели). Другая зона, представленная гнейсами выветрелыми до песков, выделена в толще слабопроницаемых пород коры выветривания (в 3-м слое) по двум скважинам (5/1ц и 50а). Контуры зоны песков в коре выветривания, а также их мощность подбирались в ходе калибровки направления потока грунтовых вод на участке распространения ореола.

Результаты калибровки геофильтрационной модели после детализации геологического строения представлены на рис. 5б (сценарий с ИП 75 мм/год). Стандартное отклонение модельных уровней от натурных составило 1.55 м, что практически идентично величине, полученной на геофильтрационной модели без детализации, но при этом локальные различия в структуре потока грунтовых вод в двух вариантах являются значимыми по отношению к направлению миграции хлорид-иона (см. рис. 5).

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

Для оценки чувствительности модельных параметров проведены расчеты по нескольким сценариям геофильтрационной и геомиграционной моделей (табл. 2). В качестве параметров рассмотрены ИП, продольная дисперсивность (δ_L) и активная пористость отложений (n_a). Поперечная дисперсивность задана на порядок меньше продольной.

Увеличение ИП с 75 до 110 мм/год вызвало рост К_Ф отложений в связи с необходимостью сохранения сходимости модельных и натурных уровней в скважинах, а также рост массового потока хлорид-иона от источника загрязнения (с 25 до 30 кг/сут). В соответствии с этим увеличился и максимальный массовый поток контаминанта (с 45 до 55 кг/сут), который поступает в зону разгрузки (ручей 2 и р. Енисей) после стабилизации ореола.

Кроме того, увеличение ИП привело к небольшому ускорению миграции загрязнения в грунтовых водах (см. табл. 2, сценарии 1 и 2), с чем связано уменьшение временных интервалов, выделяемых на графиках изменения концентрации контаминанта в скважинах ($t_{1}^{{Н~}}$–$~t_{4}^{Н}$). Для каждого интервала по наблюдательным скважинам даны два значения: для ближайшей к источнику скважины и наиболее удаленной от него.

В табл. 2 также приведены интервалы $t_{1}^{Р}$–$t_{4}^{Р}$, выделяемые на графиках массового потока хлорид-иона, поступающего в поверхностные воды. Пример такого графика для сценария 1 показан на рис. 6. Критерии выделения интервалов определяются экспертным путем. В представленном примере момент стабилизации ореола в грунтовых водах, соответствующий концу интервала $t_{2}^{Р}$, выделен, когда массовый поток контаминанта достигает максимальных значений (Q_{M max}) с относительной погрешностью не более 5%. Окончание деградации ореола (конец интервала $t_{4}^{Р}$) соответствует снижению массового потока хлорид-иона до минимальных значений (Q_{M min}), дающих прирост концентрации контаминанта в ручье 2 на два порядка ниже фона.

Рис. 6.

Модельный график поступления хлорид-иона в поверхностные воды (модель насыщенной зоны, сценарий 1).

Оценка влияния δ_L проводилась для двух сценариев (см. табл. 2). В сценарии 1 параметр задан осредненным для всех отложений (1 м), а в сценарии 3 он варьируется от 0.5 до 10 м в зависимости от типа отложений. Указанный диапазон параметра принят согласно [12] с учетом размеров ореола хлорид-иона от СХ до областей разгрузки, который составляет около 600–1100 м. При этом большие значения δ_L относятся к отложениям с большей глинистостью, для которых характерны более сложные пути миграции загрязнения.

Сравнение двух сценариев показало, что параметр δ_L влияет на результаты незначительно. Для сценария 3 заметно чуть более активное рассеивание ореола хлорид-иона в потоке грунтовых вод, что привело к увеличению на 1–2 года временных интервалов, отвечающих за рост и снижение его концентраций в скважинах ($t_{2}^{Н}$, $t_{4}^{Н}$) и массового потока в зоне разгрузки ($t_{2}^{Р}$, $t_{4}^{Р}$).

Вариации δ_L и ИП не оказывают значительного влияния на откалиброванные модельные концентрации хлорид-иона в наблюдательных скважинах (см. табл. 1).

Для оценки влияния n_a выбраны значения, характерные для соответствующих типов отложений [1]: максимальные (сценарий 1) и минимальные (сценарий 4) (см. табл. 2). Результаты расчетов показали, что изменение n_a существенно влияет на скорость распространения ореола в грунтовых водах (см. табл. 2). При минимальной n_a ореол хлорид-иона быстрее достигает области разгрузки, при ликвидации источника – быстрее деградирует. В результате в первом случае сценарий 4 можно считать консервативным, а во втором – мягким.

Важно отметить, что при продолжении мониторинга в скважинах поочередно будет зафиксирован спад концентраций хлорид-иона за счет деградации ореола, что позволит откалибровать n_a. Уточненные параметры можно будет использовать, как для расчетов распространения загрязнения от нового СХ, расположенного в 100 м к юго-западу от ликвидированного (см. рис. 1), так и при моделировании распространения других видов загрязнения в пределах рассматриваемой области.

РАЗРАБОТКА МОДЕЛЕЙ ВЛАГОПЕРЕНОСА И МИГРАЦИИ ЗАГРЯЗНЕНИЯ В ЗОНЕ АЭРАЦИИ

Наблюдения за содержанием хлорид-иона в грунтовых водах начались только с 2016 г., когда ореол контаминанта уже стабилизировался. Вследствие этого при расчетах принят консервативный сценарий, в котором моделирование миграции хлорид-иона в ЗА при активации источника не проводилось, а временные интервалы $t_{1}^{A}$ и $t_{2}^{A}$ равны нулю (см. рис. 3а).

После ликвидации СХ в 2007 г. из сформировавшейся зоны засоления под первичным источником началось вымывание хлорид-иона атмосферными осадками. Таким образом, для моделей ЗА целью расчетов была оценка времени действия вторичного источника (см. рис. 3б). При этом консервативным сценарием принимался такой, в котором временные интервалы $t_{3}^{A}$ и $t_{4}^{A}$ максимальны.

По данным бурения на участке размещения СХ (см. рис. 2) ЗА представлена чередованием суглинков, супесей и песков. При моделировании для упрощения разрез принят однородным, но рассмотрены сценарии с разным составом отложений. Мощность ЗА принята равной 14 м по результатам геофильтрационного моделирования НЗ.

При расчетах влагопереноса для связи высоты всасывания, объемной влажности пород и коэффициента влагопереноса использованы зависимости Ван Генухтена-Муалема [13, 15]. Для получения параметров уравнения Ван Генухтена, к которым относятся остаточная и полная объемная влажность пород (${{\theta }_{r}}$ и ${{\theta }_{s}}$), а также эмпирические параметры (${{\alpha }}$ и $n$), на основе известных характеристик пород (табл. 3) использована программа Rosetta [16], интегрированная в программу HYDRUS-1D.

На модели влагопереноса на верхней границе задано среднемноголетнее ИП, на нижней – постоянная высота всасывания (0 м), что соответствует УГВ. По результатам моделирования влагопереноса получены стационарные распределения влажности отложений, которые затем использованы для геомиграционных расчетов. Моделирование миграции проводилось с учетом рассеивания вещества за счет диффузии и дисперсии. При расчетах использована относительная концентрация хлорид-иона. Вторичный источник (зона засоления) задан в виде начальных концентраций (1 д.ед.) контаминанта в поровой воде отложений ЗА. Процессы растворения соли при просачивании атмосферных осадков не рассматривались.

При схематизации вторичного источника учтено, что СХ в ходе ликвидации засыпано чистым грунтом. В результате принято, что зона засоления распространяется под ликвидированным СХ на всю глубину ЗА. Принятое допущение обосновано долгим сроком работы СХ (40 лет) и высокими коррозионными свойствами NaCl по отношению к его металлической внутренней облицовке. Кроме того, оно проверено в ходе предварительных расчетов, которые показали, что для засоления всей мощности ЗА за время эксплуатации СХ минимальные фильтрационные потери рассола из него должны были составлять около 0.1 л/сут с каждого квадратного метра площади хранилища.

На геомиграционной модели рассмотрено несколько сценариев (табл. 4), в которых оценено влияние состава отложений ЗА, ее мощности, ИП, К_Ф отложений и δ_L. Влияние состава отложений ЗА учитывалось при изменении параметров уравнения Ван Генухтена согласно табл. 3.

Также принято во внимание, что после ликвидации источника в 2007 г. по данным мониторинга снижения концентраций хлорид-иона в скважинах по состоянию на 2019 г. не зафиксировано. Из этого, согласно рис. 3е, можно записать неравенство

При $t_{3}^{Н}$ для ближайшей наблюдательной скв. 50а равном от 1 до 3 лет (см. табл. 2) можно оценить минимальное значение $t_{3}^{А}$, которое составит 9 лет.

В ходе оценки чувствительности параметров выявлено, что наиболее важным из них является δ_L. Он оценен согласно [12] в диапазоне от 0.1 до 1 м с учетом мощности ЗА, отвечающей длине пути миграции хлорид-иона. Сценарии 1 и 2 при δ_L равном 1 м оказались нереалистичными, поскольку интервал $t_{3}^{А}$ был менее 9 лет (см. табл. 4). При минимальном δ_L (0.1 м) для ЗА, сложенной песками, снижение концентраций хлорид-иона в скважинах должно было уже начаться (см. сценарий 3, табл. 4). На основе этого можно заключить, что в ЗА под ликвидированным СХ преобладают суглинки и супеси, а δ_L близка к 0.1 м.

Следующие сценарии (4–7) показали, что ИП менее чувствительный параметр и может меняться от 75 до 110 мм/год. Возможные вариации К_Ф отложений не оказали существенного влияния на результаты (в табл. 4 не показаны). В качестве консервативного рассмотрен сценарий 8 с суглинистым составом ЗА, ее максимальной мощностью (16.5 м) согласно разовым замерам УГВ в скв. 58/89 (см. рис. 2) и минимальным ИП. Результаты показали, что в таком случае интервал $t_{3}^{А}$ составит 21 год.

Таким образом, по результатам рассмотренных сценариев интервалы $t_{3}^{А}$ и $t_{4}^{А}$ для миграции загрязнения в ЗА оценены в диапазонах 10–21 и 17–27 лет. Окончание интервала $t_{4}^{А}$ (см. рис. 3б) принято по снижению относительной концентрации, поступающей на УГВ, до 0.1% от начальной в зоне засоления.

СОВМЕСТНЫЙ АНАЛИЗ РЕЗУЛЬТАТОВ МОДЕЛИРОВАНИЯ МИГРАЦИИ ХЛОРИД-ИОНА В НАСЫЩЕННОЙ ЗОНЕ И ЗОНЕ АЭРАЦИИ

Для первого этапа миграции хлорид-иона (от активации источника в 1967 г. до стабилизации его ореола в грунтовых водах) задерживающие свойства ЗА не учитывались. При возможной вариативности параметров модели НЗ стабилизация ореола произошла в период с 1975 г. (консервативный сценарий) по 1988 г. (мягкий сценарий). При стабилизации ореола хлорид-ион поступал в поверхностные воды в значениях близких к максимальным (45–55 кг/сут).

Для второго этапа миграции контаминанта (от ликвидации источника в 2007 г. до окончания деградации ореола в грунтовых водах) учитывались результаты расчетов моделей обеих зон. Консервативным в этом случае можно считать сценарий с минимальным ИП, максимальной n_a (сценарий 1, табл. 2) и ЗА максимальной мощности, сложенной суглинками (сценарий 8, табл. 4). В этом случае деградация ореола в грунтовых водах начнется с 2028 г., а снижение концентраций хлорид-иона в наблюдательных скважинах – с 2031–2034 гг. Окончание деградации ореола, соответствующее снижению поступления хлорид-иона в поверхностные воды до минимальных значений (Q_{M min}) (см. рис. 6), произойдет после ликвидации источника через время $t_{3}^{А} + t_{3}^{Р} + t_{4}^{А} + t_{4}^{Р}$, т.е. ориентировочно к 2076 г. По мягкому сценарию деградация ореола происходит с 2017 г. С 2020 г. она проявится началом снижения концентраций в скважинах, а окончательно завершится к 2042 г.

ЗАКЛЮЧЕНИЕ

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

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

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

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