А3: переключатели схемы среднего поля — выбор по физике

01.10.2026. Ветка air/a3 от feature/air-model (0788cf7). Литература и чтение кода; код не менялся, прогонов не было. Все числа — из файлов проекта (путь рядом) или из статей (ссылка). Где число — моя арифметика по формулам кода, это помечено «оценка».

Оговорка по источникам: полные тексты большинства статей закрыты, я проверял выходные данные и DOI, а суть — по аннотациям и описаниям схем. Коэффициенты из статей, которые я привожу по памяти и не сверял с текстом, помечены «(по памяти)».

Итог

Смены в игре не требуется: все три выбора совпадают с AirCase.p (hb, local_k вкл, cbl). Матрицу сходимости А2 можно строить без surface. Одна необязательная проверка для local_k — в разделе 2.

Общее: что именно решает наш решатель

Это не LES. Решается установившаяся система (итерации Пикара) — средняя за ~час циркуляция (docs/plan/air_model.md, «Физика и границы модели»). Установившееся поле по смыслу — осреднённое по ансамблю (RANS): вихрей и термиков, которые приходят и уходят, в нём нет ни при какой клетке. Поэтому весь турбулентный перенос, включая крупные конвективные вихри, задаётся замыканием. Размер клетки (50–400 м) это не меняет.

Для сравнения: в неустойчивом слое толщиной z_i ≈ 1–3 км клетка 50–400 м — это Δx/z_i ≈ 0,02–0,4. Вихреразрешающий (нестационарный) счёт там мог бы разрешить часть термиков, «серая зона» по Honnert et al. 2011 — Δx от долей до ~2 z_i (Honnert, Masson, Couvreux 2011, doi:10.1175/JAS-D-11-061.1; Wyngaard 2004). Но наш решатель стационарный и сглаживает именно то, что LES разрешает по времени. Выводы «в серой зоне нелокальный перенос нужен не везде» (Shin & Hong 2013, doi:10.1175/JAS-D-12-0290.1) сюда прямо не переносятся.

1. closure: hb против const

Что в коде (tools/research/air3d/air.py).

• Толщина слоя h, w*, L — _bl_depth (строки 748–778): неустойчиво h = max(z_i − h_s, zi_min, 0,3 u*/f), w* = (g/θ0·H·h)^(1/3); устойчиво h = min(0,3 u*/f, 0,4 √(u* L/f)). Считается и при const, потому что от h зависят нагрев cbl и λ.
• hb — _closure (780–812): K = max(k_fa, κ w_m z (1 − z/h)²) при z < h; w_m = (u³ + 7κ (z_s/h) w³)^(1/3), z_s = min(z, 0,1 h) (неустойчиво; строки 797–802), u*/(1 + 5z/L) — устойчиво.
• const — тот же _closure, K = nu_const = 30 м²/с всюду (ветка «const» в начале). Комментарий в коде называет его «прикидкой».

Физическое допущение. K-профиль первого порядка: турбулентный обмен — диффузия с коэффициентом, форма которого задана свойствами слоя (масштаб скорости w_m, толщина h), а не локальным градиентом. Форма κ w_m z (1 − z/h)² — из Troen & Mahrt 1986 и Holtslag & Boville 1993; это главный вклад подходов «K-профиля», на которых стоят YSU и его предшественники (Hong, Noh, Dudhia 2006; улучшение по данным LES — Noh et al. 2003, Boundary-Layer Meteorol. 107, 401–427).

Что говорит литература для нашей постановки.

