Географический факультет Московского государственного университета им

географических наук, профессор

ФИЗИКО-МАТЕМАТИЧЕСКИЕ ОСНОВЫ ЛАНДШАФТОВЕДЕНИЯ
(Часть II)

моделей функционирования

Радиальные процессы формирования элементарных геосистем
Латеральное сопряжение геосистем элементарного

водосбора

Формирование структуры геосистем высоких иерархических порядков

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

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

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

Процессы латерального переноса на водосборах, барьеры и др.

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

вещества и энергии высокой интенсивности обладают способностью формировать специфический рельеф (флювиальный, гляциальный, эоловый и т.д.), а также прямо или косвенно обусловливать распределение и численность растений и животных, особенности почвенного покрова, воздействовать на другие потоки. Таким образом, они образуют сферу влияния, которая и есть геосистема . В бореальных условиях главным структурообразующим потоком является водный сток.
Следовательно:
геосистемы различных порядков могут быть выделены в соответствии со схемой Стралера-Философова на множестве элементов рельефа по значениям морфометрических величин, описывающих распределения воды в поле гравитации:
водоразделы любого порядка одновременно соответствуют локальным: максимумам высоты h, минимумам величины удельной площади водосбора (SCA), а также локальным максимумам kh (положительная величина);
тальвеги соответствуют локальным: минимумам высоты h, локальным максимумам SCA, а также локальным минимумам kh (отрицательная величина).

2003] (справа, цифры - порядок водотоков, автоматическое выделение программы TOPAZ).

Цветом показан тип растительного покрова, интенсивностью цвета - характер рельефа по градиенту увлажнения

(Y) от его порядка (X) имеет вид
Y=b0*X**b1;
значения параметров:

b0=0.42, b1=2.52; достоверность модели R2=0.99977.

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

морфометрические величины – параметрами моделей

(частицы)
абсолютно твердое тело
состояние физической системы
( масса

частиц mi, координаты х, у, z
и компоненты скорости vx , vy , vz
в заданный момент времени t)


параметры сплошной среды (в переменных Эйлера)













уравнение неразрывности (закон сохранения)











закон движения системы частиц (закон Ньютона)





уравнения движения вязкой жидкости (система уравнений Навье-Стокса для несжимаемой жидкости)

осреднение турбулентного потока в смысле



система уравнений Рейнольдса (для несжимаемой жидкости )




модель двумерного стекания воды по поверхности
водосбора (интегрировании системы уравнений Рейнольдса [Кучмент, 1983])













система уравнений Сен-Венана





уравнение кинематической волны

dm=ρdx
(qm+x qm)dt
Приращение массы на отрезке среды от х

до х+dx равна разности потоков:

Делим на dt·dx обе части уравнения и переносим в одну сторону

Применяя оператор пространственного дифференцирования для трехмерного случая:

следующий момент при равновероятном перемещении частиц
К выводу уравнения диффузии (градиентные

законы)

частиц через точку х=(j+1/2)h в момент t =(п+1/2) в положительном

направлении переместится за единицу времени (Nj, n-Nj+1, n) 2 частиц , если масса каждой точки m0? то переместиться вещества:

многомерном случае:
Уравнения диффузии и дисперсии
В случае дисперсии с потоком:

рyx, рyy, pyz в направлении у
и pzy, pzz,

pzx. в направлении z
Тогда на элементарную площадку dsx = dydz
в направлении х на элемент поверхности ds действует сила dsxpxx + dsypyx+ dszpzx

dfk = dsi рik

Тензор рik

симметричный

dfk =-pdsi

В общем случае

Для давления Паскаля

Уравнение Ньютона

Тензор напряжений вязкого течения

Идеальная жидкость уравнение Эйлера

несжимаемая, то divv =0
Принимая

получим:

Уравнения Навье-Стокса: движения вязкой жидкости:

Для несжимаемой жидкости:

Вывод уравнений движения вязкой жидкости Навье-Стокса

Рейнольдса в тензорном виде
Упрощения: 1) вертикальные размеры потоков малы

