Bulk models of sheared boundary layer convection

Full Text

1. ОБОЗНАЧЕНИЯ

A – коэффициент вовлечения $A=-\frac{B_h}{B_s}$

A₁ – константа отношения механической генерации ТКЭ к генерации плавучестью $A_1=\frac{C_s}{C_1}$

A₅ – константа интеграла диссипации механической части генерации для конвекции со сдвигом $A_5=\frac{V_{\text{*}}^3}{2}\left(1-h^{-1}{\int }_0^h\epsilon dz\right)$

b – плавучесть

Db – скачок плавучести в слое вовлечения КПС

B – поток плавучести

C₁ – константа интеграла диссипации свободной конвекции $C_1=1-2W_{*}^{-3}\underset{0}{\overset{h}{\int }}\epsilon dz$

c₂ – константа пропорциональности между средней ТКЭ в КПС и масштабом генерации для случая свободной конвекции $c_2=\frac{{\displaystyle {\int }_0^h{E}_{k}dz}}{h{W}_{\text{*}}^2}$ $C_{TZ}=\frac{10}{3}{c}_2$

C₃, C₄ – константы пропорциональности турбулентного потока ТКЭ из слоя вовлечения в свободную атмосферу для свободной конвекции и конвекции со сдвигом соответственно

c₆ – константа пропорциональности между средней ТКЭ в КПС и масштабом механической генерации $c_2=\frac{{\displaystyle {\int }_0^h{E}_{k}dz-c_2{W}_{\text{*}}^2}}{h{V}_{\text{*}}^2}$

C_D – аэродинамический коэффициент обмена с поверхностью

E_K – турбулентная кинетическая энергия (ТКЭ)

Df – ширина фильтра LES-модели

f – параметр Кориолиса

