Прикладная математика и механика, 2022, T. 86, № 1, стр. 105-120

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

В. К. Гончаров *

Санкт-Петербургский государственный морской технический университет
Санкт-Петербург, Россия

* E-mail: vkgonch@mail.ru

Поступила в редакцию 13.03.2020
После доработки 02.11.2020
Принята к публикации 26.11.2020

Полный текст (PDF)

Аннотация

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

Ключевые слова: атмосфера, водная масса, ледовый покров, температура, прочность, поток тепла, кристаллизация

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

Температура внутри пресноводного льда определяется процессами, происходящими на границах и в его толще, в том числе: вариацией температуры атмосферы, поглощением солнечной радиации в толще льда, потоком тепла от нижележащей водной среды и донного грунта, а также от происходящего на границе с водой намерзания льда [13].

Параметризация указанных процессов приводит к формулированию системы уравнений теплопроводности параболического типа для льда, водной массы и донного грунта, связанных условиями равенства потоков тепла и температуры на границах раздела [3]. Следует также учитывать осложняющий решение фактор: толщина льда является функцией времени, так как меняется вследствие намерзания. Это обстоятельство практически исключает возможность использовать для решения рассматриваемой задачи известные инженерные программы расчета тепловых процессов, например, ANSYS.

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

Следует отметить, что подобные задачи: решение систем уравнений, описывающих потоки тепла в соприкасающихся средах, ставились и решались ранее применительно к ледовым режимам озер и рек [5, 6]. Процессы переноса тепла и динамика поля температуры в толще материкового ледника рассмотрены [7] в связи с процессом его развития (нарастания толщины из выпадающего снега) и движения по каменному склону.

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

Выполненные ранее исследования [8, 9] позволили связать скорость намерзания пресноводного льда и поток тепла его кристаллизации с температурой атмосферного воздуха и толщиной льда. В качестве основы была использована аналогичная формуле Зубова [6] зависимость, аппроксимирующая материалы наблюдений за толщиной льда на реках Сибири и связывающую толщину льда с накопленной суммой среднесуточных отрицательных температур воздуха, а также учитывающую толщину снежного покрова [1, 10].

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

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

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

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

В общем виде система уравнений, параметризующая указанные выше процессы и характеризующих поле температуры в тех средах: лед–вода–грунт, следуя [3], в системе координат: 0 ≤ z < ∞ (начало на поверхности льда), 0 ≤ t < ∞, (при условии постоянства всех параметров в горизонтальной плоскости) имеет следующий вид

(2.1)
$\begin{gathered} {{c}_{1}}{{\gamma }_{1}}\frac{{\partial {{T}_{1}}}}{{\partial t}} = {{\lambda }_{1}}\frac{{{{\partial }^{2}}{{T}_{1}}}}{{\partial {{z}^{2}}}} + {{q}_{0}}\frac{\alpha }{{2\sqrt z }}{{e}^{{ - \alpha \sqrt z }}};\quad 0 < z < {{h}_{1}} \\ {{c}_{2}}{{\gamma }_{2}}\frac{{\partial {{T}_{2}}}}{{\partial t}} = \frac{\partial }{{\partial z}}\left( {{{\lambda }_{2}}\frac{{\partial {{T}_{2}}}}{{\partial z}}} \right) + {{q}_{0}}\frac{{{{\alpha }_{1}}}}{{2\sqrt z }}{{e}^{{ - {{\alpha }_{1}}\sqrt z }}}\beta {{e}^{{ - \beta \left( {z - {{h}_{1}}} \right)}}};\quad {{h}_{1}} < z < {{h}_{2}} \\ {{c}_{3}}{{\gamma }_{3}}\frac{{\partial {{T}_{3}}}}{{\partial t}} = {{\lambda }_{3}}\frac{{{{\partial }^{2}}{{T}_{3}}}}{{\partial {{z}^{2}}}};\quad {{h}_{2}} < z < \infty \\ \end{gathered} $

В этих формулах использованы следующие обозначения: T – температура, сi – удельная теплоемкость, γi – плотность, λi – коэффициент теплопроводности, q0 – интенсивность потока солнечной радиации, проникающей в лед, α – показатель ослабления солнечной радиации во льду, β – показатель ослабления солнечной радиации в воде, hi – толщина слоя среды. Индексы 1, 2 и 3 относятся ко льду, воде и грунту, соответственно. Коэффициент теплопроводности в водной среде λ2 является функцией глубины, так как плотность воды, связанная с температурой, и скорость течения, определяющая интенсивность турбулентности, меняются с глубиной.

Граничные условия для этой системы уравнений определяются характером энергетического взаимодействия граничащих сред. На поверхности льда – на границе с атмосферой условие имеет следующий вид

(2.2)
$z = 0;\quad - {\kern 1pt} {{\lambda }_{1}}\frac{{\partial {{T}_{1}}}}{{\partial z}} = {{\mu }_{a}}\left( {{{T}_{a}} - {{T}_{1}}} \right)$

В этой формуле Та – температура атмосферного воздуха, μа – коэффициент турбулентного теплообмена с воздухом на поверхности льда. Граничное условие на границе льда с водой – условие Стефана

