On parameterization of dissipative processes in turbulent transport models for description of thermohydrodynamics and biogeochemistry of stratified internal water bodies

Texto integral

ВВЕДЕНИЕ

Внутренние водоемы суши, к которым относятся озера и водохранилища, являются важнейшими элементами многих природных ландшафтов, а характеристики этих водных объектов (термический режим, распределение концентраций биохимических субстанций, потоки тепла и примесей в атмосферу) представляют интерес как в рамках задач гидрологии и экологии [Двуреченская и др., 2012; Онищенко, 2016; Петросян и др., 2019; Хоружая, Минина, 2017; Щепетова, Толстова, 2011], так и с точки зрения их влияния на климат Земли и его изменения [Krinner, 2003; Samuelsson et al., 2010; Thiery et al., 2015; Tranvik et al., 2009; Wang et al., 2022; Zhu et al., 2018]. При этом следует отметить, что формирование вышеназванных характеристик во многом определяется процессами турбулентного переноса (импульса, тепла и растворенных газов). Таким образом, встает вопрос о разработке корректного описания турбулентности в замкнутых водных объектах.

В моделях озер и водохранилищ, как правило, рассматривается RANS (Reynolds–Averaged Navier–Stokes) система уравнений термогидродинамики – система уравнений Навье–Стокса в приближении Буссинеска, осредненная по Рейнольдсу. Примерами подобных моделей являются как трехмерные [Гладских и др., 2021; Abbasi et al., 2016], так и одномерные [Fang and Stefan, 2009; Jöhnk, 2009; Stepanenko et al., 2016], разрабатываемые, в частности, для параметризации водоемов в глобальных и региональных моделях атмосферы. Осредненные уравнения включают дополнительные неизвестные корреляции мелкомасштабных составляющих потока – компоненты тензора напряжений Рейнольдса uʹ_i uʹ_j, где uʹ_i есть турбулентные пульсации i-той компоненты вектора скорости, а угловые скобки означают осреднение, в общем случае, по ансамблю реализаций турбулентного течения. Для определения этих напряжений привлекаются модели турбулентности, которые выражают неизвестные корреляции через известные (точнее – искомые) осредненные значения. Одним из наиболее распространенных классов турбулентных замыканий в рамках RANS подхода, применяемых в океанологии и лимнологии, являются двухпараметрические модели, основанные на гипотезе, связывающей турбулентные потоки с градиентами средних величин [Монин, Яглом, 1965]. Такие модели включают уравнения для двух прогностических переменных: первая – кинетическая энергия турбулентности (КЭТ) k, а вторая – скорость диссипации КЭТ ε [Лыкосов, 1992; Мортиков и др., 2019; Burchard, 2002; Mellor, Yamada, 1982], или иная размерная переменная, например, произведение k на турбулентный масштаб длины L, или частота турбулентных пульсаций ω [Umlauf et al., 2003], необходимая для задания коэффициентов турбулентной вязкости K_m и диффузии K_h из соотношений подобия. Приведенные замыкания, несмотря на очевидные упрощения, являются вычислительно эффективными, по сравнению с более сложными замыканиями второго порядка, и позволяют верно описать некоторые важные особенности динамики пограничного слоя, например, заглубление перемешанного слоя в идеализированных постановках [Burchard, 2002].

Отметим, что в так называемой “стандартной” k-ε модели (двухпараметрическая схема, где k и ε являются прогностическими величинами) полагается, что безразмерные функции устойчивости, турбулентные числа Прандтля Pr_T = K_m / K_h и Шмидта Sc_T – отношение коэффициентов турбулентного переноса для импульса и скалярных субстанций, – являются постоянными и не зависящими от стратификации среды. Полуэмпирические замыкания турбулентности второго порядка (см., например, [Canuto et al., 2001; Kantha, Clayson, 1994; Mellor, Yamada, 1982]) могут быть использованы для уточнения зависимости Pr_T и Sc_T от стратификации. Однако вопрос о влиянии таких параметризаций на описание процессов, в частности, биохимических, протекающих во внутренних водоемах, по-прежнему остается открытым. Например, во многих моделях турбулентное число Прандтля используется лишь как калибровочный параметр для воспроизведения термического режима конкретного водоема [Sun et al., 2020].

