Математическая модель влагопереноса в ландшафтных катенах

УДК 631.6

МАТЕМАТИЧЕСКАЯ МОДЕЛЬ ВЛАГОПЕРЕНОСА

В ЛАНДШАФТНЫХ КАТЕНАХ

,

ФГОУ ВПО МГУП, г. Москва, Россия

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

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

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

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

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

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

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

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

, (1)

где - вертикальное расчленение рельефа, то есть разность отметок поверхности элювиальной и супераквальной фаций; φ - относительная крутизна склона трансэлювиальной фации; , - эмпирические коэффициенты, вместе с относительной крутизной склона они позволяют регулировать площади составляющих катены; B – ширина катены.

Были рассмотрены примеры форм поверхности земли для катен с различным расчленением рельефа Δ0, равным 0,5; 2 и 3 м, соответственно. Для рассмотренных катен принято м; φ = 10; a1 = 1,25; a2 = 2,5.

Естественно, в данном случае пришлось рассматривать потоки влаги истинно двумерными, плоскими в сечении. Поэтому вся рассматриваемая толща разбивалась на элементарные слои переменной толщины, от 0,1 м вблизи поверхности до 1 м вблизи водоупора, . Рассматриваемый пласт разбивался вертикальными плоскостями для образования столбцов и расчетных блоков шириной этих блоков . Центры блоков находились на расстоянии от левой границы – водотока. Площадь блока равна , а его объем .

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

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

,,(2)

где – напор, м на расчетный момент времени ; при отсчете напоров от поверхности земли в самой высокой точке профиля и оси, направленной вниз

, (3)

- напор, м, эквивалентный каркасно-капиллярному давлению в зоне неполного насыщения () и эквивалентный гидростатическому давлению в зоне полного насыщения; – коэффициент влагоемкости, м3в/м4

, (4)

где - объемная влажность почвы, м3в/м3 (м3в – кубический метр почвенной влаги).

При полном влагонасыщении Сw = 0.

Связь между каркасно-капиллярным потенциалом и влажностью почвы принята в виде //

, (5)

где p – пористость; - максимальная гигроскопичность; hk – высота капиллярного подъема, м; m и n – безразмерные эмпирические коэффициенты; - расчетный шаг по времени, сут.; - вертикальное сопротивление потоку влаги между центрами и блоков, сут.

, (6)

-- горизонтальное сопротивление потоку влаги между центрами и блоков, сут

, (7)

Kω – коэффициент влагопроводности м3в/м2/сут, зависящий от объемной влажности почвы w, по [3, 5]

, (8)

при полном влагонасыщении Kω = Kф; Kф – коэффициент фильтрации.

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

Расходование влаги на испарение принято зависящим от погодных условий, от влажности почвы, оно разделялось на испарение с поверхности почвы, которое учитывалось как граничное условие, и на транспирацию, последняя распределялась по корнеобитаемому слою пропорционально влажности почвы и плотности корней и входила в уравнение в виде интенсивности влагоотбора корнями растений из единичного объема почвы ei,j, м3в/м3/сут. С этой целью для каждой декады теплого периода по известной средней температуре воздуха T,оС и относительной влажности воздуха a, % подсчитывалось потенциальное (при оптимальной влагообеспеченности) суммарное испарение (эвапотранспирация) Epot (мм/сут.) по формуле

, (9)

где Kб – биологический коэффициент, учитывающий особенности конкретного ценоза.

Потенциальная эвапотранспирация разделялась на потенциальное испарение с поверхности почвы и потенциальную транспирацию пропорциональной затененности почвы растительным покровом fр, которая изменялась по декадам: = (1-fр) и = fр. Эти потенциальные виды испарения редуцировались на каждом временном шаге:

; ; , (10)

при влажности поверхностного 2-…5 см слоя почвы wп > 0,8p e' = 1; эти зависимости согласуются, например, с исследованиями [1] и др.

, (11) где ; ; ew – коэффициент, учитывающий уменьшение транспирации при отклонении влажности почвы от оптимальной, вид этой зависимости соответствует исследованиям [4]; wk, – средняя влажность корнеобитаемого слоя почвы, переменная во времени и оптимальная в данную декаду, соответственно; ВЗ – влажность завядания.

Скорректированная величина транспирации распределялась по глубине каждого столбца в заданном корнеобитаемом слое пропорционально влажности почвы и массы корней в виде интенсивности влагоотбора корнями растений из единичного объема почвы ei,j, м3в/м3/сут.