• Над землёй K должен стремиться к κ u* z (логарифмический слой). При u* ≈ 0,3 м/с на z = 10 м это ≈ 1,2 м²/с (оценка). Постоянная ν = 30 м²/с на этой высоте завышена в ~25 раз и гасит приземный сдвиг, то есть разгон на вершине (именно его калибруем по Askervein).
• В середине развитого слоя K — сотни м²/с: при w* ≈ 2,8 м/с, h = 2 км, z = 0,3 h формула даёт K ≈ 210 м²/с (оценка по формулам кода при H = 400 Вт/м², h = 2 км). Постоянная 30 м²/с — на порядок ниже. В docs/plan/air_model.md («три обязательные правки») уже записано, что при постоянной ν в штиле с нагревом блуждают сеточные струи.
• Сравнение K-профильных (нелокальных) и локальных схем с наблюдениями и LES — исходно Holtslag & Boville 1993: нелокальная схема воспроизводит профили слоя перемешивания лучше локальной; сравнение с LES — Noh et al. 2003 (Boundary-Layer Meteorol. 107, 401–427) (уточнить DOI) и обзор схем ПС в Holtslag et al. 2013, BAMS / Cuxart et al. 2006, GABLS1 (устойчивый слой: длинные хвосты функции устойчивости перемешивают слишком глубоко; это к local_k, раздел 2).
• Граница применимости: профиль выведен для плоской однородной местности. В горах и толщина h, и то, что считать «высотой над землёй» в сдвиговых зонах за гребнем, — допущение. У нас z = высота над местным рельефом, h = z_i − сглаженный рельеф (σ = k_smooth_m = 1500 м, строки 764–768). Как реально ведёт себя верх слоя над горами — не единая картина (Serafin et al. 2018, doi:10.3390/atmos9030102); это граница модели, записывается в «Границы».

Что уже посчитано (Моррис, docs/research/air-model-sensitivity.md). Замена hb → const: Δχ² Askervein +218 (из 85,8 в номинале); подъём у старта сдвигается на 0,5 м/с (штиль) и 0,39 м/с (3 м/с); на седловине ×1,83 → сдвиг 0,59. Номинал (hb) даёт 0,55 м/с в штиль и 1,41 м/с при 3 м/с.

Выбор: hb. Уверенность высокая. const — не конкурент по физике, а упрощение, оставленное для прикидки. Параметр nu_const не нужен и его можно убрать вместе с веткой.

Что должно получиться в числах. Ничего нового: hb — текущий номинал; ожидания те же, что в Моррисе (подъём у старта при 3 м/с 1,3–1,45 м/с, штиль 0,4–0,55 м/с; docs/archive/plan/air-model-a2pre.md на ветке air/a2-pre).

Что остаётся неизвестным / фиксируется.

• Фиксируются литературой: κ = 0,4; форма профиля и множитель 7 в w_m (TM86/HB93; по памяти); z_s = 0,1 h; f = 1,13·10⁻⁴; k_fa = 1 м²/с (Моррис: слабая ручка; свободная атмосфера 0,1–1).
• zi_min Моррис показал как «не действует» (z_i из погоды выше). Фиксировать 300 м.
• Неизвестным для волны Б closure ничего не добавляет.

2. local_k: добавка по местному сдвигу

Что в коде.

• Ядро kloc (air.py, строки 313–365), вызывается каждую итерацию из update_k (1033–1040). K* = max(K_b, l²|S|F(Ri)); 1/l = 1/(κz) + 1/λ; λ = max(lam, lam_frac·h) (строка 699), lam = 40 м, lam_frac = 0,25; F = 1/(1 + 5 Ri)² при Ri > 0, √(1 − 16 Ri) при Ri ≤ 0 (строка 353); релаксация K на k_relax = 0,5 (численная).
• Горизонтальное K_h: Смагоринский c_s = 0,25 (WRF km_opt 4) — вне этого переключателя, фиксируется.

Физическое допущение. В сдвиговых зонах, которые профиль hb «не видит» (след и слой смешения за гребнем, верх слоя, сдвиг на инверсии), турбулентность производится местным сдвигом: K = l²|S| F(Ri) (Прандтль; длина перемешивания Blackadar 1962). Функция F — это φ_m^(−2) Бизингера–Дайера при Ri ≈ z/L (там 5 и 16, Dyer 1974 — по памяти; в Businger et al. 1971 коэффициенты 4,7 и 15): устойчиво 1/(1+5Ri)², неустойчиво √(1−16Ri) — точное соответствие формулам. Соответствие Ri ≈ z/L — только в приземном слое и при Pr_t = 1; в остальной области это допущение.

Что говорит литература.