по сравнению с горизонтальными размерами – тогда и ускорения в осреднением движении жидкости вдоль оси z малы по сравнению с ускорениями вдоль осей х и у; 2) вертикальные размеров потоков малы - значит касательные напряжения в жидкости меняются по оси z гораздо сильнее, чем по осям х, и у; 3) потоки воды в равнинных условиях имеют малые скорости течения - поэтому турбулентное давление обусловленное пульсациями скоростей в них мало по сравнению с осредненным гидродинамическим давлением; 4) молекулярные вязкостные напряжения малы по сравнению с турбулентными напряжениями и могут быть отброшены; 5) для природных потоков справедливо приближение Буссинеска - линейная связь между турбулентными напряжениями и градиентами осредненных скоростей – Тогда:

систему уравнений планового движения воды в речном русле:
z' (х, у,

t) - свободная поверхность воды
z0=(х, у) - поверхность водосбора
-=h - глубина потока
qx , qy - расходы воды по направлениям
Сx, Сy - коэффициенты Шези
R - интенсивность осадков
I - инфильтрация воды в почву

Интегрирование уравнений поверхностного стока на водосборе

Уравнение Сен-Венана

Уравнение кинематической волны

q - функция источника (боковой приток в единицу времени на

единицу длины водотока),
uq- относительная скорость бокового притока, if - уклон трения по формуле Шези, i0 - уклон дна, положительный в сторону уменьшения отметок дна водотока

поверхности склона;
h = h(x, t) - глубина потока воды

на поверхности почвы;
i(x) - уклон поверхности склона; n(x) – коэф-т шероховатости
Шези- Маннинга; R(t) - интенсивность осадков;
I(t) - интенсивность инфильтрации; t - время.
Начальные условия:

т.е. в период от t > 0 до t < tS, слой воды на поверхности почвы отсутствует.
Граничные условия:

ХИМИЧЕСКИЙ СОСТАВ ПОВЕРХНОСТНОГО СТОКА

C, CS - концентрация в потоке и в растворе зоны эффективного взаимодействия (мг/л);
D– коэф-т дисперсии; S - концентрация сорбированного вещества в почве (мг/г); R, I – интенсивность
атмосферных осадков и инфильтрации в почву (мм/мин); CR – концентрация атмосферных осадков (мг/л);
h - глубина потока и hs - "эффективная глубина взаимодействия";
и s - масса доступной почвы на единицу объема поверхностного и внутрипочвенного потоков (г/м3 );
 – коэф-т обмена между потоком и зоной взаимодействия (1/с);  и Kd – параметры сорбции-десорбции;
F - площадь поверхности растворяющегося твердого вещества (соли); K =D/ - коэф-т растворения соли
 – толщина переходного слоя, C0- - концентрации насыщения соли у твердой поверхности

УДОБРЕНИЙ С ДОЖДЕВАЛЬНОЙ ПЛОЩАДКИ
Экспериментальные (1) и рассчитанные (2)
гидрографы стока: поверхностный

сток с поля (А),
в замыкающем лесном створе (Б), внутрипочвенный
сток в замыкающем створе (В).
Профили влажности почвы: исходные (1) и
рассчитанные после дождевания (1) на поле (Г)
и в лесу (Д).

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

количества (В); А -распределение поля осадков под пологом леса

для отдельных деревьев (Мк -между крон, Кк - край кроны, Ск - середина кроны, Ус - приствольная часть кронового пространства, А - сосняк чернично-брусничный, Б - сосняк лишайниковый, В - ельник черничный, Г - ельник разнотравный) [Волокитина, 1979, и др.]

3.2. ПЕРЕХВАТ ОСАДКОВ ПОЛОГОМ ЛЕСА

А

В

условием:
D(z,0) = D0(z) ,
С(z,t) – концентрация компонента в атмосферных

осадках, мг/см3 ;
S(z,t) - концентрация компонентов на поверхности клеток фитоэлементов в единице объема полога леса, мг/см3 ; С0 (z) – начальная концентрация солей и аэрозолей на поверхности фитоэлементов в единице объема полога, мг/см3 ; k – эмпирический коэффициент; К(z) - коэффициент скорости растворения солей, см/мин

МОДЕЛИРОВАНИЕ ПЕРЕХВАТА ДОЖДЕВЫХ ОСАДКОВ ПОЛОГОМ ЛЕСА