F – турбулентный поток ТКЭ $F=\overline{u_i^{\text{'}}{u}_i^{\text{'}}w\text{'}}+\frac{1}{\rho }\overline{p\text{'}w\text{'}}$

h – высота КПС

l – половина толщины слоя вовлечения в интегральных моделях КПС нулевого порядка

N – частота Вяйсяля–Брента

S – механическая генерация ТКЭ, $S=\tau \cdot \frac{\partial u}{\partial z}$

u_* – скорость трения, $u_{\text{*}}^2={\tau }_s{\rho }^{-1}$

U_m – модуль скорости ветра в перемешанном подслое КПС

U_g – модуль геострофической скорости ветра

DU – скачок модуля скорости ветра в слое вовлечения

V_* – масштаб скорости, определяемый механической генерацией ТКЭ, $V_{\text{*}}^3=2u_{\text{*}}{U}_m$

W_* – конвективный масштаб скорости Дирдорффа, $W_{\text{*}}^3=B_{s}h$

W_e – комбинированный масштаб скорости, $W_e^3=W_{\text{*}}^3+{\eta }^3{u}_{\text{*}}^3$

W_new – комбинированный масштаб скорости, $W_{new}^3={\alpha }_1{W}_{\text{*}}^3+{\alpha }_5{V}_{\text{*}}^3$

W_m – комбинированный масштаб скорости, $W_m^3=c_2{W}_{\text{*}}^3+c_6{V}_{\text{*}}^3$

$\delta$ – толщина слоя вовлечения

$\epsilon$ – скорость диссипации ТКЭ

$\zeta$ – безразмерная высота в КПС, $\zeta =\frac{z}{h}$

$\tau$ – вектор вертикального потока импульса

$\Phi$ – универсальная безразмерная функция безразмерной высоты КПС

2. ВВЕДЕНИЕ

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

Данные наблюдений и лабораторных экспериментов показывают, что в целом эволюция высоты КПС хорошо описывается интегральными («балковыми») соотношениями [6–8], поэтому целесообразно уточнение существующих параметризаций конвекции на их основе. Интегральные модели КПС основаны на системе уравнений, неизвестными переменными которой являются высота КПС h, перепады метеовеличин в слое вовлечения и коэффициент вовлечения. Коэффициент вовлечения определяется как отношение экстремума потока плавучести в слое вовлечения к значению потока на поверхности $A=\frac{-B_h}{B_s}$.

В последние десятилетия наибольшее распространение получили два вида интегральных моделей КПС — модели «нулевого» [9] и «первого» порядка [10]. Они различаются в том, каким образом представлен перепад плавучести и скорости на границе КПС в соответствующих идеализированных профилях. На рис. 1 представлены эти профили для плавучести и горизонтальной скорости в моделях «нулевого» порядка. В них градиенты плавучести и скорости в слое вовлечения аппроксимируются разрывами этих величин на высоте h.

 

Рис. 1. Профили плавучести (a), потока плавучести (b), горизонтальной скорости ветра (c) и вертикального потока импульса в конвективном пограничном слое (КПС), развивающегося на фоне устойчиво стратифицированной свободной атмосферы, при наличии фонового ветра. Сплошной черной линией представлены схематичные профили, используемые в интегральных моделях КПС, пунктирной линией — осредненные по времени профили по данным LES-расчетов данной работы (см. разделы 4 и 5).

 

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

В целом точность моделей нулевого порядка практически совпадает с точностью моделей первого порядка в случаях свободной конвекции и конвекции со сдвигом при сильно устойчивой стратификации в свободной атмосфере [12, 13], ухудшаясь лишь при слабой устойчивости в свободной атмосфере или при наличии значительного градиента геострофического ветра [14]. Модели обоих типов также приводят к большим ошибкам в определении высоты КПС h и перепада скорости в слое вовлечения DU над сильно шероховатой поверхностью [15].

В предложенных к настоящему моменту моделях нулевого порядка не возникает связи между скоростью трения на поверхности u_* и сдвигом скорости в слое вовлечения DU, что, в частности, и обусловливает невысокую точность расчетов над сильно шероховатой поверхностью. Модели первого порядка эту зависимость учитывают, однако неопределенности параметризации, привносимые с новой переменной d и дополнительные математические сложности в решении получающейся системы уравнений в значительной степени сводят на нет это достоинство. Поэтому представляется перспективным предложенный в работах [16, 17] подход, при котором в модели нулевого порядка реализуется связь u_* с DU (подробно излагается в разделе 3.3).

В связи со сказанным выше целью данной работы является проверка интегральной модели проникающей конвекции со сдвигом скорости [16, 17] по данным численных экспериментов с вихреразрешающей моделью (LES — Large Eddy Simulation) Института вычислительной математики (ИВМ) РАН, в которой энергетически преобладающая часть турбулентных вихрей воспроизводится явно; эта задача подразумевает, среди прочего, определение эмпирических констант интегральной модели.

3. АНАЛИТИЧЕСКИЕ ИНТЕГРАЛЬНЫЕ МОДЕЛИ КПС

3.1. Уравнение притока тепла

Основу интегральных моделей КПС составляет интегрированное по высоте КПС уравнение притока тепла в приближении горизонтальной однородности и в предположении, что поток плавучести линейно меняется с высотой, а частота Брента–Вяйсяля N над КПС постоянная по высоте. Полученное таким образом уравнение имеет вид:

$\frac{d\left(\frac{1}{2}{N}^2{h}^2\right)}{dt}=B_s-B_h=\left(1+A\right)B_s,$ (1)

где t — время, h — высота КПС, B_s — приземный поток плавучести (принятый для простоты постоянным по времени), B_h — поток плавучести на границе КПС. Так как внутри КПС плавучесть постоянна, поток плавучести на верхней границе выражается через скачок плавучести, Db, по формуле $B_h=-\Delta b\frac{dh}{dt}$ [18], где $\Delta b\frac{dh}{dt}$ — скорость роста КПС или скорость вовлечения. При этом разность B_s–B_hв правой части (1), выражается через коэффициент вовлечения:

$A=\frac{-B_h}{B_s}=\frac{l}{h-l}.$ (2)

Если считать этот коэффициент постоянным и известным из данных наблюдений и ограничиваться режимом свободной конвекции, толщина КПС полностью описывается уравнением (1). Такой подход применялся в начале развития интегральных моделей КПС (см. [18–20]). При отсутствии вовлечения A(t) = 0 решение уравнения (1) сводится к классической формуле Н.Н. Зубова $h=\sqrt{\frac{2B_{s}t}{N^2}}$ [19]. При A = const > 0 уравнение (1) принимает вид

$\frac{dh}{dt}=\left(1+A\right)\frac{B_s}{N^{2}h}.$ (3)

Наблюдаемые при измерениях в атмосфере и в лабораторных экспериментах значения коэффициента вовлечения лежат в примерном диапазоне (см. табл. 1 в [21]), то есть практически во всем интервале между предельными режимами B_h = 0 [18] и B_h = –B_s [22]. Естественно, что при таком разбросе уравнение (3) с постоянным A, вообще говоря, неудовлетворительно, хотя наблюдаемые значения A в атмосферных КПС обычно близко к значению 0.2 [6, 8, 21, 23]. Ниже будет показано, что значение получается из анализа баланса турбулентной кинетической энергии для хорошо развитого КПС при свободной конвекции, если предположить, что интегральная по КПС генерация ТКЭ силами плавучести полностью компенсируется диссипацией. Таким образом, такие простые интегральные модели КПС, как формула (3), во многих случаях хорошо воспроизводят высоту КПС.

 

Таблица 1. Значения констант C_1, C_TZ, C₃ под данным LES-экспериментов и оценки из литературы

Константы	C1	CTZ	C3
CBL1	0.24	0.65	0.02
CBL2	0.23	0.66	0.01
CBL3	0.22	0.68	0.005
[16]	0.2	0.8	0.1
[14]	0.2	0.7	0.02

 

3.2. Уравнение интегрального баланса ТКЭ: режим свободной конвекции

В общем случае коэффициент вовлечения является переменной во времени величиной, определение которой требует добавления еще одного уравнения в систему, образующую интегральную модель КПС. Для его вывода выписывается уравнение баланса ТКЭ на границе перемешанного слоя и слоя вовлечения z = h – l [24], или уравнение баланса (ТКЭ) интегрируется по высоте для горизонтально однородного пограничного слоя [13, 16, 25]:

$\begin{array}{c}{\int }_0^h\frac{\partial E_k}{\partial t}dz={\int }_0^{h}Bdz+{\int }_0^h\tau \cdot \frac{\partial u}{\partial z}dz-\\ -{\int }_0^h\frac{\partial F}{\partial z}dz-{\int }_0^h\epsilon dz,\end{array}$ (4)

где $\epsilon$ — диссипация турбулентной кинетической энергии, $F=\overline{u{\text{'}}_{i}u{\text{'}}_{i}w\text{'}}+\frac{1}{\rho }\overline{p\text{'}w\text{'}}$ — вертикальный турбулентный поток ТКЭ, B — поток плавучести, $\tau =\left(\overline{u\text{'}w\text{'}},\overline{v\text{'}w\text{'}},\overline{w\text{'}w\text{'}}\right)$ — вертикальный поток импульса, u = (u, v, w) — вектор скорости ветра. При свободной конвекции средняя генерация ТКЭ за счет сдвига скорости ветра равна нулю. При этом условии $S=\tau \cdot \frac{\partial u}{\partial z}=0$, и, используя теорию самоподобия конвективной турбулентности Дирдорфа [20], профили ТКЭ и диссипации можно представить через универсальные функции от безразмерной высоты $\zeta =z/h$:

$E_k=W_{\text{*}}^2{\Phi }_{EK}\left(\zeta \right),\epsilon =\left(W_{\text{*}}^3/h\right){\Phi }_{DK}\left(\zeta \right),$ (5)

где $W_{\text{*}}={\left(B_{s}h\right)}^{1/3}$ — дирдорфовский масштаб скорости, характеризующий типичное значение вертикальной и горизонтальной скорости в самоорганизующихся термиках КПС. Третий член в правой части (4) обычно считается малым, так как на верхней и нижней границе КПС турбулентные потоки ${\int }_0^h\frac{\partial F}{\partial z}dz={\left.\left(\overline{u{\text{'}}_{i}u{\text{'}}_{i}w\text{'}}+\frac{1}{\rho }\overline{p\text{'}w\text{'}}\right)\right|}_0^h$ противоположны по знаку и близки по модулю. Зилитинкевич [16, 17] учел возможное излучение в свободную атмосферу гравитационных волн, возбуждаемых «ударами» термиков в устойчиво-стратифицированный слой вовлечения, и определил как вертикальный поток энергии за счет этих волн. Тогда, согласно линейной теории [26],

$F_{z=h}\propto {\lambda }^2\Lambda N^3,$ (6)

где l — длина, а — амплитуда волны. Принимая, согласно [27], $\Lambda ~\lambda ~\frac{A}{1+A}h$, получаем:

${\int }_0^h\frac{\partial F}{\partial z}dz=F_h-F_0=F_h=C_3{N}^3{h}^3{\left(\frac{A}{1+A}\right)}^3.$ (7)

Учитывая линейность профиля потока плавучести (см. рис. 1), ${\int }_0^{h}Bdz=B_{s}h\left(\frac{1-A}{2}\right)$, получаем интегральное уравнение баланса ТКЭ для режима свободной конвекции:

$C_{TZ}{W}_{\text{*}}^{-1}\frac{dh}{dt}=C_1-A-C_3\frac{N^3}{B_s^{{}^3{/}_2}}{\left(\frac{A}{1+A}\right)}^3,$ (8)

где $C_{TZ}=\frac{10}{3}{\int }_0^1{\Phi }_{EK}\left(\zeta \right)d\zeta =\frac{10/3}{c_2}, C_1=1-2{\alpha }_1=1-2{\int }_0^1{\Phi }_{DK}\left(\zeta \right)d\zeta$ — безразмерные энергетические константы. Если упростить (4), предположив что запас ТКЭ в КПС не меняется ${\int }_0^h\frac{\partial E_k}{\partial t}dz=0$ и F(h) = 0, то есть суммарная генерация ТКЭ силами плавучести и суммарная диссипация уравновешивают друг друга, то A = C₁ = const. Согласно известным в литературе оценкам (см. табл. 1), а также по данным наших LES-экспериментов (раздел 5), $0.15\le C_1\le 0.25$, чему соответствуют значения $0.15\le A\le 0.25$, неплохо согласующиеся с упомянутыми выше прямыми энергетическими оценками.

3.3. Уравнение баланса ТКЭ: режим конвекции со сдвигом скорости

При наличии средней скорости ветра член ${\int }_0^{h}Sdz$ существен, и гипотетические соотношения (5) становятся не удовлетворительными. За последние двадцать лет были предложены различные способы параметризации членов интегрального уравнения баланса ТКЭ для конвекции со сдвигом скорости. Часть авторов [9, 28, 29] пренебрегает левой частью (4), аргументируя это тем, что в развитом КПС запас ТКЭ мало меняется. Для верификации таких интегральных моделей используют результаты LES-экспериментов значительной продолжительности, исключая при этом из анализа первые 3–4 ч модельного времени. Параметризации, учитывающие тенденцию запаса ТКЭ [11, 15, 25, 30, 31], обычно используют в качестве масштаба скорости величину W_e, определяемую как $W_e^3=W_{\text{*}}^3+{\eta }^3{u}_{\text{*}}^3$, где $u_{\text{*}}=\sqrt{\frac{{\tau }_s}{\rho }}$ — скорость трения в приземном слое, ${\eta }^3=\frac{{\displaystyle {\int }_0^{h}Sdz-{\int }_0^h{\epsilon }_{s}dz}}{u_{\text{*}}^3}$ — доля сдвиговой генерации, нескомпенсированная локальной диссипацией. Привлекая дополнительное предположение о том, что производная по времени от W_e ведет себя так же, как и производная от W_* в случае свободной конвекции, получается соотношение:

${\int }_0^h\frac{\partial E_k}{\partial t}dz=C_T\frac{W_e^2}{h}\frac{dh}{dt}$ (9)

не предусматривающее самоподобного профиля $E_k=W_e^2{\Phi }_E\left(\zeta \right)$. Зилинтинкевич [16, 17] предложил, что в случае конвекции со сдвигом скорости ветра вертикальный профиль ТКЭ можно представить в виде суммы членов, описывающих вклад турбулентности механического (сдвигового) и конвективного происхождения:

$E_k=W_{\text{*}}^2{\Phi }_{EK}\left(\zeta \right)+V_{\text{*}}^2{\Phi }_{SK}\left(\zeta \right),$ (10)

где $V_{\text{*}}={\left(2u_{\text{*}}^2{U}_m\right)}^{{}^1{/}_3}$ — сдвиговый масштаб скорости, U_m — модуль средней скорости внутри КПС. Для вычисления левой части (4) в [16, 17] принималось, что $\frac{d{V}_{\text{*}}}{dt}=0$, что, однако, не согласуется с данными LES-экспериментов.

Существующие интегральные модели сдвигового КПС можно также разделить на две группы по способу аппроксимации интеграла диссипации по толщине КПС, т. е. последнего члена в правой части (4). В одной группе используется предположение, что каждой из составляющих баланса ТКЭ в КПС — плавучести, сдвигу и турбулентному переносу — соответствует своя доля диссипации, а интегральная диссипация ТКЭ представляется как линейная комбинация соответствующих интегралов. Так, согласно [16, 17]:

${\int }_0^h\epsilon dz={\alpha }_1{W}_{\text{*}}^3+{\alpha }_5{V}_{\text{*}}^3,$ (11)

где ${\alpha }_1=\frac{1-C_1}{2}$ — константа, известная из экспериментов для режима свободной конвекции, а a₅ — дополнительная константа. В другой группе моделей постулируется, что интегральная диссипация в КПС пропорциональна интегральной продукции, в сочетании с приближением сильно перемешанного КПС. Это предположение сводит интегральное уравнение ТКЭ к балансу между плавучестью и сдвигом скорости на границе слоя вовлечения. Тогда, используя в качестве альтернативных масштабов скорости W_* и u_* [32], можно выписать следующее соотношение:

${\int }_0^{h}Bdz+{\int }_0^{h}Sdz-{\int }_0^h\epsilon dz=C_1\left(W_{\text{*}}^3+A_1{u}_{\text{*}}^3\right),$ (12)

где, A_{1 =} C_s / C_1, C_s — доля сдвиговой генерации ТКЭ, нескомпенсированная диссипацией. Некоторые авторы предполагают при этом, что $A_1={\eta }^3$. Тогда правую часть (12) можно записать в виде:

${\int }_0^{h}Bdz+{\int }_0^{h}Sdz-{\int }_0^h\epsilon dz=C_1{W}_e^3.$ (13)

Перейдем к вычислению второго члена в правой части (4). Разбивая КПС на три подслоя: приземный, перемешанный и слой вовлечения — получаем:

${\int }_0^{h}Sdz={\int }_0^{h_s}Sdz+{\int }_{h_s}^{h-l}Sdz+{\int }_{h-l}^{h}Sdz.$ (14)

Первое слагаемое в правой части пропорционально кубу скорости трения ${\int }_0^{h_s}\tau \cdot \frac{\partial u}{\partial z}dz\propto u_{\text{*}}^3$ [24, 33]. Второе слагаемое мало или равно нулю из-за практического отсутствия градиента скорости в перемешанном слое. В интегральных моделях КПС «нулевого» порядка для определения сдвиговой генерации в слое вовлечения используется выражение для потока импульса $\overline{u^{\text{'}}{w}^{\text{'}}}=\Delta U\frac{dh}{dt}$, а градиент скорости оценивается как $\frac{\partial U}{\partial z}\simeq \frac{\Delta U}{l},(l=\frac{A}{1+A}h$, хотя некоторые модели используют в качестве масштаба длины высоту КПС [24]). При использовании масштаба $l=\frac{A}{1+A}h$, интеграл сдвиговой генерации выражается в виде:

${\int }_0^h\tau \cdot \frac{\partial u}{\partial z}dz\simeq C_S{u}_{\text{*}}^3+C_M\frac{{\left(\Delta U\right)}^2}{l}\frac{dh}{dt},$ (15)

где, по данным LES-моделирования и лабораторных экспериментов, $\eta =\sqrt[3]{C_s/C_1}=2$ и С_M = 0.7. Формула (15) имеет существенный недостаток, поскольку ведет к выражению $A~\frac{C_1}{1+C_{TZ}\Delta b-C_M\Delta U}$(подробный вывод см. в [34]), и при принятых значениях С_TZ и С_M и скачках Db, DU типичных для хорошо развитых КПС сильно завышает коэффициент вовлечения или даже меняет его знак. Интегральные модели КПС первого порядка менее подвержены подобному эффекту из-за интерполяции профиля потока импульса в слое вовлечения и появлению в связи с этим члена $~u_{\text{*}}^2\Delta U$ в (15), где DU = U_g – U_m. Этот член может фигурировать и в моделях нулевого порядка; например, в работе [16] такое же слагаемое появляется в результате аппроксимации Оценка правомерности этой аппроксимации ${\int }_0^{h}Sdz\simeq 2u_{*}^2{U}_m.$ также будет проверена по данным LES-экспериментов, обсуждаемых ниже.

В работе [16] уравнение интегрального баланса ТКЭ включает поток ТКЭ из слоя вовлечения за счет уходящих вверх гравитационных волн F_h. При этом в качестве масштаба длины волны принимается $\Lambda ~h$ [27], в отличие от режима свободной конвекции, где $\Lambda$ считается пропорциональной толщине слоя вовлечения $\Lambda ~\frac{A}{1+A}h$. Отсюда следует:

$F_h=C_4{\lambda }^2\Lambda N^3=C_4{h}^3{N}^3{\left(\frac{A}{1+A}\right)}^2.$ (16)

В результате, для случая конвекции со сдвигом интегральное уравнение (4), согласно [16], принимает вид:

$C_{TZ}{W}_{\text{*}}^{-1}\frac{dh}{dt}=C_1-A-C_4\frac{N^3}{B_s^{{}^3{/}_2}}{\left(\frac{A}{1+A}\right)}^2+A_5{\left(\frac{V_{\text{*}}}{W_{\text{*}}}\right)}^3,$ (17)

где введена константа $A_5=\frac{1-{\alpha }_5}{2}$.

3.4. Оптимальная модель КПС

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

Интегральный баланс тепла (раздел 3.1) выражается общепризнанным уравнением:

$\frac{d}{dt}\left[\frac{1}{2}{N}^2{h}^2-h{B}_{s}A{\left(\frac{dh}{dt}\right)}^{-1}\right]=B_s,$ (18)

которое эквивалентно (1), (2).

Интегральный бюджет ТКЭ (раздел 3.3) согласно [16, 17] с поправкой на уточненную формулу (32) в разделе 5 имеет вид:

$\begin{array}{c}{C}_{TZ}{W}_{\text{*}}^{-1}\frac{dh}{dt}+C_6{W}_{\text{*}}^{-3}\frac{d\left(V_{\text{*}}^{2}h\right)}{dt}=\\ =C_1-A-C_4\frac{N^3}{B_s^{{}^3{/}_2}}{\left(\frac{A}{1+A}\right)}^2+A_5{\left(\frac{V_{\text{*}}}{W_{\text{*}}}\right)}^3,\end{array}$ (19)

где $W_{\text{*}}={\left(B_{s}h\right)}^{{}^1{/}_3}$ — конвективный масштаб скорости Дирдорфа [20], $V_{\text{*}}={\left(2u_{\text{*}}^2{U}_m\right)}^{{}^1{/}_3}$ — масштаб скорости, характеризующий механическую (сдвиговую) генерацию турбулентности [16, 17]. Интегральный бюджет количества движения выражается парой уравнений, получаемых путем почленного интегрирования уравнений движения:

$\begin{array}{c}h\frac{d{v}_m}{dt}-\Delta v\frac{dh}{dt}=\\ =-C_D{u}_m\sqrt{u_m^2+v_m^2}+fh\left(u_m-u_g\right),\end{array}$ (20)

$\begin{array}{c}h\frac{d{u}_m}{dt}-\Delta u\frac{dh}{dt}=\\ =-C_D{v}_m\sqrt{u_m^2+v_m^2}-fh\left(v_m-v_g\right),\end{array}$ (21)

где u_m и v_m — средние значения составляющих скорости ветра внутри КПС, Du и Dv — скачки скорости на верхней границе КПС, f — параметр Кориолиса, u_g и v_g — составляющие скорости геострофического ветра. При выводе уравнений (20), (21), профили составляющих вертикального потока импульса t_x и t_y принимаются линейными, а приземные значения этих составляющих определены по аэродинамическим формулам:

${\tau }_{xs}=C_D{u}_m\sqrt{u_m^2+v_m^2},$ (22)

${\tau }_{ys}=C_D{v}_m\sqrt{u_m^2+v_m^2},$ (23)

где C_D — коэффициент сопротивления. Ниже эта модель проверяется и калибруется по данным вихреразрешающего моделирования.

4. ПОСТАНОВКА ЭКСПЕРИМЕНТОВ НА ОСНОВЕ LES-МОДЕЛИ

Детальное описание вихреразрешающей модели ИВМ РАН (далее LES INM RAS), используемой в данной работе, а также результаты ее сравнения с другими моделями для различных видов течений можно найти в [35–37]. Численное решение системы уравнений гидродинамики в приближении Буссинеска реализовано с помощью метода конечных разностей, для аппроксимации по пространству используется консервативная по импульсу и энергии схема 4-го порядка точности [38], для которой используется разнесенная С-сетка Аракавы, для аппроксимации по времени используется явная схема «предиктор–корректор» Адамса–Бэшфорта второго порядка. Программный код модели реализован для многопроцессорных ЭВМ с распределенной памятью с помощью технологии MPI.

Для оценки применимости различных аппроксимаций членов интегрального баланса ТКЭ и других предположений из [16], приводящих к системе уравнений интегральной модели КПС, описанной выше, был проведен ряд численных экспериментов, в которых КПС развивался под действием постоянного по времени и горизонтали положительного потока тепла на подстилающей поверхности H_{S =}0.35 Кмc^–1 что соответствует потоку плавучести B_{s =} 0.22 мc^–1. На верхней границе расчетной области поток тепла был задан равным нулю. Граничные условия для скорости на нижней и верхней границах задавались в терминах потока импульса, рассчитанного с использованием логарифмического профиля между границей и ближайшим модельным уровнем с шероховатостью поверхности равной $z_{0\text{BOT}}=0.1$и $z_{0\text{TOP}}={10}^{-9}$ соответственно. Во всех экспериментах использовались периодические граничные условия по горизонтали. Чтобы исключить отражение от верхней границы гравитационных волн, распространяющихся вверх от слоя вовлечения, в слое $z_{TOP}-30\Delta z<z<z_{TOP}$ к решению применялось рэлеевское демпфирование со временем релаксации $t_d=1.3\cdot 2\cdot \Delta t$, где $\Delta t$ — шаг модели по времени, $\Delta z$ — разрешение модели по вертикали. Количество расчетных узлов модели составляло N_{x =} 512, N_{y =} 512, N_{z =} 394 по координатам x, y, z соответственно.

Для корректного воспроизведения процессов в слое вовлечения и независимости статистик первого и второго порядка от разрешения модели необходимо соблюсти условие $\frac{h}{C_s\Delta f}>\mathrm{56,}$ где C_s — константа Смагоринского, Df — ширина фильтра [39]. При типичных значениях A ₌ 0.2, C_{S =} 1.2 и оценке толщины слоя вовлечения $\delta =\frac{Ah}{\left(1+A\right)}$ это условие соответствует $11.2\Delta f=\delta .$  Для исключения влияния периодических граничных условий на характеристики КПС достаточным считается выбор минимальных размеров расчетной области $\frac{\left(L_x,L_y,L_z\right)}{h_0}=\left(\mathrm{5,5,2}\right),$ где h₀ — начальная высота КПС. Динамическое замыкание для подсеточных напряжений в модели LES INM RAS, предполагащее изотропность неявного фильтра, налагало ограничение на шаг сетки модели $\frac{\Delta x}{\Delta y}=\frac{\Delta y}{\Delta z}\simeq 1$, где Dx, Dy, Dz, — шаги сетки модели по координатам x, y, z. Использовалось разрешение по пространству Dx = Dy = Dz = 5. Размеры расчетной области составили (L_x, L_y, L_z) = (2.56, 2.56, 1.92 км), и начальная высота КПС h₀ = 250 м.

Для приведенных выше параметров экспериментов $\frac{h}{C_s\Delta f}\simeq 30$. В отличие от LES-модели, использованной в [39], коэффициент Смагоринского в LES INM RAS рассчитывается динамически, поэтому условие, налагаемое на вертикальное разрешение, можно несколько ослабить. Напомним, что задачей LES-экспериментов была проверка гипотетических соотношений (5), (7), (10), (11), (16) и определение соответствующих эмпирических констант ${\alpha }_1,C_{TZ},C_3$ для случая свободной конвекции и C₄, A₄ для случая конвекции со сдвигом скорости:

${\alpha }_1=\frac{{\displaystyle {\int }_0^h\epsilon dz}}{W_{\text{*}}^3}=\frac{\left(1-C_1\right)}{2},$ (24)

$c_2=\frac{{\displaystyle {\int }_0^h{E}_{k}dz}}{h{W}_{\text{*}}^2}=\frac{3}{10}{C}_{TZ},$ (25)

$C_3=\frac{F_h}{N^3{h}^3{\left(\frac{A}{1+A}\right)}^3},$ (26)

$C_4=\frac{F_h}{N^3{h}^3{\left(\frac{A}{1+A}\right)}^2},$ (27)

${\alpha }_5=\frac{{\displaystyle {\int }_0^h\epsilon dz-{\alpha }_1{W}^3}}{V_{\text{*}}^3}=\frac{\left(1-A_5\right)}{2}.$ (28)

Определение этих констант осуществлялось следующим образом. Выполнялась серия экспериментов для безветренной конвекции при стратификации в свободной атмосфере, характеризующийся различным градиентом потенциальной температуры выше КПС: g_{0 = 0.01, 0.02, 0.03} jr^–1, обозначаемые далее как CBL1, CBL2 и CBL3. Для воспроизведения конвекции со сдвигом скорости проведены два эксперимента, в которых варьировались шероховатость поверхности z₀ и скорость геострофического ветра U_g : z₀ = 0.001 м — в эксперименте CBLu5, и z₀ = 0.1 м — в эксперименте CBLu7. Cтратификация в свободной атмосфере принималась как в эксперименте CBL3, а при вычислении констант C₄ и A₅ по данным экспериментов CBLu5 и CBLu7, использовались значения C₁, C_TZ полученные ранее в экспериментах со свободной конвекцией.

В анализе результатов LES-экспериментов ключевую роль играет определение высоты КПС. В данной работе высота КПС определялась как высота минимального значения потока плавучести, вычисленного осреднением по горизонтали внутри расчетной области вихреразрешающей модели.

5. РЕЗУЛЬТАТЫ ЭКСПЕРИМЕНТОВ

Основные результаты численных экспериментов при безветрии иллюстрируются рисунками 1 и 2. На рис. 2 показаны временные ряды значений констант C₁, C_TZ и C₃, рассчитанных по формулам (24)–(26), а также их значения из работы [16]. Дисперсия временных рядов относительно средних значений существенна, однако выраженных трендов не наблюдается, что подтверждает гипотезы подобия и дает основание принимать величины C₁, C_TZ, C₃ константами. Значения C₁, C_TZ, C₃ по вихреразрешающим расчетам заметно отличаются от предыдущих оценок: так, C₁ в среднем по трем экспериментам выше, а C_TZ — ниже значений, предлагаемых в [16, 17]. Оценка константы C₃ по данным экспериментов CBL1, CBL2, CBL3 на порядок меньше, чем в [16], но согласуется с результатами других LES-экспериментов [14]. Значения констант по данным экспериментов CBL1, CBL2, CBL3 и их значения из предыдущих работ приведены в табл. 1.

 

Рис. 2. Изменения во времени высоты КПС по данным LES-моделирования свободной конвекции (эксперименты CBL1, CBL2, CBL3 – обозначены маркерами). Линиями представлены решения уравнений интегральной модели [16] с двумя различными наборами констант C1, CTZ, C3: суффикс old обозначает набор констант из оригинальной работы [16], суффикс new – оценки констант из LES-расчетов данной работы по формулам (24)–(26) (подробнее см. табл. 1).

 

На рис. 2 показаны результатам расчета высоты КПС по модели LES и по интегральной модели [16]. Как видно, во всех трех экспериментах высота КПС растет с течением времени t пропорционально , и при этом наблюдаются колебания h с периодом $~h/W_{\text{*}}~{10}^2·c$.

Напомним, что при безветренной конвекции справедлива оценка суммарной энергии ТКЭ:

${\int }_0^h{E}_{k}dz=C_{TZ}h{W}_{\text{*}}^2.$ (29)

Потому логично ожидать, что при конвекции со сдвигом скорости разность $h^{-1}{\int }_0^h{E}_{k}dz-C_{TZ}{W}_{\text{*}}^2$ окажется пропорциональной $V_{\text{*}}^2$. Как показано на рис. 3, отношение этих величин как функция времени по данным эксперимента CBLu7 действительно практически постоянно:

$\frac{h^{-1}{{\displaystyle \int }}_0^h{E}_{k}dz-C_{TZ}{W}_{\text{*}}^2}{V_{\text{*}}^2}=C_6\approx \mathrm{0.4,}$ (30)

 

Рис. 3. Изменения во времени констант C1, CTZ, C3 по данным экспериментов CBL1, CBL2, CBL3.

 

где C₆ — константа пропорциональности. На рис. 4 показана временная эволюция средней ТКЭ $h^{-1}{\int }_0^h{E}_{k}dz$ в КПС со сдвигом скорости по мере роста КПС, для нормировки ТКЭ используются: масштаб Дирдорфа $W_{\text{*}}^2$, характеризующий ТКЭ конвективного происхождения, и масштаб, введенный Зилитинкевичем [16, 17], $V_{\text{*}}={\left(2u_{\text{*}}^2{U}_m\right)}^{{}^1{/}_3}$, характеризующий ТКЭ, порождаемую сдвигом скорости, и их линейная комбинация $W_{+}=C_{TZ}{W}_{\text{*}}^2+C_6{V}_{\text{*}}^2.$ Как видно из рисунка, первое и второе отношение существенно меняются, а нормировка с помощью обобщенного масштаба делает безразмерную среднюю ТКЭ в КПС почти постоянной. Отсюда следует полезное соотношение:

${\int }_0^h{E}_{k}dz=h\left(C_{TZ}{W}_{\text{*}}^2+C_6{V}_{\text{*}}^2\right).$(31)

 

Рис. 4. Изменения во времени отношения интеграла ТКЭ в КПС и различных интегральных масштабов, по данным эксперимента CBLu7.

 

Согласно (31), левая часть уравнения интегрального баланса ТКЭ записывается в виде:

$\frac{d}{dt}\left({\int }_0^h,E_k,d,z\right)=C_{TZ}\frac{d}{dt}\left(h{W}_{\text{*}}^2\right)+C_6\frac{d}{dt}\left(h{V}_{\text{*}}^2\right).$ (32)

Приведенные выше результаты связаны с разнонаправленным изменением во времени масштабов W_* и V_*. Масштаб Дирдорфа W_* растет вместе с h, а V_* снижается из-за уменьшения трения в приземном слое и средней скорости в КПС (не показано). Существенное уменьшение масштаба сдвиговой части генерации ТКЭ V_* во времени нарушает предположение из [16] о неизменности этой величины $\frac{d{V}_{\text{*}}^2}{dt}=0$. Таким образом, левая часть интегрального уравнения баланса ТКЭ должна иметь вид:

$\begin{array}{c}\frac{d}{dt}\left({\int }_0^h,E_k,d,z\right)=\frac{d}{dt}h\left(c_2{W}_{\text{*}}^2+C_6{V}_{\text{*}}^2\right)=\\ =C_{TZ}{W}_{\text{*}}^2\frac{dh}{dt}+C_6\frac{d\left(h{V}_{\text{*}}^2\right)}{dt},\end{array}$ (33)

где справа появляется член $C_6\frac{d\left(h{V}_{\text{*}}^2\right)}{dt}$, отсутствующий в оригинальном уравнении (17). Выражение ТКЭ через масштаб W_e и вынесение соответствующей константы пропорциональности C_T за знак производной по времени, как в (10), не отражает универсального свойства интегрального баланса TKЭ, поскольку соответствующая кривая на рис. 3 падает со временем. Это вероятнее всего является следствием того, что для упрощения масштаба W_e в уравнении (13) из оценки сдвиговой генерации (15) отбрасывается член с DU, и получаемый масштаб скорости W_e растет быстрее, чем средняя ТКЭ в КПС.

Рассмотрим теперь, насколько приведенные выше масштабы подходят для оценки интеграла диссипации ТКЭ. Для этого обратимся к рис. 5, на котором кривые отношения интеграла диссипации к различным масштабам построены таким же способом, как и для рис. 4, описанным выше, при этом ${W_{new}^{\text{'}}}^3=W_{\text{*}}^3+A_5{V}_{\text{*}}^3$. Из него видно, что значительный тренд имеет лишь отношение $\left({\int }_0^h,\epsilon ,d,z\right)/\left(\overline{c_1^{\text{'}\text{'}}}{W}_e^3\right)$, что также может быть объяснено отсутствием члена с DU в определении W_{e $\frac{{\displaystyle {\int }_0^h{E}_{k}dz}}{C_{2}h{W}_{*}^2}\frac{{\displaystyle {\int }_0^h{E}_{k}dz}}{h{V}_{*}^3}\frac{{\displaystyle {\int }_0^h{E}_{k}dz}}{C_{2new}h{W}_{new}^2}\frac{{\displaystyle {\int }_0^h{E}_{k}dz}}{h{\left(B_{int}+S_{int}\right)}^{2/3}}\frac{{\displaystyle {\int }_0^h{E}_{k}dz}}{C_{T}h{W}_{*}^2}$}.

 

Рис. 5. Изменения во времени отношения интеграла диссипации в КПС и различных интегральных масштабов, по данным эксперимента CBLu7.

 

Оценки констант ${\alpha }_5,C_6,C_4$, а также величины A по данным эксперимента CBLu7 составили 0.12, 0.22, 0.005 соответственно. Поскольку $\frac{{\alpha }_1}{{\alpha }_5}\ne \frac{c_2}{C_6}$, уравнение (17) нельзя переписать в виде (8), заменяя W_* на W_new. Как видно из значения параметра C₄, характеризующего поток энергии гравитационных волн из КПС в вышележащую инверсию, при характерных значениях A ₌ 0.2, B_{s =} 10^–2, N ₌ 10^–2 этот поток (третий член справа в уравнении (17)) $C_4\frac{N^3}{B_s^{{}^3{/}_2}}{\left(\frac{A}{1+A}\right)}^2\simeq {10}^{-9}$, т. е. значительно меньше остальных слагаемых баланса ТКЭ в конвективном слое (например, A), и им можно пренебречь.

6. ЗАКЛЮЧЕНИЕ

Был проведен ряд численных экспериментов по воспроизведению свободной конвекции и конвекции со сдвигом над однородной поверхностью с помощью LES-модели ИВМ РАН. Данные расчетов подтверждают, что профили ТКЭ и диссипации ТКЭ хорошо масштабируются с помощью высоты КПС и масштаба Дирдорфа, соответствующие безразмерные константы оказались близкими к предыдущим оценкам из лабораторных экспериментов. Анализ экспериментов со сдвиговой конвекцией показал, что в этом случае масштабирование интеграла ТКЭ и интеграла диссипации по высоте КПС с помощью линейной комбинации W_* и V_* является оправданным предположением и позволяет строить интегральные модели КПС в приближении «нулевого порядка» — мгновенного скачка температуры и скорости на верхней границе. Также показано, что поток ТКЭ, выносимый из слоя вовлечения гравитационными волнами, мал и им можно пренебречь в уравнении интегрального баланса ТКЭ.

Источники финансирования. Настоящая работа поддержана Российским научным фондом (грант 17-17-01210) — в части проведения численных экспериментов с вихреразрешающей моделью и грантом РФФИ 17-05-41095 РГО_а — в части обзора состояния проблемы и верификации интегральных моделей КПС. С.С. Зилитинкевич благодарит за поддержку Российский научный фонд (гранты 15-17-20009 и 15-17-30009) и Академию Финляндии (грант ABBA № 280700 и ClimEco № 314 798/799). Коллектив авторов благодарит Суперкомпьютерный центр МГУ имени М.В. Ломоносова за предоставленные вычислительные ресурсы.

About the authors

A. V. Debolskiy

Lomonosov Moscow State University; Obukhov Institute of Atmospheric Physics, RAS 

Author for correspondence.
Email: and.debol@srcc.msu.ru
Russian Federation, Leninskiye Gory 1/4, Moscow 119234; Pyzhyevsky per., 3, Moscow 119017

V. M. Stepanenko

Lomonosov Moscow State University

Email: stepanen@srcc.msu.ru
Russian Federation, Leninskiye Gory 1/4, Moscow 119234

A. V. Glazunov

Lomonosov Moscow State University; Institute of Numerical Mathematics, RAS

Email: glazunov@inm.ras.ru
Russian Federation, Leninskiye Gory 1/4, Moscow 119234; Pyzhyevsky per., 3, Moscow 119017

S. S. Zilitinkevich

Lomonosov Moscow State University; Finnish Meteorological Institute; Lobachevsky State University of Nizhni Novgorod; Institute of Atmospheric and Earth System Research (INAR), Umiversity of Helsinki

Email: sergej.zilitinkevich@fmi.fi
Russian Federation, Leninskiye Gory 1/4, Moscow 119234; Helsinki; 23, Gagarin avenue, Nizhni Novgorod, 603600; Helsinki

References

Supplementary files