• Сочетание «нелокальный профиль в слое + местный Ri-коэффициент» — стандарт. В HB93 и YSU местный K = l² |S| F(Ri) стоит в свободной атмосфере над слоем и там, где слой устойчив (HB93, Hong et al. 2006; функции устойчивости типа Louis 1979). Наш max(K_b, K_loc) — тот же смысл в одном выражении.
• Над сложным рельефом на клетке 100–250 м выбор схемы турбулентности вторичен по сравнению с рельефом и поверхностным обменом (Goger & Dipankar 2023/24, arXiv:2311.05528); это довод не вводить тут лишней сложности, а не против local_k.
• Без местного K в следе нет производства турбулентности, и отрыв «дышит» (Аньези H = L = 300 м: 3000 итераций без сходимости, с local_k — 711; tools/research/air3d/reference.md, раздел «Турбулентная вязкость K»).

Что посчитано. Моррис: вкл → Δχ² Askervein −43; у Онгудая подъём у старта сдвигает на 0,46 м/с (штиль) и 0,26 м/с (3 м/с); ротор за гребнем — сильнейший фактор (S до 7,7). Из разведки А2 (air/a2-pre): при малом λ/h (0,031) λ почти везде равна порогу 40 м, фоновый K в слое мал, и «жёсткой» становится связь K(Ri) ↔ θ′ в полосе верха слоя перемешивания (0,75–1,1 z_i) — это предельный цикл, а не ошибка физики.

Две оговорки физика.

1. Неустойчивая ветвь внутри слоя. Внутри развитого неустойчивого слоя K_b уже задаёт конвективное перемешивание (hb — для этого и сделан). Добавка √(1−16Ri) > 1 при Ri < 0 частично считает то же второй раз. Литература (YSU, HB93) применяет локальный Ri-коэффициент в неустойчивом слое не на этом участке. Кандидат на проверку, не меняя физики заявленного переключателя: ограничить F ≤ 1 при Ri ≤ 0 (внутри слоя). Проверить как вариант в А2/волне Б: есть ли влияние на ротор и на сходимость. Я этого не считал.
2. Устойчивая ветвь — длинный хвост. 1/(1+5Ri)² не обрывается при конечном Ri; такие функции в GABLS1 перемешивают устойчивый слой глубже LES (Cuxart et al. 2006, Holtslag et al. 2013). Это важно ночью и на инверсии (верх слоя термиков — там у нас цикл из А2), для дневных расчётов второстепенно. Решение пользователя по ночи — вне А3.

Выбор: local_k вкл. Уверенность: средняя–высокая (вкл — да; конкретная форма F в неустойчивой ветви — средняя).

Что должно получиться в числах. То же, что при номинале Морриса: Askervein χ² на 43 лучше, ротор за гребнем (минимум u на 20 м) слабее на 0,30 м/с при включении. Точнее, чем Моррис, числа из этих файлов не определить.

Что остаётся для калибровки (волна Б).

• Длина перемешивания. Блэкадар: λ_∞ = 0,00027 G/f ≈ 27 м при G = 10 м/с (docs/research/air-model-sensitivity.md, таблица факторов), HB93 — 30 м. Это литературное значение для нейтрального слоя. При h = 0,3 u*/f оно даёт λ/h ≈ 0,02–0,03 (то же, что A2-pre назвал λ/h = 0,031 и «новый» набор). Значение λ/h = 0,25 из AM-09 (λ ≈ 330 м на нейтральном h Askervein) — верхняя граница, которую данные Askervein тянули, а не то, что следует из Блэкадара (docs/research/air-model-tune.md: «данные тянут выше»; Моррис: вырождение с α, cos 0,94).
• В конвективном слое (h ≈ 2 км) λ = max(40 м, λ/h·h) при λ/h = 0,25 даёт 500 м. Единой формулы λ ~ z_i в литературе для нейтрального и конвективного слоя не нашёл. Это физически неизвестный параметр, калибровать на следе за гребнем (Perdigão, Bolund), как уже записано в рекомендации Морриса.
• lam (порог): диапазон 15–150 м (HB93 30; ECMWF 150, Beljaars & Viterbo 1998). Неизвестен в этом диапазоне; действует на ротор и подъём в штиль (S ≈ 1,1–1,3).
• Фиксировать: c_s = 0,25; коэффициенты 5 и 16 (Бизингер–Дайер).

3. heat_mode: cbl против surface

Что в коде.

