Determination of Thermophysical Parameters of the Soil According to Dynamic Data on its Temperature

Full Text

ВВЕДЕНИЕ

Управление тепловым режимом почв является важнейшим условием интенсивного сельскохозяйственного производства в условиях открытого и закрытого грунта. Оно основано на использовании тепловых свойств почвы, которые в значительной степени зависят от содержания твердой фазы, включая содержание органических веществ, минеральный состав почвы, содержание воды в почве и пористость. Тепловые свойства почвы, такие как теплопроводность, теплоемкость и температуропроводность, в основном контролируют температуру и тепловой поток почвы. Определение этих параметров важно для понимания поведения теплового режима почвы и управления температурой почвы в масштабе поля [1, 7, 10, 23].

Существует несколько методов моделирования тепловых свойств почвы по наблюдаемой температуре почвы [1–4, 7, 9, 18, 20, 21]. Большинство моделей основано на решении одномерного уравнения теплопроводности с постоянной диффузией.

Цель работы – определение тепловых свойств почв на основании экспериментальных послойных данных о температуре почвы в летний период и выбор расчетных зависимостей для вычисления коэффициента температуропроводности и потока тепла в почву.

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

Известно, что одномерное распространение тепла в почве описывается классическим уравнением теплопроводности, которое имеет вид [3, 6, 7, 20, 23]:

$\frac{дT}{дt}=\mathrm{\kappa }\frac{{\partial }^2\mathrm{T}}{\partial {\mathrm{ɀ}}^2}\left(\mathrm{\kappa } =\frac{\mathrm{\lambda }}{{\mathrm{C}}_{\mathrm{v}}}\right), \left\{0<\mathrm{ɀ} < L \& \infty ; t >0\right\}$                                                 (1)

где T(ɀ, t) – температура почвы в точке ɀ в момент времени t; Ƙ и λ – соответственно коэффициенты температуропроводности (м²/с) и теплопроводности почвы, Вт/(м °C); C_v – объемная теплоемкость, Дж/(м³ °C). L – глубина почвы (м) начиная с которой T(ɀ, t) = const или ∂T(L, t)/∂ɀ = 0.

Имеется ряд работ [3, 7, 21, 23], в которых рассматривается решения уравнения теплопроводности при различных краевых условиях.

Для исследования переноса тепла в почве необходимо поставить начальные и граничные условия для уравнение (1). Известно, что влияние начального условия не сказывается на распределения температуры почвы в момент наблюдения. Особенно если при периодической постановке задачи начальное условие отсутствует (так называемые “задачи без начального условия”) [8, 23].

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

$T(0,t)$ $= T_0+\underset{}{\overset{}{\sum _j^m{T}_j·\cos (j\omega t + {\epsilon }_j)}}$                                                (2)

где T₀ – среднесуточная (или годовая) температура деятельной поверхности почвы; m – число гармоники; T_j – амплитуды колебаний температуры поверхности почвы; ω = 2π/τ₀ – круговая суточная частота; τ₀ – период (длина) волны, выраженный в сутках или годах; ε_j – сдвиги фаз, зависящие от начала отсчета времени; j – номер гармоники.

Обычно почва рассматриваются как полуограниченный массив. Тогда, с учетом того, что температура почвы на бесконечности постоянна, нижнее граничное условие 1-го рода имеет вид:

при $ɀ\stackrel{}{\to }\infty : T(\infty ,t) =T_0$                                                              (3)

Если температурные колебания быстро затухают с глубиной, и, начиная с некоторой глубины ɀ ≥ L, температура почвы на расчетном интервале времени практически не меняется, то вместо уравнения (3) на нижней границе следует использовать условия [5, 6, 21, 23]:

при $ɀ \stackrel{}{\to }L : \partial T(L,t)/\partial ɀ = 0$                                                       (4)

Решение уравнения (1) при краевых условиях (2) и (3) в безразмерных переменных имеет вид [3, 6]:

$T(y,\tau ) =$$T_{0 }+ \sum _{j=1}^m{\mathrm{\Phi }}_{\mathrm{j}}(\mathrm{y},{\mathrm{b}}_{\mathrm{j}}) · \cos \left[\mathrm{j}\overline{\mathrm{\omega }}\mathrm{\tau } + {\mathrm{\alpha }}_{\mathrm{j}}(\mathrm{y},{\mathrm{b}}_{\mathrm{j}})\right]$                                    (5)  

где

$y = ɀ /L,\tau = \mathrm{\kappa }t /L^2, b_j =\sqrt{j\overline{\omega }/2} , \overline{\mathrm{\omega }} = \mathrm{\omega }{L}^2/\mathrm{\kappa }$     и

$\begin{array}{c}\begin{array}{c}\begin{array}{c}{\mathrm{\Phi }}_j (y, b_j) = T_j · e^{-b_{j}y} , {\alpha }_j (y, b_j) = {\epsilon }_j - {\psi }_j (y ,b_j)’\\ {\omega }_j(у ,b ) = b_j-y.\end{array}\end{array}\end{array}$                     (6)

Эта задача изучалась Фурье; впервые она была применена Кельвином для определения хода температуры в почве Эдинбурга [3, с. 86].

Как отмечено в работах [5, 23] при выполнении практических расчетов нельзя задавать в качестве исходных данных значения температуры почвы на бесконечности, так как они неизвестны. Измерить его невозможно, поэтому вместо Τ(∞, t) следует задать температуру на некоторой глубине L, начиная с которой при ɀ ≥ L величина T(ɀ, t) = const или ∂T(L, t)/∂ɀ = 0.

Таким образом, условие (4) более соответствует реальным условиям, чем условие (3). Поэтому следует также рассмотреть краевую задачу (1), (2) и (4).

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

$y = \frac{ɀ}{L}, \tau = \frac{\mathrm{\kappa }}{L^2}t, \overline{\mathrm{\omega }} =\mathrm{\omega }\frac{L^2}{\mathrm{o}},$ u(y, τ) = T(y, τ) – T₀                                (7)

из уравнений (1), (2) и (4) получим следующую краевую задачу с безразмерными параметрам, новыми переменными и функцией:

$\frac{\partial u}{\partial \mathrm{\tau }} = \frac{{\partial }^{2}u}{\partial y^2},$ {0 < y < 1 ;  τ  > 0},                                                          (8)

$u (0, \tau ) = \sum _{j=1}^m{T}_j · \cos (j\overline{\omega }\tau + {\epsilon }_j)$,                                                          (9)