(2.3)
$z = {{h}_{1}}{\text{:}}\quad {{T}_{1}} = {{T}_{2}} = 0,\quad - {\kern 1pt} {{\lambda }_{i}}\frac{{\partial {{T}_{i}}}}{{\partial z}} = {{\gamma }_{1}}{{L}_{m}}\frac{{\partial {{h}_{1}}}}{{\partial t}}$

В этой формуле Lm – теплота фазового перехода вода–лед. Граничное условие на границе: вода–грунт и для грунта – геотермический поток QE

(2.4)
$z = {{h}_{2}}{\kern 1pt} :\quad {{T}_{2}} = {{T}_{3}},\quad {{\lambda }_{2}}\frac{{\partial {{T}_{2}}}}{{\partial z}} = {{\lambda }_{3}}\frac{{\partial {{T}_{3}}}}{{\partial z}};\quad z \to \infty {\kern 1pt} :\quad {{\lambda }_{3}}\frac{{\partial {{T}_{3}}}}{{\partial z}} = {{Q}_{E}}$

В качестве начальных условий следует определить начальную температуру воздуха и профили температуры во всех трех средах, а также начальную толщину льда, то есть

(2.5)
$t = 0{\text{:}}\quad {{T}_{a}} = {{T}_{{a0}}},\quad {{T}_{i}}\left( {z,0} \right) = {{T}_{{i0}}}{{f}_{i}}\left( z \right),\quad {{h}_{{10}}} = {{h}_{0}}$

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

(2.6)
$t > 0{\text{:}}\quad {{T}_{a}}\left( t \right) = {{T}_{{a0}}}{{F}_{T}}\left( t \right),\quad {{q}_{0}}\left( t \right) = {{q}_{s}}{{F}_{s}}\left( t \right)$

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

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

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

Результаты исследований соотношения температуры атмосферного воздуха Tair и температуры поверхности льда Tice на реках, приведенные в [1], показали, что эти величины связаны следующим линейным соотношением

(3.1)
${{T}_{{\operatorname{ice} }}} = \beta {{T}_{{\operatorname{air} }}}$

Для больших толщин льда по многолетним наблюдениям на реках Сибири среднее значение коэффициента пропорциональности оказалось равным β = 0.88. Эти результаты позволяют значительно упростить выражение для граничного условия на верхней поверхности льда, пренебрегая эффектами, связанными со снежным покровом и термическим пограничным слоем. Допустимо задать изменение во времени температуры поверхности льда, считая, что она “следит” за температурой атмосферного воздуха, то есть температура поверхности льда меняется пропорционально температуре воздуха.

На границе раздела льда с водной средой происходит кристаллизация и намерзание льда. Интенсивность этого процесса по результатам наблюдений на реках [1, 10], и в морских акваториях (формула Зубова) [4, 11] зависит главным образом от температуры атмосферного воздуха, и может быть аппроксимирована зависимостью от “накопленной” суммы отрицательных среднесуточных температур воздуха. Это наблюдение позволяет параметризовать поток тепла кристаллизации льда функцией температуры атмосферного воздуха и толщины льда в рассматриваемый момент, или, соответственно, функцией времени.

Учитывая изложенное, система уравнений (2.1) переходит в одно уравнение, граничные условия (2.2) и (2.3) значительно упрощаются, и рассматриваемая проблема сводится к следующей известной [12, 13] краевой задаче математической физики:

(3.2)
$\begin{gathered} c\gamma \frac{{\partial T}}{{\partial t}} = \lambda \frac{{{{\partial }^{2}}T}}{{\partial {{z}^{2}}}} + {{q}_{s}}\frac{{\alpha {{f}_{s}}(t)}}{{2\sqrt z }}\exp \left( { - \alpha \sqrt z } \right) \\ T\left( {0,t} \right) = {{T}_{{a0}}}{{f}_{a}}\left( t \right),\quad T\left( {H,t} \right) = 0 \\ \lambda \frac{{\partial T\left( {H,t} \right)}}{{\partial z}} = {{Q}_{{\text{f}}}}{{f}_{q}}\left( t \right) \\ T\left( {z,0} \right) = {{T}_{{a0}}}{{f}_{z}}\left( z \right),\quad H(t) = {{H}_{0}}{{f}_{h}}(t),\quad 0 < z < H,\quad 0 < t < \infty \\ \end{gathered} $

В этих формулах использованы следующие обозначения: T – температура льда, с – удельная теплоемкость льда, γ – плотность, λ – коэффициент теплопроводности льда, qs – интенсивность потока солнечной радиации, проникающей в лед, Qf  – поток тепла кристаллизации льда, Tа0 – начальная температура поверхности льда, Н – толщина льда, α – показатель ослабления солнечной радиации во льду,  fs(t), fa(t) и fq(t) – функции, характеризующие изменение во времени интенсивности солнечной радиации, температуры поверхности льда и потока тепла кристаллизации, соответственно,  fz(z) – начальное распределение температуры (профиль) по толщине льда, функция fh(t) – характеризует изменение толщины льда вследствие его намерзания.

Используя обозначения: KT  = λ/с⋅γ, S  = q0/сγ и Rf  = Qf/λ, преобразуем выражения (3.2) к обычному для уравнения теплопроводности виду