R(z,t)- интенсивность осадков, см/мин;
Е0 (z,t) - испарение с увлажненной поверхности фитоэлементов (испаряемость), см/мин;
U(z) – площадь поверхности фитоэлементов в единице объема полога на высоте z, 1/см;
G(z) – проекция площади фитоэлементов на единицу поверхности горизонтального сечения на высоте z;
D(z,t) – безразмерная доля сухих листьев и ветвей в момент t на уровне z;.
(z)- средняя сорбируемость воды на единицу площади фитоэлемента.
z=0 на вершине полога леса, положительна вниз, z=Н – нижняя граница полога

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

леса во времени (R0 - над пологом, D = 1

сухие кроны, ZH - низ полога)

Интерфейс расчета проникновения
дождевых осадков под крону ели
с учетом разбрызгивания

- объемная влажность почвы; z - положительно вниз; t -

время;
S(х, t) - расход влаги корнями растений на единицу объема почвы;
К() - коэффициент гидравлической проводимости почвы:

Кf - коэффициент фильтрации в вертикальном направлении;
ВЗ влажность, соответствующая прочно- и рыхлосвязанной воде (ВЗ); ПВ - влажность насыщения (ПВ);
DW() - коэффициент диффузии влаги:

,

РB - капиллярно-сорбционный потенциал насыщения почвы с защемленным воздухом.
(z) - концентрация корней

.

E(t) – испарение с поверхности почвы; d(t) - дефицит влажности воздуха,
r1 ,r2- эмпирические коэффициенты.

Начальные условия: Граничные условия

ВЕРТИКАЛЬНЫЙ МАССОПЕРЕНОС В ПОЧВАХ

С - концентрация подвижного растворенного компонента, мг/л;
S- концентрация адсорбированного компонента, мг/г;
D - коэффициент гидродинамической дисперсии;
v - скорость течения воды в порах;  - объемное содержание воды в порах;
 - объемный вес пористой среды; t - время; z - расстояние;
 - коэффициент массопереноса в 1/день; K, N - коэффициенты изотермы сорбции
Индексы m, im -обозначение мобильной и иммобильной зон.

Для почвенной колонки мощностью L граничные условия: начальные условия:

ДЕРНОВО-ПОДЗОЛИСТЫХ ПОЧВАХ
Ход осадков (а), рассчитанные интенсивности транспирации (б) и испарения

(в), а также профили влажности почв 12 июня (г), 22 июня (д) и измеренная концентрация корней (е). 1 - рассчитанные профили; 2 - измеренные профили влажности почв

Интенсивность фильтрации растворов и промывных вод
в экспериментах на монолитах дерново-подзолистых почв.
А - намытая среднесуглинистая, монолит 3;
Б - среднесуглинистая на покровных суглинках, монолит 4;
1 - температура воды; 2 - подача раствора удобрений;
3 - расходы фильтрата, л/ч.

Экспериментальные данные (точки) и рассчитанные в
численном эксперименте выходные кривые для
дерново-подзолистой среднесуглинистой почвы
под ельником-кисличником

= gl
P1–P2+G sin =  или p1-p2+ glsin  =

, Так как

lsin  =z1-z2

p1-p2+g(z1-z2) =.

H =p/g+z – пьезометрический (гидростатический) напор

Q=-k(H2-H1)/l=-k(H/l)

Равновесие элементарного столбика жидкости

совершая предельный переход:

3.4. ГРУНТОВЫЙ И ПОДЗЕМНЫЙ СТОК В ГЕОСИСТЕМАХ

x=x1;
H=h2 при x=x2;
Граничные условия:
разделим переменные

и проинтегрируем дифференциальное уравнение Дарси с учетом пределов изменения переменных величин h и х:

Кривая депрессии потока грунтовых вод:

(Окский заповедник). А –

расчет установившегося среднегодового УГВ (1) и измеренные УГВ в летнюю межень (2 – август 1986 г, 3 - август 1987 г).
Б – расчет инфильтрационного питания грунтовых вод в августе 1987 года

тела в точке х в момент t, тогда плотность тепловой