$u^{\text{'}}(1, \tau ) = 0$,                                                                                    (10)

Решение задачи (8)–(10) может быть представлено в виде суммы

$u (y,т) = u_1 (y,t) + u_2 (y,t) + ... + u_m(y, \tau ) = \sum _{j=1}^m{u}_j(y, \tau )$,                  (11)

где u_j (y, τ) решение следующих частных краевых задач:

$\frac{\partial u}{\partial \mathrm{\tau }} = \frac{{\partial }^{2}u}{\partial y^2},$  {0 < y < 1 ;  τ  > 0},  j = 1,2,..., m,                                (12)

$\left\{\begin{array}{l}{u}_j(0,\tau ) =T_j · \cos (j\overline{\mathrm{\omega }}\tau +{\epsilon }_j)\\ u_j^{\text{'}}(1,\tau ) = 0\end{array}\right.$                                                            (13)

Из линейности уравнение теплопроводности следует, что действительная и мнимая части некоторого комплексного решения уравнения теплопроводности каждая в отдельности удовлетворяют уравнению (12) [8, с. 239].

Если комплексное решение уравнения теплопроводности, т. е.

$\hat{u}j(y,\tau ) =\mathrm{Re}\left[\hat{u}j(y,\tau )\right]+i · \mathrm{Im}\left[\hat{\mathit{u}}j(y,\mathrm{\tau })\right] ={ u}_{1j}(y,\mathrm{\tau }) +\mathit{ }i ·u_{2j}(y\mathit{,}\tau )$                (14)

удовлетворяет уравнению (12) с граничными условиями (13), то его действительная Re[û_j (y, τ)] = = u_1j (y, τ) и мнимая Im[û_j (y, τ)] = u_2j (y, τ) части удовлетворяют следующим условиям:

$\left\{\begin{array}{l}{u}_{1j} (y, \tau )= T_j · \cos (j\overline{\mathrm{\omega }}\tau + {\epsilon }_j)\\ u_{1j}^{\text{'}}(1,\tau ) = 0\end{array}\right.$   и

$\left\{\begin{array}{l}{u}_{2j} (y, \tau )= T_j · \sin (j\overline{\mathrm{\omega }}\tau + {\epsilon }_j)\\ u_{2j}^{\text{'}}(1,\tau ) = 0\end{array}\right.$                      $(i = 1,2, ....., m)$             (15)

Рассмотрим задачу:

$\begin{array}{l}\frac{\partial \hat{u}j}{\partial \tau }=\frac{{\partial }^2\hat{u}j}{\partial y^2}, \{0<y<1; \tau >0\}\\ \left\{\begin{array}{l}{\hat{u}}_j (y, \tau )= T_j · \exp \left[-i(j\overline{\mathrm{\omega }}\tau + {\epsilon }_j)\right]\\ {\hat{u}}_j^{\text{'}}(1,\tau ) = 0\end{array}\right.\end{array}$                                                (16)

Будем искать ее решение в форме

$\begin{array}{c}\hat{u}j (y, \tau ) = Y_j \left(y\right) · \exp [-i (j\omega \tau + {\epsilon }_j\left)\right],\\ \left\{0<y<1; \tau >0\right\}.\end{array}$                                               (17)

где Y_j(y) – новые функции относительно переменной у, которых следует определить (j = 1,2,…, m).

Можно показать, что это решение имеет вид:

$\hat{u}j (y, \tau ) =T_j$ $\frac{ch \left[\right(1 - i ) · b_j (1 - у\left)\right]}{\mathrm{ch}{\left[(1-i) · {\mathrm{b}}_j\right]}_{{{}_{}}_{}}}{e}^{-i(j\overline{\omega }\tau +{\epsilon }_j)}$   

$= T_j · \left[Y_{1j}(b_j,y) +i{Y}_{2j}(b_j,y)\right]e^{-i(j\overline{\omega }\tau +{\epsilon }_j)}$                                              (18)

Выделяя вещественную часть функции û_j (y, τ), т. е. Re[û_j (y, τ)] = u_j (y, τ), находим решение исходной задачи (16) без начальных условий в виде:

$u_j (y, \tau ) = T_j · \left[Y_{1j}(b_j,y) \cos (j\overline{\omega }\tau + {\epsilon }_j) +Y_{2j}(b_j,y) \sin j\overline{\omega }\tau + {\epsilon }_j)\right]$              (19)

Учитывая решение (11) и подстановку u(y,τ) = T(y, τ) – T₀ в выражение (19), окончательно имеем решение исходной задачи (1), (2) и (4) с безразмерными параметрами и переменными в следующем виде:

$T(y, \tau ) = T_0 +\sum _{j=1}^m{T}_j · \left[Y_{1j}(b_j,y) \cos (j\overline{\omega }\tau + {\epsilon }_j) +Y_{2j}(b_j,y) \sin j\overline{\omega }\tau + {\epsilon }_j)\right]$    (20)

После несложных преобразований в решении (18) можно найти явные выражения функции  Y_1j (b_j, y) и Y_2j (b_j,v)  в следующем виде:

$\begin{array}{c}\begin{array}{c}{Y}_{1j} (b_j, y) =\frac{\mathrm{ch}\left[b_j\right(2-y\left)\right] \cos \left(b_{j}y\right) + \mathrm{ch} \left(b_{j}y\right) \cos \left[b_j\right(2-y\left)\right]}{\mathrm{ch}\left(2b_j\right) + \cos \left(2b_j\right)}\\ Y_{2j} (b_j, y) =\frac{\mathrm{ch}\left[b_j\right(2-y\left)\right] \sin \left(b_{j}y\right) + \mathrm{ch} \left(b_{j}y\right) \sin \left[b_j\right(2-y\left)\right]}{\mathrm{ch}\left(2b_j\right) + \cos \left(2b_j\right)}\end{array}\end{array}$                  (21)

Решение (20) запишем в более компактном виде

$T(y, \tau ) = T_0 +\sum _{j=1}^m{\mathrm{\Phi }}_{aj} ·(b_j,y) · \cos \left[j\overline{\omega }\tau + {\alpha }_j(b_j,y)\right]$                                  (22)

где

${\mathrm{\Phi }}_{aj}(b_j,y) = T_j · \sqrt{\frac{\mathrm{ch} \left(d_j\right) + \cos \left(d_j\right)}{\mathrm{ch} \left(2b_j\right) + \cos \left(2b_j\right)}}$                                                 

$b_j=L\sqrt{j\frac{\mathrm{\omega }}{2\mathrm{\kappa }}}, d_j = 2b_j(1 - y),$

${\alpha }_j (b_j,y )= {\epsilon }_j — arc\tan$ $\left[\frac{\mathrm{sh} \left(q_j\right) \sin \left(b_{j}y\right) + \mathrm{sh} (b_{j}y ) \sin \left(q_j\right)}{\mathrm{ch} \left(q_j\right) \cos \left(b_{j}y\right) + \mathrm{ch} (b_{j}y ) \cos \left(q_j\right)}\right]$             (23)

$q_j = b_j (2 - y),$

и $\mathrm{ch} \left(ɀ\right) = \left[\exp \right(ɀ) + \exp (-ɀ \left)\right]/2,$  $\mathrm{sh} \left(ɀ\right) = \left[\exp \right(ɀ) + \exp (-ɀ \left)\right]/2,$ – гиперболический косинус и синус соответственно.

Решение обратной задачи теплопереноса в почве. На основе аналитических решений одномерного уравнения теплообмена в почве разными авторами выведены формулы для определения коэффициента температуропроводности [1–4, 6–9, 19, 21, 22]. Эти методы основаны на измерении температуры дневной поверхности и с глубиной.

Классические методы. Ниже приведены классические алгоритмы, основанные на решении обратных задач уравнения теплопереноса, полученного для полуограниченной мощности почвы, когда суточный ход температуры на поверхности почвы представлен одной (М1, М4) и двумя (М2, М3) гармониками.

М1 – амплитудный метод [3]:

${\mathrm{\kappa }}_i = \frac{{\mathrm{\pi }}^{{{{}^{}}^{}}^{}}}{{\mathrm{\tau }}_0} \frac{{({\mathrm{ɀ}}_{\mathrm{i}+1}- {\mathrm{ɀ}}_{\mathrm{i}})}^2,}{{\ln }^2\left[{\displaystyle \frac{T_{\max }\mathit{ }\left({\mathrm{ɀ}}_{\mathrm{i}}\right) - T_{\min }\mathit{ }\left({\mathrm{ɀ}}_{\mathrm{i}}\right)}{{\mathrm{T}}_{\max }\mathit{ }\left({\mathrm{ɀ}}_{\mathrm{i}+1}\right) - {\mathrm{T}}_{\min }\mathit{ }\left({\mathrm{ɀ}}_{\mathrm{i}+1}\right)}}\right]}$                                        (24)

М2 – арктангенсный метод [2]:

${\mathrm{\kappa }}_i = \frac{{\mathrm{\pi }}^{{{{}^{}}^{}}^{}}}{{\mathrm{\tau }}_0} \frac{{({\mathrm{ɀ}}_{\mathrm{i}+1}- {\mathrm{ɀ}}_{\mathrm{i}})}^2,}{arc{\tan }^2{\displaystyle \frac{\left[{\displaystyle T_1\mathit{ }\left({\mathrm{ɀ}}_{\mathrm{i}}\right) - T_3\mathit{ }\left({\mathrm{ɀ}}_{\mathrm{i}}\right)}\right]\left[T_2\mathit{ }\left({\mathrm{ɀ}}_{\mathrm{i}+1}\right) - T_4\mathit{ }\left({\mathrm{ɀ}}_{\mathrm{i}+1}\right)\right] - \left[T_2\mathit{ }\left({\mathrm{ɀ}}_{\mathrm{i}}\right) - T_4\mathit{ }\left({\mathrm{ɀ}}_{\mathrm{i}}\right)\right]{\left[T_1\mathit{ }\left({\mathrm{ɀ}}_{\mathrm{i}+1}\right) - T_3\mathit{ }\left({\mathrm{ɀ}}_{\mathrm{i}+1}\right)\right]}^{{{}^{}}^{}}}{\left[T_1\mathit{ }\left({\mathrm{ɀ}}_i\right) - {\mathrm{T}}_3\mathit{ }\left({\mathrm{ɀ}}_{\mathrm{i}}\right)\right]\left[T_1\mathit{ }\left({\mathrm{ɀ}}_{\mathrm{i}+1}\right) - T_3\mathit{ }\left({\mathrm{ɀ}}_{\mathrm{i}+1}\right)\right] + \left[T_2\mathit{ }\left({\mathrm{ɀ}}_i\right) - {\mathrm{T}}_4\mathit{ }\left({\mathrm{ɀ}}_{\mathrm{i}}\right)\right]{\left[T_2\mathit{ }\left({\mathrm{ɀ}}_{\mathrm{i}+1}\right) - T_4\mathit{ }\left({\mathrm{ɀ}}_{\mathrm{i}+1}\right)\right]}^{{{{}^{}}^{}}^{}} }}}$ (25)

М3 – логарифмический метод [4]:

${\mathrm{\kappa }}_i =\frac{4\pi ({ɀ}_{i+1}- {ɀ}_i)}{{\tau }_0}{\ln }^{-2}\times \frac{{\left[T_1{ɀ}_i) - T_3({ɀ}_i){]}^2 + [T_2\left({ɀ}_i\right) - T_4\left({ɀ}_i\right)\right]}^2}{{\left[T_1{ɀ}_{i+1}) - T_3({ɀ}_{i+1}){]}^2 + [T_2\left({ɀ}_{i+1}\right) - T_4\left({ɀ}_{i+1}\right)\right]}^2}$            (26)

М4 – метод фазового сдвига [8]:

${\mathrm{\kappa }}_i =\frac{\pi }{{\tau }_0}{\left(\frac{{ɀ}_{i+1}- {ɀ}_i}{{\epsilon }_i - {\epsilon }_{i+1}}\right)}^2$                                                                         (27)

где τ₀ в алгоритмах (24)–(27) период тепловой волны (например, 24 ч для ежедневных наблюдений); T_min (ɀ) и T_max (ɀ) – минимальное и максимальное значения температуры во время измерения на глубинах ɀ = ɀ_i и ɀ = ɀ_i+1 соответственно; T_i(ɀ) и T_i+1(ɀ) в алгоритмах (25) и (26) являются значениями температуры почвы на глубинах ɀ = ɀ_i и ɀ = ɀ_i+1 соответственно в момент времени t_i = i·τ₀/4(i = 1, 2, 3, 4) (для настоящего примера τ₀ = 24 ч и t₁ = 6, t₂ = 12, t₃ = 18 и t₄ = 24 ч); ε_i – в алгоритме (27) – начальная фаза на глубине ɀ_i.

Усовершенствованные методы. В отличие от вышеперечисленных алгоритмов продолжили разработку новых методов, учитывающих реальный процесс теплообмена в почве. Для этого использовали решение уравнения теплопереноса, полученное для ограниченной и полуограниченной мощности почвы [6, 18, 21, 23]. Предлагаемые методы расчета температуропроводности почвы представлены ниже.

М5 – предлагаемый метод 1 [23]:

${\mathrm{\kappa }}_i =\frac{\pi }{{\tau }_0}$ $\frac{{\left(2{ɀ}_i\right)}^2}{{\ln }^2\left\{\sum _{j=1}^4 {\left[T({ɀ}_i,t_i) - T({ɀ}_i,t_{i+2})\right]}^2/4T_1^2\right\}}$                                  (28)

М6 – предлагаемый метод 2 [23]:

$M_{\exp }^{m=1}(b;y) =$$\frac{{\displaystyle \sum _{i=1}^2}{\left[T(y_i,t_i) - T(y_i,t_{j+2})\right]}^2}{4{T_1^2}_{{{{{}_{}}_{}}_{}}_{}}}$ $=\frac{\mathrm{ch} [2b (1 — y\left)\right] + \cos \left[2b\right(1 — y\left)\right]}{\mathrm{ch}\left(2b\right) + \cos \left(2b\right)}$$=M_{cal}^{m=1}(b;y)$     (29)

где T(y_i, t_j) в алгоритме (29) – значения температуры почвы на глубине ɀ = ɀ_i, в момент времени t_j = j·τ₀/4 (j = 1, 2, 3, 4) (для настоящего примера τ₀ = 24 ч и t₁ = 6, t₂ = 12, t₃ = 18 и t₄ = 24 ч).

Определение Ƙ по формуле (29) осуществляется методом подбора на ЭВМ значения параметра b^* из условия совпадения значения левой части $M_{\exp }^{m=1}(b;y)$ и правой части $M_{\exp }^{m=1}(b;y)$, рассчитанных по исходным данным.

Из следующего равенства:

${\mathrm{\kappa }}_i =\frac{\pi }{{\tau }_0} · {\left(\frac{L}{b_i^{*}}\right)}^2$,                                                                               (30)

находим значение коэффициента температуропроводности Ƙ_i на глубине ɀ = ɀ_i, где b* – значения параметра b, при которых

$\left|M_{\exp }^{m=1}(b;y) - M_{\mathrm{cal}}^{m=1}(b;y)\right|\le {10}^{-4}$.

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

Используя решения (5) с (6) и (22) с (23) при m = 2 можно вывести формулу для определения коэффициента температуропроводности κ для произвольного периода τ₀ и безразмерной глубины y.

Для этого необходимо знать распределение температуры в почвенном слое [0, L] для восьми моментов времени на расчетном интервале времени τ₀. Предлагаемые методы расчета температуропроводности почвы представлены ниже.

Используя решение (22) для двух гармоник, т. е. m = 2, сначала для произвольной безразмерной глубины y и времени t_j = j·τ₀/8 (j = 1, 2, 3,…, 8), записываются следующие восемь уравнений:

$T(y,t_i) = T_0 + {\mathrm{Ф}}_1(y,b ) · \cos$$\left(\frac{\pi }{4}i + {\alpha }_1\right)$$+ {\mathrm{Ф}}_2(y,b_2 ) · \cos$$\left(\frac{\pi }{2}i + {\alpha }_2\right),$ $(i =\overline{1,8})$         (31)

На самом деле:

для j = 1 имеем:

$T(y,t_1) = T_0 + {\mathrm{Ф}}_1(y,b_1 ) · \cos$$\left(\frac{\pi }{4}·1 + {\alpha }_1\right)$$+ {\mathrm{Ф}}_2(y,b_2 ) · \cos$$\left(\frac{\pi }{2}·1 + {\alpha }_2\right),$

для j = 2 имеем:

$T(y,t_2) = T_0 + {\mathrm{Ф}}_1(y,b_1 ) · \cos$$\left(\frac{\pi }{4}·2 + {\alpha }_1\right)$$+ {\mathrm{Ф}}_2(y,b_2 ) · \cos$$\left(\frac{\pi }{2}·2 + {\alpha }_2\right),$

для j = 8 имеем:

$T(y,t_8) = T_0 + {\mathrm{Ф}}_1(y,b_1 ) · \cos$$\left(\frac{\pi }{4}·2 + {\alpha }_1\right)$$+ {\mathrm{Ф}}_2(y,b_2 ) · \cos$$\left(\frac{\pi }{2}·8 + {\alpha }_2\right),$

так как имеет место

$j\overline{\mathrm{\omega }}{\tau }_i=j · \frac{\mathrm{\omega }{L}^2}{\mathrm{\kappa }} · \frac{\mathrm{\kappa }}{L^2}{t}_i = j · \mathrm{\omega }{t}_i=j · \frac{2\pi }{{\tau }_0}·i·\frac{{\tau }_0}{8} =j·\frac{2\pi }{8}·i$ и $\overline{\mathrm{\omega }}{\tau }_i=\frac{2\pi }{8}· i ·=\frac{\pi }{4} · i$, $2\overline{\mathrm{\omega }}{\tau }_i=\frac{\pi }{2} · i$.

После некоторых преобразований с решениями (31) имеем [7]:

$\sum _{i=1}^4{\left[T(y,t_i) - T(y,t_{i+4})\right]}^2 =8{\mathrm{\Phi }}_1^2({\mathrm{b}}_1, \mathrm{y})$                                                  (32)

С учетом обозначений (6) и (23) для функции Φ₁(b₁, y) в равенстве (32) имеем следующие выражения, которые соответствуют граничным условиям (3) и (4):

$\frac{{\displaystyle \sum _{i=1}^4}{\left[T(y,t_i) - T(y,t_{i+4})\right]}^2}{8{T_1^2}_{{{{}_{}}_{}}_{}}}=e^{-2by}$                                                       (33)

$\frac{{\displaystyle \sum _{j=1}^4 }{\left[T(y,t_i) - T(y,t_{i+4})\right]}^2}{8{T_1^2}_{{{{}_{}}_{}}_{}}}=\frac{\mathrm{ch} [2b (1 — y\left)\right] + \cos \left[2b\right(1 — y\left)\right]}{\mathrm{ch}\left(2b\right) + \cos \left(2b\right)}$             (34)

Так как $b = L\sqrt{\mathrm{\pi }/{\mathrm{\kappa \tau }}_0}$, из выражений (33) и (34) получим более улучшенные методы для вычисления параметра $\mathrm{\kappa }$:

М7 – предлагаемый алгоритм 3 [6, 21, 22, 23]:

${\mathrm{\kappa }}_i =\frac{\pi }{{\tau }_0}$$\frac{{\left(2{ɀ}_i\right)}^2}{{\ln }^2\left\{\sum _{j=1}^4{\left[T({ɀ}_i,t_j) - T({ɀ}_i,t_{j+4})\right]}^4/8T_1^2\right\}}$                                    (35)

М8 – предлагаемый алгоритм 4:

$M_{\exp }^{m=2}(b;y) =$$\frac{\sum _{j=1}^4{\left[T(y_i,t_j) - T(y_i,t_{j+4})\right]}^2}{8T_a^2}=$$\frac{\mathrm{ch} [2b (1 — y\left)\right] + \cos \left[2b\right(1 — y\left)\right]}{\mathrm{ch}\left(2b\right) + \cos \left(2b\right)}=$$M_{\mathrm{cal}}^{m=2}(b;y)$  (36)

где T(y, t_j) – значения температуры при измерении на глубине ɀ_i почвы и в моменты времени t_j и t_j+4, τ₀ – период (длина) температурной волны выраженный в сутках или в годах (например, τ₀ = 24 ч и t₁ = 3, t₂ = 6, …, t₇ = 21, и t₈ = 24 ч).

Определение Ƙ по формуле (36) осуществляется также методом подбора на ЭВМ, как это было описано выше для формулы (29). Сначала находится значения параметра b^* из условия совпадения значения левой $M_{\mathrm{cal}}^{m=2}(b;y)$ и правой частей, рассчитанных по исходным данным, т. е. $M_{\exp }^{m=2}(b;y)$ Далее из формулы (30) находим значение коэффициента температуропроводности κ_i на глубине ɀ = ɀ_i, где b* – значения параметра b, при которых

$\left|M_{\exp }^{m=2}(b*;y_i) - M_{\mathrm{cal}}^{m=2}(b*;y_i)\right|\le {10}^{-4}$ .

В отличие от ранее разработанных методов для определения ϰ требуется знать заранее распределение температуры по времени в почвенном слое [0, L] на произвольной безразмерной глубине y* = ɀ/L и T(y*, t_j) для восьми моментов времени, которое позволяет определить параметр κ с большей точностью по формулам (35) и (36).

Для определения коэффициента температуропроводности κ (с использованием формул (30) и (36)) необходимо знать: Т₁ – амплитуду колебаний температуры деятельной поверхности почвы; τ₀ – период (длину) суточной (годовой) волны, выраженный в сутках или в годах; $T(y^{*},t_j^{*})$ – значения температуры почвенного слоя [0, L] на произвольной глубине y_i = ɀ_i/L для .

Например, если τ^*₀ = 24 ч, то t* = 3, 6, 9,…, 24 ч. Имея эти данные, сначала подсчитываем разницу $\left[T(y_i,t_j) - T(y_i,t_{j+4})\right]$ для всех  $i = \overline{1.8}$. Далее из формулы (36) находим значение коэффициента температуропроводности на глубине y_i* = ɀ*_i /L соответственно по формуле:

${\mathrm{\kappa }}_{{\mathrm{i}}_{{{}_{}}_{}}}^{*}= \frac{\mathrm{\pi }}{{\mathrm{\tau }}_0}\frac{{\left(2{\mathrm{ɀ}}_i^{*}\right)}^2}{{\ln }^2\left\{\sum _{j=1}^4{\left[T({\mathrm{ɀ}}_i^{*}, t_{\mathit{j}}^{\mathit{*}}) - T({\mathrm{ɀ}}_i^{*}, t_{j+4}^{*}\right]}^2/8T_1^2\right\}}$.

Расчет потока тепла в почве. Расчет потока тепла в почве основан на использовании данных об изменении температуры почвы с глубиной и во времени при известных теплофизических характеристиках. Если известно, как изменилась температура почвы за некоторый период времени и определены основные ТФС-теплофизические свойств (объемная теплоемкость – C_v, температуропроводность – κ и параметры дневной деятельности поверхности почв такие, как среднесуточная температура поверхности почвы (T₀), амплитуды (T_i) и фаза (ε_i), то можно рассчитать количество тепла, которое прошло через поверхность почвы и вызвало данное изменение температуры. Важность теплового потока почвы (q) для поверхностного энергетического баланса и исследований испарения способствовала разработке методов оценки и прогнозирования, когда измеренные значения q недоступны. Другой подход для расчета Ƙ и прогнозирования q использует измерения температуры почвы в сочетании с различными решениями уравнений теплового потока почвы. Поскольку температуру почвы легко измерить и часто доступны непрерывные записи данных для нескольких глубин, эти методы дают возможность оценить q, когда прямые измерения отсутствуют [25].

Прямых методов оценки расчетного теплового потока почвы на поверхности почвы (q) не существует. Для расчета q обычно пользуются традиционном калориметрическом методом. Однако в самом калориметрическом методе есть неопределенность [16].

Поток тепла в почву q(ɀ, t), рассчитывается путем подстановки решения уравнения (1) в закон теплопроводности Фурье [3, 7, 10, 15, 16, 19, 25]:

$q(ɀ, t) = - \lambda \frac{\partial T(ɀ, t)}{\partial ɀ}$,                                                                        (37)

а суммарный поток за время от первого наблюдения до данного момента t есть

$Q\left(t\right) = {\int }_0^{t}q\left(t\right)dt$,                                                                                 (38)

Используя (5) для для m = 1, в ɀ = 0 имеем

$q (ɀ = 0,t, m = 1) = q_{ɀ=0}^{m=1}\left(t\right) =$${\overline{)(-\lambda )\frac{\partial T}{\partial ɀ}}}_{ɀ=0}=B_1 · \cos \left[\frac{\pi }{4}+(\mathrm{\omega }\tau +{\epsilon }_1)\right]$(39)

а для m = 2, можно легко показать, что выражение для теплового потока на поверхности почвы ɀ = 0 в момент времени t получим следующее:

$q (ɀ = 0,t, m = 1) = q_{ɀ=0}^{m=1}\left(t\right) =B_1 · \cos \left[\frac{\pi }{4}+(\mathrm{\omega }\tau +{\epsilon }_1)\right] + =B_2 · \cos \left[\frac{\pi }{4}+(\mathrm{\omega }\tau +{\epsilon }_2)\right]$(40)

где $B_1 = T_{a1}{C}_v\sqrt{\mathrm{\omega \kappa }}, B_2 = T_{a2}{C}_v\sqrt{2\mathrm{\omega \kappa }}$.

Необходимо отметить, что при вычислении теплового потока по формулам (39) и (40) следует использовать значения коэффициентов температуропроводности, которые определены соответственно для m = 1 и 2.

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

Можно показать, что для вычисления теплового потока на глубине ɀ = h почвы в момент времени t имеет место:

$\begin{array}{l}q (ɀ = h,t, m) ={\overline{)(-\lambda )\frac{\partial T(\mathrm{ɀ},t)}{\partial ɀ}}}_{ɀ=h}=-C_v\mathrm{\kappa }{\overline{)\frac{\partial \mathrm{T}(\mathrm{ɀ},t)}{\partial \mathrm{ɀ}}}}_{\mathrm{ɀ}=\mathrm{h}}=\\ =-\lambda \frac{1}{L}\sum _{j=1}^m\frac{T_j · b_j}{∆\left(b_j\right)} · \begin{array}{c}\left\{\left[\prod _{1j}(b_j,y_h)\cos (j\mathrm{\omega }t + {\epsilon }_j)+\prod _{2j}(b_j,y_h)\sin (j\mathrm{\omega }t + {\epsilon }_j)\right]\right\}\end{array}\end{array}$(41)

где

$\begin{array}{l}\prod _{1j}(b_j,y_h) = \mathrm{sh}\left(b_j{y}_h\right) · \cos \left[b_j\mathit{(}\mathit{2}\mathit{-}{y}_h\mathit{)}\right]-\mathrm{sh}\left[b_j(2-y_h)\right]\cos (b_j-y_h)+\mathrm{ch}\left(b{y}_h\right) · \sin \left[b_j (2 - y_h)\right]-\mathrm{ch}\left[b_j(2 - y_h)\right] · \sin \left(b{y}_h\right),\\ \prod _{2j}(b_j,y_h) = \mathrm{sh}\left[b_{\mathit{j}}-(2-y_h)\right]\cos \left(b_j{y}_h\right)-\mathrm{ch}\left[b_j(2 - y_h)\right] \sin \left(b_j{y}_h\right)+\mathrm{ch}\left(b_j{y}_h\right) · \sin \left[b_j(2 - y_h)\right]- \mathrm{sh}\left(b_j{y}_h\right)\cos \left[b_j(2 - y_h)\right],\end{array}$(42)

$∆(b_j ) = \mathrm{ch}\left(2b_j\right) + \cos \left(2b_j\right),$

$b_j=L\sqrt{j\frac{\mathrm{\omega }}{2\mathrm{\kappa }}}, y_h=\frac{h}{L}, 0\le h \le L$.

Используя равенство (41)–(42), можно получить выражение для теплового потока на поверхности почвы (ɀ = 0) в момент времени t:

$q (ɀ = 0,t,m ) = -\lambda$$\frac{\partial T\left(ɀ=0,t\right)}{\partial ɀ}=-\lambda \frac{1}{L}\sum _{j=1}^m\frac{T_j ·b_j}{∆\left(b_j\right)}\left\{\left[\mathrm{sh}\left(2b_j\right) + \sin \left(2b_j\right)\right] \sin (j\omega t + {\epsilon }_j) +\left[\sin \left(2{\mathrm{b}}_{\mathrm{j}}\right)-\mathrm{sh}\left(2b_j\right)\right]\cos (j\mathrm{\omega }t+{\epsilon }_j)\right\}$   (43)

Можно показать, что имеет место:

$\left[\mathrm{sh}\left(2b_j\right) + \sin \left(2b_j\right)\right] \sin (2b_j+{\epsilon }_j)+\left[\sin \left(2b_j\right) - \mathrm{sh}\left(2b_j\right)\right] \cos (j\omega t + {\epsilon }_j)=$

$-\sqrt{2\left[\mathrm{sh}\left(2b_j\right)\cos \left(\frac{\pi }{4}+j\omega t+{\epsilon }_j\right)-\sin \left(2b\right)\cos \left(\frac{\pi }{4}-j\omega t-{\epsilon }_j\right)\right]}$

Тогда окончательно имеем

$q (ɀ = 0,t,m ) = -\lambda$$\frac{\partial T\left(ɀ=0,t\right)}{\partial ɀ}=-\lambda \frac{\sqrt{2\lambda }}{L}\sum _{j=1}^m\frac{T_j ·b_j}{∆\left(b_j\right)}$$\left[\mathrm{sh}\left(2b_j\right)\cos \left(\frac{\pi }{4}+j\omega t+{\epsilon }_j\right)-\sin \left(2b\right)\cos \left(\frac{\pi }{4}-j\omega t-{\epsilon }_j\right)\right]$(44)

Для m = 1 из решения (44) имеем:

$q_{z=0}^{m=1}\left(t\right)={\mathrm{\Phi }}_1\left\{\mathrm{sh}\left(2\mathrm{b}\right)\cos \left(\frac{\mathrm{\pi }}{4}+\mathrm{\omega t}+\mathrm{\epsilon }\right)-\sin \left(2\mathrm{b}\right)\cos \left(\frac{\mathrm{\pi }}{4}-\mathrm{j\omega t}-\mathrm{\epsilon }\right)\right\}$

${\mathrm{\Phi }}_1=\mathrm{B} · \frac{{\mathrm{T}}_1}{∆\left(b\right)}, \mathrm{B}=\mathrm{\lambda }\sqrt{\frac{\mathrm{\omega }}{\mathrm{\kappa }}}$,                                                                  (45)

Для m = 2 из решения (44) имеем:

$q_{z=0}^{m=2}\left(t\right)={\mathrm{\Phi }}_1\left\{\mathrm{sh}\left(2{\mathrm{b}}_1\right)\cos \left(\frac{\mathrm{\pi }}{4}+\mathrm{\omega t}+{\mathrm{\epsilon }}_1\right)-\sin \left(2{\mathrm{b}}_1\right)\cos \left(\frac{\mathrm{\pi }}{4}-\mathrm{j\omega t}-{\mathrm{\epsilon }}_1\right)\right\}+$

${\mathrm{\Phi }}_2\left\{\mathrm{sh}\left(2{\mathrm{b}}_2\right)\cos \left(\frac{\mathrm{\pi }}{4}+\mathrm{\omega t}+{\mathrm{\epsilon }}_2\right)-\sin \left(2{\mathrm{b}}_2\right)\cos \left(\frac{\mathrm{\pi }}{4}-\mathrm{j\omega t}-{\mathrm{\epsilon }}_2\right)\right\},$                     (46)

${\mathrm{\Phi }}_1=\frac{{\mathrm{T}}_1}{∆\left(b_1\right)}\left(\mathrm{\lambda }\sqrt{\frac{\mathrm{\omega }}{\mathrm{\kappa }}}\right), {\mathrm{\Phi }}_2=\frac{{\mathrm{T}}_2}{∆\left({\mathrm{b}}_2\right)}\left(\mathrm{\lambda }\sqrt{\frac{2\mathrm{\omega }}{\mathrm{\kappa }}}\right),$

$b_1=L\sqrt{\frac{\mathrm{\omega }}{2\mathrm{\kappa }}}, b_2=L\sqrt{2\frac{\mathrm{\omega }}{2\mathrm{\kappa }}}$

$∆\left(b_1\right) = ch\left(2b_1\right) + \cos \left(2b_1\right),$

$∆\left(b_2\right) = ch\left(2b2\right) + \cos \left(2b_2\right), \omega =\frac{2\pi }{{\tau }_0}.$

При вычислении теплового потока по формулам (45) и (46) следует использовать значения коэффициентов температуропроводности (κ) и параметры дневной деятельности поверхности почв (такие как среднесуточная температура поверхности почвы T₀, амплитуды T_i и фазы ε_i), которые определены соответственно для гармоники m = 1 и 2.

ОБЪЕКТЫ И МЕТОДЫ

Исследование проводили в Центре сельскохозяйственных исследований и прикладных разработок Игдырского университета. В Игдире лето жаркое, сухое и ясное, а зима короткая, очень холодная, снежная и местами облачная. В течение года температура обычно колеблется от –6 до 35°C, редко бывает ниже –11°C или выше 39°C. Наибольшее количество осадков наблюдается в мае, а наименьшее – в августе, среднегодовое количество осадков составляет 254.2 мм, испарение 1094.9 мм [11]. Некоторые свойства глеевой пойменной почвы (Calcaric Gleyic Pantofluvic Fluvisol) опытного участка приведены в табл. 1.

 

Таблица 1. Некоторые физические и химические свойства почвы

Глубина, см	ρb, кг/м³	θ, м³/м³	EC, дС/м	OM, %	Глина	Ил	Песок	Грануломет- рический состав
					%
0–10	974.3	0.1449	0.3877	1.40	28.58	32.07	39.75	CL*
10‒20	1023.3	0.1521	0.4126	1.62	27.36	31.48	42.82	CL
20–25	1105.4	0.1717	0.5417	2.35	25.55	29.52	45.13	L
25–30	1245.7	0.1806	0.5524	2.81	20.41	23.22	56.38	SCL
30–35	1158.1	0.1907	0.5507	3.07	17.44	19.63	66.37	SL
35–40	1216.5	0.1925	0.5121	2.66	14.28	27.11	58.68	SL
0 –40	1120.6	0.1721	0.4928	2.32	22.27	27.17	51.52

Примечание. ρ_b – объемная плотность; θ – объемная влажность; EC – электропроводность почвы; ОМ –органическое вещество почвы; CL – глинистый суглинок; L – суглинок; SCl – супесчаный глинистый суглинок; SL – супесчаный суглинок.

 

В исследовании использовались датчики измерения и регистрации температуры Elitech RC-4 с инструментальной погрешностью ±0.5°C в диапазоне измерения температуры –20…+40°C [26]. Для определения тепловых свойств почв измеряли температуру профиля с 01.06.2020 по 31.08.2020. Использовали средние значения температуры в эти даты. Датчики располагались в почвенном профиле на глубинах 0, 5, 10, 15, 20 и 40 см.

Почвенные анализы. Для определения свойств почвы были отбирали образцы почвы с земель, на которых размещены датчики. В образцах нарушенной почвы определяли состав почвы, содержание влаги в почве и органическое вещество, а в образцах ненарушенной почвы определяли объемную плотность почвы. Текстура почвы, содержание влаги, органического вещества и плотность почвы определяли согласно [13, 14, 17, 24].

Датчики были запрограммированы на получение почасовых данных о температуре, а данные о температуре собирали с датчиков через компьютер в течение летного сезона. Тепловые свойства почв рассчитывали по данным, полученным с датчиков, размещенных на разной глубине в летный период.

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

Объемная теплоемкость. Объемную теплоемкость почвы рассчитывали по общепринятой формуле [10]:

$\begin{array}{l}{C}_v — C_{m,s} {\rho }_b + C_{v,w}· \mathrm{\theta } =\\ \left[C_{m,org}· \frac{m_{org}}{m}+C_{m,\min } · \left(1-\frac{m_{org}}{m}\right)\right]·{\rho }_b + C_{v,w}· \mathrm{\theta },\end{array}$                           (47)

где C_m,s – удельная теплоемкость твердой части почвы, Дж/(кг °C); ρ_b – плотность почвы, кг/м³; C_v,w – объемная теплоемкость почвенной влаги, равная 4186.6 кДж/(м³ °C); θ – объемное содержание влаги, м³/м³; C_m,org и C_m,min – удельная теплоемкость органических и минеральных компонентов твердой фазы почвы, соответственно, Дж/(кг °C); m_org – масса органического вещества почвы, кг; m –масса почвы, кг; $\frac{m_{org}}{m}$ – содержание органического вещества в почве, %.

Другие теплофизические параметры: теплопроводность (λ), глубину затухания (d) и теплоусвояемости тепла (e), рассчитывали по формулам соответственно [3, 7, 10, 22]:

$\lambda = \mathrm{\kappa }{С}_v, d =\sqrt{{\tau }_{0 }\frac{\mathrm{\kappa }}{\pi }}, e = C_v\sqrt{\mathrm{\kappa }}$                                                (48)

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

Среднесуточную температуру поверхности почвы (T₀), амплитуды (T_i) и фаза (ε_i) рассчитывали по уравнениям:

$T_0=\frac{1}{N}\sum _{i=1}^{N}T(0.t_i),$

$A_j=\frac{2}{N}\sum _{i=1}^{N}T(0.t_i)\cos \left(j\omega t_i\right),$                                                       (49)

$B_j=\frac{2}{N}\sum _{i=1}^{N}T(0.t_i)\sin \left(j\omega t_i\right),$ 

$T_j=\sqrt{A_j^2} + B_j^2,$

${\epsilon }_j=\left\{\begin{array}{l}arc\tan (-\frac{B_j}{A_j} ), A_j > 0\\ \frac{\pi }{2}, A_j= 0,(j=1,2,..., m)\\ \pi -arc\tan \left(\frac{B_j}{A_j}\right), A_j < 0\end{array}\right.$

где T(0, t_i) – температура поверхности почвы в момент времени t_i; N – число данных о температуре в сутке, т. е. N = 24.

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

$r=\frac{\sum _{i=n}^{i=n}(T_i-{\overline{T}}_i)({\tilde{T}}_i-{\overline{\tilde{T}}}_i)}{\sqrt{\sum _{i=n }^{i=n}{(T_i-{\overline{T}}_i)}^2 \sum _{i=n }^{i=n}({\stackrel{~ }{T}}_i-{\overline{\tilde{T}}}_{})}}$                                                               

${\sigma }_{\frac{T}{t}}=\sqrt{\frac{\sum _{i=n }^{i=n}{(T_i-{\tilde{T}}_i)}^2}{n-p-1}}$                                                              (50)

где $T_i$ и $\widetilde{T_i}$ – соответственно наблюдаемое и расчетное значения температуры, °C;  $\overline{T}$ – средние значения $T_i$ ; p – количество оцениваемых параметров в моделях.

РЕЗУЛЬТАТЫ И ОБСУЖДЕНИЕ

Значения удельных теплоемкостей органической и минеральной составляющих исследуемой почвы равны C_m,org = 1925.928 Дж/(кг °C) и C_m,min = 753.624 Дж/(кг °C) соответственно (табл. 2).

 

Таблица 2. Некоторые тепловые свойства исследуемых почв

Глубина, ɀ см	morg /m	Cm,s , Дж/(кг °C)	Cv , кДж/(м³ °C)
0–10	0.014	770.04	1356.91
10–20	0.016	772.62	1427.43
20–25	0.023	781.17	1582.38
25–30	0.028	786.57	1735.96
30–35	0.031	789.61	1712.87
35–40	0.027	784.81	1760.68
0–40	0.0232	780.80	1596.04

Примечание. m_org/m – содержание органического вещества в почве; C_m,s – удельная теплоемкость твердой части почвы; C_v – объемная теплоемкость почвы

 

В результате проведенных анализов проб почв определена средняя плотность почвы опытного участка: ρ_b = 1120.6 кг/м³, содержание органического вещества $\frac{m_{org}}{m} = 0.0232$ и объемная влажность θ = 0.1721 м³/м³.

В этом случае значение удельной теплоемкости (C_m,s) для слоя 0–10 см твердой части почвы, рассчитанное по формуле (2), будет следующим:

$C_{m,s} = C_{m, org}\frac{m_{org}}{m}+ C_{m.min} · \left(1-\frac{m_{org}}{m}\right)=1925.928\frac{\mathrm{Дж}}{\mathrm{кг}°\mathrm{C}}0.014+753.624\frac{\mathrm{Дж}}{\mathrm{кг}°\mathrm{C}}(1-0.014)=$

$= 26.96299\frac{\mathrm{Дж}}{\mathrm{кг}°\mathrm{C}}+ 743.0733\frac{\mathrm{Дж}}{\mathrm{кг}°\mathrm{C}}=770.0363 \frac{\mathrm{Дж}}{\mathrm{кг}°\mathrm{C}}=0.770036\frac{\mathrm{Дж}}{\mathrm{кг}°\mathrm{C}}$.

Используя формулы для вычисления C_v и учитывая, что для слоя 0–10 см почвы θ = 0.1449 м³/м³, C_m,s = 770.0363 Дж/(кг °C), а удельная объемная теплоемкость воды C_v,w = 4186.8 кДж/(м³ °C), рассчитываем значение объемной теплоемкости для слоя почвы 0–10 см следующим образом:

$\begin{array}{l}{C}_v = C_{ms} · {\rho }_b + C_{v,w} · \mathrm{\theta } = 770.0363\frac{\mathrm{Дж}}{\mathrm{кг}°\mathrm{С}} ·974.3\frac{\mathrm{кг}}{{\mathrm{м}}^3}+4186.8\frac{\mathrm{кДж}}{{\mathrm{м}}^3°\mathrm{С}} · 0.1449\frac{{\mathrm{м}}^3}{{\mathrm{м}}^3}=\\ =750246.324 \frac{\mathrm{Дж}}{{\mathrm{м}}^3°\mathrm{С}} + 606.6673 \frac{\mathrm{кДж}}{{\mathrm{м}}^3°\mathrm{С}} = 1356.913644\frac{\mathrm{кДж}}{{\mathrm{м}}^3°\mathrm{С}}\end{array}$

Для всех слоев почвы расчетные значения объемной теплоемкости почвы приведены в табл. 2.

Среднее значения распределения температур за период с 01.06.2020 по 31.08.2020 приведены на рис. 1. С увеличением глубины почвы и времени температура почвы не теряла своей синусоидальности (рис. 1).

 

Рис. 1. Среднесуточные значения температуры почвы T(ɀi, ti) на различных глубинах за период с 01.06.2020 по 31.08.2020.

 

В табл. 3 приведены результаты расчета параметров, а также статистические характеристики аппроксимации исходных данных, рассчитанных по уравнению (49) для m = 1 и 2. Введение второй гармоники позволяет с высокой точностью определять параметры распределения температуры на поверхности почвы.

 

Таблица 3. Параметры (T₀, T_i и ε_i) поверхности почвы

Параметр на поверхности почвы	Число гармоник
	m = 1		m = 2
Среднесуточная температура, °С	T0	24.5478	T0	24.5478
Амплитуда колебаний температуры, °С	T1	5.4270	T2	1.5067
Сдвиг фазы	ε1	2.3607	ε2	–0.4932
Статистические параметры аппроксимации
Коэффициент корреляции Пирсона	r	0.961	r	0.997
Среднеквадратическая ошибка (RMSE)	σ	1.19	σ	0.07

 

Средние значения температуропроводности почвы (Ƙ), теплопроводности (λ), глубины затухания (d) и теплоусвояемости (е), рассчитанные классическими (послойно-амплитудным, арктангенсным, логарифмическим и фазовым) и предложенными (точечными) методам, приведены в табл. 4.

 

Таблица 4. Средние значения тепловых свойств почвы

Формула	Число гармоник	10 –6 · Ƙ, м²/с	λ, Вт/м °С	d, см	e, Вт ч0.5/м² °C
Классические послойные методы
(24)	1	0.7422	1.4882	14.29	22.9163
(25)	2	0.9287	1.5995	15.98	25.6347
(26)	2	0.7643	1.6338	14.50	23.2549
(27)	1	0.8174	1.7612	14.99	24.0490
Предложенные точечные методы
(28)	1	0.9325	1.4882	16.01	25.6860
(29)–(30)	2	1.0022	1.5995	16.60	26.6294
(35)	1	1.0237	1.6338	16.78	26.9134
(36)–(30)	2	1.1035	1.7612	17.42	27.9431

 

Тепловые свойства почв (Ƙ, λ, d и e) различались в каждой используемой модели. Видно, что наиболее адекватной моделью, отражающей реальность по глубине затухания, являются точечная (17.42 см). Таким образом, введение второй гармоники позволяет уточнить результат и сказать о том, что затухание температурных колебаний происходит примерно на глубине 3 см по сравнению с классическими методами.

Значения температуры почвы T(ɀ, t) на глубине ɀ = 5, 10, 15, 20 и 40 см были рассчитаны по формулам (5), (6) и (22), (23). Для значений Ƙ использовали из табл. 4.

Измеренные T_изм(ɀ, t) и прогнозируемые T_прог (ɀ, t) значения температуры сравнивали для оценки эффективности методов. Результаты приведены в табл. 5.

 

Таблица 5. Эффективность моделей (М1–М8) для прогнозирования температуры почвы на пяти глубинах ɀ = 5, 10, 15, 20 и 40 см

Параметр	M1	М2	M3	M4	M5	M6	M7	M8
	ɀ = 5 см
r, %	95.62	96.08	95.69	95.84	96.09	98.67	96.28	99.44
σT/t	1.10	1.02	1.09	1.06	1.02	0.62	0.98	0.33
	ɀ = 10 см
r, %	80.22	75.91	79.65	78.36	75.83	67.04	73.84	80.97
σT/t	1.46	1.63	1.48	1.53	1.63	1.97	1.66	0.51
	ɀ = 15 см
r, %	97.31	97.97	97.47	97.75	97.97	91.69	97.94	99.68
σT/t	0.93	0.91	0.92	0.91	0.91	1.17	0.89	1.77
	ɀ = 20 см
r, %	97.87	87.94	86.07	86.95	98.43	83.07	98.49	99.64
σT/t	1.29	1.28	1.27	1.27	1.29	1.58	1.26	0.84
	ɀ = 40 см
r, %	87.86	70.88	85.96	81.20	70.54	26.34	62.89	56.31
σT/t	3.05	3.05	3.05	3.05	3.05	3.19	3.05	1.26

Примечание. M1 – амплитуда, M2 – арктангенс, M3 – логарифм, M4 – сдвиг фаз, M5–M8 – улучшенные методы

 

Результаты расчетов показали, что предложенный точечный метод (М8) (формула (36)) является наиболее эффективной моделью, поскольку она дает более точные прогнозы для T(ɀ, t), чем другие алгоритмы.

В соответствии с критериями выбора моделей установлено, что адекватной моделью является предложенная точечная модель (36), поэтому этот модель была учтена при вычисления теплового потока с поверхности почвы.

 

Рис. 2. Сезонный ход ежечасных температур T(0, t_i) и плотность теплового потока q(0, t_i) на поверхности почвы, ɀ = 0, вычисленные различными методами (М1–М8).

 

Используя данные табл. 2–4, а также формулы (24)–(30) и (36) определили тепловой поток q, на поверхности (ɀ = 0) почвы в момент времени t, рис. 2. Согласно точечной модели (М8), наибольший тепловой поток был определен в 12:00 дня (q = 106.85 Вт м⁻²), а наименьший – в 3:00 ночи (q = –64.62 Вт м⁻²) на поверхности почвы.

ЗАКЛЮЧЕНИЕ

На основе исследований модели переноса тепла в почве при учете динамики граничных условий на поверхности, описываемых двумя гармониками, получены аналитические решения при условии второго граничного условии конечной глубине. На основе этого решения предложен точечный метод для определения коэффициента температуропроводности почвы. При использовании этого метода используется распределение температуры в почвенном слое [0, L] для восьми моментов времени в расчетном интервале времени τ₀ для каждой глубины.

При моделировании теплообмена в почвах необходимо учитывать условия второго рода на нижней границе. Так как на определенных глубинах почвы конвективный теплообмен отсутствует и тепловой поток равен 0, поэтому необходимо принять граничные условия второго рода. По критериям выбора моделей установлено, что более адекватной моделью является предложенная точечная модель (М8). Глубина затухания температурных волн, вычисленная по этой методике (d = 17.45 см), соответствуют реальному распределению волн по профилю почв. Следует отметит, что точечные методы, в отличие от классических, учитывает важный параметр – амплитуду колебаний температуры поверхности почвы.

КОНФЛИКТ ИНТЕРЕСОВ

Авторы заявляют об отсутствии у них конфликта интересов.

About the authors

R. Mikail

Igdir University

Email: fariz.mikailsoy@igdir.edu.tr

Department of Mathematics

Turkey, Igdır, 76000

E. Hazar

Igdir University

Email: fariz.mikailsoy@igdir.edu.tr

Department of Mathematics

Turkey, Igdır, 76000

E. Shein

LomonosovMoscow State University

Email: fariz.mikailsoy@igdir.edu.tr
Russian Federation, Moscow, 119991

F. Mikailsoy

Igdir University

Author for correspondence.
Email: fariz.mikailsoy@igdir.edu.tr

Department of Mathematics

Turkey, Igdır, 76000

References

Supplementary files