I. Основные уравнения

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

I) в почвенном растворе –

2) в почвенно-поглощающем комплексе -;

З) в кристаллической форме возможно существование гипса

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

В одномерной постановке уравнения модели имеют следующий вид:

Уравнение /1/ описывает одномерное движение влаги в зовах неполного и полного водонасыщения почвогрунта. Уравнение /2/ определяет скорость влагопереноса. Уравнения /З/ - /7/ списывают конвективно-диффузионный перенос рассматриваемых ионов с учетов ионообменной сорбции и растворения - кристаллизации солей. Уравне­ний /8/, /9/ - изотермы катионного обмена Б. П. Никольского. Условия /10/ задает постоянство емкости поглощения во времени. Уравнение /11/ описывает кинетику растворения -кристаллизации гипса. Здесь приняты следующие обозначения:

H - обобщенный потенциал почвенной влаги, м;

                                                H = P + x,

P - капиллярный потенциал, м;

X - вертикальная координата (x = 0 на поверхности почвы, положительное направление вниз);

W -  объемная влажность, дм/дм³;

- капиллярная влагоемкость, м;

k(W) - коэффициент влагопроводности, м/сут;

e - функция отбора, влаги корнями растений, сут;

V - скорость влагопереноса, м/сут;                  

- концентрации  в почвенном растворе, мг.экв/л;

- содержание  в ППК, мг.экв/100 г. почвы;

- концентрация кристаллического гипса в почве, мг.экв/100г почвы;

 - концентрация гипотетической соли  в  почвенном растворе, мг.экв/л (в программе предусмотрена связка ионов в гипотетические соли);

  - концентрация предельного насыщения раствора гипсом, мг.экв/л. Значение   .вычисляется в программе по эмпирической формуле Ю.Н. Чиркина в зависимости от содер­жания в растворе хлоридов и сульфатов ( - содержание в растворе гипотетических солей

):

 

D* - коэффициент конвективной диффузии ионов, м2/сут;

g - объемная масса скелета почвы, г/см³;

- константы изотерм ионного обмена Na - Ca и Ca – Mg

  - сумма обменных оснований, мг.экв/100 г почвы;

- коэффициенты скорости растворения гипса, сут;

    t  - время, сут.

2. Начальные и краевые условия

Начальные условия задают исходное распределение потенциа­лов влаги /или влажности/, концентрации ионов   и гипса по координате

где:

 -начальные значения потенциала (м) или влажности (дм3/дм3);

 - начальное содержание растворенных ионов  (мг.экв/л или %)

 - исходное содержание гипса в почвогрунте  (мг.экв/100 г или %).

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

Для уравнения влагопереноса на верхней границе предусмотре­но задание двух видов краевых условий:

1)  условие 1-го рода (поверхностный полив)

 

 

 

где:

 - слой воды на поверхности почвы (м) при поливе затопле­нием и по полосам; при поливе по бороздам следует принимать  = -0,2 м.

2) условие 2-го рода (осадки, испарение, полив дождеванием)

где:

 - величина потока через поверхность (м/сут/)

В программе WAS052  величина   на заданный интер­вал времени    вычисляется по формуле

где   - соответственно осадки, поливная норма и суммарное испарение, транспирация, (м3/га) за период  (сут) Если  <0  , то испарение превышает осадки, при  > 0 осадки и поливная норма превышают испарение.

Для уравнения солепереноса на верхней: границе условно  3-го рода Данквертса - Бреннера

 - концентрация соответствующих ионов в поливной воде (мг.экв) при поливе; в межполивной период = 0.

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

1) условие 1-го рода - потенциал                 

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

 2) условие 2-го рода - поток

где  - заданный поток через нижнюю границу; < 0 отток за пределы области расчета (переток через нижележащие слои, боковое растекание);

> 0  - приток в области расчета (напорное питание, приток с соседней территории/:

= 0 - отсутствие потока (водоупор или равенство протека и оттока);

З) условие 3 - го рода - потек, зависящий от H;  в выражении (21)

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

4) условие 2-го рода - свободное стекание влаги

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

т.е. полагается, что диффузионный поток ионов мал.

3. Параметры массопереноса

3.I. Параметры уравнения влагопереноса

Для уравнения влагопереноса необходимо задавать два основных параметра: зависимости капиллярного потенциала от влажности

(основная гидрофизическая характеристика – ОГХ) и влагопроводности от влажности.       

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

P(W)  должна быть аппроксимирована функцией.

где Wп - полная влагоемкость (меньше значения пористости на величину содержания защемленного воздуха); - нижний предел иссушения почвогрунта (максимальная  гигроскопичность);  - условная высота капиллярного поднятия. Параметры ,  должны находиться на основе аппроксимации экспери­ментальных кривых P(W). Если таковых не имеется, то априори можно принять  - максимальная гигроскопичность;

