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

Полный текст

Введение

В последние десятилетия в геоморфологии склонов активно применяются математические методы, практическая цель использования которых в этой области науки о рельефе – решение задач сохранения почв и их плодородия [1–10]. Противоэрозионные мероприятия, препятствующие оврагообразованию, разрабатываются с учетом морфодинамики профиля делювиальных склонов. Процессы, протекающие в естественных ландшафтах, продолжительны во времени, поэтому, пожалуй, только математическое моделирование позволяет спрогнозировать характер развития конкретного склона на длительном временном интервале. Изучение динамики формирования склонов методами математического моделирования позволяет определить характерные, наиболее типичные пути развития склонового профиля. Такая информация важна как с точки зрения сохранения плодородия почв, так и для понимания процессов рельефообразования. Модель плоскостного смыва, несмотря на свою ограниченность, дает адекватную картину эволюции эрозионных процессов на пологих склонах: потеря верхнего слоя, где интенсивность смыва является функцией скорости водного потока, приводит к деградации склоновых экосистем.

Физико-математические концепции моделирования склоновых процессов

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

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

Физические концепции теории массопереноса в ламинарном приближении использовались в работах [7, 9, 10], где авторы предполагают различные механизмы отрыва и перемещения частиц грунта с учетом вязкости, гравитации, сил Архимеда, трения и т.д. Формирование взвеси начинается с некоторой критической скорости (v_кр) отрыва частиц от поверхности движущимся потоком воды. Изменение мутности потока (ρ) при этом пропорционально кубу скорости (v³) [1, 9, 11]. Процесс смыва моделируется уравнением с коэффициентом диффузии, отражающим интенсивность денудации [7].

Смысл коэффициента диффузии  (мутности) при моделировании смыва с помощью диффузионных уравнений до конца не выяснен. Как правило, он характеризует интенсивность взаимодействия потока с поверхностью, нарастающую вниз по склону. В такой модели отсутствует важнейшая характеристика – уклон. До настоящего времени не существует единого обоснованного выбора вида функции для коэффициента диффузии (или, в контексте плоскостного смыва, денудации), а его зависимость только от длины склона не представляется физически обоснованной. Действительно, при эволюции делювиального склона должны меняться скорость смыва и осаждения материала в разных точках профиля. Это обусловливает саморазвитие склоновых систем за счет обратной связи, когда зависящий от уклона коэффициент диффузии ответственен за саму форму профиля. Отметим, что система склонов и их эволюция могут быть рассмотрены с позиций синергетики как процесс самоорганизации систем с обратной связью [1, 9, 12].

Результаты и обсуждение

Указанные соображения привели нас к необходимости более полного раскрытия физической картины процессов, протекающих на делювиальных склонах, находящихся под антропогенной нагрузкой. Подойдем к описанию процесса смыва с иной точки зрения и попытаемся приблизить математическую модель к физике явления. Будем считать, что поток по склону плоский и ламинарный, а смыв верхнего слоя пропорционален кубу скорости потока (v³) [1, 7, 9]. Смывом в точке будем называть количество вещества, перенесенное за единицу времени через поперечное сечение потока единичной площади. Рассмотрим геоморфологические процессы на плоском склоне в приближении однородного (гомогенного) плоского потока. Пусть имеется некоторый склон с профилем в виде функции z(х), по которому осуществляется эрозионный перенос вещества потоком воды от снеготаяния или после дождя. Форма профиля при этом будет изменяться со временем. Необходима такая математическая модель, решение которой описывало бы эволюцию профиля адекватно.

В связи с этим в основу нашей модели были заложены балансовые и диффузионные уравнения. Такой выбор оправдан двумя обстоятельствами. Во-первых, несмотря на “скрытость” физических параметров, уравнение диффузии позволяет получать весьма важные результаты о коэффициенте денудации (подвижности верхнего слоя грунта на склоне), необходимого для дальнейшего моделирования с привлечением динамических дифференциальных уравнений [1, 13]. Во-вторых, диффузионная модель позволяет описать и спрогнозировать развитие склона, что также важно для изучения устойчивости рельефа и его самоорганизации.

Основные идеи разработки нашей модели состояли в использовании системы уравнений:

$\frac{\partial z}{\partial t}+\frac{\partial q}{\partial x}=0,q=-k\left(x\right)\cdot \frac{\partial z}{\partial x},\frac{dq}{dt}=0,$ (1)

где z – высота, х – координата по горизонтальной оси или расстояние от начала профиля склона, q – твердый расход материала. Таким образом, процесс денудации моделируется балансовыми диффузионными уравнениями, в которых коэффициент пропорциональности k(x) отражает взаимодействие потока с поверхностью.

Для постоянного коэффициента решение системы (1) можно получить в аналитическом (замкнутом) виде [11, 14, 15]. Исследования с различными зависимостями k(х) показали, что параллельное отступание профиля в процессе его развития описывается с помощью экспоненциального нарастания интенсивности смыва от водораздела к тальвегу. Таким образом, если считать, что денудационный расход материала изменяется от точки к точке, то следует искать некоторую нелинейную зависимость k(х), адекватно описывающую наблюдаемые профили. В принципе, получив экспериментально значения k(х) для пород различного литологического и механического состава можно, при прочих равных условиях, разделить этот коэффициент на две составляющие, одна из которых описывает физические свойства, что существенно увеличит информативность модели. Ранее [16] было установлено эмпирическое соотношение между расходом (q), точкой на склоне (х) и уклоном:

$q=a{x}^m{(-\partial z/\partial x)}^n$, (2)

где $\partial z/\partial x$ – уклон, m, n, a – эмпирические константы. Уравнение (2) указывает на нелинейность процесса (m, n не равны единице).

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

