diff --git a/apps/class-solid/src/components/Analysis.tsx b/apps/class-solid/src/components/Analysis.tsx index 9235c05..2729a44 100644 --- a/apps/class-solid/src/components/Analysis.tsx +++ b/apps/class-solid/src/components/Analysis.tsx @@ -1,5 +1,5 @@ import type { Config } from "@classmodel/class/config"; -import { type Parcel, calculatePlume } from "@classmodel/class/fire"; +import { calculatePlume, transposePlumeData } from "@classmodel/class/fire"; import { type ClassOutput, type OutputVariableKey, @@ -45,7 +45,7 @@ import { MdiCamera, MdiDelete, MdiImageFilterCenterFocus } from "./icons"; import { AxisBottom, AxisLeft, getNiceAxisLimits } from "./plots/Axes"; import { Chart, ChartContainer, type ChartData } from "./plots/ChartContainer"; import { Legend } from "./plots/Legend"; -import { Line, Plume, type Point } from "./plots/Line"; +import { Line, type Point } from "./plots/Line"; import { SkewTPlot, type SoundingRecord } from "./plots/skewTlogP"; import { Button } from "./ui/button"; import { Card, CardContent, CardHeader, CardTitle } from "./ui/card"; @@ -239,16 +239,24 @@ export function VerticalProfilePlot({ "Specific humidity [kg/kg]": "qt", "u-wind component [m/s]": "u", "v-wind component [m/s]": "v", + "vertical velocity [m/s]": "w", "Pressure [Pa]": "p", "Exner function [-]": "exner", "Temperature [K]": "T", "Dew point temperature [K]": "Td", "Density [kg/m³]": "rho", + "Relative humidity [%]": "rh", } as const satisfies Record; const classVariable = () => variableOptions[analysis.variable as keyof typeof variableOptions]; + const showPlume = createMemo(() => { + return ["theta", "qt", "thetav", "T", "Td", "rh", "w"].includes( + classVariable(), + ); + }); + const observations = () => flatObservations().map((o) => observationsForProfile(o, classVariable())); @@ -266,10 +274,17 @@ export function VerticalProfilePlot({ const firePlumes = () => flatExperiments().map((e, i) => { const { config, output, ...formatting } = e; - if (config.sw_fire) { + if (config.sw_fire && showPlume()) { + const plume = transposePlumeData( + calculatePlume(config, profileData()[i].data), + ); return { ...formatting, - data: calculatePlume(config, profileData()[i].data), + linestyle: "4", + data: plume.z.map((z, i) => ({ + x: plume[classVariable()][i], + y: z, + })), }; } return { ...formatting, data: [] }; @@ -288,10 +303,12 @@ export function VerticalProfilePlot({ })) as ChartData[]; const allX = () => [ + ...firePlumes().flatMap((p) => p.data.map((d) => d.x)), ...profileDataForPlot().flatMap((p) => p.data.map((d) => d.x)), ...observations().flatMap((obs) => obs.data.map((d) => d.x)), ]; const allY = () => [ + ...firePlumes().flatMap((p) => p.data.map((d) => d.y)), ...profileDataForPlot().flatMap((p) => p.data.map((d) => d.y)), ...observations().flatMap((obs) => obs.data.map((d) => d.y)), ]; @@ -322,10 +339,6 @@ export function VerticalProfilePlot({ setResetPlot(analysis.id); } - const showPlume = createMemo(() => { - return ["theta", "qt", "thetav", "T", "Td"].includes(classVariable()); - }); - return ( <>
@@ -356,10 +369,7 @@ export function VerticalProfilePlot({ {(d) => ( - keyof Parcel} - /> + )} diff --git a/apps/class-solid/src/components/plots/Line.tsx b/apps/class-solid/src/components/plots/Line.tsx index aa5aa53..8ff9c1d 100644 --- a/apps/class-solid/src/components/plots/Line.tsx +++ b/apps/class-solid/src/components/plots/Line.tsx @@ -1,4 +1,3 @@ -import type { Parcel } from "@classmodel/class/fire"; import * as d3 from "d3"; import { createSignal } from "solid-js"; import type { ChartData } from "./ChartContainer"; @@ -36,34 +35,3 @@ export function Line(d: ChartData) { ); } - -export function Plume({ - d, - variable, -}: { d: ChartData; variable: () => keyof Parcel }) { - const [chart, _updateChart] = useChartContext(); - const [hovered, setHovered] = createSignal(false); - - const l = d3.line( - (d) => chart.scaleX(d[variable()]), - (d) => chart.scaleY(d.z), - ); - - const stroke = () => (hovered() ? highlight("#ff0000") : "#ff0000"); - - return ( - setHovered(true)} - onMouseLeave={() => setHovered(false)} - fill="none" - stroke={stroke()} - stroke-dasharray={"4"} - stroke-width="2" - d={l(d.data) || ""} - class="cursor-pointer" - > - {`Fire plume for ${d.label}`} - - ); -} diff --git a/packages/class/src/config.ts b/packages/class/src/config.ts index 5f18b1a..11a945d 100644 --- a/packages/class/src/config.ts +++ b/packages/class/src/config.ts @@ -408,7 +408,7 @@ const untypedSchema = { symbol: "Lfire", type: "number", unit: "m", - default: 10000, + default: 1000, title: "Length of the fire", "ui:group": "Fire", }, @@ -416,7 +416,7 @@ const untypedSchema = { symbol: "dfire", type: "number", unit: "m", - default: 300, + default: 10, title: "Depth of the fire", "ui:group": "Fire", }, @@ -432,16 +432,24 @@ const untypedSchema = { symbol: "Cfire", type: "number", unit: "J kg⁻¹", - default: 17.781e6, + default: 18.6208e6, title: "Heat stored in fuel", "ui:group": "Fire", }, + Cq: { + symbol: "Cq,fire", + type: "number", + unit: "kg kg⁻¹", + default: 0, + title: "Moisture per kg fuel released into plume", + "ui:group": "Fire", + }, omega: { symbol: "ωfire", type: "number", unit: "kg m⁻²", - default: 7.6, - title: "Fuel mass per area", + default: 0.1, + title: "Consumed fuel mass per area", "ui:group": "Fire", }, spread: { @@ -458,7 +466,7 @@ const untypedSchema = { unit: "-", default: 0.7, title: - "Fraction of F converted into radiative heating, and not into diffused into the atmosphere", + "Fraction of F converted into radiative heating, and not diffused into the atmosphere", "ui:group": "Fire", }, }, @@ -522,6 +530,7 @@ export type FireConfig = { d: number; h0: number; C: number; + Cq: number; omega: number; spread: number; radiativeLoss: number; diff --git a/packages/class/src/fire.ts b/packages/class/src/fire.ts index 02d911d..5faa872 100644 --- a/packages/class/src/fire.ts +++ b/packages/class/src/fire.ts @@ -4,7 +4,12 @@ import type { FireConfig } from "./config.js"; import type { ClassProfile } from "./profiles.js"; -import { calcThetav, dewpoint } from "./thermodynamics.js"; +import { + dewpoint, + qsatLiq, + saturationAdjustment, + virtualTemperature, +} from "./thermodynamics.js"; // Constants const g = 9.81; // Gravitational acceleration [m/s²] @@ -15,20 +20,24 @@ const Lv = 2.5e6; // Latent heat of vaporization [J/kg] * Configuration parameters for plume model */ interface PlumeConfig { - zSl: number; // Surface layer height [m] - lambdaMix: number; // Mixing length in surface layer [m] + fac_ent: number; // factor for fractional entrainment (0.2 in Morton model) beta: number; // Fractional detrainment above surface layer + aW: number; // Factor scaling acceleration due to buoyancy + bW: number; // Factor scaling deceleration due to entrainment dz: number; // Grid spacing [m] + fac_area: number; // prescribed surface-layer area growth factor } /** * Default plume configuration */ const defaultPlumeConfig: PlumeConfig = { - zSl: 100.0, - lambdaMix: 30.0, + fac_ent: 0.8, beta: 1.0, + aW: 1.0, + bW: 0.2, dz: 1.0, + fac_area: 10, }; /** @@ -49,6 +58,7 @@ export interface Parcel { T: number; // Temperature [K] Td: number; // Dewpoint temperature [K] p: number; // Pressure [hPa] + rh: number; // Relative humidity [%] } /** @@ -57,6 +67,7 @@ export interface Parcel { function initializeFireParcel( background: ClassProfile, fire: FireConfig, + plumeConfig: PlumeConfig = defaultPlumeConfig, ): Parcel { // Start with parcel props from ambient air const z = background.z[0]; @@ -71,26 +82,38 @@ function initializeFireParcel( const area = fire.L * fire.d; const FFire = ((fire.omega * fire.C * fire.spread) / fire.d) * (1 - fire.radiativeLoss); - const FqFire = 0.0 * FFire; // Dry plume for now + const FqFire = (fire.omega * fire.Cq * fire.spread) / fire.d; + const FvFire = FFire * (1 + 0.61 * theta * FqFire); // Use cube root as the root may be negative and js will yield NaN for a complex number result - const w = Math.cbrt( - (3 * g * FFire * fire.h0) / (2 * rho * cp * thetavAmbient), - ); + const fac_w = + (3 * g * plumeConfig.aW * FvFire) / + (2 * rho * cp * thetavAmbient * (1 + plumeConfig.bW)); + const w = Math.cbrt(fac_w * fire.h0); // Add excess temperature/humidity and update thetav/qsat accordingly const dtheta = FFire / (rho * cp * w); - const dqv = FqFire / (rho * Lv * w); + const dqv = FqFire / (rho * w); theta += dtheta; qt += dqv; - const [thetav, qsat] = calcThetav(theta, qt, p, exner); + // Thermodynamics + const T = saturationAdjustment(theta, qt, p, exner); + const qsat = qsatLiq(p, T); + const ql = Math.max(qt - qsat, 0); + const thetav = virtualTemperature(theta, qt, ql); + const rh = ((qt - ql) / qsat) * 100; + const Td = dewpoint(qt, p / 100); // Calculate parcel buoyancy - const b = (g / background.thetav[0]) * (thetav - thetavAmbient); + const b = (g / thetavAmbient) * (thetav - thetavAmbient); + + // Calculate initial entrainment/detrainment + const m = rho * area * w; + let e = ((rho * area) / (2 * w)) * b; // Entrainment assuming constant area with height in surface layer + e = e + (rho * w * area * (1 + plumeConfig.fac_area)) / fire.h0; // Additional, prescribed plume growth over surface layer + const d = 0; - const T = background.exner[0] * theta; - const Td = dewpoint(qt, p / 100); // Store parcel props return { z, @@ -100,13 +123,14 @@ function initializeFireParcel( thetav, qsat, b, - area: fire.L * fire.d, - m: rho * area * w, - e: ((rho * area) / (2 * w)) * b, - d: 0, + area, + m, + e, + d, T, Td, p: background.p[0] / 100, + rh, }; } @@ -122,10 +146,9 @@ export function calculatePlume( let parcel = initializeFireParcel(bg, fire); const plume: Parcel[] = [parcel]; - const detrainmentRate0 = plumeConfig.lambdaMix ** 0.5 / parcel.area ** 0.5; - let crossedSl = false; - let epsi = 0; - let delt = 0; + // Constant fractional entrainment and detrainment with height above surface layer + const epsi = plumeConfig.fac_ent / Math.sqrt(parcel.area); + const delt = epsi / plumeConfig.beta; for (let i = 1; i < bg.z.length; i++) { const z = bg.z[i]; @@ -136,48 +159,33 @@ export function calculatePlume( const theta = parcel.theta - emz * (parcel.theta - bg.theta[i - 1]); const qt = parcel.qt - emz * (parcel.qt - bg.qt[i - 1]); - // Calculate virtual potential temperature and buoyancy - const [thetav, qsat] = calcThetav(theta, qt, bg.p[i], bg.exner[i]); + // Thermodynamics and buoyancy + const T = saturationAdjustment(theta, qt, bg.p[i], bg.exner[i]); + const qsat = qsatLiq(bg.p[i], T); + const ql = Math.max(qt - qsat, 0); + const thetav = virtualTemperature(theta, qt, ql); + const rh = ((qt - ql) / qsat) * 100; + const Td = dewpoint(qt, bg.p[i] / 100); const b = (g / bg.thetav[i]) * (thetav - bg.thetav[i]); // Solve vertical velocity equation - const aW = 1; - const bW = 0; - const w = - parcel.w + - ((-bW * parcel.e * parcel.w + - aW * parcel.area * bg.rho[i - 1] * parcel.b) / - parcel.m) * + const w2 = + parcel.w * parcel.w + + 2 * + (plumeConfig.aW * b - plumeConfig.bW * epsi * parcel.w * parcel.w) * dz; - - // Calculate entrainment and detrainment - let e: number; - let d: number; - - if (z < plumeConfig.zSl) { - // Surface layer formulation - e = ((parcel.area * bg.rho[i - 1]) / (2 * parcel.w)) * parcel.b; - d = - parcel.area * - bg.rho[i - 1] * - detrainmentRate0 * - ((z ** 0.5 * (w - parcel.w)) / dz + parcel.w / (2 * z ** 0.5)); + let w: number; + if (w2 < 0) { + w = 0; } else { - // Above surface layer - if (!crossedSl) { - epsi = parcel.e / parcel.m; - delt = epsi / plumeConfig.beta; - crossedSl = true; - } - - e = epsi * m; - d = delt * m; + w = Math.sqrt(w2); } - const area = m / (bg.rho[i] * w); + // Calculate entrainment and detrainment + const e = epsi * m; + const d = delt * m; - const T = bg.exner[i] * theta; - const Td = dewpoint(qt, bg.p[i] / 100); + const area = m / (bg.rho[i] * w); // Update parcel parcel = { @@ -195,9 +203,10 @@ export function calculatePlume( T, Td, p: bg.p[i] / 100, + rh, }; - if (w < 0 || area <= 0) { + if (w <= 0 || area <= 0) { break; } diff --git a/packages/class/src/profiles.ts b/packages/class/src/profiles.ts index 63410b0..06b8de3 100644 --- a/packages/class/src/profiles.ts +++ b/packages/class/src/profiles.ts @@ -2,7 +2,12 @@ import type { MixedLayerConfig, NoWindConfig, WindConfig } from "./config.js"; import type { ClassOutputAtSingleTime } from "./output.js"; -import { dewpoint, virtualTemperature } from "./thermodynamics.js"; +import { + dewpoint, + qsatLiq, + saturationAdjustment, + virtualTemperature, +} from "./thermodynamics.js"; const CONSTANTS = { g: 9.81, // Gravity [m/s²] @@ -20,11 +25,13 @@ export interface ClassProfile { qt: number[]; // Total specific humidity [kg/kg] u: number[]; // U-component of wind [m/s] v: number[]; // V-component of wind [m/s] + w: number[]; // W-component of wind [m/s] p: number[]; // Pressure [Pa] exner: number[]; // Exner function [-] T: number[]; // Temperature [K] Td: number[]; // Dew point temperature [K] rho: number[]; // Density [kg/m³] + rh: number[]; // Relative humidity [%] } export const NoProfile: ClassProfile = { @@ -34,11 +41,13 @@ export const NoProfile: ClassProfile = { qt: [], u: [], v: [], + w: [], p: [], exner: [], T: [], Td: [], rho: [], + rh: [], }; /** @@ -51,7 +60,7 @@ export function generateProfiles( ): ClassProfile { const { Rd, cp, g } = CONSTANTS; const { h, theta, qt, u, v, dtheta, dqt, du, dv } = output; - const { z_theta, z_qt, gamma_theta, gamma_qt } = config; + const { z_theta, z_qt, gamma_theta, gamma_qt, divU } = config; const { p0 } = config; // Determine top of profile based on the lowest z value across all variables @@ -68,18 +77,23 @@ export function generateProfiles( const thetah = piecewiseProfile(zh, h, theta, dtheta, z_theta, gamma_theta); const qth = piecewiseProfile(zh, h, qt, dqt, z_qt, gamma_qt); - // Calculate virtual temperature, asssume base state is dry, so no saturation adjustment - const thetav = thetaProf.map((t, i) => - virtualTemperature(t, qtProfile[i], 0), - ); + // Calculate virtual temperature on half levels (no saturation) for pressure calc. const thetavh = thetah.map((t, i) => virtualTemperature(t, qth[i], 0)); - // Pressure and other thermodynamic variables + // Pressure and other thermodynamic variables, incl. saturation adjustment (taking theta = theta_l) const p = calculatePressureProfile(zh, p0, Rd, cp, g, thetavh, dz); const exner = p.map((pressure) => (pressure / p0) ** (Rd / cp)); - const T = exner.map((ex, i) => ex * thetaProf[i]); + const T = thetaProf.map((t, i) => + saturationAdjustment(t, qtProfile[i], p[i], exner[i]), + ); + const qsat = T.map((t, i) => qsatLiq(p[i], t)); + const ql = qtProfile.map((q, i) => Math.max(q - qsat[i], 0)); + const thetav = thetaProf.map((t, i) => + virtualTemperature(t, qtProfile[i], ql[i]), + ); const Td = p.map((p, i) => dewpoint(qtProfile[i], p / 100)); const rho = p.map((pressure, i) => pressure / (Rd * exner[i] * thetav[i])); + const rh = qtProfile.map((q, i) => ((q - ql[i]) / qsat[i]) * 100); // Include wind let uProfile: number[]; @@ -93,18 +107,23 @@ export function generateProfiles( vProfile = new Array(z.length).fill(999); } + // Vertical velocity from constant divergence + const w = z.map((zi) => -divU * zi); + return { z, theta: thetaProf, qt: qtProfile, u: uProfile, v: vProfile, + w, thetav, p, exner, T, Td, rho, + rh, }; } diff --git a/packages/class/src/thermodynamics.ts b/packages/class/src/thermodynamics.ts index 9d79659..2413dae 100644 --- a/packages/class/src/thermodynamics.ts +++ b/packages/class/src/thermodynamics.ts @@ -29,18 +29,17 @@ export function dqsatdTLiq(p: number, t: number): number { ); } -export function calcThetav( +export function saturationAdjustment( thl: number, qt: number, p: number, exner: number, -): [number, number] { +): number { const tl = exner * thl; let qsat = qsatLiq(p, tl); - let ql = 0.0; if (qt <= qsat) { - return [virtualTemperature(thl, qt, 0.0), qsat]; + return tl; } // Newton-Raphson iteration @@ -57,9 +56,7 @@ export function calcThetav( tnr -= f / f_prime; iter++; } - - ql = qt - qsat; - return [virtualTemperature(tnr / exner, qt, ql), qsat]; + return tnr; } /** diff --git a/packages/class/src/validate.test.ts b/packages/class/src/validate.test.ts index 98348eb..40fcde2 100644 --- a/packages/class/src/validate.test.ts +++ b/packages/class/src/validate.test.ts @@ -77,11 +77,12 @@ describe("parse", () => { gamma_v: [0], z_v: [5000], ustar: 0.3, - L: 10000, - d: 300, + L: 1000, + d: 10, h0: 20, - C: 17781000, - omega: 7.6, + C: 18.6208e6, + Cq: 0, + omega: 0.1, spread: 1.5, radiativeLoss: 0.7, name: "",