• surface: источник тепла Q = H/(ρc_p Δz) в первую клетку воздуха столба (строка 666).
• cbl: при H > 0 то же количество тепла раскладывается равномерно по клеткам слоя 0 < z < h над рельефом: Q = H/(ρc_p · n · Δz), n — число таких клеток (строки 667–677). Это эквивалент профиля потока H(z) = H0(1 − z/h): дивергенция линейного потока — постоянная скорость нагрева H0/(ρc_p h). Охлаждение (H < 0) — всегда в первую клетку. Уход потока на вовлечение у верха слоя не задан.
• ρc_p = 1,2·1005 (строка 38); z_i, zi_min как выше.

Физическое допущение. В среднем по времени слой перемешивания прогревается почти равномерно по высоте, потому что тепло от земли поднимают термики (нелокальный перенос), а не молекулярно-турбулентная диффузия по градиенту. Поток линейно падает с высотой от H0 у земли до ~0 (строго: −0,2 H0 из-за вовлечения) у z_i. Это классика Deardorff 1966: в середине слоя подъём тепла идёт при нулевом или даже «неправильном» знаком градиента θ.

Что говорит литература.

1. Градиентная диффузия не может перенести тепло при нейтральном θ̄. В слое перемешивания средний градиент θ ≈ 0, значит, по градиентной схеме поток ≈ 0. K-профильные схемы (TM86, HB93, YSU) добавляют «противоградиентный» член γ_c ≈ b H0/(w_s h), b ≈ 8 (TM86: 8,5; HB93: 7,8 — по памяти), именно для этого (Troen & Mahrt 1986; HB93; Hong et al. 2006). Схема EDMF добавляет к диффузии массовый поток термиков — он несёт основную часть потока в середине слоя (Siebesma, Soares, Teixeira 2007). Доли нелокального и локального переноса в серой зоне по LES разобраны в Shin & Hong 2013.
2. Как это делают мезомасштабные модели. WRF с YSU/MM5-подобными схемами кладёт поверхностный поток граничным условием в первый слой, а вверх его несёт нелокальный член схемы. То есть «весь поток в первую клетку и только градиентная диффузия» — не практика ни при какой клетке. Режим без схемы ПС (LES-режим WRF) кладёт поток в первую клетку, но там термики разрешаются нестационарно и потом осредняются. У нас осредняется сразу (см. «Общее»).
3. Оценка, что даёт surface в нашей системе (моя арифметика по формулам кода, H = 400 Вт/м² ≈ 0,33 К·м/с, h = 2 км, z = 0,3 h: K ≈ 210 м²/с, Pr_t = 0,85): чтобы одной диффузией унести 0,7 H0, нужен градиент ≈ 0,9 К/км сверхадиабатический по всей толщине слоя (≈ 1,8 К между низом и верхом). В жизни слой перемешан (Deardorff 1966; γ_c по HB93 с b ≈ 8 даёт тот же порядок — ≈ 0,7 К/км, то есть ровно столько «ложной» неустойчивости пришлось бы убрать противоградиентным членом). В решателе с полной полунеявной плавучестью эта неустойчивость превращается в разрешённое течение. Отсюда и подъём: surface завышает.
4. Что уже видно в счёте. Подъём у старта, штиль (old, окно 50 м): cbl 0,41 м/с против surface 1,97 (docs/archive/plan/air-model-a2pre.md на ветке air/a2-pre); Моррис: surface +1,48 м/с (штиль) и +0,95 м/с (3 м/с), S = 14,8 и 9,5; разведка А2 при 3 м/с: +0,25 (old) … +0,65 (new). Порядок w* для тех же условий ≈ 2,8 м/с (оценка): подъём 2,0–2,3 м/с, усреднённый за час по клетке 50 м в штиль, был бы ~0,8 w* — для осреднённого поля это неправдоподобно много. Парапланы Франции (1,47 млн наборов высоты): набор 1,03–1,41 м/с (docs/research/experimental_data.md, строка H) — это набор аппарата, не w термика, сравнивать по модулю нельзя, но порядок cbl (0,4–1,5) согласуется, surface — нет. Данных, по которым можно было бы прямо отличить режимы, у нас нет.
5. Сходимость. cbl: 111–161 итерация; surface: «сеточные термики» на 400 м, 3000 итераций, w гуляет (tools/research/air3d/reference.md, «Нагрев»; а2-pre: surface не сходится уже в областях).