=  где  - высота капиллярного поднятия; = 1,2…(чем тяжелее грунт, тем больше )

Экспериментальная зависимость влагопроводности K(W) должна быть   аппроксимирована формулой С.Ф. Аверьянова

где     - коэффициент фильтрации,   m - пористость. Значение  принимается тем же, что и в зависимоcти P(W) . Пара­метр  находится путем аппроксимации экспериментальных кри­вых K(W) (рис.16).

Для расчета в программу закладываются значения  m, Wп, ,,,. Если экспериментально определенных гидрофизи­ческих характеристик не имеется, то для ориентировочных расчетов можно использовать данные табл. 1.

Таблица 1

Механичес­кий состав	M                Wп
Песчаный	0.41-0,44	0.37-0.40	0.05-0,10	0,15-0,17	0.5-1.2	3,5-4.5
Супесчаный	0.44-0,48	0,40-0.44	0.08-0,13	0,19-0.24	0,7-1.5	4,0-5,С
Легкоcу-
глинистый	0,48-0,52	0,44-0.48	0.10-0,15	0,24-0.30	1.5-2,5	4.5-6,5
Среднесуг-
линистый	0.44-0,48	0.40-0.44	0.13-0.18	0,30-0,35      2.5-3,5		5.5-9.0
Тяжелосуг- линиотыв и
глинистый	0,47-0.54	0.43-0.49   | 0,I6-0,25		0,35-0,39    | 4.0-7.0		6,0-12,0

 

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

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

где Wз - влажность завядания (Wз~1.4Wмг), Wппв – предельно-полевая влагоемкость;f(x) - функция распределения корневой системы /рис. 2/

, x,  - параметры распределения, определяемые в программе по коду культуры; fкор - мощность корне обитаемой зоны.

В качестве исходных данных вводятся значения Wз, Wппв, fтр (fтр – дол транспирации от суммарного испарения).

 

3.2. Параметры солепереноса

 

Для расчета задаются следующие параметры солепереноса: доля «сквозных» пор - ; приведенный параметр гидродисперсии - ; константы изотерм катионного обмена , коэффициент скорости растворения .

Если расчет ведется при полном водонасыщении, т.е. уравнения (1), (2)  в расчетах не участвуют, то доля «сквозных» пор определятся как  =  /m- активная пористость, m - полная пористость. Если уравнения (1) и (2) в расчетах участвуют, то = 1. Параметр гидродисперсии (приведенный)    и "активная" пористость  определяются по экспериментальным данным промывки монолитов для каждого литологического слоя (по иону ) Коэффициент конвективной диффузии в программе вычисляется по формуле

где V - скорость фильтрации. Для ориентировочных расчетов можно определять    по графику  приведенному на рис.3. Значение , используемое в программе, связано с     соотношением = /.

Koэфициенты изотерм катионного обмена рассчитывают на основе экспериментальных данных совместного определения катионов  в почвенном растворе и ППК. Навески почвы заливаются растворами солей (обычно  или ) различной концент­рации (1,5,10 г/л) и выдерживаются 1-2 сут до установления катионного равновесия. Затем раствор отфильтровывается и выполняется анализ на содержание в нем и ППК почвы ионов . По результатам эксперимента для каждого горизонта составляется таб­лица, в которую записывают значения (мг.экв/л), (мг.экв/100г),  Полученные значения наносятся на графики: на одном по оси абсцисс откладывают значения  , а по оси ординат  другом - соответственно и  . Подученные обла­сти точек, аппроксимируются прямыми линиями. По тангенсу угла на­клона прямой на первом графике определяется  , а на втором

Примечание: Значение , определяемое для расчета по

программе WAS052, меньше в  ~ 3,16 раз значения  определяемого для расчета по программе WAS051. Это связано о различием в размерности концентраций перового раствора (мг.экв/л – в WAS052, мг.экв/100 мл – в WAS051). Вели­чина  не изменяется.

Для чернозема ориентировочные значения коэффициентов изотерм катионного обмена составляют: ~0.04-0.08; ~1-3.

Величина коэффициента интенсивности растворения гипса  варьирует в довольно широких пределах (0,01...10 сут). Она должна подбираться по данным промывки монолитов, содержащих гипс, путем сопоставления расчетных и экспериментальных выходных кривых ионов  

                  4. Численная реализация модели на ЭВМ

Система уравнений (1)-(12) решается численным методом ко­нечных разностей. Область расчета [0,L] разбивается на узлы равномерной или неравномерной сеткой /рис.4а/. Задается число узлов сетки   Kx    и их координаты . Первы узел рас­полагается на поверхности почвы= 0, последний - на ниж­ней границе  = L. Внутри области расчета строение почвогрунта может быть неоднородным. Задается число литологических слоев  Ks    и верхняя граница каждого слоя , тогда