Ранее при участии авторов настоящей работы было предложено турбулентное замыкание [Соустова и др., 2020; Gladskikh et al., 2023], полученное на основе модели нестационарных турбулентных течений в стратифицированной жидкости [Островский, Троицкая, 1987]. В этой модели с помощью методов кинетической теории газов [Власов, 1966; Lundgren, 1967] получены выражения для турбулентных потоков импульса, кинетической и потенциальной энергий турбулентности в стратифицированной среде. Показано, в частности, что учет двухсторонней трансформации кинетической и потенциальной энергии турбулентных пульсаций позволяет объяснить поддержание турбулентности небольшими сдвигами скорости при любых значениях градиентного числа Ричардсона $Ri= \frac{N^2}{S^2}>0$ (где N – частота плавучести, S – сдвиговая частота среднего течения), характеризующего соотношение между сдвиговой неустойчивостью и действием плавучести. Процедура получения соответствующих уравнений для средних аналогична той, которая используется в кинетической теории газов – замыкается и решается уравнение для одноточечной функции распределения – так называемое кинетическое уравнение, и по известной функции распределения вычисляются напряжения Рейнольдса. При этом учитываются некоторые важные, но обычно пренебрегаемые эффекты, например, зависимость вертикальной анизотропии турбулентности от стратификации, а также неградиентная поправка к турбулентному потоку массы. С привлечением модели [Островский, Троицкая, 1987] получено выражение, связывающее Pr_T c градиентным числом Ричардсона Ri, и тем самым было снято ограничение в “стандартном” k-ε замыкании, где турбулентное число Прандтля полагается постоянным.

В работах [Соустова и др., 2020; Gladskikh et al., 2023] показано, что турбулентное замыкание учитывает вклад сдвига скорости и плотностной стратификации в мелкомасштабную турбулентность внутреннего водоема и, в том числе, обеспечивает отсутствие явления “вырождения” турбулентного обмена (ламинаризации течения) при значениях градиентного числа Ричардсона, превышающих некоторое эмпирическое пороговое значение. Учет данных эффектов влияет на вертикальные распределения кинетической энергии турбулентности, температуры и пассивных примесей.

Основные выводы модели [Соустова и др., 2020; Gladskikh et al., 2023] согласуются с результатами теории, развиваемой в работах группы С. С. Зилитинкевича [Zilitinkevich et al., 2007, 2013], в рамках которой на основе анализа уравнений для одноточечных моментов гидродинамических полей получено замыкание EFB (Energy and Flux Budget), также допускающее поддержание перемешивания в пограничных слоях при любой устойчивости среды. В отличие от EFB подхода, модель [Соустова и др., 2020; Gladskikh et al., 2023] позволяет получить параметризацию турбулентного числа Прандтля, зависящую лишь от одного параметра (который может быть определен по данным численного моделирования или лабораторных экспериментов), описывающего анизотропию в вертикальном и горизонтальном масштабах корреляции поля плотности. Однако, как будет показано в настоящей работе, такая формулировка приводит к дополнительным ограничениям на максимальное потоковое число Ричардсона, которые не согласуются с известными данными прямого и вихреразрешающего моделирования турбулентных течений, данными измерений в устойчиво-стратифицированных пограничных слоях [Mauritsen, Svensson, 2007; Mauritsen et al., 2007; Stretch et al., 2009; Stroscio, 1982; Venayagamoorthy, Stretch, 2010; Yamada, 1985; Zilitinkevich et al., 2010].

В настоящей работе представлено обобщение модели, частный случай которой приведен в [Соустова и др., 2020; Gladskikh et al., 2023], для разрешения данной проблемы. Модифицированное замыкание позволило рассмотреть задачу об оценке чувствительности результатов моделирования биохимических характеристик водоема к параметризации двух режимов перемешивания – режима, связывающего нейтральную и устойчивую стратификацию, и режима перемешивания при сильно устойчивой стратификации. В разделе 2 статьи обсуждаются параметризация скорости диссипации кинетической и потенциальной энергий турбулентности в полной модели второго порядка, позволяющая снять ограничения в представленном ранее замыкании, и приведен вывод выражения для турбулентного числа Прандтля. Раздел 3 посвящен краткому описанию численной модели термодинамики и биогеохимии внутреннего водоема, а раздел 4 – результатам численного моделирования и анализу влияния выбора параметров замыкания (параметра анизотропии и максимального потокового числа Ричардсона) на распределение температуры, биохимических концентраций и потоков на границе с атмосферой. В заключении приведены основные выводы работы и возможные направления дальнейших исследований.