(3.3)
$\begin{gathered} \frac{{\partial T}}{{\partial t}} = {{K}_{T}}\frac{{{{\partial }^{2}}T}}{{\partial {{z}^{2}}}} + S\frac{{\alpha {{f}_{s}}(t)}}{{2\sqrt z }}\exp \left( { - \alpha \sqrt z } \right){\text{;}}\quad 0 < z < H,\quad 0 < t < \infty \\ T\left( {0,t} \right) = {{T}_{{a0}}}{{f}_{a}}\left( t \right),\quad T\left( {H,t} \right) = 0{\text{;}}\quad 0 < t < \infty \\ \frac{{\partial T\left( {H,t} \right)}}{{\partial z}} = {{R}_{{\text{f}}}}{{f}_{q}}\left( t \right){\text{;}}\quad 0 < t < \infty \\ T\left( {z,0} \right) = {{T}_{{a0}}}{{f}_{z}}\left( z \right),\quad H(t) = {{H}_{0}}{{f}_{h}}(t) \\ \end{gathered} $

Перейдем к безразмерным переменным. Для этого вертикальную координату нормируем по толщине льда: ζ = z/H. Время нормируем на некоторый временной масштаб D, равный, например, 1 суткам (τ = t/D). Тогда выражения (3.3) примут следующий вид

(3.4)
$\begin{gathered} \frac{{\partial T}}{{\partial \tau }} = {{R}_{T}}\frac{{{{\partial }^{2}}T}}{{\partial {{\zeta }^{2}}}} + {{Q}_{s}}\frac{{{{\alpha }_{h}}{{f}_{{st}}}(\tau )}}{{2\sqrt \zeta }}\exp \left( { - {{\alpha }_{h}}\sqrt \zeta } \right) \\ T\left( {0,\tau } \right) = {{T}_{{a0}}}{{f}_{{at}}}\left( \tau \right),\quad T\left( {1,\tau } \right) = 0 \\ \frac{{\partial T\left( {1,\tau } \right)}}{{\partial \zeta }} = {{R}_{{{\text{fh}}}}}{{f}_{{qt}}}\left( \tau \right) \\ T\left( {\zeta ,0} \right) = {{T}_{{a0}}}{{f}_{{zh}}}\left( \zeta \right) \\ {{R}_{T}} = \frac{{{{K}_{T}}D}}{{{{H}^{2}}}},\quad {{Q}_{s}} = \frac{{{{q}_{s}}D}}{{c\lambda }},\quad {{R}_{{{\text{fh}}}}} = {{R}_{{\text{f}}}}H,\quad {{\alpha }_{h}} = \frac{\alpha }{{\sqrt H }},\quad 0 < \zeta < 1,\quad 0 < \tau < \infty \\ \end{gathered} $

В этих выражениях функции fst(τ), fat(τ) и fqt(τ), а также fzh(ζ) записаны в безразмерном виде.

Так как толщина льда является функцией времени, то, строго говоря, зависящими от времени являются величины RT, Rfr и αh. Намерзание льда – достаточно медленный процесс: по данным полевых измерений [1] при обычном зимнем температурном режиме толщина льда увеличивается в среднем на 1% в сутки. Поэтому во втором слагаемом в правой части уравнения (3.4) этим эффектом можно пренебречь в сравнении с суточной вариацией солнечной радиации. Зависимость величины Rfr от времени будет учтена при параметризации скорости намерзания льда и соответствующего потока тепла кристаллизации.

Зависимость коэффициента температуропроводности RT от времени представляется возможным устранить, если нормирование этого коэффициента и времени процесса представить в следующем виде

(3.5)
${{\bar {R}}_{T}} = \frac{{{{K}_{T}}D}}{{H_{0}^{2}}},\quad \bar {\tau } = \int {\frac{{dt}}{{{{{\left[ {{{f}_{h}}\left( t \right)} \right]}}^{2}}}}} $

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

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

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

4. Эффект поглощения льдом солнечной радиации. Для исследования влияния поглощения льдом солнечной радиации на его температуру, уравнение и граничные условия (3.4) представим следующим образом

(4.1)
$\begin{gathered} \frac{{\partial {{T}_{{sn}}}}}{{\partial \tau }} = {{R}_{T}}\frac{{{{\partial }^{2}}{{T}_{{sn}}}}}{{\partial {{\zeta }^{2}}}} + {{Q}_{s}}\frac{{{{\alpha }_{h}}{{f}_{{st}}}(\tau )}}{{2\sqrt \zeta }}\exp \left( { - {{\alpha }_{h}}\sqrt \zeta } \right);\quad 0 < \zeta < 1,\quad 0 < \tau < \infty \\ {{T}_{{sn}}}\left( {0,\tau } \right) = 0,\quad {{T}_{{sn}}}\left( {1,\tau } \right) = 0;\quad 0 < \tau < \infty \\ \end{gathered} $

Решение этой задачи известно [12], оно имеет следующий вид

(4.2)
${{T}_{{sn}}}(\zeta ,\tau ) = {{\alpha }_{h}}{{Q}_{{hd}}}\int\limits_0^\tau {{{f}_{{st}}}\left( \vartheta \right)d\vartheta } \int\limits_0^1 {\frac{{{{e}^{{ - {{\alpha }_{h}}\sqrt \eta }}}}}{{\sqrt \eta }}} \sum\limits_{n = 1}^{ + \infty } {{{e}^{{ - {{n}^{2}}{{\pi }^{2}}{{R}_{T}}(\tau - \vartheta )}}}\sin \left( {n\pi \eta } \right)} \sin \left( {n\pi \zeta } \right)d\eta $