$\rho =\underset{\Delta V\to 0}{\lim }\frac{\Delta q}{\Delta V}$,

где  – количество вынесенного материала в выделенном объеме. Полагая, что мутность потока пропорциональна кубу скорости, а сам процесс смыва начинается с некоторой (критической) скорости (v_кр), можно записать:

$\rho =\alpha \left(\frac{v-v_{кр}}{v_{кр}}\right){\left(\frac{v}{v_{кр}}\right)}^3$, (3)

где a – коэффициент размерности [3, 9]. Будем считать, что толщина потока аппроксимируется степенной функцией (h~x^m), тогда скорость выноса субстрата равна:

$v=\sqrt{gh(-\partial z/\partial x)}=\sqrt{g{x}^m(-\partial z/\partial x),}$ (4)

где g – ускорение свободного падения.

Высота профиля (z) связана с мутностью потока ρ соотношением z=αρ, где α – коэффициент пропорциональности для данного типа поверхности.  Отсюда можно получить связь между расходом наносов q и высотой z(x) в точке наблюдения. Действительно, на элементарном участке склона за одну секунду уносится объем воды равный ΔzΔxDy с массой грунта ΔQ, где Dz – изменение высоты склона на концах отрезка Δх. Для плоских склонов можно считать, что процессы денудации по латеральной оси (вдоль направления протяжения склона) изотропны. Тогда, полагая, что Dy = 1, получим:

$\Delta z\Delta x=\rho (x_0+\Delta x)-\rho \left(x_0\right)=\Delta Q$, (5)

где ρ(x₀) – мутность потока в точке x₀, ρ(x₀+Δх) – мутность потока в точке, отстоящей от x_о на расстоянии Δх. Балансовое уравнение для смыва можно теперь записать в виде:

$\Delta Q=v{\left(\frac{d\rho }{dx}\right)}_{x=x_0}\cdot \Delta x\cdot \Delta t\text{ ,}$ (6)

где v – скорость потока, $\Delta Q$ – смыв за время Δt с участка Δх, ${\left(\frac{d\rho }{dx}\right)}_{x=x_0}$ – градиент мутности в точке х=х₀. Переходя к пределу, получим уравнение баланса в виде:

$\frac{\partial z}{\partial t}=Cv\frac{\partial \rho }{\partial x}=C\sqrt{x^m(-\partial z/\partial х)}\frac{\partial \rho }{\partial x}$, (7)

где С – константа размерности. Это уравнение связывает деформацию высоты склона с пространственным изменением мутности потока на нем. Если учесть, что уклон определяет мутность потока через его скорость, а коэффициент денудации – функция мутности, то равенство (7) описывает саморазвитие склона. С учетом тождеств (3–5) можно выразить мутность потока через уклон поверхности в явном виде и получить окончательную математическую модель денудационного процесса.

Полученное в итоге нелинейное дифференциальное уравнение в частных производных решалось численно методом конечных разностей по неявной схеме [1–17]. Область склона аппроксимировалась конечным множеством точек (сеткой) с различными начально-краевыми формами: в виде ступенчатой, кусочно-линейной, логистической кривых. Краевые условия задавали поведение границ: неподвижные на водоразделе и свободные у подножья, что соответствовало фиксированному водоразделу и полному выносу субстрата через границу расчетной области (т.е. на рассматриваемом интервале профиля аккумуляция отсутствует).

Результаты численного решения уравнения (7) свидетельствуют о том, что исходное слабонаклонное плато с начальным обрывом постепенно размывается и принимает форму логистоподобной (или S-образной) кривой, показанный на рисунке кривой  2. Обращает на себя внимание то, что логистический характер профиля сохраняется вплоть до 4800 лет при средней скорости выноса 1 м в 100 лет. Вместе с тем совершенно очевидно вырождение формы кривой на длительных этапах эволюции: как видно из рисунка (кривая 4), происходит “выполаживание” кривой. Рассмотренный пример для m=1/2 демонстрирует, на наш взгляд, адекватность качественного описания формы профиля в процессе развития. Показатель степени m описывает ход нарастания высоты слоя и скорости потока при движении вниз по склону. Этот параметр может отличаться от 0.5, как это утверждается в работе [1]. Экспериментальное определение значения m для различных типов склонов может явиться важным показателем для разработки теории массопереноса в плоскостном смыве при его дальнейшем переходе к ручейковому.

Эволюция развития склона. Математическая модель (10) для m=1/2. Исходная форма склона выбрана в виде ступенчатой функции.

Численная реализация на различных шагах счета: 1 – исходный профиль, 2 – 240, 3 – 4800, 4 – 5600 лет эволюции

Заключение

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

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

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

Отметим, что параметры m и n допустимо определять эмпирически путем сравнения результатов модели с реальными склонами. Кроме этого, коэффициент денудации, зависящий от положения на склоне, может служить своеобразным индикатором «интуитивной» адекватности модели, одновременно проявляя свою физическую сущность. Рассмотренные зависимости позволяют, на наш взгляд, ближе подойти к пониманию эволюции рельефа, стимулируют развитие математической парадигмы в геоморфологии склоновых систем делювиального генезиса. Очевидная практическая ценность полученных результатов состоит в возможности прогноза темпов и пространственной изменчивости процессов разрушения плодородного слоя почвы.

Об авторах

А. Н. Салугин

Автор, ответственный за переписку.
Email: anastasiya-kulik@yandex.ru
Россия, Волгоград

А. В. Кулик

Email: anastasiya-kulik@yandex.ru
Россия, Волгоград

Список литературы

Дополнительные файлы

Доп. файлы

Действие

1. JATS XML

Скачать

2. Эволюция развития склона. Математическая модель (10) для m=1/2. Исходная форма склона выбрана в виде ступенчатой функции.

Скачать (60KB)

Метаданные