ОСНОВНЫЕ ПОЛОЖЕНИЯ МОДЕЛИ ТУРБУЛЕНТНОГО ПЕРЕНОСА И ОБОБЩЕНИЕ ЗАМЫКАНИЯ

Кратко остановимся на результатах работ [Островский, Троицкая, 1987; Соустова и др., 2020; Gladskikh et al., 2023]. Приведем выражения для турбулентных потоков импульса и массы:

$u_i^{\text{'}}{u}_j^{\text{'}}=U^2{\delta }_{ij}-LU\left(\frac{\partial <u_i>}{\partial x_j}+\frac{\partial <u_j>}{\partial x_i}\right)$, (1)

$<{\rho }^{\text{'}}{u}_i^{\text{'}}>=-LU\left(\frac{\partial <\rho >}{\partial x_i}+\frac{1}{{\rho }_0{U}^2}\left[g_i<{{\rho }^{\text{'}}}^2>-g{\beta }_i\right]\right)$. (2)

Здесь $U\simeq \sqrt{<{u_i^{\text{'}}}^2>}$ – характерный масштаб турбулентных пульсаций скорости, $\frac{1}{{\rho }_0{U}^2}\left[g_i<{{\rho }^{\text{'}}}^2>-g{\beta }_i\right]$ – компоненты вектора ускорения силы тяжести g, g_i – компоненты вектора  характеризующего пульсации давления в стратифицированной жидкости.

Заметим, что в выражение для потока массы (2) входит дополнительное слагаемое $\frac{1}{{\rho }_0{U}^2}\left[g_i<{{\rho }^{\text{'}}}^2>-g{\beta }_i\right]$, что приводит к ряду существенных отличий от результатов, полученных в рамках обычных градиентных моделей [Монин, Озмидов, 1981]. Физический смысл приведенного выше дополнительного слагаемого можно объяснить следующим образом. Рассмотрим среду с равновесным распределением плотности ρ₀ и флуктуациями поля ρ'. Предположим, что в этой среде возник жидкий объем, плотность которого превышает равновесную, то есть ρ' > 0. Он движется вниз по оси z, направленной вверх, со скоростью w' < 0. Таким образом, поток массы ρ'w' отрицательный. В противоположном случае он также будет отрицательным: жидкий объем с плотностью, меньшей плотности среды ρ' < 0, будет двигаться вверх по z, со скоростью w' > 0. Величина возвращающей силы при этом будет зависеть от формы жидкого объема, а именно, от соотношения его масштабов, что и описывается параметром анизотропии.

Параметр анизотропии для статистически однородного поля флуктуаций плотности (где компоненты вектора β имеют вид β_x = β_y = 0, β_z = <ρ'²> × R) определяется как:

$R=\left\{\begin{array}{c}\mathrm{1,} L_z\ll L_r, \\ \approx {\left(\frac{L_r}{L_z}\right)}^2, L_z\gg L_r,\end{array}\right.$

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

В [Островский, Троицкая, 1987] получена замкнутая модель турбулентного течения в стратифицированной жидкости в виде уравнений для средних – полей скорости u и плотности <ρ>, а также для кинетической энергии турбулентности k и дисперсии пульсаций плотности <ρ'²>, через которые можно выразить потенциальную энергию течения как

$П=\frac{1}{2}{\left(\frac{g}{{\rho }_0}\right)}^2\frac{<{{\rho }^{\text{'}}}^2>}{N^2}$.

Приведем последние два уравнения:

$\begin{array}{c}\frac{\partial k}{\partial t}+〈u_i〉\frac{\partial k}{\partial x_i}-L\sqrt{k}{\left(\frac{\partial 〈u_i〉}{\partial x_j}+\frac{\partial 〈u_j〉}{\partial x_i}\right)}^2-\frac{g}{{\rho }_0}L\sqrt{k}\times \\ \times \left(\frac{\partial 〈\rho 〉}{\partial z}+\frac{3}{2{\rho }_{0}k}\left(g_i〈{{\rho }^{\text{'}}}^2〉-g{\beta }_i\right)\right)+C\frac{k^{3/2}}{L}=\frac{5}{3}\frac{\partial }{\partial x_i}\left(L\sqrt{k}\frac{\partial k}{\partial x_i}\right),\end{array}$ (3)

$\begin{array}{c}\frac{\partial 〈{{\rho }^{\text{'}}}^2〉}{\partial t}+〈u_i〉\frac{\partial 〈{{\rho }^{\text{'}}}^2〉}{\partial x_i}-2\frac{\partial 〈\rho 〉}{\partial x_i}L\sqrt{k}\left(\frac{\partial \rho }{\partial x_i}+\frac{3}{2{\rho }_{0}k}\left(g_i〈{{\rho }^{\text{'}}}^2〉-g{\beta }_i\right)\right)+\\ + D\frac{k^{1/2}}{L}〈{{\rho }^{\text{'}}}^2〉=\frac{\partial }{\partial x_i}L\sqrt{k}\frac{\partial 〈{{\rho }^{\text{'}}}^2〉}{\partial x_i}.\end{array}$ (4)

Заметим, что в уравнениях (3) и (4) приняты колмогоровские гипотезы для учета диссипативных процессов (см. [Колмогоров, 1942]). При этом скорость диссипации кинетической энергии турбулентности ε и дисперсии поля плотности ε_П определяются выражениями

$\epsilon =C\frac{k^{3/2}}{L}$

и

${\epsilon }_{П}=D\frac{k^{1/2}}{L}<\rho {\text{'}}^2>$

соответственно [Rodi, 1980]. В работах [Соустова и др., 2020; Gladskikh et al., 2023] предполагается подобие временных масштабов динамического турбулентного поля t_k = k / ε и турбулентного скалярного поля t_П = <ρ'²> / ε_П, что позволяет считать $C\equiv D$ и получить выражение для турбулентного числа Прандтля, зависящего, помимо Ri, лишь от одного параметра – параметра анизотропии R:

${\mathrm{Pr}}_T\left(Ri, R\right)= \frac{\left(4-3R \right)Ri+1+{\left[{\left(\left(4-3R\right)Ri+1\right)}^2-Ri\right]}^{1/2}}{2}$. (5)

Для случая сильной устойчивости, Ri >> 1, Pr_T является линейной функцией градиентного числа Ричардсона:

${\mathrm{Pr}}_T\left(Ri, R\right)=\left(4-3R\right)Ri$. (6)

С другой стороны, согласно определению Pr_T:

${\mathrm{Pr}}_T=\frac{Ri}{R{i}_f}\to \frac{Ri}{R{i}_{f\infty }}$ при Ri >> 1, (7)

где Ri_f∞ – критическое значение потокового числа Ричардсона

$R{i}_f=-B/P=\frac{-g_i<{\rho }^{\text{'}}{u}_i^{\text{'}}>}{{\rho }_0\left(d<u_i>/dz\right)u_i^{\text{'}}{u}_j^{\text{'}}}$

величины, характеризующей отношение сдвиговой генерации КЭТ P и преобразования B КЭТ в потенциальную энергию турбулентности. Из балансового уравнения для полной энергии турбулентности в случае статистически стационарного потока следует, что Ri_f ≤ 1 и потому существует некоторое максимальное значение Ri_f∞. Например, согласно результатам численного моделирования и данным измерений в атмосферном пограничном слое [Mauritsen, Svensson, 2007; Mauritsen et al., 2007; Stretch et al., 2009; Stroscio, 1982; Venayagamoorthy, Stretch, 2010; Yamada, 1985; Zilitinkevich et al., 2010], . Можно показать (см. [Zilitinkevich et al., 2010; 2019]), что значение максимального потокового числа Ричардсона определяет зависимость градиента средней скорости течения и, в том числе, скорости диссипации КЭТ от устойчивости в теории подобия Монина–Обухова.

Предложенное выше замыкание накладывает ограничения на допустимые значения $R{i}_{f\infty }=1/\left(4-3R\right)$. Так, если положить Ri_f∞ = 0.2, то из (6) получим R = –1 / 3, что противоречит определению параметра анизотропии, который должен принимать значения 0 ≤ R ≤ 1.

С использованием данных прямого численного моделирования (DNS, Direct Numerical Simulation) турбулентного течения Куэтта при устойчивой стратификации [Mortikov et al., 2019; Zasko et al., 2022] были рассчитаны члены балансовых уравнений (3) и (4). В DNS расчетах использовалась численная модель, развиваемая в НИВЦ МГУ [Mortikov, 2016; Mortikov et al., 2019], а эксперименты проводились при числах Рейнольдса Re до 80000 при последовательном увеличении интегрального числа Ричардсона Ri_b (определяемого по разнице температуры и скорости на верхней и нижней стенках канала) до ламинаризации потока. На рис. 1 показана зависимость соотношения временных масштабов t_k и t_П от градиентного числа Ri в двух экспериментах при Re = 4 × 10⁴ и Re = 8 × 10⁴ и фиксированном Ri_b. Выполнение предположений [Островский, Троицкая, 1987; Соустова и др., 2020; Gladskikh et al., 2023] предполагает, что соотношение t_k / t_П = 1, что противоречит данным DNS расчетов, где величина

$\frac{C}{D}\ne 1$

и, более того, зависит от Ri.

 

Рис. 1. Соотношение временных масштабов ${\mathit{t}}_{\mathbf{k}}\mathbf{=}\frac{\mathbf{k}}{\mathbf{\epsilon }}$ и ${\mathit{t}}_{\mathbf{П}}\mathbf{=}\frac{\mathbf{П}}{{\mathbf{\epsilon }}_{\mathbf{П}}}$ в зависимости от градиентного числа Ричардсона по данным прямого численного моделирования

 

Введем безразмерный параметр $\alpha \equiv C/D$, характеризующий соотношение временных масштабов динамического турбулентного поля t_k и турбулентного скалярного поля t_П. Тогда решение системы уравнений (3) и (4) для случая однородного течения позволяет получить выражение для турбулентного числа Прандтля:

${\mathrm{Pr}}_T\left(Ri, R, \alpha \right)= \frac{\left(4-3R \right)\alpha Ri+1+{\left[{\left(\left(4-3R\right)\alpha Ri+1\right)}^2-4\alpha Ri\right]}^{1/2}}{2}$. (8)

Отметим, что параметры α и R также зависят от устойчивости среды. Рассмотрим случай Ri >> 1, тогда:

${\mathrm{Pr}}_T=\left(4-3R_{\infty }\right){\alpha }_{\infty }Ri$, (9)

где R_° и α_∞ предельные значения в сильно устойчивой стратификации. Тем самым из (8) с учетом (7) получим:

${\alpha }_{\infty }=\frac{1}{R{i}_{f\infty }\left(4-3R_{\infty }\right)}$.

Определение непосредственно предельных значений Ri_f° и R_° (и, вообще говоря, зависимостей α и R от Ri) потребует специализированных расчетов статистически стационарной сдвиговой турбулентности с помощью DNS и LES (Large-Eddy Simulation) моделей крайне высокого разрешения и выходит за рамки настоящей работы. Однако результаты прямого численного моделирования устойчиво-стратифицированного течения Куэтта, представленные выше, показывают, что отношение t_k / t_П относительно слабо изменяется с ростом устойчивости – в данных лабораторных экспериментов [Beguir et al., 1978] для нейтрально стратифицированного турбулентного потока отмечается близкое к постоянной значение α ≈ 1. По этим причинам далее ограничимся случаем, где параметры α и R являются фиксированными константами модели. Тогда перепишем (5) с учетом выражения для параметра α:

${\mathrm{Pr}}_T\left(Ri, R{i}_{f\infty }\right)= \frac{\frac{Ri}{R{i}_{f\infty }}+1+{\left[{\left(\frac{Ri}{R{i}_{f\infty }}+1\right)}^2-4\frac{Ri}{R{i}_{f\infty }\left(4-3R\right)}\right]}^{1/2}}{2}$. (10)

Таким образом, полученная параметризация турбулентного числа Прандтля (10) зависит от двух параметров: Ri_f° и R. Максимальное потоковое число Ричардсона Ri_f∞ задает асимптотическую линейную зависимость Pr_T от градиентного числа Ричардсона при сильной устойчивости, а параметр анизотропии R задает переход от режима нейтральной стратификации к режиму сильной устойчивости. Выражение (10) очевидно согласуется с любым значением Ri_f∞ ≤ 1.

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

ЧИСЛЕННАЯ МОДЕЛЬ ТЕРМОГИДРОДИНАМИКИ И БИОГЕОХИМИИ ВНУТРЕННЕГО ВОДОЕМА

Для исследования влияния предложенного замыкания на описание турбулентных процессов во внутренних водоемах и формирование определяемых ими характеристик использовалась разработанная авторами настоящей работы трехмерная модель термогидродинамики и биогеохимии внутреннего водоема на основе единого гидродинамического кода, объединяющего DNS-, LES- и RANS подходы для расчета геофизических турбулентных течений при высоком пространственном и временном разрешении (см., например, [Mortikov, 2016; Mortikov et al., 2019; Kadantsev et al., 2021; Gladskikh et al., 2023]). Численная модель включает уравнения гидродинамики в стратифицированном турбулентном вращающемся слое жидкости в приближении мелкой воды, а также уравнения для переноса тепла с учетом горизонтальной и вертикальной диффузии [Гладских и др., 2021]:

$$\begin{gathered} \frac{\partial u}{\partial t}=-A\left(u\right)+D_H\left(u, {\lambda }_m\right)+D_z\left(u, K_m+\nu \right)-g\frac{\partial \eta }{\partial x}-\frac{g}{{\rho }_0}\frac{\partial }{\partial x}{\int }_z^{\eta }\rho dz\text{'}+fv, \\ \frac{\partial \nu }{\partial t}=-A\left(\nu \right)+D_H\left(\nu , {\lambda }_m\right)+D_z\left(\nu , K_m+\nu \right)-g\frac{\partial \eta }{\partial y}-\frac{g}{{\rho }_0}\frac{\partial }{\partial y}{\int }_z^{\eta }\rho dz\text{'}+fu, \\ \nabla u=\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}+\frac{\partial w}{\partial z}=\mathrm{0,} \\ \frac{\partial T}{\partial t}=-A\left(T\right)+D_H\left(T, {\lambda }_h\right)+D_z\left(T, K_h+\chi \text{'}\right), \\ \rho =\rho \left(T\right), \\ \frac{\partial \eta }{\partial t}=\left(u_h\times {\nabla }_h\right)\eta =w. \end{gathered}$$

Здесь u = (u, v, w) – вектор скорости, u_h – горизонтальные компоненты, K_m(λ_m) и K_h(λ_h) – коэффициенты вертикальной (горизонтальной) турбулентной вязкости и температуропроводности соответственно; ν, χ׳ – коэффициенты молекулярной вязкости и температуропроводности, η – отклонение свободной поверхности от равновесного состояния, z – вертикальная координата, проходящая от дна водоема. Также в приведенной система A(q) – оператор адвекции:

$A\left(q\right)=\frac{\partial uq}{\partial x}+\frac{\partial vq}{\partial y}+\frac{\partial wq}{\partial z}$,

а D_H(q, λ) и D_z(q, K) – операторы горизонтальной и вертикальной диффузии с коэффициентами λ и К соответственно:

$$\begin{gathered} D_H(q, \lambda )=\frac{\partial }{\partial x}\lambda \frac{\partial q}{\partial x}+\frac{\partial }{\partial y}\lambda \frac{\partial q}{\partial y}, \\ {D}_z\left(q, K\right)=\frac{\partial }{\partial z}K\frac{\partial q}{\partial z}. \end{gathered}$$

Для описания процессов вертикального турбулентного перемешивания в модели использовалась двухпараметрическое k-ε замыкание с параметризацией турбулентного числа Прандтля (10), а турбулентные числа Шмидта Sc_T полагались равными Pr_T.

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

Описание биогеохимических процессов в модели включает в себя уравнения для расчета концентраций и потоков биохимических веществ. Уравнения описывают перенос, диффузию и реакции для таких веществ как: метан (CH₄), кислород (O₂), углекислый газ (CO₂), живые и отмершие частицы фито- и зоопланктона и т. д., и имеют вид:

$\begin{array}{c} \frac{\partial C_{\left(.\right)}}{\partial t}+\left[\frac{\partial u_i{C}_{\left(.\right)}}{\partial x_i}+\frac{\partial w{C}_{\left(.\right)}}{\partial z}+ w_s\frac{\partial C_{\left(.\right)}}{\partial z}\right]=\\ =\left[\frac{\partial }{\partial x_i}\left(\lambda +\chi \right)\frac{\partial C_{\left(.\right)}}{\partial x_i}+\frac{\partial }{\partial z}\left(K_h+\chi \right)\frac{\partial C_{\left(.\right)}}{\partial z}\right]+R\left(.\right),\end{array}$ (11)

где C_(.) – концентрации веществ, K_h(λ) и χ – коэффициенты вертикальной (горизонтальной) турбулентной и молекулярной диффузии соответственно, w_s – скорость оседания, а член R(.) описывает сколько молекул того или иного вещества было добавлено к раствору или извлечено из него в результате реакций.

Учитываются следующие процессы взаимодействия веществ – аэробное окисление метана в водной толще; фотосинтез; дыхание; биохимическое потребление кислорода в водной толще – окисление органических соединений в теле водоема; биохимическое потребление кислорода в донных отложениях – сток кислорода на окисление органических соединений в донных отложениях; отмирание фито- и зоопланктона.

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

Для описания газообмена использована т. н. “модель обновления поверхности” [MacIntyre et al., 2010], учитывающая турбулентные процессы, которые влияют на перемешивание около раздела воздух–вода:

$k_{600}=\frac{C_{\mathrm{1,}SR}\sqrt[4]{{\left.\epsilon \right|}_{z=0}\cdot \nu }}{\sqrt{600}}$.

Здесь k₆₀₀ – опорная величина для расчёта коэффициента газообмена конкретного газа, представляющая собой такой коэффициент для углекислого газа при температуре среды 20 °C: число Шмидта в этом случае равно 600. Константа C_1,SR = 0.5, ν – кинематическая вязкость, а скорость диссипации ε рассчитывается согласно модели турбулентности.

На нижней границе поток метана рассчитывается через его производство донными отложениями по формуле, приведенной в [Stepanenko et al., 2016], а поток кислорода – из логарифмического закона.

ПРИМЕНЕНИЕ ОБОБЩЕННОГО ЗАМЫКАНИЯ В МОДЕЛИ ТЕРМОГИДРОДИНАМИКИ И БИОГЕОХИМИИ ВНУТРЕННЕГО ВОДОЕМА

Для исследования вклада параметров Ri_f° и R в описание процессов турбулентного обмена и, тем самым, формирование концентраций примесей и их потоков в атмосферу были рассмотрена конфигурации модели, соответствующие описаниям озера Куйваярви и Рыбинского водохранилища. Куйваярви – небольшой водоем, расположенный в Южной Финляндии. В качестве данных атмосферного форсинга в течение всего времени расчета использовались результаты измерений в 2013 г. [Heiskanen et al., 2015; Mammarella et al., 2015] (составляющие скорости ветра, потоки коротковолновой радиации, явного и скрытого тепла). Из данных измерений также были взяты начальные условия: вертикальное распределение температуры, концентрации газов и растворенного углерода. Расчет проводился приблизительно за весь период открытой воды с мая по ноябрь. Глубина водоема: 12.5 м, согласно оценке в точке проведения измерений. Что касается Рыбинского водохранилища, то для этого крупного (площадь – 4580 км², наибольная глубина – 30 м) искусственного водоема, расположенного на севере Центральной России, нами использовались для задания атмосферного форсинга данные глобального метеорологического реанализа NCEP/NCAR в период с мая по ноябрь 2020 г.

Для задания турбулентного числа Прандтля использовалась предложенная в настоящей статье параметризация (10) с различными значениями Ri_f° и R. Рассматривались три набора параметров:

• (а) Ri_f° = 0.2; R = 0.2;
• (б) Ri_f° = 0.7; R = 0.2;
• (в) Ri_f° = 0.2; R = 0.7.

Были сопоставлены временные ряды вертикального распределения квадрата частоты плавучести N². Результаты расчетов для комбинаций (б) и (в) оказались визуально идентичны, поэтому комбинацию (в) не приводим. Результаты для (а) и (б) приведены на рис. 2 и 3.

 

Рис. 2. Временные ряды N², полученные для параметров: Ri_f∞ = 0.2 и R = 0.2 (слева), Ri_f∞ = 0.7 и R = 0.2 (справа) в конфигурации модели, соответствующей озеру Куйваярви

 

Рис. 3. Временные ряды N², полученные для параметров: Ri_f∞ = 0.2 и R = 0.2 (слева), Ri_f∞ = 0.7 и R = 0.2 (справа) в конфигурации модели, соответствующей Рыбинскому водохранилищу

 

Таким образом, выбор значения Ri_f° оказывает наибольшее влияние на распределение температуры и максимума градиента плавучести. Можно ожидать, что при увеличении Ri_f° процессы переноса вещества через термоклин будут подавляться.

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

 

Рис. 4. Вертикальное распределение температуры (вверху), кислорода (в центре), метана (внизу) на седьмой расчетный месяц при различных параметрах замыкания для озера Куйваярви

 

Показано, что описание перемешивания и тем самым формирования вертикальных распределений различных характеристик намного более чувствительно к значению критического числа Ричардсона, чем к параметру анизотропии. В частности, отметим накопление концентрации метана на глубине при Ri_f∞ = 0.7. Параметризация переходного режима между нейтральной и сильно-устойчивой стратификацией напротив слабо влияет на распределения концентраций в водоеме.

Наконец, были проанализированы значения потоков газов в атмосферу, поскольку эти характеристики представляют наибольший интерес с точки зрения вклада озер и водохранилищ в глобальные климатические процессы. Были рассмотрены потоки, полученные в конфигурации параметров (а) – flux_a и (б) – flux_b. На рис. 5 и 6 приведены соотношения flux_a /flux_b для метана и кислорода.

 

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

 

Рис. 6. Соотношения потоков метана (слева) и кислорода (справа) для Рыбинского водохранилища, полученных в конфигурациях параметров (а) и (б)

 

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

ВЫВОДЫ

В настоящей работе представлено обобщение турбулентного замыкания [26, 27], которое учитывает неградиентную поправку на поток массы и допускает возможность поддержания турбулентных пульсаций при сильно устойчивой стратификации. Обобщение модели основано на учете различий во временных масштабах динамического турбулентного поля t_k и турбулентного скалярного поля t_п. Данное предположение согласуется с результатами прямого численного моделирования турбулентного течения Куэтта при устойчивой стратификации и позволило снять ограничения в модели [Соустова и др., 2020; Gladskikh et al., 2023] на величину максимального потокового числа Ричардсона. Важным следствием обобщенного замыкания является то, что Ri_f° в общем случае определяется предельными значениями как отношения временных масштабов α = t_k / t_П, так и параметра анизотропии R.

Предполагая, что  и  лишь слабо зависят от устойчивости, получено выражение для турбулентного числа Прандтля, которое включает лишь два параметра – параметр анизотропии и максимальное потоковое число Ричардсона. Первый определяет различия в вертикальном и горизонтальном масштабах корреляции поля плотности, описывает изменения  в переходном режиме от нейтральной к сильно устойчиво стратификации. Второй параметр, максимальное потоковое число Ричардсона, напрямую связан с условием неотрицательности диссипативных членов в статистически стационарной турбулентности и определяет, в частности, асимптотический рост Pr_T с градиентным числом Ричардсона при Ri >> 1, а также форму выражений для градиента скорости и скорости диссипации в теории подобия Монина–Обухова и ее обобщениях [Zilitinkevich et al, 2019]. Введенное предположение очевидно требует дальнейшей проверки по данным DNS и LES расчетов сдвиговой турбулентности, в особенности в режимах сильной устойчивости.

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

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

Вывод обобщенного замыкания, а также численное исследование вклада его отдельных характеристик в формирование турбулентности выполнены в рамках проекта РНФ № 23–27–00002.

Sobre autores

D. Gladskikh

Institute of Applied Physics of the Russian Academy of Sciences; Lomonosov Moscow State University

Autor responsável pela correspondência
Email: daria.gladskikh@gmail.com
Rússia, Nizhny Novgorod; Moscow

E. Mortikov

Lomonosov Moscow State University; Institute of Numerical Mathematics of the Russian Academy of Sciences

Email: daria.gladskikh@gmail.com
Rússia, Moscow; Moscow

Bibliografia

Arquivos suplementares