Для использования решения следует определить величину Qs и вид функции fst(t), которая характеризует зависимость от времени интенсивности солнечной радиации, проникающей в лед.

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

(4.3)
${{f}_{{st}}}\left( \tau \right) = \left\{ \begin{gathered} {{s}_{0}} - {{k}_{s}}\sin \left( {2\pi \tau + {{\sigma }_{s}}} \right)\quad при\quad \sin \left( {2\pi \tau + {{\sigma }_{s}}} \right) < 0 \hfill \\ 0\quad при\quad \sin \left( {2\pi \tau + {{\sigma }_{s}}} \right) < 0 \hfill \\ \end{gathered} \right.$

В этой формуле σs – сдвиг по фазе, который позволяет совместить должным образом начало отсчета времени с появлением солнечной радиации (с восходом солнца). Вариацией величины s0 можно регулировать соотношение длительности светлого и темного времени суток, а величина ks позволяет нормировать до “1” максимальное значение рассматриваемой функции  fst(τ).

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

(5.1)
$\begin{gathered} \frac{{\partial T}}{{\partial \tau }} = {{R}_{T}}\frac{{{{\partial }^{2}}T}}{{\partial {{\zeta }^{2}}}};\quad 0 < \zeta < 1,\quad 0 < \tau < \infty \\ T\left( {0,\tau } \right) = {{T}_{{a0}}}{{f}_{{at}}}\left( \tau \right),\quad T\left( {1,\tau } \right) = 0;\quad 0 < \tau < \infty \\ T\left( {\zeta ,0} \right) = {{T}_{0}}{{f}_{{hh}}}\left( \zeta \right);\quad 0 < \zeta < 1 \\ \end{gathered} $

Решение краевой задачи такого типа известно [12], и в рассматриваемом случае решение имеет следующий вид

(5.2)
$\begin{gathered} {{T}_{{am}}}\left( {\zeta ,\tau } \right) = {{T}_{{a0}}}\left\{ {{{f}_{{hh}}}\left( \zeta \right) + \left. {{{f}_{{at}}}\left( \tau \right)} \right|_{{\zeta = 0}}^{{^{{^{{^{{^{{}}}}}}}}}} + } \right. \\ \left. { + \;\frac{1}{{2\sqrt {\pi {{R}_{T}}} }}\int\limits_0^\tau {\frac{{{{f}_{{at}}}\left( \vartheta \right)}}{{{{{\left( {\tau - \vartheta } \right)}}^{{3/2}}}}}\sum\limits_{ - \infty }^{ + \infty } {\left( {\zeta + 2n} \right)\exp \left( { - \tfrac{{{{{\left( {\zeta + 2n} \right)}}^{2}}}}{{4{{R}_{T}}\left( {\tau - \vartheta } \right)}}} \right)d\vartheta } } } \right\} \\ \end{gathered} $

Начальное распределение температуры по толщине льда для целей упрощения решения целесообразно принять в виде линейной зависимости, предполагая, что в предшествующий моделированию период времени температурный режим не менялся, и в результате внутри льда установился некоторый стационарный термический режим, при котором температура поверхности льда равна температуре воздуха βТа0 в соответствии с (3.1), а температура на границе с водной средой постоянна и равна температуре замерзания воды, то есть Т = 0°С. Соответствующая зависимость имеет такой вид

(5.3)
${{f}_{{hh}}}\left( \zeta \right) = 1 - \zeta $

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

(5.4)
${{f}_{{at}}}(\tau ) = 1 + {{k}_{T}}\tau + \delta T\sin \left( {2\pi \tau + {{\sigma }_{T}}} \right)$

В этой формуле kT  характеризует понижение температуры воздуха относительно начальной Tа0 каждые сутки, δT – вариация температуры на протяжении суток (некоторая доля начальной температуры Tа0), постоянная на протяжении всего срока, σТ  – сдвиг по фазе для совмещения начала и характера вариации температуры (потепление/похолодание) с солнечной радиацией (2.3).

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

(6.1)
$\begin{gathered} \frac{{\partial {{T}_{{{\text{cr}}}}}}}{{\partial \tau }} = {{R}_{T}}\frac{{{{\partial }^{2}}{{T}_{{{\text{cr}}}}}}}{{\partial {{\zeta }^{2}}}};\quad 0 < \zeta < 1,\quad 0 < \tau < \infty \\ {{T}_{{{\text{cr}}}}}\left( {1,\tau } \right) = 0,\quad \frac{{\partial {{T}_{{{\text{cr}}}}}\left( {1,\tau } \right)}}{{\partial \zeta }} = {{R}_{{{\text{fh}}}}}{{f}_{{qt}}}\left( \tau \right);\quad 0 < \tau < \infty \\ \end{gathered} $

Условия на границе с водой (ζ = 1) соответствуют условию Стефана: постоянство температуры, равной температуре замерзания воды, и градиент температуры, определяющийся потоком тепла кристаллизации. Поток тепла кристаллизации является функцией времени, так как интенсивность намерзания льда, очевидно, связана с вариацией температуры на границе льда с атмосферой (ζ = 0).

В канонической форме краевая задача (6.1) имеет следующий вид

(6.2)
$\begin{gathered} {{u}_{t}} = {{a}^{2}}{{u}_{{\zeta \zeta }}};\quad 0 < \zeta < 1,\quad 0 < t < \infty \\ {{u}_{z}}(1,t) = - \varphi (t),\quad u(1,t) = 0;\quad 0 < t < \infty \\ \end{gathered} $

Задача относится к классу параболических уравнений на отрезке конечной длины l = 1, и ее решение можно получить с использованием функции Грина [12, 13]. Решение в данном случае имеет такую форму

(6.3)
$\begin{gathered} u(\zeta ,t) = - {{a}^{2}}\int\limits_0^t {\varphi (t)} G\left( {\zeta ,t - \tau } \right)d\tau \\ {{a}^{2}} = {{R}_{T}},\quad \varphi \left( t \right) = - {{R}_{{{\text{fh}}}}}{{f}_{{qt}}}\left( t \right) \\ \end{gathered} $

В рассматриваемой задаче (отрезок конечной длины) функцию Грина следует представить в следующем виде

(6.4)
$G\left( {\zeta ,\tau - \vartheta } \right) = \frac{1}{H} + \frac{2}{H}\sum\limits_{n = 1}^{ + \infty } {{{e}^{{ - {{n}^{2}}{{\pi }^{2}}{{R}_{T}}\left( {\tau - \vartheta } \right)}}}\cos \left[ {n\pi \left( {1 - \zeta } \right)} \right]} $

Выполняя необходимые преобразования, получаем следующее решение, которое удовлетворяет граничному условию при ζ = 1 (z = H), представляющему собой поток тепла, порождаемый кристаллизацией льда на границе с водной средой

(6.5)
${{\tilde {T}}_{{{\text{cr}}}}}(\zeta ,\tau ) = - {{R}_{T}}{{R}_{{{\text{fh}}}}}\int\limits_0^t {{{f}_{{qt}}}(\vartheta )} \left[ {1 + 2\sum\limits_{n = 1}^{ + \infty } {{{e}^{{ - {{n}^{2}}{{\pi }^{2}}{{R}_{T}}(\tau - \vartheta )}}}} \cos \left[ {n\pi \left( {1 - \zeta } \right)} \right]} \right]d\vartheta $

Для того, чтобы обеспечить равенство T = 0 на границе раздела с водой, следует вычесть из него частное решение уравнения (6.1) при ζ = 1

(6.6)
$\begin{gathered} {{T}_{{{\text{cr}}}}}(\zeta ,\tau ) = - {{R}_{T}}{{R}_{{{\text{fh}}}}}\left\{ {\int\limits_0^t {{{f}_{{qt}}}(\vartheta )} \left[ {1 + 2\sum\limits_{n = 1}^{ + \infty } {{{e}^{{ - {{n}^{2}}{{\pi }^{2}}{{R}_{T}}(\tau - \vartheta )}}}} \cos \left[ {n\pi \left( {1 - \zeta } \right)} \right]} \right]d\vartheta } \right. - \\ \left. {{{{\left. { - \;\int\limits_0^t {{{f}_{{qt}}}(\vartheta )} \left[ {1 + 2\sum\limits_{n = 1}^{ + \infty } {{{e}^{{ - {{n}^{2}}{{\pi }^{2}}{{R}_{T}}(\tau - \vartheta )}}}} \cos \left[ {n\pi \left( {1 - \zeta } \right)} \right]} \right]d\vartheta } \right|}}_{{_{{\zeta = 1}}}}}} \right\} \\ \end{gathered} $

В решении (6.6) остается определить вид функции φ(t) = Rfh fqt(t), которая характеризует поток тепла от намораживания льда. С этой целью воспользуемся эмпирической формулой [1, 10], основанной на материалах контроля толщины льда и температуры воздуха на сибирских реках. Для толщин снега менее 20 см имеем

(6.7)
${{h}_{{\operatorname{ice} }}} = 2.4{{\left( {\sum {\left| { - {{T}_{a}}} \right|} } \right)}^{{1/2}}}$

В этой формуле Ta – среднесуточная отрицательная температура атмосферного воздуха, $\sum {\left| {--{{Т}_{а}}} \right|} $ – абсолютная величина суммы среднесуточных температур воздуха к моменту регистрации толщины льда hice. Нужно отметить, что эта эмпирическая зависимость практически учитывает, и усредняет эффекты всех возможных факторов, влияющих на образование и рост толщины пресноводного речного льда.

Формула (6.7) позволяет численно смоделировать скорость увеличения толщины льда при различных режимах вариации температуры воздуха и начальных величинах толщины льда, а также для задаваемых приращений толщины льда. В работе [9] на основе моделирования этих параметров получены следующие формулы для скорости намерзания и толщины льда

(6.8)
${{W}_{{{\text{fr}}}}}({{h}_{{\operatorname{ice} 0}}},{{T}_{a}}) = 1.66 \times {{10}^{{ - 4}}}(2.73h_{{\operatorname{ice} 0}}^{{ - 0.68}} - 1)\left| {{{T}_{a}}} \right|$
(6.9)
${{h}_{{\operatorname{ice} }}}({{h}_{{\operatorname{ice} 0}}},{{T}_{a}}) = 1.66 \times {{10}^{{ - 4}}}(2.73h_{{\operatorname{ice} 0}}^{{ - 0.68}} - 1)N\sum\limits_1^N {\left| {{{T}_{a}}} \right|} $

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

(6.10)
${{Q}_{{{\text{fr}}}}} = {{W}_{{{\text{fr}}}}}{{L}_{{\operatorname{ice} }}}{{\gamma }_{{\operatorname{ice} }}}$

В этой формуле Lice – удельная теплота кристаллизации пресноводного льда, γice – плотность пресноводного льда. Поток тепла кристаллизации льда распределяется между льдом и водной средой пропорционально теплопроводности этих сред [14]. В толщу льда тепло кристаллизации распространяется молекулярной теплопроводностью, а в водную среду за счет турбулентного перемешивания. Долю тепла кристаллизации, поступающую в лед можно определить коэффициентом

(6.11)
${{\delta }_{{i/w}}} = \frac{{{{\lambda }_{{\operatorname{ice} }}}}}{{{{\lambda }_{{\operatorname{ice} }}} + {{\lambda }_{w}}}}$

В этой формуле λice – теплопроводность льда. Теплопроводность водной среды можно оценить на основании исследований динамики речных течений. Для равнинной реки Амур целесообразно принять среднее по данным [15] значение λw = 60 Вт/м °С.

Для использования аппроксимации (6.10) и (6.11) в качестве граничного условия в задаче (6.1) этот поток следует отнести к теплопроводности пресноводного льда λice. Кроме того, следует учесть, что температура воздуха является функцией времени. После подстановки значений констант получаем выражение

(6.12)
${{R}_{{{\text{fr}}}}}\left( {{{h}_{0}},{{T}_{a}}} \right) = {{\delta }_{{i/w}}}\frac{{{{Q}_{{{\text{fr}}}}}}}{{{{\lambda }_{{\operatorname{ice} }}}}} = {{\delta }_{{i/w}}}{{W}_{{{\text{fr}}}}}\frac{{{{L}_{{\operatorname{ice} }}}{{\gamma }_{{\operatorname{ice} }}}}}{{{{\lambda }_{{\operatorname{ice} }}}}} = 0.276{{\delta }_{{i/w}}}(2.73h_{0}^{{ - 0.68}} - 1)\left| {{{T}_{a}}} \right|$

7. Модель вариации температуры льда. Общее решение краевой задачи (3.4), которое определяет вариацию температурного поля в толще пресноводного льда под действием солнечного прогрева и вариации температуры атмосферного воздуха записывается в виде линейной суперпозиции решений (4.2), (5.2) и (6.6)

$T(\zeta ,\tau ) = {{T}_{{{\text{sn}}}}}(\zeta ,\tau ) + {{T}_{{{\text{am}}}}}(\zeta ,\tau ) + {{T}_{{{\text{cr}}}}}(\zeta ,\tau )$
$T(\zeta ,\tau ) = {{\alpha }_{h}}{{Q}_{{hd}}}\int\limits_0^\tau {{{f}_{{{\text{st}}}}}\left( \vartheta \right)d\vartheta } \int\limits_0^1 {\frac{{{{e}^{{ - {{\alpha }_{h}}\sqrt \eta }}}}}{{\sqrt \eta }}\left\{ {\sum\limits_{n = 1}^{ + \infty } {{{e}^{{ - {{n}^{2}}{{\pi }^{2}}{{R}_{T}}(\tau - \vartheta )}}}\sin \left( {n\pi \eta } \right)} \sin \left( {n\pi \zeta } \right)} \right\}d\eta } + $
$ + \;{{T}_{0}}\left\{ {{{f}_{{hh}}}\left( \zeta \right) + {{{\left. {{{f}_{{{\text{at}}}}}\left( \tau \right)} \right|}}_{{\zeta = 0}}} + \frac{1}{{2\sqrt {\pi {{R}_{T}}} }}\int\limits_0^\tau {\frac{{{{f}_{{{\text{at}}}}}\left( \vartheta \right)}}{{{{{\left( {\tau - \vartheta } \right)}}^{{\tfrac{3}{2}}}}}}\sum\limits_{ - \infty }^{ + \infty } {\left( {\zeta + 2n} \right)\exp \left( { - \tfrac{{{{{\left( {\zeta + 2n} \right)}}^{2}}}}{{4{{R}_{T}}\left( {\tau - \vartheta } \right)}}} \right)d\vartheta } } } \right\} + $
$ - \;{{R}_{T}}{{R}_{{{\text{fh}}}}}\left\{ {\int\limits_0^t {{{f}_{{qt}}}(\vartheta )} \left[ {1 + 2\sum\limits_{n = 1}^{ + \infty } {{{e}^{{ - {{n}^{2}}{{\pi }^{2}}{{R}_{T}}(\tau - \vartheta )}}}} \cos \left[ {n\pi \left( {1 - \zeta } \right)} \right]} \right]d\vartheta } \right. - $
(7.1)
$\left. {{{{\left. { - \;\int\limits_0^t {{{f}_{{qt}}}(\vartheta )} \left[ {1 + 2\sum\limits_{n = 1}^{ + \infty } {{{e}^{{ - {{n}^{2}}{{\pi }^{2}}{{R}_{T}}(\tau - \vartheta )}}}} \cos \left[ {n\pi \left( {1 - \zeta } \right)} \right]} \right]d\vartheta } \right|}}_{{_{{\zeta = 1}}}}}} \right\}$

Входящие в это выражение функции определены в следующем виде

${{f}_{{{\text{at}}}}}(\vartheta ) = 1 + {{k}_{T}}\vartheta + \delta T\sin \left( {2\pi \vartheta + {{\sigma }_{T}}} \right)$
(7.2)
$\begin{gathered} {{f}_{{hh}}}\left( \eta \right) = 1 - \eta \\ {{f}_{{qt}}}\left( \vartheta \right) = {{\delta }_{{i/w}}}(2.73H_{0}^{{ - 0.68}} - 1)\left| {{{T}_{0}}} \right|{{f}_{{{\text{at}}}}}\left( \vartheta \right) \\ \end{gathered} $
${{f}_{{{\text{st}}}}}\left( \vartheta \right) = \left\{ \begin{gathered} {{s}_{0}} - {{k}_{s}}\sin \left( {2\pi \vartheta + {{\sigma }_{s}}} \right)\quad при\quad \sin \left( {2\pi \vartheta + {{\sigma }_{s}}} \right) < 0 \hfill \\ 0\quad при\quad \sin \left( {2\pi \vartheta + {{\sigma }_{s}}} \right) > 0 \hfill \\ \end{gathered} \right.$

Для оценки толщины льда в моделируемый период времени использовалась основанная на выражении (6.9) формула

(7.3)
${{H}_{j}} = {{H}_{0}} + 0.276\left| {{{T}_{0}}} \right|\sum\limits_{j = 1}^M {(2.73H_{{j0}}^{{ - 0.68}} - 1)\int\limits_{{{t}_{{j - 1}}}}^{{{t}_{j}}} {{{f}_{{{\text{at}}}}}\left( t \right)dt} } $

Здесь М – количество интервалов времени, по прошествии которых оценивается приращение и новая толщина льда вследствие намораживания в этот период. Этот интервал времени, очевидно, должен равняться или превышать интервал численного интегрирования по времени в формулах (7.1). По результатам расчета толщины льда на каждом  j-интервале уточняются значения коэффициентов, которые были нормированы по толщине льда при переходе от уравнения (3.3) к безразмерной его форме (3.4), т.е.

(7.4)
${{R}_{{Tj}}} = \frac{{{{K}_{T}}D}}{{H_{j}^{2}}},\quad {{R}_{{{\text{fh}}j}}} = {{R}_{{\text{f}}}}{{H}_{j}},\quad {{\alpha }_{{hj}}} = \frac{\alpha }{{\sqrt {{{H}_{j}}} }}$

Проверка адекватности решения сопоставлением с материалами полевых измерений.

Для проверки адекватности поученных результатов было выполнено сопоставление расчетов по полученным формулам с материалами полевых измерений толщины льда и вертикального распределения температуры от воды до атмосферы на реке Амур (Хэйлунцзян) в январе 2014 и 2015 г., выполненных соисполнителями Проекта из Тайюаньского технологического университета (КНР). Описаны измерительные средства и методологии измерений [16]. Для моделирования приняты физические характеристики пресноводного льда [3].

На рис. 1 показана вариация температуры воздуха в январе 2014 г. и сопоставлены измеренные толщины льда (точки) и рассчитанные по формуле (7.3) величины (сплошная линия). Согласование измеренных и рассчитанных толщин льда можно признать вполне удовлетворительным как по величине, так и по характеру изменения, который соответствует характеру вариации температуры воздуха (проявляется в изменении кривизны расчетной кривой).

Рис. 1.

Сравнение измеренных толщин льда (точки) и рассчитанных (сплошная кривая) толщин льда по среднесуточным температурам атмосферного воздуха (точки и ломаная линия).

В период 6–10 января 2015 г. выполнена регистрация вертикального распределения температуры от водной среды до атмосферы (каждые 8 ч) и измерения толщины льда. На рис. 2 показаны характерные вертикальные профили температуры: 1 – 16:00 (период потепления) и 2 – 08:00 (похолодание). Стрелками показаны границы слоев: вода–лед, лед–снег и снег–атмосфера, которые определены по резкому изменению градиента температуры [17].

Рис. 2.

Измеренные вертикальные профили температуры воздуха, снега, льда и водной массы, и границы их раздела.

На рис. 3 сопоставлены измеренные в январе 2015 г. (точки) и рассчитанные по формуле (7.3) (сплошная кривая) толщины льда, а также зарегистрированные температуры воздуха. Росту толщины льда по данным измерений достаточно хорошо соответствуют результаты расчета, при этом резким изменениям температуры воздуха соответствует изменение интенсивности намерзания льда.

Рис. 3.

Сравнение измеренных толщин льда (точки) и рассчитанных (сплошная кривая) толщин льда по среднесуточным температурам атмосферного воздуха (точки и ломаная линия).

Для моделирования вертикальных профилей температуры в ледовом покрове реки Амур в январе 2015 г. с использованием модели (7.1) и (7.2) требовалось задать вариацию температуры атмосферного воздуха, интенсивность солнечной радиации, начальную толщину ледового покрова и начальный профиль температуры.

Измеренные температуры воздуха были “приведены” к условиям отсутствия теплоизолирующего снежного покрова в соответствии с формулой (3.1), и результаты пересчета показаны на рис. 4.

Рис. 4.

Измеренная температура воздуха в период регистрации вертикальных профилей температуры в ледовом покрове и температура, “приведенная” к условиям отсутствия снежного покрова.

Интенсивность солнечной радиации определена применительно к географической широте места полевых измерений формулой (4.3). Начальный профиль температуры внутри льда принят линейным, соответствующим начальному перепаду между водной средой и зарегистрированной температурой на поверхности льда. Результаты моделирования сопоставлены с данными полевых измерений на рис. 5 и 6. На рисунках в соответствующем масштабе показан слой снега, разделяющий лед и атмосферу. Вертикальная координата нормирована по толщине ледового покрова. Эти материалы показывают достаточно хорошее соответствие результатов расчетов по разработанной модели (7.1) и (7.2) результатам полевых измерений.

Рис. 5.

Сопоставление моделирования температуры внутри ледового покрова (сплошная линия) с данными полевых измерений (точки).

Рис. 6.

Сопоставление моделирования температуры внутри ледового покрова (сплошная линия) с данными полевых измерений (точки).

Исключение составляют слои, примыкающие к границе раздела лед–снег, где имеют место наиболее высокие градиенты температуры во времени. Это обстоятельство следует связывать с тем, что локальная теплопроводность пресноводного льда и снега зависят от их температуры, а процесс вариации их температуры обладает определенной “инерцией”. Этот эффект в наибольшей степени проявляется в результатах измерений, выполненных после повышения температуры атмосферного воздуха до максимальных значений в конце светового дня (рис. 7).

Рис. 7.

Сопоставление моделирования температуры внутри ледового покрова (сплошная линия) с данными полевых измерений (точки), выполненных после резкого повышения температуры атмосферного воздуха.

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

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

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

Исследования выполнены при поддержке Российского фонда фундаментальных исследований, проект № 15-58-53013-ГФЕН_а.

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

  1. Донченко Р.В. Ледовый режим рек СССР. Л.: Гидрометеоиздат, 1987. 242 с.

  2. Лепихин А.П. К расчету ледяного покрова на пресноводных водных объектах // Геогр. вестн. 2008. № 1. С. 1–13.

  3. Красс М.С., Мерзликин В.Г. Радиационная теплофизика снега и льда. Л.: Гидрометеоиздат, 1990. 262 с.

  4. Доронин Ю.П. Физика океана. СПб.: РГГМУ, 2000. 340 с.

  5. Leppäranta Matti. Freezing of Lakes and the Evolution of their Ice Cover. Berlin; Heidelberg: Springer, 2015. 301 p.

  6. Voyevodin A.F., Grankina T.B. Mathematical Modeling of the Ice-Thermal regime of Water Body // Proceedings of the 19th IAHR International Symposium on Ice, Vancouver, Canada, 2008. V. 1. P. 727–735.

  7. Cuffey K.M., Paterson W.S.B. The Physics of Glaciers. Elsevier, 2010. 716 p.

  8. Гончаров В.К., Клементьева Н.Ю., Jianmin Qin и др. Прогнозирование роста толщины пресноводного ледового покрова // Труды международной конференции по судостроению и океанотехнике NAOE2016, 2016. Сборник статей. С. 219–222.

  9. Goncharov V.K., Klementeva N.Yu., Jianmin Qin et al. Investigation of correlation between the temperature on air-snow and snow-ice interfaces and the atmospheric air temperature // J.Earth Sci. and Eng. 2016. V. 6. № 5. P. 245–253.

  10. Чижов А.Н. Формирование ледяного покрова и пространственное распределение его толщины. Л.: Гидрометеоиздат, 1990. 128 с.

  11. Жуков Л.А. Общая океанология. Л.: Гидрометеоиздат, 1976. 376 с.

  12. Будак Б.М., Самарский А.А., Тихонов А.Н. Сборник задач по математической физике. М.: Наука, 1972. 688 с.

  13. Тихонов А.Н., Самарский А.А. Уравнения математической физики. М.: Наука, 1977. 736 с.

  14. Kreith F., Black W.Z. Basic Heat Transfer. Harper and Row Publishers, 1980. Крейт Ф., Блэк У. Основы теплопередачи. М.: Мир, 1983. 512 с.

  15. Великанов М.А. Динамика русловых потоков. Т. 1. Структура потока. М.: Гостехиздат, 1954. 323 с.

  16. Qin Jianmin, Wang Lijuan, Goncharov V.K. et al. Realization of the Auto-Measurement of ice in riverway based on the difference of the air, ice and water // Proc. 23th IAHR Symp. on Ice, Lahti, Finland, 2010. V. 1. P. 571–580.

  17. Hammonds K., Lieb-Lappen R., Baker I., Wang X. Investigating the thermophysical properties of the ice-snow interface under a controlled temperature and gradient. Part I: Experiments and observations // Cold Regions Sci.&Technol., 2015. V. 120. P. 157–167.

Дополнительные материалы отсутствуют.