CC BY
5
УДК 556.3.04; 556.3.07
doi: 10.55959/MSU0579-9406-4-2023-63-3-77-87
МОДЕЛЬНАЯ ВЕРИФИКАЦИЯ ДАННЫХ МОНИТОРИНГА УРОВНЕЙ ГРУНТОВЫХ ВОД
Елена Владимировна Кононченко11, Александр Борисович Петраш2, Сергей Олегович Гриневский3
1 ФГБУ «Гидроспецгеология», Москва, Россия; lena.konon@mail.ruH
2 ФГБУ «Гидроспецгеология», Москва, Россия; A.B.Petrash@yandex.ru
3 Московский государственный университет имени М.В. Ломоносова, Москва, Россия; sogrin@geol.msu.ru
Аннотация. Рассмотрен пример нестандартного использования математического моделирования, когда данные режима уровней грунтовых вод выступают в качестве объекта верификации. Для рассматриваемой площадки предприятия, где размещены источники загрязнения грунтовых вод, разработан комплекс моделей, описывающих процессы трансформации осадков на поверхности земли, влагопереноса в зоне аэрации и геофильтрации. На основе модельных расчетов обоснована корреляция между натурными значениями амплитуды сезонных колебаний уровней и годовой инфильтрацией, а также выявлены некорректные данные мониторинга.
Ключевые слова: математическое моделирование, зона аэрации, режим уровней подземных вод, ин-фильтрационное питание подземных вод
Для цитирования: Кононченко Е.В., Петраш А.Б., Гриневский С.О. Модельная верификация данных мониторинга уровней грунтовых вод // Вестн. Моск. ун-та. Сер. 4. Геология. 2023. № 3. С. 77-87.
MODEL VERIFICATION OF GROUNDWATER TABLE MONITORING DATA
1 l^l 2 3
Elena V. Kononchenko , Aleksandr B. Petrash , Sergey O. Grinevskiy
1 «Hydrospetzgeologiya», Moscow, Russia; lena.konon@mail.ru1
2 «Hydrospetzgeologiya», Moscow, Russia; A.B.Petrash@yandex.ru
3 Lomonosov Moscow State University, Moscow, Russia; sogrin@geol.msu.ru
Abstract. An example of unconventional modelling using when groundwater table regime data is considered as a verification object has been presented in the article. A models set has been developed for the site where sources of pollution are situated. The models set describes such processes as precipitation transformation on the earth surface, unsaturated flow in vadose zone and groundwater flow. Based on modelling results correlation between seasonal variation amplitudes of groundwater levels measured at monitoring wells and annual groundwater recharge has been proved, inconsistent monitoring data have been revealed.
Key words: mathematical modelling, vadose zone, groundwater table regime, groundwater recharge
For citation: Kononchenko E.V., Petrash A.B., Grinevskiy S.O. Model verification of groundwater table monitoring data. Moscow University Geol. Bull. 2023; 3: 77-87. (In Russ.).
Введение. Математическое моделирование в настоящее время широко используется для решения различных задач, связанных с проведением гидродинамических расчетов и оценкой воздействия источников загрязнения на подземные воды. Разработка математических моделей геофильтрации и миграции загрязнения в подземных водах требует сбора и анализа разнородных фактических данных о природно-техногенных условиях изучаемой территории. Особое место среди них занимают результаты мониторинга состояния подземных вод, к которым относятся режим уровней, показатели радиохимического и химического состава и др. Указанные данные используются в качестве основы при разработке и калибровке моделей, с ними сравниваются модель-
ные уровни и концентрация отдельных компонентов [Веселов и др., 2004].
В связи с этим актуальность приобретает верификация самих исходных данных, как на этапе их комплексного анализа, так и в процессе создания математических моделей. Проверка данных позволяет отделить те, которые можно использовать в качестве опорных при моделировании, от сомнительных, которые в дальнейшем либо отсеиваются, либо используются как вспомогательные, с более низкой достоверностью.
В работе в качестве данных для верификации нами выбраны режимные наблюдения за уровнями грунтовых вод (УГВ) в сети скважин на площадке предприятия, где размещены объекты — потен-
Рис. 1. Схема расположения площадки предприятия (А) и карта гидроизогипс грунтового водоносного горизонта уфимских отложений (осредненные данные за 2010-2012 гг.) (Б): 1 — изолинии рельефа с шагом 5 м; 2 — водоемы и площадные водотоки; 3 — временные водотоки; 4-5 — границы области моделирования: 4 — внешние граничные условия III рода по контурам ВГ и ГД (временные водотоки); 5 — внешние граничные условия II рода по контурам АБ, БВ и АД (непроницаемые границы); 6 — предполагаемая зона разгрузки потока грунтовых вод в виде рассредоточенных высачиваний и родников, граничное условие III рода; 7 — границы площадки; 8 — линия гидрогеологического разреза; 9 — гидроизогипсы с шагом 0,5 м; 10-11 — скважины и их номер:
10 — наблюдательные, 11 — инженерно-геологические
циальные источники радиоактивного загрязнения грунтовых вод. Эти объекты относятся к категории приповерхностных источников загрязнения, поскольку они расположены на поверхности земли или заглублены на несколько метров ниже рельефа.
Основная цель работы — обоснование достоверности данных режима УГВ. Для этого использованы специально разработанные математические модели трансформации осадков на поверхности земли, влагопереноса в зоне аэрации (ЗА) и геофильтрации. Дополнительно проведена модельная оценка параметров, необходимых для дальнейших расчетов миграции загрязнения.
Геологическое строение и гидрогеологические условия изучаемой территории. Рассматриваемый район находится на юго-восточном склоне Восточно-Европейской платформы, где ее осадочный чехол разделен на две части гипсами и ангидритами кунгурского яруса нижней перми, которые служат региональным водоупором. Верхняя часть геологического разреза сложена карбонатно-терригенны-ми отложениями уфимского яруса нижней перми, а также терригенными четвертичными отложениями (рис. 1; 2).
Площадка предприятия расположена на коренном берегу долины крупной реки. Уфимские отло-
жения, к которым приурочен горизонт грунтовых вод, представлены толщей глинистых пород с многочисленными прослоями выветрелых песчаников, трещиноватых известняков и мергелей мощностью до 1-1,5 м. Для уфимской толщи характерны неравномерная трещиноватость, частые фациальные замещения и невыдержанность прослоев по мощности и простиранию. Существенная гетерогенность пород обусловливает неоднородность их фильтра-ционно-емкостных свойств.
Материалы и методы исследований. Анализ режима уровней грунтовых вод по данным мониторинга. На изучаемой территории сеть мониторинга состояния грунтовых вод расположена в северной части площадки, где находятся потенциальные источники загрязнения. Мониторинг ведется с 2010 г. в 7-ми скважинах глубиной 20 м (рис. 1; 2). Замеры УГВ проводятся в среднем с частотой раз в месяц, во время половодья и паводков — раз в декаду или чаще. Режим УГВ за период наблюдений 2010-2019 гг. представлен на рис. 3, А и Б на примере двух скважин. Режим УГВ проанализирован для всех наблюдательных скважин. Весенний подъем уровней фиксируется с конца марта — середины апреля и длится до мая. Затем в целом наблюдается спад до следующей весны, на фоне которого осенью
Рис. 2. Схематический гидрогеологический разрез по линии 1-1': 1 — техногенные отложения голоценового возраста; 2 — верхнечетвертичные и голоценовые отложения (пески, галечники, суглинки); 3 — нижнечетвертичные делювиальные отложения (глины, суглинки); 4-5 — отложения пермского возраста: 4 — уфимского яруса (глины, суглинки с прослоями песчаников, известняков, мергелей); 5 — кунгурского яруса (гипсы с прослоями глин и известняков); 6 — трещиноватые породы; 7 — нисходящий родник; 8 — УГВ; 9 — наблюдательная скважина, снесенная на разрез, с указанием положения фильтра, вверху — номер скважины; внизу — ее глубина, м; справа — абсолютная отметка УГВ, м; 10 — инженерно-геологическая скважина
периодически фиксируется временная стабилизация уровней, реже — их небольшой рост.
За время наблюдений амплитуды сезонных колебаний и среднегодовые уровни (СГУ) существенно менялись (рис. 3, А, Б). Первые три года (2010-2012) амплитуды были минимальны и составляли 1-2,5 м, а СГУ был практически одинаков. В следующие два года (2013-2014) произошел рост амплитуд до 2-4 м, а СГУ возросли на 0,5-1 м по сравнению с первым периодом. В течение 2015-2017 гг. амплитуды возросли до 3-6 м, а затем за период 2018-2019 гг. стали снижаться. При этом СГУ также выросли на 0,5-1 м с последующим снижением до прежних значений.
Между режимом УГВ и величиной годовых осадков намечается некоторая зависимость. В годы с количеством осадков ниже нормы (555 мм за период 2010-2019 гг.) СГУ и амплитуды колебаний были обычно меньше по сравнению с таковыми в годы, когда осадков было больше нормы (рис. 3, А, Б). Вследствие этого можно предполагать, что СГУ и амплитуды зависят от вариаций инфильтрацион-ного питания (ИП) в течение периода наблюдений. Однако для выяснения вопроса о том, является ли ИП основным режимообразующим фактором, необходима его оценка с учетом внутригодовых и многолетних вариаций.
Помимо изменчивости ИП на режим УГВ, вероятно, влияли и другие факторы. Различия амплитуд
сезонных колебаний в скважинах могут быть в определенной мере обусловлены фильтрационно-емкост-ной неоднородностью водовмещающих пород. Так, например, ускорение подъема уровня и увеличение амплитуды колебаний будет происходить в случае, когда он достигнет зоны с меньшей гравитационной емкостью по сравнению с емкостью нижележащей зоны.
Необходимо отметить практически мгновенную реакцию уровней в скважинах на начало весеннего половодья при значительной мощности ЗА (1318 м). Указанная закономерность может быть связана с тем, что в условиях, близких к полному насыщению, при наличии макропор и неравномерной трещино-ватости пород происходит ускорение влагопереноса в ЗА. Кроме того, быстрому весеннему повышению уровней, вероятно, способствует конструкция скважин, затрубное пространство которых по всей длине заполнено гравием. В результате в него может поступать вода из систем трещин и макропор ЗА, что дополнительно увеличивает амплитуду колебаний.
На изменение УГВ в скважинах могут влиять и факторы искусственного происхождения, такие, как техногенные источники питания подземных вод, однако по имеющимся данным в северной части площадки предприятия таковые отсутствуют.
При анализе режима УГВ обращает на себя внимание их скачок, пришедшийся на рубеж 2014-
Рис. 3. Графики сопоставления модельных и натурных УГВ по наблюдательным скважинам 4 и 5 (А и Б) и изменения модельной инфильтрации на основе моделей водного баланса на поверхности земли и в зоне аэрации (В и Г) за период режимных наблюдений
2015 гг. Разница между уровнями в конце 2014 г. и начале 2015 г. составляет 1,5-2,75 м в зависимости от скважины. Характерно, что после скачка уровней во всех скважинах далее фиксировалось обычное для зимней межени стабильное их положение (скв. 4), либо спад (скв. 5) вплоть до весеннего подъема. Анализ температуры воздуха по данным ближайшей метеостанции г. Уфа, расположенной в 20 км от
площадки предприятия (рис. 4), показал отсутствие положительных температурных аномалий, которые могли бы вызвать активное таяние снега и увеличение ИП в зимнюю межень 2014-2015 гг. В связи с этим причина скачка уровней вместе с ростом амплитуд сезонных колебаний, начиная с 2013 г., возможно, заключается в технической неисправности скважин, обусловленной разгерметизацией устьев
2010 2011 2012 2013 2014 2015 2016 2017 2018 2019
месячные суммы осадков, мм/месяц средняя температура воздуха, °С скользящее среднее значений месячных осадков - скользящее среднее температур воздуха
Рис. 4. Графики изменения средних температур воздуха и месячных сумм осадков по метеостанции г. Уфа за период режимных наблюдений. Месячные суммы осадков с октября по март (осадки в виде снега) показаны голубым, с мая по август (осадки в виде дождя) — темно-синим, сентябрь и апрель (переходный период) — светло-синим; для средних скользящих осреднение проведено по 3-м (для месячных сумм осадков) и 29-ти (для температур воздуха) значениям
и поступлением воды с поверхности по затрубному пространству. Также нельзя исключать некорректность данных об УГВ, связанную с ошибками при измерениях и/или обработке.
Таким образом, комплексный анализ фактического материала показал, что особенности режима УГВ, фиксируемого наблюдательными скважинами, могут быть обусловлены несколькими различными факторами, к основному из которых, предположительно, относятся вариации естественного ИП.
Оценка инфильтрационного питания на основе моделей водного баланса на поверхности земли и в зоне аэрации. Для определения степени влияния ИП на режим УГВ разработаны модели водного баланса на поверхности земли и в ЗА, с помощью которых оценено среднемноголетнее ИП, а также его сезонные и годовые вариации за период режимных наблюдений.
Оценка среднемноголетнего ИП состояла из двух этапов. На первом этапе проведено моделирование трансформации осадков на поверхности земли в программе 8иг£Ба1 [Гриневский, 2017] с учетом их задержки растительностью, накопления в виде снега и снеготаяния, испарения с поверхности растительности и снега, промерзания и оттаивания почвы, стокообразования. В качестве исходных данных использованы метеоданные с суточным разрешением по ближайшей метеостанции за период 1970-2019 гг. (осадки, минимальная, максимальная и средняя
температуры воздуха, относительная влажность воздуха, средняя скорость ветра).
Поверхностный сток (ПС) в период положительных температур рассчитан с использованием метода нумерованных кривых стока [USDA, 1985], который основан на обобщении результатов многочисленных эмпирических наблюдений на малых водосборах и связывает осадки и ПС с учетом состава почв и ландшафтных условий. В зимний период ПС определен, исходя из проницаемости и впитывающей способности мерзлой почвы. В соответствии с преимущественно полевым типом ландшафта территории, травянистым растительным покровом и суглинистым составом почв номер кривой поверхностного стока (НКПС) варьировал от 80 до 95, а индекс развития листьев — от 0 до 3 [Гриневский, 2010]. Потенциальная эвапотранспирация определена по методу Шуттлевота-Валласа ^ЬиШе'^гШ, 1985]. Параметры модели, характеризующие таяние и уплотнение снега, откалиброваны по данным о его среднемноголетней мощности, рассчитанной для периода 1970-2019 гг. Мощность снега в зимнее время в среднем составляла 30 см, достигая 45 см к концу февраля.
Результатами первого этапа расчетов стали суточные значения поступления влаги на поверхность почвы, потенциального испарения с поверхности почвы и потенциальной транспирации корнями растений. Они использованы на втором этапе модели-
рования влагопереноса в ЗА в программе Hydrus-1D [Simtoek, 2013].
Отметим, что модели влагопереноса разработаны на основе весьма ограниченных данных о характеристиках пород ЗА изучаемой площадки, что вызвало необходимость упрощения схематизации и минимизации количества модельных параметров. В результате влагоперенос принят одномерным и рассмотрен для вертикальной модельной колонки мощностью 16 м, что равно средней мощности ЗА в северной части площадки. Расчеты проведены для периода 1990-2019 гг.: первые 20 лет использованы для получения начальных условий по высоте всасывания и влажности отложений модельной колонки, следующие 10 лет, соответствующие периоду режимных наблюдений за УГВ, — для оценки ИП.
Строение ЗА схематизировано двумя квазиоднородными слоями (рис. 2). Первый слой представлен нижнечетвертичными делювиальными суглинками и глинами мощностью 5 м, второй — глинами уфимского яруса мощностью 11 м. С целью упрощения схематизации почвенный слой не учитывался, как и гетерогенность пород.
Для описания связи высоты всасывания, объемной влажности отложений и коэффициента влагопереноса использованы зависимости Ван Генух-тена-Муалема [Mualem, 1976; Van Genuchten, 1980]. Входящие в них параметры (табл. 1), к которым относятся остаточная и полная объемная влажность пород (0r и 9S), их коэффициент фильтрации при полном водонасыщении (к), а также эмпирические параметры (а и n), для каждого слоя определены с помощью программы Rosetta [Schaap et al., 2001]. В последней используются педотрансферные функции, представляющие собой эмпирические зависимости указанных выше параметров от свойств отложений, таких, как гранулометрический состав, плотность и др.
Данные о плотности и полевой влагоемкости отложений обоих слоев взяты по результатам лабо-
раторных измерений. Гранулометрический состав отложений задан приближенно исходя из их лито-логического состава. Содержание глинистых частиц принято равным 35%, что соответствует глинам по классификации В.В. Охотина или грунту «clay loam» по классификации, принятой в США [Грунтоведение, 2005]. Содержание песчаных и алевритовых частиц принято одинаковым (32,5%). Для второго слоя модели коэффициент фильтрации принят 0,05 м/сут по результатам предварительной калибровки геофильтрационной модели (ГФМ) при стационарном режиме уровней и значении ИП 110 мм/год.
На верхней границе модели задано поступление влаги в ЗА, на нижней — постоянная высота всасывания (0 м), что соответствует УГВ. Для описания зависимости транспирационной способности растений от высоты всасывания использована модель Ван Генухтена [Simtoek, 2013]. Значение высоты всасывания, при котором интенсивность транспира-ции падает в 2 раза по сравнению с потенциальной величиной, рассчитано по методике, описанной в [Гриневский, 2011], и составляет 72 м. В качестве модели распределения корней по глубине, от плотности которых зависит интенсивность поглощения ими влаги, использована модель Хофмана-Ван Генухтена [Simtoek, 2013]. В этой модели основная плотность корней сосредоточена в приповерхностном слое почвы и резко снижается с глубиной, что характерно для травянистой растительности. Модель определяется одним параметром — мощностью распространения корневой системы, которая для суглинистых почв площадки принята равной 0,3 м.
Расчеты показали, что наиболее чувствительный параметр моделей трансформации осадков на поверхности земли и влагопереноса в ЗА — НКПС. В связи с этим рассмотрено несколько сценариев с различными НКПС в диапазоне от минимального до максимального (табл. 2). Совместный анализ водного баланса двух этапов моделирования за 20102019 гг. показал, что при увеличении НКПС с 80 до
Таблица 1
Тип пород зоны аэрации Делювиальные суглинки и глины Глины уфимского яруса
Гранулометрический состав, %; размер частиц, мм пески (>0,05) 32,5*
алевриты (0,002-0,05) 32,5*
глины (<0,002) 35*
Плотность сухого грунта, г/см3 1,54 1,71
Полевая влагоемкость, % 25 19
Параметры уравнения Ван Генухтена Влажность, % Остаточная, 0r 7,3 7
Полная, 9S 39 33
Коэффициент фильтрации при полном водонасыщении (к), м/сут 0,16 0,05
Эмпирические параметры 0, 1/м 2,48 4
n, (-) 1,34 1,39
Характеристики пород зоны аэрации и параметры уравнения Ван Генухтена
Примечание: * — гранулометрический состав принят на основе литологического состава.
Рис. 5. График зависимости среднемноголет-ней инфильтрации от глубины залегания УГВ по сценарию 1.2 (номер кривой стока 85)
Таблица 2
Составляющие водного баланса на поверхности Земли и в зоне аэрации по различным сценариям (осреднение за 2010-2019 гг.)
Элемент баланса, мм/год Сценарий расчета (номер кривой стока)
1.1 (80) 1.2 (85) 1.3 (90) 1.4 (95)
Атмосферные осадки 555
Потенциальная эва-потранспирация (испаряемость) по методу Шуттлевота-Валласа 701
Испарение с поверхности листьев и снега 13
Транспирация корнями растений 141 136 133 119
Испарение из почвы 203 206 199 189
Эвапотранспирация (суммарное испарение) 357 355 345 321
Поверхностный сток 73 86 107 158
Инфильтрация 120 109 97 71
Дебаланс 5 5 6 5
95 на значительную величину (85 мм/год) возрастает ПС, что составляет 15% от количества осадков. При этом суммарное испарение снижается на 35 мм/год. Среднемноголетнее ИП, таким образом, варьирует от 71 до 120 мм/год в зависимости от выбранного НКПС. По имеющимся фондовым материалам для соседнего района, расположенного в сходных ландшафтных условиях, ИП, оцененное по колебаниям УГВ, составляет 65-110 мм/год при принятой для расчетов гравитационной емкости уфимских пород (ц), равной 0,025-0,05, что не противоречит полученным модельным оценкам.
Дополнительно определена зависимость ИП от глубины залегания УГВ для различных сценариев на
основе моделирования влагопереноса при разной мощности ЗА. В качестве примера на рис. 5 приведен график для сценария с кривой стока 85. Расчеты показали, что среднемноголетнее ИП достигает максимальных значений при минимальной глубине УГВ, равной 4 м. При уменьшении глубины залегания УГВ инфильтрация сначала снижается, затем на глубине около 1,5-2 м нисходящий поток влаги сменяется восходящим, а инфильтрация — испарением. По результатам расчетов значение эвапотранспирации (ЭТ) со свободной поверхности УГВ при нулевой глубине его залегания составляет от 223 до 310 мм/ год в зависимости от принятого сценария. В балансовом отношении эта величина характеризует дополнительный подток влаги с зеркала грунтовых вод, в результате чего суммарное испарение становится равным потенциальному (701 мм/год).
Для оценки сезонных и годовых вариаций ИП за период режимных наблюдений в качестве входных данных использованы результаты моделирования трансформации осадков на поверхности земли по сценариям с различными НКПС. Основные сезонные изменения влажности пород ЗА, возникающие за счет неравномерного поступления влаги на верхнюю границу модели, а также испарения и транс-пирации влаги корнями растений, происходят в ее верхней части. Ниже расположена подзона транзита, где влажность пород близка к их полевой влагоемко-сти и отвечает среднемноголетнему ИП [Мироненко, 1998]. В связи с этим мощность модельной колонки сокращена до 4 м, что соответствует глубине, на которой модельное ИП становится максимальным (рис. 5). Поскольку значимые для водного баланса процессы происходят выше этой глубины, ИП, поступающее через нижнюю границу короткой модельной колонки, может быть перенесено на нижнюю границу ЗА, отвечающую положению УГВ [Гавич, 1988]. Данный подход с учетом принятых упрощений в модельной схематизации позволяет
Таблица 3
Модельная годовая инфильтрация (мм/год) по различным сценариям
Год Модель водного баланса на поверхности земли и в зоне аэрации, сценарий расчета (номер кривой стока) Модель геофильтрации, сценарий расчета
1.1 (80) 1.2 (85) 1.3 (90) 1.4 (95) 2.1 2.2
2010 32,5 37,7 40,4 39,3 89,1 123,8
2011 99,4 75 87,3 19,8 91,5 142,6
2012 54,2 58 62,1 63,2 84,4 105
2013 203,4 149,4 107,8 31,3 116,2 163,8
2014 197,3 166,2 165,5 123,4 96,4 119
2015 113,7 108,3 102,6 85,5 97,6 128,1
2016 164,8 167,1 156,3 139,1 110,2 156,2
2017 193,2 180,5 84,5 45,7 122 176,9
2018 101,8 101,7 96,4 80 85,4 97,6
2019 35,7 46,4 64,7 82,7 85,4 125,7
среднее 119,6 109 96,8 71,0 97,8 133,9
избежать отставания по времени начала модельной инфильтрации от весеннего подъема УГВ (рис. 3).
Результаты расчетов для каждого года по различным сценариям (табл. 3) показали, что уменьшение НКПС закономерно ведет к снижению ПС и росту ИП. В некоторых случаях (2010 и 2019 г.) отмечается противоположная зависимость, которая, вероятно, связана с особенностями расчета ПС в зимний период, когда проницаемость мерзлой почвы зависит от содержания в ней льда.
При сравнении ИП по отдельным годам для разных сценариев видно, что оно отличается в основном на 10-50 мм. Исключение составляют два года (2013 и 2017 г.), для которых уменьшение НКПС вызвало увеличение годового ИП в 4-6 раз. Сравнение графиков изменения модельного ИП с режимом УГВ показывает, что сценарии с меньшим НКПС (рис. 3, В) лучше согласуются с относительно высокими амплитудами сезонных колебаний, характерными для двух указанных лет. При рассмотрении всего периода режимных наблюдений можно сделать вывод, что максимальные значения коэффициента корреляция между натурными амплитудами колебаний и модельным годовым ИП получены для сценария 1.2 с кривой стока 85 (табл. 4), что позволяет принять его в качестве наиболее реалистичного. Для сцена-
рия 1.4 (кривая стока 95) корреляция практически отсутствует.
Таким образом, в качестве возможного диапазона изменения ИП при разработке ГФМ принят интервал от 97 до 120 мм/год, соответствующий сценариям с НКПС от 90 до 80. Результаты моделирования показали, что изменения амплитуды колебаний УГВ за период режимных наблюдений в значительной мере обусловлены вариациями ИП. Другие факторы (гетерогенность пород, техническое состояние скважин), вероятно, имеют подчиненное значение.
Обращает на себя внимание тот факт, что по модельным оценкам практически каждый год до начала весеннего снеготаянья в зимнее время имело место ИП, доходящее до 15-20% от годового, а зимой 2015/2016 гг. — до 40% (рис. 3, В, Г). Этот результат не согласуется с режимом УГВ. В частности, зимой 2015/2016 гг., наоборот, фиксировался спад уровней. В то же время в конце осени — начале зимы 2015 г. практически во всех скважинах наблюдался их увеличение. Можно предположить, что модельное ИП для зимнего периода отстает по времени от натурного и должно приходиться на осень, а не на зимнюю межень. Полученное несоответствие, по-видимому, связано с принятыми упрощениями в модельной
Таблица 4
Коэффициенты корреляции между натурными амплитудами сезонных колебаний УГВ и модельной годовой инфильтрацией по различным сценариям модели водного баланса на поверхности земли и в зоне аэрации
Сценарий расчета (номер кривой стока) Коэффициент корреляции
Скв. 1 Скв. 2 Скв. 3 Скв. 4 Скв. 5 Скв. 6 Скв. 7 Средний
1.1 (80) 0,58 0,64 0,58 0,71 0,73 0,56 0,66 0,64
1.2 (85) 0,73 0,78 0,73 0,80 0,82 0,68 0,78 0,76
1.3 (90) 0,59 0,66 0,47 0,48 0,47 0,33 0,45 0,49
1.4 (95) 0,53 0,59 0,41 0,29 0,20 0,28 0,30 0,37
схематизации, это требует проведения отдельной работы по созданию моделей, основанных на более обширных фактических материалах и учитывающих гетерогенность пород ЗА.
Разработка моделей геофильтрации при стационарном и нестационарном положении уровня грунтовых вод. Разработка геофильтрационных моделей проведена с целью воспроизведения динамики УГВ в течение периода наблюдений, анализа соответствия натурных и модельных УГВ для различных временных интервалов, а также оценки вариативности ИП и других модельных параметров, необходимых для дальнейших расчетов миграции загрязнения. Для моделирования использован программный комплекс GMS [GMS..., 2022], модельные расчеты выполнены с помощью программного кода MODFLOW-2005 [Chiang Wen-Hsing, 2005].
Моделирование геофильтрации состояло из двух этапов. На первом этапе разработана ГФМ при стационарном положении УГВ, на основе которой на втором этапе создана нестационарная ГФМ. Границы области моделирования показаны на рис. 1. Южная граница задана непроницаемой по зоне выклинивания водовмещающих пород на коренном склоне долины реки (рис. 2). Нижняя непроницаемая граница модели принята по кровле регионального кунгурского водоупора.
Водовмещающие уфимские породы характеризуются невыдержанностью более проницаемых прослоев в основной глинистой толще, что препятствует их корреляции в межскважинном пространстве. В результате уфимские породы рассмотрены на модели в качестве квазиоднородного модельного слоя, который разбит на 10 расчетных слоев с мощностью модельных ячеек от 1 до 10 м и размером в плане 10x10м.
Калибровка стационарной ГФМ проведена на основе осредненных УГВ за период режимных наблюдений в виде серии расчетов при значениях ИП, изменявшихся в диапазоне 97-120 мм/год. С ее помощью оценен диапазон вариативности значений коэффициента фильтрации (Кф) уфимских пород, который составляет 0,03-0,055 м/сут. Указанный диапазон Кф далее использован в качестве ориентира на втором этапе моделирования при разработке нестационарной ГФМ.
Типичный график сопоставления модельных и натурных УГВ для стационарной ГФМ (рис. 6) показывает, что наилучшая сходимость получена для двух скважин (скв. 4 и 5). Вследствие этого для облегчения процесса калибровки ГФМ при нестационарном положении УГВ выбраны данные режимных наблюдений за период 2010-2019 гг. в указанных двух скважинах (450 замеров). Калибровка осуществлена в автоматическом режиме с помощью модуля PEST [Doherty, 2004], а на завершающих этапах — ручным подбором. Основными калибруемыми параметрами были значения ИП, Кф и ц.
Параметр ц калибровали в диапазоне 0,01-0,1 с учетом того, что его ориентировочные значения для основной глинистой части пород составляют 0,005-0,05 [Гавич, 1988], а для прослоев песчаников и известняков они могут быть существенно выше.
ИП калибровали с учетом его сезонных и годовых вариаций. Для оптимизации расчетов ИП задано постоянным в пределах трех временных интервалов, которые выделены по результатам анализа режима УГВ. Для интервала с марта-апреля по май, на который приходится основной весенний подъем уровней, ИП принято максимальным, для периода с июня по сентябрь — относительно невысоким. Для интервала с октября по март предположено, что ИП отсутствует. При калибровке учтено, что среднемноголетнее ИП за период режимных наблюдений не должно существенно отличаться от ранее определенного диапазона (97-120 мм/год).
Процессы испарения и транспирации моделировали с помощью задания значения ЭТ со свободной поверхности УГВ при нулевой глубине их залегания. ЭТ линейно уменьшалась от максимальных значений у поверхности земли до нуля на глубине, равной 4 м. Указанное значение ЭТ задано меняющимся в течение года. Для интервала с апреля по сентябрь оно принято постоянным, для остального времени — нулевым. При этом учтено, что среднегодовое значение ЭТ должно соответствовать 235 мм/год (рис. 5).
Период расчета для нестационарной ГФМ выбран с 2009 по 2019 г. Начало расчета с 2009 г. связано с необходимостью задать начальные уровни, от положения которых существенно зависят результирующие модельные уровни и амплитуда их сезонных колебаний.
Калибровка нестационарной ГФМ показала, что положение СГУ в основном определяется соотношением годового ИП и Кф пород, а амплитуда сезонных колебаний уровней — величиной ИП за весенний период и ц. В ходе калибровки рассмотрена сходимость модельных и натурных УГВ, а также амплитуды их колебаний для разных интервалов наблюдений. Воспроизвести удовлетворительную сходимость УГВ для всего периода наблюдений одновременно не удалось из-за их скачка на рубеже 2014/2015 гг.
Отдельно рассмотрен вариант увеличения ИП с 2015 г., для которого проведена калибровка ГФМ по натурным УГВ за период 2015-2019 гг. В этом случае удовлетворительная сходимость достигнута при увеличении среднемноголетнего ИП до 250-300 мм/год, что в 2,5 раза выше определенных ранее значений. При этом модельные значения амплитуды сезонных колебаний выросли, а максимальные модельные уровни в весеннее время стали близки к натурным. Однако для зимней межени характерен значительный спад модельных уровней, что противоречит натурным данным.
Таким образом, для окончательной калибровки параметров модели использованы натурные УГВ
124 125 126 127 128 129
Наблюденные уровни, м
Рис. 6. График сопоставления модельных и натурных УГВ по наблюдательным скважинам для стационарной геофильтрационной модели
и амплитуда их колебания за 2010-2014 гг. Калибровка годового ИП за 2015-2019 гг. проведена только на основе натурных значений амплитуды колебания УГВ. По итогам калибровки получено два сценария расчетов: с минимальным (98 мм/год) и максимальным (134 мм/год) среднемноголетним ИП (табл. 3). Для них воспроизведена динамика УГВ в течение года, а также получены модельные значения амплитуды колебания уровней, максимально близкие к натурным (рис. 3, А, Б). Откалиброванный диапазон Кф уфимской толщи составляет 0,037-0,05 м/сут, для ц — 0,06-0,08. Сценарию с минимальным ИП соответствуют минимальные значения параметров, и наоборот.
Отметим, что полученные на модели значения Кф пород отражают осредненные фильтрационные свойства уфимской толщи. При этом прослои песчаников, известняков и мергелей, вероятно, характеризуются более высокими значениями Кф, а основная глинистая часть толщи — более низкими. Модельные значения ц, вероятно, относятся к емкости прослоев (поровой — для песчаников и трещинной — для известняков и мергелей), а также трещинной емкости глин. При этом основная глинистая масса (за исключением трещин), по всей видимости, содержит в разной степени связанную воду, т. е. занимаемая ею емкость не осушается и не заполняется при сезонных колебаниях УГВ.
Для сравнения с полученными модельными значениями для пород уфимской толщи ц дополнительно оценена по данным режима УГВ с использованием модельного среднемноголетнего ИП (табл. 3).
Оценка проведена по методике, описанной в работе [Шестаков, 2003], для периодов, когда, предположительно, происходил независимый спад УГВ. Для периодов после весеннего половодья (рис. 3, А, Б), в которые скорость спада уровней была максимальна, результаты оценки дают величину ц от 0,004 до 0,015-0,025, что существенно ниже модельных значений. Можно предположить, что эта величина соответствует по большей части трещинной емкости пород, которая быстро заполняется весной и также быстро осушается после прекращения обильного ИП. Для периодов зимней межени, когда происходит относительно медленное снижение уровней, получены значения ц от 0,035 до 0,12, которые, вероятно, отражают суммарную гравитационную емкость пород, как трещинную, так и поровую. Таким образом, можно предположить, что полученные по результатам калибровки значения ц (0,06-0,08) соответствуют осредненной гравитационной емкости пород уфимской толщи.
Скорость миграции загрязнения в трещинах и наиболее проницаемых прослоях может быть существенно выше, чем в остальной части пород. В связи с этим откалиброванное значение ц может служить в качестве консервативной (минимальной) оценки активной пористости водовмещающих пород при дальнейших геомиграционных расчетах.
Заключение. Показан пример нестандартного применения математического моделирования, при котором данные о режиме уровней грунтовых вод выступают в качестве объекта верификации. Для этих целей специально разработан комплекс
моделей, описывающих процессы трансформации осадков на поверхности земли, влагопереноса в зоне аэрации и геофильтрации.
По результатам модельных расчетов проведена проверка корректности данных мониторинга, которая показала, что абсолютные отметки уровней грунтовых вод за 2010-2014 гг. можно использовать в качестве опорных при дальнейшем моделировании миграции загрязнения. В то же время данные за 2015-2019 гг., вероятно, некорректны и должны быть исключены из модельных расчетов.
На основе комплексного анализа фактического материала и проведенного моделирования обо-
снована корреляция между значениями амплитуды сезонных колебаний уровней грунтовых вод, фиксируемых в скважинах, и естественной годовой инфильтрацией. Поэтому значения амплитуды сезонных колебаний за весь период наблюдений, в том числе и за период 2015-2019 гг., могут быть использованы для дальнейших модельных расчетов.
На этом примере также продемонстрирована возможность получения новой информации об объекте исследований посредством разработки комплекса математических моделей для случая, когда объем фактических данных весьма ограничен.
СПИСОК ЛИТЕРАТУРЫ
1. Веселое В.В., Паничкин В.Ю., Мирошниченко О.Л. Решение обратных задач геофильтрации на математических моделях гидрогеологических объектов. Геология Казахстана // Сб., посвященный XXXII сессии междунар. геол. конгресса (Флоренция, Италия, 2004). Алматы: МОН РК, 2004. С. 417-426.
2. Гавич И.К. Гидрогеодинамика. М.: Недра, 1988. 349 с.
3. Гриневский С.О. Моделирование поглощения влаги корнями растений при расчетах влагопереноса в зоне аэрации и инфильтрационного питания подземных вод // Вестн. Моск. ун-та. Сер. 4. Геология. 2011. № 3. С. 39-50.
4. Гриневский С.О., Поздняков С.П. Принципы региональной оценки инфильтрационного питания подземных вод на основе геогидрологических моделей // Водн. ресурсы. 2010. № 5. С. 543-557.
5. Гриневский С.О., Поздняков С.П. Ретроспективный анализ влияния климатических изменений на формирование ресурсов подземных вод // Вестн. Моск. ун-та. Сер. 4. Геология. 2017. № 2. С. 42-50.
6. Грунтоведение / Под ред. В.Т. Трофимова. M.: Изд-во Моск. ун-та, 2005. 1024 с.
7. Мироненко В.А., Румынин В.Г. Проблемы гидроэкологии: В 3 т. Т. 1. Теоретическое изучение и моделирование геомиграционных процессов. М.: Изд-во МГГУ 1998. 611 с.
8. Шестаков В.М., Поздняков С.П. Геогидрология. М.: Изд. ИКЦ «Академкнига», 2003. 176 с.
9. Chiang Wen-Hsing. 3D-groundwater modeling with PMWIN: a simulation system for modeling groundwater flow
and transport processes. Berlin, Heidelberg: Springer-Verlag, 2005. 397 р.
10. Doherty J. PEST. Model-Independent Parameter Estimation. User Manual: 5th Edit. Australia, Watermark Computing. 2004. 336 р.
11. Genuchten M.T. van. A closed form equation for predicting the hydraulic conductivity of unsaturated soils // Soil Sci. Soc. Amer. J., 1980. Vol. 44. P. 892-898.
12. GMS: User Manuals. URL: https://www.aquaveo.com (дата обращения: 12.03.2023).
13. Mualem Y. A new model for predicting the hydraulic conductivity of unsaturated porous media // Water Res. Res. 1976. Vol. 12. P. 513-522.
14. Schaap M.G., Leij F.L., van Genuchten M.T. Rosetta. A computer program for estimating soil hydraulic parameters with hierarchical pedotransfer functions // J. Hydrology. 2001. Vol. 251. P. 163-176.
15. Shuttleworth W.J., Wallace J.S. Evaporation from sparse crops-an energy combination theory // Quart. J. Royal Meteorol. Soc. 1985. Vol. 3. P. 839-855.
16. Simimek J., Sejna M., van Genuchten M.T. et al. The HYDRUS-1D software package for simulating the one-dimensional movement of water, heat, and multiple solutes in variably-saturated media. Ver. 3.0. // Prepr. Depart. of Environ. Sci. University of California Riverside, California, Riverside, 2013. 340 р.
17. USDA, Soil Conservation Service. National engineering handbook. Section 4. Hydrology. US Government Printing Office. Washington, D.C. 1985.
Статья поступила в редакцию 20.03.2023, одобрена после рецензирования 22.03.2023, принята к публикации 22.07.2023