Определение напоров почвенной влаги с помощью системы алгебраических уравнений (2) представляет собой громоздкую вычислительную задачу, так как сводится к нахождению порядка 600 неизвестных (при принятой разбивке на блоки) с шагом около 1 суток на протяжении нескольких десятков лет. Следует также отметить существенную нелинейность этой системы уравнений, в которой емкостной коэффициент и проводимость существенно зависят от напоров почвенной влаги, следовательно, и от влажности почвы, что требует 3…7 итераций на каждом временном шаге. Поэтому алгоритм решения этой системы должен быть наиболее эффективным. В настоящее время таковым является метод матричной прогонки [6], который введением вектора напоров по всем i-м столбцам для каждого слоя j позволяет понизить размерность задачи до одномерной

Uj = | ; ; ; . . . | при этом j = 0,1,2,3...Nx. (12)

С помощью этого вектора система уравнений (2) запишется в матричном виде

, (13)

где и - квадратные диагональные матрицы размером , учитывающие вертикальные потоки влаги между и блоками и между и блоками

и т. д., где ;

(14)

и т. д.,

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

и т. д., где , (15)

и т. д., где ;

; (16)

Левое граничное условие, то есть отсутствие потока в центре рассматриваемого пласта, учитывается особыми правилами вычисления элементов этой матрицы D11 и D12; D11 = - D12; . Аналогичное правое граничное условие учитывается при вычислении последних элементов этой матрицы ; .

В случае, если в каком-то слое j имеется источник или сток, они учитываются при вычислении соответствующих элементов матрицы .

Входящий в систему (13) вектор Fj объединяет все свободные члены:

Fj =|; ;

;…|. (17)

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

Решение системы матричных уравнений (13) ищется в виде рекуррентной формулы

, (18)

для этого при прямой прогонке вычисляют матрицы прогоночных коэффициентов PPj и прогоночные векторы-столбцы QQj :

, (19)

. (20)

В этих формулах обозначает обращенную матрицу. Матрицы PPj квадратные, размером , их общее количество равно Nx, при принятой разбивке области фильтрации они содержат элементов. Длина вектора-столбца QQj 21 элемент, всего таких векторов Nx.

При обратной прогонке вычисляют искомый вектор напоров почвенной влаги Uj на конец временного шага по формуле (18).

Верхнее граничное условие учитывают, особым образом вычисляя первую матрицу прогоночных коэффициентов PP0 и первый вектор-столбец QQ0 . Так, если через верхнюю границу (х = 0) нет потока влаги, то , то как это следует из (18), элементы диагонали матрицы PP0 равны , а остальные – нулевые. Все элементы вектора-столбца QQ0 равны нулю. При физическом испарении через поверхность почвы поток влаги равен

.

Поэтому диагональные элементы матрицы PP0 равны , а остальные – нулевые. Элементы вектора-столбца QQ0 равны: . Если поверхность почвы увлажняется поливом, то в приведенных выражениях величина физического испарения заменяется интенсивностью водоподачи (м/сут.), взятой с обратным знаком. Если поливается часть поверхности, то эту замену осуществляют для соответствующих номеров столбцов i. При поливе затоплением, когда на поверхности почвы создается слой воды hув, для соответствующих i , Qi=Hув=hув - Δi.

Условия на нижней границе () реализуются при особом вычислении последнего вектора напоров UNx:

1) при отсутствии потока влаги через нижнюю границу (водоупор) , поэтому в соответствии с (18) напоры вычисляются по формуле

, (21)

EE – квадратная диагональная единичная матрица, то есть у которой элементы Eii =1;

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

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

, (22)

где ddi - вектор-столбец с элементами di;

3) при заданной постоянной во времени глубине грунтовых вод, точном уровне грунтовых вод, отсчитываемом от поверхности земли на возвышенности

; (23)

4) при очень глубоких грунтовых водах, не участвующих в круговороте почвенных вод,

, (24)

где |0,5hNx-1| - вектор-столбец с одинаковыми элементами, равными 0,5hNx-1.

После вычисления вектора по формуле (18) определяют все остальные векторы напоров, включая и U0, то есть получают матрицу напоров почвенной влаги, при этом напоры в крайних левом и правом столбцах приравнивают к соседним, исходя из отсутствия потока влаги на этих границах: ; . Если потоки влаги имеются, напоры вычисляют по тем же формулам, которые закладывают в граничные условия, например, отток в дрену, находящуюся на определенной глубине. Матрицу напоров переводят в матрицу влажностей почвы с помощью выражений (3) и (5). Так как коэффициенты влагопроводности и влагоемкости зависят от заранее неизвестной влажности, то приходится выполнять ряд итераций (около 3…7), пока ответы по влажности не будут отличаться более чем на 0,00001 объема. Знание напоров и сопротивлений позволяет подсчитать потоки влаги в любых сечениях как на границах области, так и внутри нее, например, переток влаги из одной фации в другую, или вертикальные потоки, характеризующие промываемость почвенного слоя. Приведенные модели были реализованы в виде компьютерных программ.