рис. 4б. Для каждого литологического слоя задаются свои гидрофизи­ческие и гидрохимические характеристики. Исходные значения потенциалов или влажности и концентрации солей задаются в узлах (рис. 4в), число которых  Kj   может не совпадать с Kz. По числу узлов задаются координаты Hj, в которых в начальный момент времени заданы исходные эпюры. Программа их аппроксимирует по координате x - ломаной (потенциалы, влажности и концентрации Cj) и ступенчатой () функциями.

Для прогнозирования влагосолепереноса каждый расчетный год разбивается esl отдельные периоды (задается число периодов в году

 и их продолжительность , i = 1,2,…    (с определенными значениями испарения   осадков , поливных норм    и нижнего краевого условия). Каждый период програм­ма разбивает неравномерной сеткой по времени и производит расчет. За начальное распределение компонент для каждого периода принима­ется распределение, полученное на конец предыдущего периода. В программе предусмотрено несколько режимов расчета, определяемых переменными MR, MH, MC. Уравнение влагопереноса может решаться без уравнения солепереноса (при MH = 1, MC = 0), и наоборот (при MH = 0, MC = 1), а также совместное их решение (при MH = 1, MC = 1).

Схема ввода исходной информации приведена на рис.5, вводимые параметры приведены ниже - построчно.

Форматы ввода:

1) целые переменные KX, KS, KJ, MR, MHW, MPF, G2, MD, MC, MGP, MIZ, MAK,MCP, ND, NM, NPEP, KE вводятся в формате 12 (целое число, заснимающее две позиции), NP, KP, KCP  - в формате 11(целое число, занимающее одну позицию), NG - в формате 14 (целое число, занимающее четыре позиции);

2) остальные переменные вводятся в формате E6.2 - вещест­венное, число со знаком или без него, с десятичной точкой (обяза­тельно), с порядком или без него, занимающее 6 позиций, например: +0.,100, -13.57,   243. ,  -1.2-4 (последняя запись означает -0.000I2).

Общие исходные данные:

1-строка; текстовая информация о решаемой задаче  (шифр, наименование

объекта и др.), всего не более 80 символов, включая пробелы;

2 -строка; целые переменные, определяющие число узлов сетки и режимы

расчета

KX - число узлов разностной сетки по координате (Kx), Kx <50;

KS - число литологических слоев в области расчета (Ks), Ks <50;

KJ - число узлов, в которых заданы исходные эпюры (Kj), Kj <50;

MR - признак, определяющий режим расчета может принимать следующие значения:

MR = +1 ("водно-солевой" режим) - расчет перераспреде­ления компонент на ряд заданных интервалов с изменяющейся интенсивностью потоков на верхней границе;

MR = +2 ("промывка") - расчет ведется с шагом T0 по времени до тех пор, пока средняя концентрация солей в заданном слое /НКОР/ не станет меньше допустимой /СДОП/;

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

MHW - признак, определяющий вводимую эпюру для уравнения влагопереноса, может принимать 2 значения:

MHW = +1 - в качестве исходных вводятся значения потенциала Н, м;        MHW = +2 - в качестве исходных вводятся значения объемной влажности;  MPF = 0 - рабочий параметр;

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

G2 = +1 - условия 1-го рода на нижней границе, модели­руется заданный потенциал;

G2 = +2 – условие 2-го рода на нижней границе, модели­руется заданный поток;

G2 = +3 - условие 3-го рода на нижней границе, моделируется отток при работе горизонтального дренажа (при этом MD > 0);

G2 = +4 - условие 3-го рода, моделируется свободное гравитационное стекание влаги при глубоких грунтовых водах.

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

MD = 0 - дренажа нет, параметры не задаются;

MD = 1 - имеется горизонтальный дренаж в однородном пласте, подстилаемом водоупором  (схема сводится к однослойной);

MD = 2 - имеется горизонтальный дренаж в слоистой толще (схема сводится к двухслойной, дрены в верхнем пласте, нижний

пласт хорошо проницаемый.

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

MH - признак решения уравнения -влагопереноса, может принимать следующие значения:

  MH = 0 - уравнение влагопереноса в расчете не участвует, расчет солевого режима ведется при полном водонасыщении и постоянной по профилю скорости фильтрации;

  MH = 1 - уравнение влагопереноса в расчете участвует.

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

  MC = 0 - солевой режим не рассчитывается; уравнение солепереноса в   расчетах не участвует;

  MC = 1 - солевой режим рассчитывается;

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

MGP = 0   - исходные концентрации солей вводятся b мг/экв;

  MGP = 1 - исходные концентрации солей вводятся в %

MIZ = 1,  MAK = 0 - рабочие параметры.

MCP -  число задаваемых вариантов качества оросительной воды  MCP < 6).

  3 – строка: координаты узлов разностной сетки (Xi, м) и признаки печа­ти результатов в узлах (NPi) . Если  Npi = 0   , то вывод в узле на печать не производится, если Npi = 1  - производится вывод на печать в соответствующем узле Xi, за которым стоит значение Npi,  i = 1,…Kx

