CC BY
25Математические заметки СВФУ Апрель—июнь, 2020. Том 27, № 2
УДК 519.63
ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ПРОЦЕССА ПРОСАЧИВАНИЯ В ТРЕЩИНОВАТО-ПОРИСТЫЙ ГРУНТ В УСЛОВИЯХ КРАЙНЕГО СЕВЕРА
С. П. Степанов, А. В. Григорьев, Н. М. Афанасьева
Аннотация. Приводится математическое моделирование комплексной мультифи-зичной задачи, актуальной для территорий Крайнего Севера и Арктики. Актуальность данной задачи характеризируется важностью процесса просачивания при формировании и оттаивании слоя вечной мерзлоты. Современные прикладные задачи по большей части требуют учета сложных геометрий, а также большого количества разных процессов и их взаимовлияния. Мультифизичная модель состоит из уравнения Ричардса для описания процесса просачивания, модели двойной пористости для описания естественной трещиноватости грунта, задачи Стефана для описания температурного режима грунта в условиях криолитозоны. Вычислительный алгоритм базируется на конечно-элементной аппроксимации по пространству на триангулированных сетках Делоне и использовании схемы расщепления по времени с применением линеаризации с предыдущего временного слоя.
Б01: 10.255877SVFU.2020.15.67.007
Ключевые слова: уравнение Ричардса, задача Стефана, модель двойной пористости, трещиновато-пористые среды.
Введение
Для территорий Крайнего Севера и Арктики характерно наличие слоя толщи вечной мерзлоты в грунте, который напрямую влияет на все процессы жизнедеятельности человека. Огромное значение имеют процессы протаи-вания и формирования «вечной мерзлоты», подвижность грунта несет в себе опасность технических аварий, обвала зданий и сооружений, вывода из строя дорожных покрытий. Для математического моделирования процессов прота-ивания/промерзания используется модель Стефана. Корректность постановки задачи Стефана изучалась в работах О. А. Олейник, А. М. Мейрманова, Л. И. Рубинштейна и др. В работах С. Л. Каменомостской, Б. М. Мудака и А. А. Самарского были впервые предложены численные методы сквозного
Афанасьева Н. М. благодарит Мегагрант правительства РФ N14.Y26.31.0013 (за средства гранта выполнены методические расчеты), Григорьев А. В. благодарит грант РФФИ 17—01—00732А (за средства гранта выполнена реализация модели двойной пористости и модели просачивания), Степанов С. П. благодарит грант РНФ 19—11—00230 (за средства гранта выполнена реализация температурной модели).
© 2020 Степанов С. П., Григорьев А. В., Афанасьева Н. М.
счета при построении дискретных аналогов задач типа Стефана. Если вопросы образования и таяния вечной мерзлоты хорошо изучены [1—7], то процесс просачивания не изучен так детально в данном контексте в силу сложности физической модели. Процесс просачивания воды критически важен для описания формирования и таяния многолетнемерзлых грунтов. Реальные грунты, с одной стороны, содержат естественную трещиноватость — развитую сеть трещин разных масштабов, с другой стороны, трещиноватость может возникать за счет периодического просачивания талых вод. Таким образом необходимо учитывать многие факторы: наличие минимум двух разных континуумов, влияние системы трещин на потоки, конвективный перенос тепла, коэффициенты просачивания, зависящие от температуры и фазовых переходов. Для изучения и моделирования такого рода процессов необходимо изучать комплексные, а самое главное сильно связанные между собой физические процессы, необходимо уделить внимание вопросам адекватности физических моделей [8—11]. Данные вопросы требуют проработанных математических моделей, позволяющих учитывать продвинутую физику процессов [17-20], задействованных в просачивании, а также особенности строения вечномерзлых грунтов. В этих работах рассмотрены основные подходы, применяемые для такого рода задач. Данная работа является комбинацией предложенных ранее идей и методов.
В данной статье рассматривается численное моделирование процесса просачивания талой воды в трещиновато-пористый грунт с учетом температурного режима и фазовых переходов. Математическая модель конструируется путем сочетания нескольких моделей. Приведем основные предположения для построения мультифизичной модели. Процесс просачивания реализуется на основе уравнения Ричардса, где в качестве коэффициентов проницаемости и производной насыщенности по давлению принимаются эмпирические зависимости, полученные исходя из оценок скорости просачивания. Область просачивания разделена на два накладывающихся, но разделенных континуума, которые позволяют описывать процесс просачивания в трещиноватой среде и порах соответственно. Процесс теплопереноса в грунте описывается уравнением теплопроводности и учитывает фазовый переход поровой влаги в лед и наоборот. Влияние насыщенности на температуру учитывается через введение дополнительного конвективного слагаемого, а учет влияния температуры на процесс просачивания производится через коэффициент проницаемости.
Вычислительный алгоритм базируется на конечно-элементной аппроксимации математической модели по пространственным переменным [12]. Для аппроксимации по времени строится чисто неявная разностная схема с линеаризацией с предыдущего временного слоя [13]. Возможности вычислительного алгоритма иллюстрируются расчетами модельной задачи в двухмерной и трехмерной постановках.
1. Математическая модель
Рассмотрим процесс просачивания воды в трещиновато-пористый грунт в условиях криолитозоны. Для этого запишем связанную математическую модель.
Уравнение Ричардса и модель двойной пористости. С целью описания течения жидкости в трещиновато-пористой среде Баренблатт, Желтов, Кочина [14] разработали модель двойной пористости, для описания процесса просачивания Ричардс предложил уравнение, обобщающее закон Дарси. Отметим, что существует три различных формы записи уравнения Ричардса [15,16]: в терминах давления, в терминах насыщенности и в смешанной форме. Мы, в свою очередь, воспользуемся моделью двойной пористости к уравнению Ричард-са, записанного в терминах давления:
_ у(р1 + 2,)) +r(x)(p1 -p2) = 0,
dp i dt (1)
TO2^fi_ v(p2 2)) _ r{x){pi -p2) = 0,
dp2 dt
здесь индексами 1 и 2 обозначены параметры трещиноватой и пористой сред соответственно, pa = pa/pg — приведенные давления, pa — давления, ma — пористости, sa(pa) — насыщенности, Ka(pa) — гидравлические проницаемости, r(x)(p1 — p2) — обменный переток между средами. Для коэффициентов верны следующие зависимости:
Sa(pa) = 1 — exp(—YPa), Ka(pa) = Ksasa(pa)a. (2)
Ksa — проницаемость полностью насыщенной среды, y, <г — коэффициенты задачи, a = 1, 2 соответственно.
Задача Стефана для описания температурного режима. Для моделирования теплового режима грунтов используется уравнение теплопроводности с учетом фазовых переходов поровой влаги
дТ
ЫФа) + рМ*)-^ - div(A(^)VT) = 0. (3)
На практике фазовые превращения происходят в малом интервале температуры [Т* — Д,Т* + А]. Пусть
0, Т < Т* — А,
{ ' 2\ 7 А 7 7 А- (4)
Т > Т* + А,
. 0, Т < Т* — А, ^ J 1 2Д. / А / / А. (5)
0, Т > Т* + Д.
Для коэффициентов уравнения верны следующие соотношения:
ср(фд) = р-с- + фд(р+ с+ - р-с-), А(фд) = А- + фд(А+ - А-), (6)
где Ь — удельная теплота фазового перехода, р+, с+, А+, р-, с-, А- — плотность, удельная теплоемкость, теплопроводность талой и мерзлой зон соответственно.
Полная модель. Проведем адаптацию полной физической модели. Влияние насыщенности на температуру учитывается через введение дополнительного конвективного слагаемого:
с+ р+ М(К1(Р1)УР1 + К2(р2)Ур2, УТ). (7)
Учет влияния температуры на процесс просачивания производится через коэффициент проницаемости (если отметить гидравлическую проницаемость через
Ка):
К ( ) / Ka(Pa), Т> 0 , 2 (8)
К а (Ра) = < ~ , ч а = 1, 2, (8)
\еКа(ра), Т < 0,
где е = 10-6 — малое число. Таким образом, основываясь на (1)—(3), (7), (8), запишем полную систему уравнений, описывающую процесс просачивания в трещиновато-пористой среде с учетом температуры и фазовых переходов:
_ + 2,)) +г(ж)(р1 _р2) = о,
др1 Ы
_ div(^Í2(p2)V(p2 + -г)) - г{х){р1 - р2) = 0, (9)
др2 Ы дТ
{ср + - сиу(а(Фа)УТ) + с+р+р{К1{р1)Ур1 + К2{р2)Ур2, УТ) = 0.
Граничные и начальные условия. Рассмотрим двумерную область О С В? с границей Г = дО (рис. 1). Граница имеет вид Г = Ги + Г, + Ц, в свою очередь, Ги = + Г^. Дополним полную систему граничными и начальными условиями:
Для температуры есть три типа граничных условий. На верхней части области Ги:
дТ
-КФл)7Г=13(т-та1Г), хети, (10)
дп
здесь Тагг — температура воздуха. На боковых частях границы:
—МФА)тг~ = 0, (И)
дп
На нижней части области требуется выполнение условия Дирихле:
Т = -2 0°, х е Гь. (12)
Рис. 1. Область О. Для давления (насыщенности) всюду стоят условия непротекания жидкости:
-Ка(Ра)(У(ра + *), п) = 0, x е дО, (13)
кроме части Г^п, где при условии Та^г > 15 С°:
Ра = 1, x е , Тагг > 15 С°. (14)
Начальные условия для температуры:
Т = 0, х е О, г = 0. (15)
Начальные условия для давления (насыщенности):
Р1 = 0, Р2 = 0, x е О, г = 0. (16)
2. Конечно-элементная дискретизация
Численная реализация основана на методе конечных элементов, а именно на использовании вычислительной библиотеки ЕЕшСЯ [12]. Для описания аппроксимации по времени введем равномерную сетку с шагом т.
¿V = гп = пт, п = 0,1,...,^, т^ = Т, (17)
и введем обозначения рП = Р1 (гп), рП = Р2(гп), Тп = Т(гп). Для аппроксимации по пространству использованы стандартные лагранжевы конечные элементы. Дальше введем пространство конечных элементов V С Н1 (О). Здесь Н 1(О) — пространство Соболева. Для аппроксимации по времени используем аналог
неявной разностной схемы с линеаризацией с предыдущего временного слоя: А«, «"+1 _ пп
г " ^^ +К®)(рГ1 -РГ1) = 0,
д « р"+1 _
та2£(Г) г - ММ^МРГ1 +*)) -Р2 + 1) = О,
Т п+1 _Т П
(ср™ --сИу(А(0д)пУГп+1) + (18)
т
+с+р+м(К1(рп) Урп + к^рп) vтn) = о.
В дальнейшем представим систему (18) в вариационном виде. Для этого домножим уравнения на соответствующие тестовые функции и проинтегрируем по области О. Тем самым с учетом граничных условий получим вариационную формулировку:
я«, «"+1 - рп
77111 т + !{К^у^1 + г)уу1) с1х+
п п
+ 1 г(ж)(рп+1 - рп+>^х + у к1(рп)^1 - I к1(рп)«1 = 0,
У1 - ^2 У 1й-Х У 8 — J *4V.F1
п г„ г(
т2/+ + / (М^Мр^ + г)^)^
п о
-1 Л™ I („П
-I г(х)(рп+1 - рп+>2 ¿X + У к2(рп)^2 -I к2(рп)«2^8 = 0, (19) п г„ гь
/* ТП+1 _ ТП /*
/ (ср™ + ---«з с1х+ I (\(фА)п\7Тп+1, Чу3) ¿х
п п
+ 1 с+р+м(к1(рп)Урп + к2(рп)Урп, УТ">2 ¿X п
+ 1 в(Тп+1 - ТаггЬ = 0.
3. Результаты расчетов и выводы
Проведем численное моделирование рассматриваемой задачи в двумерной постановке. Расчетная область состоит из нескольких слоев грунта (рис. 2).
Грунт имеет начальную температуру То = 2°С.
Температура на дневной поверхности задавалась с учетом амплитуды колебания температуры воздуха на рассматриваемой территории, которая варьируется в пределах от Тш = -40°С зимой до Та = +20°С летом, и задавалась следующей формулой:
_Тп)-Т3 . 7г(30 • (М — 1) + т + 75) , тп)+т; агг~ 2 '8Ш 180 + 2
Рис. 2. Расчетная сетка
где М — начальный месяц.
Для численного моделирования были использованы следующие параметры: коэффициент вязкости / = 1.0 • 10-3;
коэффициенты проницаемости Кз1 = 1.0е-6, К,2 = 1.0е-9; коэффициент обменного перетока г(х) = г = 2.0 • 10-4; параметры задачи а = 2.0, 7 = 1.0, в = 14.0.
1 слой. Объемная теплоемкость ср: талый 2397.6 • 103, мерзлый 1886.4 • 103; коэффициент теплопроводности а: талый 1.37, мерзлый 1.72;
теплота фазового перехода Ь: 75330 • 103.
2 слой. Объемная теплоемкость ср: талый 2130.0 • 103, мерзлый 2090.0 • 103; коэффициент теплопроводности а: талый 2.67, мерзлый 3.37;
теплота фазового перехода Ь: 64769 • 103.
3 слой. Объемная теплоемкость ср: талый 2960.0 • 103, мерзлый 2700.0 • 103; коэффициент теплопроводности а: талый 1.4, мерзлый 1.56;
теплота фазового перехода Ь: 130544 • 103.
4 слой. Объемная теплоемкость ср: талый 2480.0 • 103, мерзлый 1890.0 • 103; коэффициент теплопроводности а: талый 1.57, мерзлый 1.86;
теплота фазового перехода Ь: 125532 • 103.
Расчеты проводились на 1 год с шагом 1 сутки (24 часа) в расчетной сетке размером 59683 треугольных элементов.
На рис. 3 представлены результаты моделирования. Расчетные данные иллюстрируют, что влагонасыщенность грунта существенно влияет на сезонную глубину протаивания. Данная работа носит модельный характер и демонстрирует применимость метода как для случая двумерной модельной задачи, так и для случая трехмерной модельной задачи. Таким образом, методические расчеты были проведены на двумерном случае и после отработки методики на простом примере задача была масштабирована до трехмерного квазиреального случая.
Представим результаты численного моделирования распространения тепла в трехмерной постановке. Расчетная сетка содержит 246508 тетраэдров (рис. 4).
з2 Т
Рис. 3. Результаты для двумерной задачи: слева насыщенность для трещин, в середине для пор, справа распределение температуры для различных временных слоев
Рис. 4. Расчетная сетка
Все коэффициенты, начальные условия и граничные условия ставятся аналогично двумерному случаю.
Результаты расчетов для трехмерной задачи представлены на рис. 5. Выполнены исследования по моделированию процесса просачивания воды
Рис. 5. Трехмерные расчеты: слева насыщенность для трещин, в середине — для пор, справа — распределение температуры для различных временных слоев.
в криолитозону. Проведено разделение процессов просачивания в трещиновато-пористый грунт от температурных процессов с учетом фазовых переходов, тем самым реализована схема расщепления по физическим процессам. В будущем планируется провести дальнейшую адаптацию предложенного метода для случая реальных данных, в том числе с применением нейронных сетей, для более качественной процедуры адаптации к данным.
ЛИТЕРАТУРА
1. Samarskii A. A, Moiseenko B. D. Through calculation scheme for the multidimensional Stefan problem // Comput. Mathematics and Math. Phys. 1965. V. 5, N 5. P. 816-827.
2. Samarskii A. A., Vabishchevich P. N., Iliev O. P., Churbanov A. G. Numerical simulation of convection/diffusion phase change problems: a review // Intern. J. Heat and Mass Transfer. 1993. V. 36, N 17. P. 4095-4106.
3. Samarskii A. A., Vabishchevich P. N., Mathematical Modelling. V. 1. Computational Heat Transfer. Wiley, 1995.
4. Cheng Q., Sun T., Jones S. B. et al. In situ measured and simulated seasonal freeze-thaw
cycle: A 2-year comparative study between layered and homogeneous field soil profiles //J. Hydrology. 2014. V. 519. P. 1466-1473.
5. Васильев В. И., Васильева М. В., Степанов С. П. и др. Математическое моделирование температурного режима грунтов фундаментов в условиях многолетнемерзлых пород // Вестник Моск. гос. техн. ун-та им. НЭ Баумана. Сер. Естественные науки. 2017. № 1. С. 142-159.
6. Vabishchevich P. N., Varlamov S. P., Vasiliev V. I. et al. Numerical simulation of the temperature dynamics of railway foundation material in permafrost // Math. Models Comput. Simulations. 2017. V. 9, N 3. P. 292-304.
7. Иванов В. А. Численное исследование влияния теплоизоляции на режим работы магистрального газопровода в условиях Крайнего Севера // Мат. заметки СВФУ. 2017. Т. 24, №4. С. 96-108. doi: https://doi.org/10.25587/SVFU.2018.4.11320.
8. Masters I., Pao W. K. S, Lewis R. W. Coupling temperature to a double-porosity model of deformable porous media // Intern. J. Numer. Meth. Engin. 2000. V. 49, N 3. P. 421-438.
9. Helmig R. et al. Multiphase flow and transport processes in the subsurface: a contribution to the modeling of hydrosystems. Berlin: Springer-Verl., 1997.
10. Bai M., Roegiers J.-C. Fluid flow and heat flow in deformable fractured porous media // Intern J. Engin. Sci. 1994. V. 32, N 10. P. 1615-1633.
11. Васильева М. В. Прокопьев Г. А. Численное решение задачи двухфазной фильтрации с неоднородными коэффициентами методом конечных элементов // Мат. заметки СВФУ. 2017. Т. 23, № 2. C. 46-62. doi: https://doi.org/10.25587/SVFU.2017.2.9245.
12. Langtangen H. P. A FEniCS tutorial // Automated Solution of Differential Equations by the Finite Element Method. Berlin; Heidelberg: Springer-Verl., 2012. P. 1-73.
13. Самарский А. А. Теория разностных схем. М.: Наука, 1989.
14. Barenblatt G. I., Zheltov I. P., Kochina I. N. Basic concepts in the theory of seepage of homogeneous liquids in fissured rocks //J. Appl. Math. Mech. 1960. V. 24, N 5. P. 12861303.
15. Rathfelder K., Abriola L. M. Mass conservative numerical solutions of the head-based Richards equation // Water Resources Research. 1994. V. 30, N 9. P. 2579-2586.
16. Ross P. J. Efficient numerical methods for infiltration using Richards equation // Water Resources Research. 1990. V. 26, N 2. P. 279-290.
17. Celia M. A., Bouloutas E. T., Zarba R. L. A general mass-conservative numerical solution for the unsaturated flow equation // Water resources research. 1990. V. 26, N 7. P. 1483-1496.
18. van Genuchten M. T. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils 1 // Soil Sci. Soc. Amer. J. 1980. V. 44, N 5. P. 892-898.
19. Васильев В. И., Максимов А. М., Петров E. E., Цыпкин Г. Г. Тепломассоперенос в промерзающих и протаивающих грунтах. M.: Наука. 1996.
20. Droniou J., Eymard R., Guichard C. Uniform-in-time convergence of numerical schemes for Richard's and Stefan's models // Finite Volumes for Complex Applications VII-Methods and
Theoretical Aspects. Cham: Springer-Verl., 2014. P. 247-254.
Поступила в редакцию 19 ноября 2019 г. После доработки 13 февраля 2020 г. Принята к публикации 30 апреля 2020 г.
Степанов Сергей Павлович
Северо-Восточный федеральный университет им. М. К. Аммосова, Институт математики и информатики, Международная научно-исследовательская лаборатория
«Многомасштабное математическое моделирование и компьютерные вычисления»,
ул. Кулаковского, 42, Якутск 677000
cepe2a@inbox.ru
Григорьев Александр Виссарионович
Северо-Восточный федеральный университет им. М. К. Аммосова, Институт математики и информатики,
Научно-исследовательская кафедра «Вычислительные технологии»,
ул. Кулаковского, 42, Якутск 677000
re5itsme@gmail.com
Афанасьева Надежда Михайловна
Северо-Восточный федеральный университет им. М. К. Аммосова, Институт математики и информатики,
Научно-исследовательская кафедра «Вычислительные технологии», ул. Кулаковского, 42, Якутск 677000 afanasieva.nm@gmail.com
Математические заметки СВФУ Апрель—июнь, 2020. Том 27, № 2
UDC 519.63
SIMULATION OF THE PROCESS OF INFILTRATION
INTO FRACTURED POROUS SOIL IN PERMAFROST
S. P. Stepanov, A. V. Grigoriev, and N. M. Afanasyeva
Abstract: The article provides mathematical modeling of the complex multiphysical problem relevant for the territories of the Far North and the Arctic. The relevance of this task is characterized by importance of the seepage process in the formation and thawing of the permafrost layer. Modern applications for the most part require consideration of complex geometries, as well as a large number of different processes and their mutual relationship. The multiphysical model consists of the Richards equation to describe the seepage process, the double porosity model to describe natural soil fracturing, the Stefan task to describe the temperature regime of the soil in permafrost zone conditions. The computational algorithm is based on finite-element space approximation on triangulated Delone meshes and using of a time splitting scheme using linearization from a previous time layer.
DOI: 10.25587/SVFU.2020.15.67.007 Keywords: the Richards equation, the Stefan problem, double porosity, fractured porous media.
REFERENCES
1. Samarskii A. A. and Moiseenko B. D., "Through calculation scheme for the multidimensional Stefan problem," Comput. Math. Math. Phys., 5, No. 5, 816-827 (1965).
2. Samarskii A. A., Vabishchevich P. N., Iliev O. P., and Churbanov A. G., "Numerical simulation of convection/diffusion phase change problems: a review," Int. J. Heat Mass Transfer, 36, No. 17, 4095-4106 (1993).
3. Samarskii A. A. and Vabishchevich P. N., Mathematical Modelling, vol. 1, Computational Heat Transfer, Wiley (1995).
4. Cheng Q., Sun Y., Jones S. B. et al., "In situ measured and simulated seasonal freeze-thaw cycle: A 2-year comparative study between layered and homogeneous field soil profiles," J. Hydrology, 519, 1466-1473 (2014).
5. Vasil'ev V. I., Vasilyeva M. V., Stepanov S. P. et al., "Mathematical modeling of the temperature regime of soils of basement foundation on permafrost [in Russian]," Vestn. Mosk. Gos. Tekhn. Univ. im. N. E. Baumana, Ser. Estestv. Nauki, No. 1, 142-159 (2017).
6. Vabishchevich P. N., Varlamov S. P., Vasiliev V. I. et al., "Numerical simulation of the temperature dynamics of railway foundation material in permafrost," Math. Models Comput. Simul., 9, No. 3, 292-304 (2017).
7. Ivanov V. A., "Numerical research on effect of thermal insulation on a gas pipeline's performance in Far North environment," Mat. Zamet. SVFU, 24, No. 4, 96-108 (2017).
doi: https://doi.org/10.25587/SVFU.2018A11320.
8. Masters I., Pao W. K. S., and Lewis R. W., "Coupling temperature to a double-porosity model of deformable porous media," Int. J. Numer. Methods Engin., 49, No. 3, 421-438 (2000).
9. Helmig R., Multiphase Flow and Transport Processes in the Subsurface, A Contribution to the Modeling of Hydrosystems, Springer-Verl., Berlin (1997).
© 2020 S. P. Stepanov, A. V. Grigoriev, N. M. Afanasyeva
10. Bai M. and Roegiers J.-C., "Fluid flow and heat flow in deformable fractured porous media," Int. J. Eng. Sci., 32, No. 10, 1615-1633 (1994).
11. Vasilyeva M. V. and Prokopiev G. A., "Numerical solution to the problem of two-phase filtration with heterogeneous coefficients by the finite element method," Mat. Zamet. SVFU, 24, No. 2, 46-62 (2017). doi: https://doi.org/10.25587/SVFU.2017.2.9245.
12. Langtangen H. P., "A FEniCS tutorial," in: Automated Solution of Differential Equations by the Finite Element Method, pp. 1-73, Springer-Verl., Berlin; Heidelberg (2012).
13. Samarskii A. A., The Theory of Difference Schemes, Marcel Dekker, New York; Basel (2001).
14. Barenblatt G. I., Zheltov I. P., and Kochina I. N., "Basic concepts in the theory of seepage of homogeneous liquids in fissured rocks," J. Appl. Math. Mech., 24, No. 5, 1286-1303 (1960).
15. Rathfelder K. and Abriola L. M., "Mass conservative numerical solutions of the head-based Richards equation," Water Resources Res., 30, No. 9, 2579-2586 (1994).
16. Ross P. J., "Efficient numerical methods for infiltration using Richards equation," Water Resources Res., 26, No. 2, 279-290 (1990).
17. Celia M. A., Bouloutas E. T., and Zarba R. L., "A general mass-conservative numerical solution for the unsaturated flow equation," Water Resources Res., 26, No. 7, 1483-1496 (1990).
18. Van Genuchten M. T., "A closed-form equation for predicting the hydraulic conductivity of unsaturated soils 1," Soil Sci. Soc. America J., 44, No. 5, 892-898 (1980).
19. Vasiliev V. I., Maksimov A. M., Petrov E. E., and Tsypkin G. G., Heat and Mass Transfer in Freezing and Thawing Soils, Nauka, Moscow (1996).
20. Droniou J., Eymard R., and Guichard C., "Uniform-in-time convergence of numerical schemes for Richards' and Stefan's models," in: Finite Volumes for Complex Applications VII: Methods and Theoretical Aspects, pp. 247-254, Springer, Cham (2014).
Submitted November 19, 2019 Revised February 13, 2020 Accepted April 30, 2020
Sergei P. Stepanov
Ammosov North-Eastern Federal University, Institute of Mathematics and Informatics, 42 Kulakovsky Street, Yakutsk 677000, Russia cepe2a@inbox.ru Aleksandr V. Grigoriev
Ammosov North-Eastern Federal University, Institute of Mathematics and Informatics, 42 Kulakovsky Street, Yakutsk 677000, Russia re5itsme@gmail•com
Nadezhda M. Afanasyeva Ammosov North-Eastern Federal University, Institute of Mathematics and Informatics, 42 Kulakovsky Street, Yakutsk 677000, Russia afanasieva.nm@gmail.com