В таблице приведены в качестве примера результаты моделирования передвижения влаги в ландшафтных катенах Московской области за многолетний период (42 года). Статьи водного баланса в таблице представлены за теплый биологически активный период года (со среднесуточной температурой воздуха t>50С). При моделировании водного режима сопряженных фаций учитывали следующие природные и антропогенные факторы: атмосферные осадки, испарение с поверхности почвы, транспирацию растениями, переменную во времени мощность корнеобитаемой зоны, орошение, осушение (дренаж в виде ловчей дрены или систематических дрен). Анализ данных таблицы позволяет сделать следующие выводы при обосновании мелиоративных мероприятий.

На северо-западе Московской области (метеостанция Клин) дефицит увлажнения в среднемноголетнем разрезе практически отсутствует. При слабовыраженном рельефе (с перепадом высоты около 0,5 м) для получения удовлетворительной продуктивности (0,77) требуется довольно интенсивный дренаж даже на элювиальных фациях – в среднем 64 мм с площади этой фации. Орошение на этих фациях практически не нужно. На трансэлювиальных фациях нужно дренирование основания склона, а супераквальные фации нуждаются в систематическом дренаже довольно высокой интенсивности. С увеличением расчленения рельефа (при перепаде высоты 3 м) элювиальные фации практически перестают нуждаться в дренаже (4 мм/год), но зато возникает заметный дефицит влаги, что требует в средний год около 2 поливов, а в сухие годы оросительная норма достигает 200 мм. Из таблицы следует закономерное изменение промываемости почвенного слоя. На элювиальных фациях она возрастает с увеличением расчленения рельефа, а на супераквальных фациях она несколько уменьшается.

На юго-востоке (метеостанция Кашира) увеличивается засушливость климата, что приводит к необходимости других мелиоративных воздействий:

элювиальные фации не нуждаются в осушении, но растет потребность в орошении в 2…3 раза. Примерно во столько же раз снижается также потребность в искусственном дренаже супераквальных фаций;

при слабо расчлененном рельефе промывной водный режим почв сменяется на выпотной (капиллярное подпитывание со стороны грунтовых вод). При более расчле-ненном рельефе в условиях орошения промываемость почв на элювиальных фациях может быть значительной (в среднем 186 мм/год), но без орошения она существенно меньше (85 мм/год). Из-за лучшей теплообеспеченности юго-востока Московской области относительная продуктивность земель на 10…15 % выше.

Водный режим ландшафтных катен Московской области

Статьи вод-	Перепад высоты 0,5 м	Перепад высоты 2 м	Перепад высоты 3 м
ного балан- са (мм) и продуктивность	Ни-зина	Склон	Возвы- шен-ность	Ни-зина	Склон	Возвы- шен-ность	Ни-зина	Склон	Возвы- шен-ность
1	2	3	4	5	6	7	8	9	10
Метеостанция Клин
Осадки 385
Испаряемость 377
Дефицит увлажнения 8
Дренажный сток:
в дрены на возвышенности			64			15			4
в ловчую дрену		108			188			234
в дрены в низине	53			104			141
в реку	71			80			85
Оросительная норма 50%			7			30			64
Промывае-мость 0,7м слоя	94	90	119	86	99	157	83	106	199
Относительная урожайность	0,84	0,74	0,77	0,78	0,68	0,91	0,75	0,65	0,95
Метеостанция Кашира
Осадки 330 Испаряемость 456
Испаряемость 456
Дефицит увлажнения 126
Дренажный сток:
в дрены на возвышенности			2			0			0
Продолжение табл.
1	2	3	4	5	6	7	8	9	10
в ловчую дрену		24			93			156
в дрены в низине	4			31			64
в реку	28			49			58
Оросительная норма 50%			52			109			157
Промываемость 0,7 м слоя	-2	-1	63	-11	12	132	-11	28	186
Относительная урожайность	0,85	0,85	0,94	0,91	0,86	0,96	0,88	0,82	0,95

Примечание. В таблице супераквальная фация обозначена низиной, трансэлювиальная и трансаккумулятивная – склоном, элювиальная – возвышенностью; их протяженность по склону, соответственно, 350, 100 и 350 м.

Библиографический список

1. Будаговский почвенной влаги. М.: Наука, 19с.

2. , , Сухарев . М.: КолосС, 2005, 216 с.

3. , Новиков модель переноса влаги и растворов солей в почвах на орошаемых землях. //Тр. МГМИ, М. 1974. 6 с.

4. Константинов в природе. Л.: Гидрометеоиздат, 1968, 532 с. 2-е изд.

5. , Голованов почвоведение. М.: Колос, 19с.

6. , Гулин методы. М.: Наука, 1989.