4 – строка – HJ(j), j = 1,…Ky - координаты узлов, в которых определены исходные эпюры, м; при этом HJ(1) = 0,  HJ(KJ) = L. В одной строке до 12 значений HJ

5 – “строка”   - заполняется KS таких “строк”, водно-физическими характеристиками профиля по слоям, т.е. сколько слов - столько строк;

HS - верхняя граница каждого литологического слоя в области расчета [0,L], м; HS(1) =0

ПРО - пористость слоя /^/, доли I;

ПВ -  полная влагоемкость слоя (Wп) , м3/М3;

   ППВ - предельно-полевая влагоемкость слоя (Wппв), м3/м3

ИЗ - влажность завядания слоя (Wз), м3/м3;

 

МГ - максимальная гигроскопичность слоя /или влажность, при

которой прекращается движение влаги - ), дм3/дм³;

КФ - коэффициент .фильтрации слоя ( , м/сут;

HK* - "условная" высота капиллярного поднятия (), м;

NKW - показатель степени, в формуле влагопроводности С.Ф. Аверьянова ();

G - объемная масса скелета (g), г/см3

NU - резервный параметр;

ТЕКСТ - буквенное наименование слоя.

6  - ”строка” - вводится, только при MC = 1, заполняетcя  KS  строк гидрохимическими характеристиками профиля по слоям:

КАППА - доля, "сквозных" пор  в общей пористости ( = 1 при MH =1);

LMD - коэффициент гидродисперсии солей (), м;

KNA-CA - константа изотермы Na - Ca  ();

KCA – MG константа изотермы Ca – Mg ();

ППК - сумма обменных оснований () слоя, мг.экв/100 г почвы;

 -.коэффициент скорости растворения гипса, сут,

7 – строка: вводится только при MH = 1 исходная -эпюра для уравнения  влагопереноса:

при MHW = 1 вводится значения потенциала H (метрах) в узлах HJ(j);

при MHW = 2 вводятся значения объемной влажности (дм3/дм3) в узлах HJ(j), j = 1,…,KJ.

8 -12 строки - вводится только при MC = 1,  исходные эпюры для    уравнений солепереноса (содержание ионов в почвенном растворе),  C1(j), C2(j), C3(j), C4(j), C5(j) - соответственно концентрации ионов   почвенном растворе в точках HJ(j) в начальный момент времени (вводятся в мг.экв/л при или % при MCP=1) содержание      в ППК рассчитывается программой автоматически по значениям ;

ТЗ - строка вводится только при MC = 1.  N15(j),  j = 1,…KJ - исходное содержание гипса в твердой фазе в точках HJ(j) в начальный момент времени, мг.экв/100г почвы при MGP =1 ИЛИ % при MGP =

I6 строка - параметры горизонтального дренажа, вводятся только при      MD = 1 или 2 и МН = 1.

 B - междренное расстояние, м;

 DDR - эффективный диаметр дрены, м;

 HDR - глубина заложения дренажа, м;

    TB - мощность слоя от поверхности земли до водоупора при одно­словной геофильтрационной схеме (MD = 1)  или до подстилаю­щего хорошо проницаемого слоя при двухслойно схеме (MD = 2, коэффициент фильтрации подстилающего слоя в 10 раз и более выше коэффициента фильтрации верхнего слоя), м;

 КФ - средний коэффициент фильтрации слоя ТВ (считается, что дре­на заложена в этом слое, т.е. ТВ > HDR ), м/сут;

 PH - проводимость подстилающего слоя (при MD = 2; PH = 0 при MD = 1), м2/сут.

17 «строка» - вводится только при MC = 1 (заполняется MCP «строк»).  - соответственно концентрации   в поливной воде, мг.экв/л.

 18 строка- HD, NM, NG - соответственно календарные число, месяц и год - дата, с которой начинается счет.

 19 строка - имеет 2 различные формы.

                    При MD = 2 ("промывка"):

 19 строка: G1 - признак способа промывки

     (G1 = 1 -       затоплением, G1 = 2  - дождеванием);

КCР - номер варианта качества воды, подаваемой на промывку;

                     из таблицы 17 строки;

T0 - интервал времени, через который периодически производится   выдача результатов о ходе промывки на печать, сут.;

G1 - слой воды на поверхности, м (при G1 = 1) или интенсив­ность дождевания, м/сут (при G1 = 2);  

  H2 - значения нижнего краевого условия;

при G1 = 2  H2 - потенциал на нижней границе, м;

  при G2 = 2  H2 - поток на нижней границе ( H2 > 0 - приток,

H2 < 0  - отток,  H2 = 0 - водоупор или равенство притока и оттока, м/сут);

при G2 =3  Н2 - дополнительные потоки на нижней границе, влияющие на УГВ при работе дренажа (фильтрационные потери из каналов, напорные питание - со знаком плюс; перетек че­рез разделяющие слои, отток - со знаком минус м/сут;

при G2 = 4; H2 = 0;

     СDОП - средняя концентрация растворенных солей (за вычетом гипса),

                 до которой надо опреснять заданный слой,  %;

   HKOR - глубина опресняемого слоя, м.

                При MR = 2    19 строка - последняя.

При MR = 1   (водно-солевой режим) вводятся:

19 строка: NПЕР - число периодов в году NПЕР < 50;

                  KE- код сельскохозяйственной культуры, по которому вычисляется распределение корневой системы (KE = 0   - корневой сис­темы нет, 1 - кукуруза, 2 - подсолнечник, 3 - озимая пшеница, 4 - яровая пшеница, 5 - свекла, 6 - люцерна, 7 - картофель, 8 - томаты, 9 - яблоня, 10 - виноград, 11 - хлопчатник).

 

20 «строка» - заполняется NПЕР «строк» при MR = 1   по одной на каж­дый расчетный интервал времени года здесь:

КР - признак печати (0 - печати эпюр на конец интервала времени нет; 1 - печатаются эпюры в узлах Xi, где NP =1);

     KСР - номер варианта качества оросительной воды: 0 - если полива нет (PNOR = 0 ; 0 < KCP £ MCP, если полив есть (PNOR > 0);

          T0 - расчетный интервал времени, сут.;

       ESM - суммарное испарение за время Т0, м3/га;

        FTR - доля транспирации от ECM  (доля 1);

          OS - осадки за период T0   за вычетом поверхностного стока, м3/га;

PNOR - поливная норма за период   T0, м3/га;

HKOR -глубина корневой системы, м;

    H2 - нижнее краевое условие (см. 19 строку при MR =2).

 

После ввода 20-ой «строки»  в количестве NПЕР штук вводится 19 «строка» на следующий год расчета, затем снова группа «строк»  20 и т.д., пока в 19 строке NПЕР > 0 . Если в 19 строке NПЕР = 0, то программа заканчивает работу.

Вывод информации на печать

Все вводимые переменные выводятся на печать.

При MR = 1 на каждый расчетный период времени выводится его номер, дата на конец периода, время от начала периода ( TV ,сут), транспирация за время TV(TRAN, м3/га), время суммарное от нача­льной даты расчетного года (STV, сут. ), транспирация за это вре­мя (STRAN , м3/га), границы зон полного и неполного водонасыщения (УГВ или верховодки, точки, где Р = 0).

При значении KP = 1 на конец каждого интервала выдаются массивы в узлах сетки, где NP = 1: координата - x, м; по­тенциал H, м; капиллярный потенциал - Р, м; влажность - W, дм3/дм3; средняя влажность в слое [0,x] – WSR,  дм3/дм3; скорость на момент времени TV  - V , м/сут; поток влаги через сечение -x   за время TV – QTV, мм; поток влаги через сечение x   за время STV  - QSTV, мм; влагозапасы в слое [0,x] – SW, мм; изменение влагозапасов в слое [0,x]  по сравнению с предыдущим периодом - DM, мм; CA, MG, NA, CL, SO4, HCO3 - концентрации соответствующих ионов в почвенном растворе, мг.экв/л и %, минерализация почвенного раствора C, г/л и %; натриевое адсорбционное отношение SAR = NA/ концентрация гипса в твердой фазе – CASO4, мг.экв/100 г и %; концентрация кальцита в твердой фазе -  CACO3, мг.экв/100. г и % (в программе WAS052 концентрация кальцита равна нулю); сумма солей в твердой фазе - N , %; общая сумма солей твердой и жидкой фазы-

C + N, %; содержание в ППК – CA, MG, NA,  мг.экв/100 г; обменный натриевый процент ESP = . Значения по составу фаз солей выводятся на печать только при MC = 1.

При MR = 2 (промывка) печать эпюр идет о заданным интервалом времени T0 ; в конце расчета печатается промывная норма нетто ( Q,м3/га) и эпюры влажности и засоления во всех узлах сетки.

Последняя запись, выводимая на печать: КОНЕЦ РАСЧЕТА.

Состав программы

Программа WAS052 требует 160 Кб памяти. Ввод исходной информации по файлу с номером 5; вывод на печать по файлу с номе­ром 6. Состоит из головной программы WAS052 и 31 подпрограммы:

1. SETX 61  - ввод исходной сетки - x, вычисление расчет­ных шагов. ввод координат - HJ;

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

3. ZASLR1 - подпрограмма аппроксимации исходных эпюр;

4. ISEP62 - ввод исходных эпюр и формирование рабочих массивов;

5. GRAND2 - ввод параметров дренажа и расчет оттока в дренаж;

6. VGUS61 - ввод в формирование краевых условий;

7. VFIL01    - вычисление скоростей влагопереноса и фильт­рации;

8. POTN5    - вычисление P(W), W(P), s(P);    

         9  PKW01  - вычисление влагопроводности;

      10. TRAN01  - вычисление источников отбора влаги корнями;

      11. FKOR01 - выбор функции корневой системы;

12.COFC52- решение уравнения солепереноса;

13. RESH01- решение уравнения влагопереноса;

14. YPROG1- метод прогонка;

      15. ZONA01 - расчет границ верховодок и УГВ;

      16. ZAPW61 - расчет влагозапасов;

      17. WSR01 - расчет средней влажности;

      18. DAT01  - расчет даты по номеру дня;

      19. CONDT1 - контроль исходной сетки;

      20. WRIT61 - печать результатов;

 21 OBNL7- пересылка данных в рабочие массивы;

      22. YMER52  -  перевод концентраций из мг.экв/л в % и обратно;

      23. PPK61 - расчет состава ППК;

 

      24. SALG52  - связка ионов в гипотетические соли;

      25. CSUM52 - расчет средних и суммарных концентраций;

      26. MATDR   - матричная прогонка;

      27. RGORD - обращение матрицы;

       28. ESAR61 - вычисление ESP и SAR;

       29. AIJC61  - расчет вспомогательных коэффициентов для

уравнений солепереноса;

        30. OSAD61 - расчет твердей фазы солей;

       31. ALF61 - расчет коэффициентов активности.

Расчет на один год (365 сут.) с разбивкой на 10 узлов по глубине и 15 интервалов времени в году занимает 5-6 мин. машинного времени.

Программные прерывания

1. Если при MR = 2 получится STV > 1000, то печатаются эти сообщения и текст "НЕ ИДЕТ ПРОМЫВКА" и программа заканчивает выполнение; подобная ситуация возникает при промывках, когда рассоление не достигается либо в результате неверно заданных преде­лов, расселения, недостаточной интенсивности дождевания, либо в результате отсутствия отточности при промывке.

2. Если значения X или HJ или HS следуют не в порядке возрастания печатается соответствующее сообщение с ука­занием параметра и номера значения и программа прекращает выпол­нение.

3. Если значения ПВ > ПОР, или ППВ > ПВ, или ВЗ > ППВ, или МГ > ВЗ, печатается сообщение и программа прекращает выполнение.

4. Если при задании параметров дренажа ТВ < ВДР, или КФ £ 0 или ДDР = 0. или при МД = 2 РН £ 0, то печатается сооб­щение и программа прекращает выполнение; если при G2 = 3 не за­даны параметры дренажа (МД = 0), печатается сообщение и выполнение прекращается.

5. Если при задании в качестве начальных значений влажности W £ МГ для соответствующего слоя, печатается сообщение и выполнение прекращается»

5. Рекомендации по проведению расчетов .на ЭВМ

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

При близком залегании УГВ и наличии дренажа расчетная глуби­на L принимается большей или равной глубине заложения дренажа, на нижней границе задается условие 3-го рода (G2 = 3). Если дренаж отсутствует, что задается условие 2-го рода (21) (G2 =2). Значение потока  (идентификатор H2 в 20 строке) опреде­ляется из гидрогеологических условий: например, интенсивность напорного питания, перетока через нижележащие слои или извест­ных боковых притоков и оттоков.

При глубоком залегании УГВ, если нас интересует вся зона аэрации или ставится задача прогноза подъема УТВ, расчетная глубина принимается до водоупора или немного ниже существующего УГВ, и

на нижней границе ставится условие 2-го рода (G2 = 2). При глу­боких грунтовых водах, если подъем их до нижней границы за расчет­ное время не происходит, можно задавать краевое условие (22), мо­делирующее свободное гравитационное отекание (G2 = 4, вода, про­шедшая вниз за границу x = L , обратно не возвращается). Для урав­нения солепереноса краевые условия жестко заданы.

В области [0,L] разбивается сетка. В однородной области реко­мендуется неравномерная сетка с координатами: 0;.0,1; 0,2; 0,3;

0,5; 0,75; 1,0; 1,25; 1,5; 1,75; 2,0; 2,25; 2,5; 3,0; 3,5; 4,0; 5,0 и т.д. На нижней границе сетку рекомендуется сгущать, напри­мер, если L = 4,0 м, то берем узлы до 3,5 м так как предложено выше и вводим дополнительный узел на глубине 3,75 м. Если область [0,L] включает несколько литологических слоев, рекоменду­ется для повышения точности узлы располагать так, чтобы в каждый слой попало не менее 2-х узлов. В особо ответственных случаях узлы на границе сдоев следует располагать на расстоянии 0,01 м по обе стороны границы. Например, граница двух слоев находится на глубине 1,0 м. Для удобства задаем HS = 0,99. Граничные узлы располагаем на глубинах 0,98 м (попадаем в верхний слой) и 1,0 м (попадает в нижний слой).

Затем готовятся сведения о водно-физических и гидрохимических свойствах почвогрунта по литологическим слоям, попадающим в об­ласть расчета (рекомендуется выделять пахотный горизонт, даже если расчетный профиль однородный). Экс­периментальные кривые P(W) и k(W) аппроксимируются о помо­щью предложенных выше аналитических формул (24) и (25) и по мето­ду наименьших квадратов находятся необходимые параметры,.  Если экспериментальных кривых P(W) и k(W) не имеется, то зависимость P(W) принимается оприори по формуле (24) о учетом имеющихся данных (высота капиллярного поднятия, пористость, полная влагоемкость, наименьшая влагоемкость максимальная гигроскопичность и т.д.), а параметр коэффициента влагопроводности  подбирается путем сопоставления расчетных и экспериментальных данных по перераспределению влаги при поливе. Для этого на ЭВМ моделируется полив при известном времени и норме полива. Расчет проводится несколько раз при различных значениях .  Полученные значения расчетных эпюр влажности сравнива­ются о экспериментальной эпюрой и по наилучшему совпадению выбирается

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

Значения полной влагоемкости можно принимать: Wп~(0,9...0,95)m, где m - пористость. О параметрах солепереноса говорилось выше.

В качестве начальных условий для уравнения влагопереноса можно задавать влажность (причем в зоне полного насыщения она полага­ется равной полной влагоемкости соответствующего слоя). При близ­ком УГВ можно задавать равновесную эпюру через потенциал: напри­мер при УГВ = 2 м задается значение потенциала H на всех глу­бинах HJ равным H = -2 м.

Начальные условия для уравнения солепереноса можно задавать по данным отбора естественного перового раствора путем вакуумирования через микропористые зонды. Можно также задавать исходные концентрации по данным водной вытяжки. Программа автоматически пересчитывает состав фаз солей на заданную влажность,  при этом рассчитывается состав ППК и выпадение "лишних" солей в твердую фазу (гипс). Если в почве содержится гипс - в твердей фазе или пред­полагается его внесение в качестве мелиоранта, то его концентра­ции должны быть определены каким-либо способом и заданы в качестве начального условия. Если расчетный состав ППК не совпадает с фактическим, то рекомендуется выполнить корректировку исходных данны ( или лучше исходных концентраций ионов). При увели­чении    содержание Na* в ППК увеличивается по сравнению с , а при увеличении    содержание  в ППК увеличивается по сравнению с

После выбора параметров уравнений весь расчетный период раз­бивается на годы (по 365 сут»), а каждый год на интервалы  (i £ 50). На каждый интервал здаются суммарное испарение (ЕСМ), осадки (ОС), поливная норма (NП), доля транспирации FTR, глубина корневой системы (НКОP), нижнее краевое (Н2). Если peжим орошения задан, то разбиение на интервалы рекомендуется тип: полив перерыв, полив, перерыв и т.д. При этом можно задавать полную норму нетто, тогда на период полива ECM = OC = 0. Во вневегетационный период можно разбивку делать  более крупную. Зимнее снеготая­ние и влагозарядочные поливы задаются отдельно. Долю транспирации рекомендуется в начале вегетации задавать 0.1...0,3, постепенно

увеличивая до 0,6...0,8 в конце вегетации. Даже если растения нет (пар), то рекомендуется фиктивно его вводить.

Значение признака печати (КР) на последний интервал времени в каждом году рекомендуется задавать равным 0, т.к. в коце каждо­го года происходит вывод на печать эпюр во всех узлах.

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

Пример: Расчет водно-солевого режима орошаемых черноземов в Каушанском районе МССР (по данным института молдгипроводхоз)

Грунтовые воды на рассматриваемом массиве расположены глубоко. Требуется выполнить прогноз водно-солевого режима на одну ротацию севооборота (8 лет). Известно, что за 8 лет орошения подъема УГВ до глубины Зм не произойдет. За расчетную принимаем толщу немного боль­ше корнеобитаемого слоя L =1,5 м. На нижней границе принимаем усло­вие свободного гравитационного отекания влаги  (G2 = 4, H2 = 0). Водно-физические и гидрохимические свойства определены в лабораторных ус­ловиях. Осреднение выполнено по слоям 0-0,5м (почва), 0,5-0,8 м (суг­линок средний), 0.8-1,5 м (суглинок средний) табл. 1

Табл. 1 Гидрофизические свойства почвогрунта

Глубина, м	m	Wп	Wппв	Wз	Wмг	Кф		g
0,0-0,5	0,502	0.460	0,340	0,124	0,092	0,40	2.90	1,26	5,0
0.5-0,8	0,480	0.440	0,300	0,134	0,099	0,50	2.56	1.31	7.0
0.8-1,5	0,470	0.430	0,260	0,122	0.090	0.40	2,59	1,16	7,5

 

Табл. 3 Гидрохимические свойства почвогрунта

Глубина, м	l*
0,0-0.5	0,5	0,055	2,00	30,7	0.05
0,5-0,8	0,5	0,045	1,90	29,6	0,05
0,8-1,5	0,5	0,0I7	1,45	28.0	0,05

 

Разностную сетку принимаем в узлах X: 0; 0.1; 0,2; 0.3; 0,4; 0,5;

0,6; 0.8; 1,0; 1,2; 1,5 м. Разбивка до 0.5, выполнена  часто для получения более подробной информации в этом слое при расчете.

Число литологических слоев Ks = 3, границы слоев приникаем на глу­бинах Нs: 0; 0.49; 0,79 м. Тогда первые 5 узлов попадают в 1-ый слой, следующие 2 - во 2-ой, и последние 4 -в третий.

Исходные влажность и засоление (по данным водной вытяжки) зада­ны на 11 апреля 1988 г (табл. 4), варианты качества оросительной во­ды, для которых нужно произвести расчет, приведены в табл. 5.

Для расчета приняты метеоданные за 1961-68 гг. Режим орошения рас­считан по дефициту водного баланса, Так, как сезонное промерзание почвы отсутствует, то на зимний период времени задаются потоки влаги, как и на остальные ( в противном случае следовало бы на зимний пери­од нулевые статьи водного баланса, а осадки за вычетом поверхностно­го стока и суммарное испарение - отнести на период весеннего снеготаяния.

Табл. 4. Исходные эпюры влажности и засоления

Глубина, м	Влажность,    W
0.0  0.2  0.4  0,6  0,8 1,0 1.5	0.309 0.315 0.316 0.321 0.260 0.271 0.274	0.010 0.010 0,012 0.014 0.014 0,013 0,013	0,002 0,002 0,003 О.00З O.004 O.004 O.004	O.003 0,003 0,003 0,003 0,007 0,0I2 0,0I2	0,003 0,003 0,U09 0,009 0,009 0,009 0,009	0,007 0,007 O.007 0,007 0,007 0,007 O.007	0,02 0.02 0.02: 0.02 0,02 0.02 0,02

 

Табл.5 Варианты качества оросительной воды /мг.экв/д/

№							С г/л	SAR
I 2 3 4  5	1,1    2,2    2.7    4,8    4,8	2.1 .        4,0      6,6      9,0      9.6	3,9    7.7     12,7    17,8    31.5	2.2   4,4   7.1   11,1   15.0	3,0      8,4      11,1      16,0      25,1	1,9     1,1     3,8     4.5     5,8	0,48        0,90        1,44        2,05        3,03	3.1 4,4    5,9    6,8  11.7

 

Данные, используемые для расчета, на 1-ый год приведены в табл.6. Доля транспирации задается приближенно, поливная норма-нетто. На ос­тальные годы заполнение производится аналогично. Распечатки исходных данных и результатов расчета на 1-ый год приведены ниже.

На рис.   приведены результаты прогноза изменения засоления и сос­тава ППК для двух значений качества оросительной воды: 3 и 5 из табл. 5.

 Орошение водой с минерализацией 1,44 г/л не вызывает засоления за 8 лет за счет значительных поливных норм и промывки осадками, хо­тя в отдельные "сухие" годы засоление превышает критическое. При оро­шении водам минерализацией 3 г/л происходит засоление.  В обоих слу­чаях значительно повышаются доли натрия и магния в ППК. При орошении водой - 3 г/л кальций из ППК, замещенный натрием и магнием, попадая в почвенной раствор, соединяется с , которого довольно много пос­тупает с оросительной водой. В результате по прогнозу на глубинах 0,4…0,8 м должна образоваться гипсовая прослойка, которая может ухудшить водно-физические свойства почвы. Ухудшение может также произойти в результате осолонцевания.