энергии Q можно записать Q = с0Т
0 - массовая линейная плотность тела, с - коэффициент теплоемкости кал/(г·град) при небольших изменениях Т можно принять постоянным.
Чтобы получить из общего уравнения переноса уравнение переноса теплоты, вводим потенциал переноса в виде  =СТ.
Согласно градиентному закону Фурье: (1)

 - коэффициент теплопроводности, кал/(см °С·сут) или Дж/(с∙м∙0С) = Вт/(м·К), можно считать постоянным. Подставляя (1) в уравнение неразрывности для тепловой энергии: (2)

вынося  и с0 за скобки и обозначая kT=/c0 коэффициент температуро-проводности, получаем уравнение теплопроводности:
(3)

В общем виде уравнение (2) переноса теплоты

-примеры описания экспериментальных данных для образцов из горизонтов Апах(1), АЕ(2),

и ЕВ(3) [Архангельская, 2012]

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

(4)

раздела фаз
с дополнительными условиями:
граничное условие:

начальное условие:
С условиями на границе

промерзания воды:

Распределения температуры при наличии фазового перехода и движения раздела фаз

Условие в почве с содержанием воды:

Совместный эффект влияния на тепло- и влагоперенос термических и влажностных градиентов в почвах:

коэффициенты теплопроводности мерзлого грунта

структура болотных ландшафтов

2012]

по вертикали от начала координат до уровня грунтовых вод,

y1 - до поверхности болота; z - уровень грунтовых вод от поверхности массива: z =f(r, у, ); zо толщина деятельного горизонта: z0 = (r, у, ); 
Р - интенсивность прихода влаги на единицу площади болота:
w и f - соответственно объемы осадков и испарения за время t.

Уклон поверхности болотного массива i ;
Коэффициент фильтрации kz в деятельном горизонте, k0 – средний в слое (zo-z), kмакс - коэффициент фильтрации при z=0 (т. е. у поверхности болота)
Проточность qz (горизонтальный фильтрационный расход через вертикальное поперечное сечение высотой (zo—z) и шириной, равной единице)
Модуль горизонтальной проточности Mz.

Фильтрационные характеристики болотной фации

В простом одномерном случае:

1984]. 1 - осоково-гипновый; 2 - вейниково-березовый; 3 - сосново-сфагново-кустарничковый,

элемент микрорельефа - повышение; 4 - комплекс группы фаций ленточногрядовой структуры, элемент микрорельефа - гряда сфагново-кустарничково-сосновая; 5 - то же, элемент микрорельефа - гряда сфагново-кустарничковая, облесенная сосной; 6 - сфагново-кустарничковый, облесенный сосной, элемент микрорельефа - повышение

болоте Обловское

верховом болоте Обловское»
Параметры зондирования:
є=45, развертка 1600.
Запасы биомассы умножены на

20 для наглядности

линии cтока.
1- поверхность болота, 2 - уровень грунтовых вод,
3

- граница деятельного горизонта и торфяная залежь,
4 - дно болота и подстилающий минеральный грунт.

s - расстояния вдоль линии стока; bs - длина дуги горизонтали на расстоянии s; b0—длина дуги горизонтали в начале участка; q0 - единичный расход в деятельном горизонте в начале участка; qs - единичный расход на расстоянии s;
Тип ПТК характеризуется: средним модулем проточности Ms, формой потока в плане - b0, bS ,.и внешним (q0) и внутренним (pS) питанием.

Гидродинамический профиль

При наличии водораздельных точек (линий), для которых s=0 и внешнее питание q0 =0 и соответственно i=0 уравнение гидродинамического профиля упрощаем:

Интегрирование уравнения при ряде упрощающих
предположений дает выражение:

n – характеризует сходимость-расходимость линий тока, на симметричном круглом болоте растекание от водо-раздела внутри массива, при n=1 и замене s на r (линии тока совпадают с радиусами) дает [Иванов, 1975 ]:

i =dy/ds

- радиус болотного массива
Если ω - площадь данного болотного

массива то r0 = rn

Нормальные профили болот при М = 2,5 и 20 см2/сек. и максимальных выпуклостях уmax= 2, 4 и 6 м.

Нормальный профиль выпуклых болотных массивов верхового типа

►

►

В соответствии с уравнением Дарси

облесенной сосной фацией в центре [Иванов, 1957].
1—расчетный нормальный профиль

при M =10 см2/сек., уmax = 4м, m=p= 200 мм/год; 2— профиль по нивелировке; 3— профиль летних уровней грунтовых вод.