Где cbl слабее. Над крутым прогретым склоном есть тонкий слой склонового течения у земли (десятки метров), за которым идёт конвективный слой. Схема cbl не отличает их: весь нагрев размазан по z_i − h_s. Литература по LES горных долин (Schmidli 2013, doi:10.1175/JAS-D-13-083.1, Serafin et al. 2018) показывает, что в горах нагрев атмосферы — сумма турбулентной дивергенции и средних течений (подъём по склону, компенсирующее опускание). Точной доли «местного» нагрева, которую надо положить в нижние клетки, по ним я не извлёк, и чисел у меня нет. Это и есть нерешённое.

Выбор: cbl. Уверенность: высокая для установившегося поля развитого слоя на склонах умеренной крутизны; средняя у крутых склонов и над скалами. surface я не рекомендую ни при какой клетке в текущем стационарном решателе: он физически неверен (пп. 1–3) и к тому же хуже сходится.

Что должно получиться в числах. Подъём у старта (max w на 200 м, 1,5 км от старта; Онгудай 12:00, 150°): 3 м/с 1,33 (old) … 1,45 (new) м/с; штиль 0,41–0,55 м/с (у new в окне 50 м решение не сошлось, число условно; a2pre). Ожидание от cbl для 1-го порядка на клетке 50–400 м: занижение на 5–25 % (docs/plan/air_model.md, сходимость по клетке).

Что остаётся неизвестным. H0 из погодной модели (surface_heating.gd, 360–430 Вт/м²); tau_cool (S = 2,1 / 1,3 на подъём; 7200 с, диабатическая часть — А1); k_smooth_m (S = 1,6 в штиль); Pr_t = 0,85 — фиксирован (Kays 1994). Если примем гибрид (вопрос ниже) — доля f_loc — новый неизвестный параметр.

4. Схема против параметров

Принцип: ошибка схемы — то, что меняется при измельчении клетки или смене порядка переноса при тех же физических параметрах; физика — то, что задаётся законами и литературой. Переключатели выше — физические решения: выбираются по модели, не по тому, что даёт меньше χ².

Предостережение, которое вытекает из Морриса: Askervein различает одну комбинацию (94 % на первом направлении SVD; cos λ/h–α = 0,94, λ/h–порядок = 0,93, λ/h–local_k = 0,86). Поэтому:

• Нельзя выбирать переключатель по χ² Askervein: Δχ² closure 218 и local_k 43 показывают, что схема «сдвигает уровень разгона», а не что она правильна.
• Нельзя считать, что 1-й порядок в игре «подгоняется» λ/h: поправка сетки — отдельный множитель (AM-09).
• Нельзя брать surface, потому что численная диффузия «размазывает» нагрев вверх и так «почти» повторяет cbl: это накладывает ошибку схемы на выбор физики и сломается при смене клетки.

5. Что остаётся физически неизвестным для волны Б (на схеме hb + local_k + cbl)

Порядок действий при калибровке: сначала структура (эта записка), затем поправка сетки, затем параметры. Волну Б надо запускать на cbl + hb + local_k.

6. Вопрос пользователю (heat_mode)

Литература однозначна для развитого слоя: нагрев распределять по толщине слоя (cbl). Неоднозначно одно: крутой склон, где у земли есть тонкое течение «вверх по склону» отдельно от конвективного слоя выше. По литературе я не смог получить, какую долю нагрева надо класть у земли. Поэтому вопрос.

Как класть тепло от земли в модель воздуха? Сейчас тепло «размазывается» по всей толщине слоя перемешивания (километр–два): так ведёт себя настоящий прогретый слой, где тепло переносят термики, и так делают погодные модели. Подъём у старта получается 0,4–0,5 м/с в штиль и ~1,4 м/с при 3 м/с (осреднённое за час поле). Варианты:

1. Как сейчас (рекомендую). Физически обоснован, сходится, не требует новых параметров. Цена: тонкое течение прямо у склона (десятки метров) модель отдельно не выделяет; его вклад в среднем поле занижен, а его вариации — в масштабе 2–3 (термики, порывы).
2. Всё тепло — в нижнюю клетку у земли. Подъём у старта в 3–4 раза больше (2 м/с в штиль), но это не физика осреднённого поля: модель решает одно установившееся поле и получает неустойчивый слой с сеточными «термиками», хуже сходится (часто не сходится вовсе). Не рекомендую.
3. Часть у земли, остальное по слою (гибрид). Даёт возможность учесть тонкое склоновое течение, но это новый параметр (доля у земли), у которого нет значения ни в литературе, ни в наших данных; подобрать его можно только по моделированию LES склона (нужна отдельная работа) либо «на вкус». Делать, только если по полётам окажется, что подъём над склоном занижен. Ответ нужен только если хотите вариант 3; иначе остаётся 1.

7. Открытые вопросы и границы

• Формула λ для конвективного слоя и форма F(Ri) в неустойчивой ветви не решены литературой; записаны как кандидаты проверки в А2/волне Б.
• Числа b (7,8/8,5), 5 и 16, l∞ = 0,00027 G/f — из памяти/вторичных документов проекта, не из прочитанных текстов статей.
• Все «что должно получиться» взяты из уже посчитанного (Моррис, А2-pre) без новых прогонов; несошедшиеся решения у new в штиле условны.
• Параметры K и h считаются в оценках только для иллюстрации порядка величин (h = 2 км — допущение для Онгудая в полдень, z_i = 3362 м над морем по a2pre).
• У Noh et al. 2003 DOI не искал: даны журнал, том, страницы.

Источники

• Troen I., Mahrt L., 1986, Boundary-Layer Meteorol. 37, 129–148, doi:10.1007/BF00122760.
• Holtslag A.A.M., Boville B.A., 1993, J. Climate 6, 1825–1842, OSTI 5644985.
• Hong S.-Y., Noh Y., Dudhia J., 2006, Mon. Wea. Rev. 134, 2318–2341, doi:10.1175/MWR3199.1.
• Noh Y., Cheon W.G., Hong S.-Y., Raasch S., 2003, Boundary-Layer Meteorol. 107, 401–427.
• Blackadar A.K., 1962, J. Geophys. Res. 67, 3095, doi:10.1029/JZ067i008p03095.
• Louis J.-F., 1979, Boundary-Layer Meteorol. 17, 187–202, doi:10.1007/BF00117978.
• Deardorff J.W., 1966, J. Atmos. Sci. 23, 503–506, doi:10.1175/1520-0469(1966)023<0503:TCGHFI>2.0.CO;2.
• Siebesma A.P., Soares P.M.M., Teixeira J., 2007, J. Atmos. Sci. 64, 1230–1248, doi:10.1175/JAS3888.1.
• Honnert R., Masson V., Couvreux F., 2011, J. Atmos. Sci. 68, 3112–3131, doi:10.1175/JAS-D-11-061.1.
• Wyngaard J.C., 2004, J. Atmos. Sci. 61, 1816–1826, doi:10.1175/1520-0469(2004)061<1816:TNMITT>2.0.CO;2.
• Shin H.H., Hong S.-Y., 2013, J. Atmos. Sci. 70, 3248–3261, doi:10.1175/JAS-D-12-0290.1.
• Cuxart J. et al., 2006, Boundary-Layer Meteorol. 118, 273–303, doi:10.1007/s10546-005-3780-1.
• Holtslag A.A.M. et al., 2013, Bull. Amer. Meteor. Soc. 94, 1691–1706, doi:10.1175/BAMS-D-11-00187.1.
• Schmidli J., 2013, J. Atmos. Sci. 70, 4041–4066, doi:10.1175/JAS-D-13-083.1.
• Serafin S. et al., 2018, Atmosphere 9, 102, doi:10.3390/atmos9030102.
• Goger B., Dipankar A., 2023/24, arXiv:2311.05528.
• Проект: tools/research/air3d/air.py, tools/research/air3d/reference.md, docs/research/air-model-sensitivity.md, docs/research/air-model-tune.md, docs/plan/air_model.md, docs/research/experimental_data.md, docs/archive/plan/air-model-a2pre.md (ветка air/a2-pre, коммит eab5c47).