поверхности болота и растительные ассоциации вдоль трансекта; Б – радарограмма

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

Сравнение измеренных профилей [Сандлерский и др. 2016] и рассчитанного по модели нормального профиля болотного массива с формой купола

пределах изменения r от 0 до rn и у от

0 до yтaх получим

уравнение связывает: средний модуль стока (т) с болотного массива, средний модуль проточности M (сочетание типов болотных фаций, составляющих болотный массив), средний линейный размер массива r (или его площадь ω) и величину максимального превышения поверхности над его окраинами ymax.

Массивы верховых сосново-сфагновых ландшафтов объединены в зону олиготрофных и евтрофных сосново-сфагновых торфяников [Пьявченко, 1985] с границами приблизительно по изолиниям годовой нормы стока 250 мм на севере и 160 мм на юге ЕТР (до 100мм/год в Сибири). Большие грядово-мочажинные комплексы реализуются в зоне выпуклых олиготрофных торфяников с нормой стока 250-300мм/год в ЕТР (до 150мм/год в Сибири).

лавины (кг); g - ускорение свободного падения (м/с2); φ -

угол склона, градусы; Кg - проекция суммарной силы сопротивления движению лавины (м/с2)

ρ – плотность потока; H – высота потока; U – скорость потока; μ - коэффициентом силы сухого кулоновского трения; ξ – коэффициент силы турбулентного сопротивления, пропорциональной квадрату скорости движения снежной массы

Система моделирования снежных лавин RAMMS

Математическое моделирование снежных лавин

ρс - плотность снега; ∆р – разность давления (р1-р2); С – сила сцепления;

Критическая толщина пласта Hn

Коэффициенты сопротивления для моделирования лавин в горах Алатау по модели RAMMS [Благовещенский и др., 2017]

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

и лавин небольшого объема. Небольшие лавины могут моделироваться путем увеличения коэффициентов μ и ξ, - нужна предварительная настройка модели.
Модель не описывает движение снега за фронтом лавины, в области небольших движущихся масс и больших коэффициентов сопротивления.
Еще меньшую точность RAMMS имеет при моделировании селевых потоков, поскольку сели являются многокомпонентными системами (жидкость, грязь, камни, ветки и др)

Моделирования лавин в горах Алатау по модели RAMMS [Благовещенский и др., 2017]

элементы i и j не связаны, то
X – вдоль склона,

Y – вертикальная

Состояние элемента снега i определяется декартовыми координатами его центра (xi ,yi) и двумя составляющими скорости (vxi, vyi); FУij и FВij – силы упругого и вязкого взаимодействия элементов снега i и j; NЭ – общее количество элементов снега в модели.

2) Если элементы i и j связаны и слабо взаимодействуют

3) Если элементы i и j связаны и сильно взаимодействуют

FУxij и FУyij – составляющие силы FУij;
cО и cС – жесткости упругого слабого и сильного взаимодействия элементов,

Ньютона
коэффициент (rij – (dЭ + dm)), характеризует взаимное проникновение элементов

снега друг в друга

FВij –вязкая силы прямо пропорциональна скорости движущегося в среде тела:

При уменьшении коэффициента вязкости снежной массы kв (а) и одновременном согласованном уменьшении параметра R0 в некоторый момент времени начинается образование снежной лавины: рост скорости снежной массы (б) [Соловьев, 2014]

Белые круги - последовательные положения снежных элементов [Соловьев, 2014]. На

склонах малой крутизны, снежные фрагменты скользят по склону или по нижележащим слоям снега. Скорость скольжения элементов, находящихся в непосредственном соприкосновении со склоном минимальна, и чем выше по отношению к склону расположен фрагмент снега, тем больше его скорость (φ ≤ 40°). При большой крутизне склона движение фрагментов снега начинается так же как скольжение, однако через незначительное время сменяется "прыжками" (φ = 50°). Фрагменты снега в лавине находятся в «прыгающем» постоянно ускоряющемся движении фрагментов с интервалом 0,5 с

Верификация модели
Распределения давления по высоте снежной массы h (а), и влияние угла склона α на длину пробега лавины L(б) [Соловьев, 2014]

и функционирования геосистем. Наименьшими материальными объектами (точками), из которых состоят

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