Fix to translate saturation adjustment to buoyancy in plume

packages/class/src/fire.ts

@@ -46,8 +46,9 @@ const defaultPlumeConfig: PlumeConfig = {
 export interface Parcel {
   z: number; // Height levels [m]
   w: number; // Vertical velocity [m/s]
+  thetal: number; // Liquid-water potential temperature (conserved) [K]
   theta: number; // Potential temperature [K]
-  qt: number; // Total specific humidity [kg/kg]
+  qt: number; // Total specific humidity (conserved) [kg/kg]
   thetav: number; // Virtual potential temperature [K]
   qsat: number; // Saturation specific humidity [kg/kg]
   b: number; // Buoyancy [m/s²]
@@ -71,7 +72,7 @@ function initializeFireParcel(
 ): Parcel {
   // Start with parcel props from ambient air
   const z = background.z[0];
-  let theta = background.theta[0];
+  let thetal = background.theta[0]; // CLASS assumes RH_near-surf < 100%
   const thetavAmbient = background.thetav[0];
   let qt = background.qt[0];
   const rho = background.rho[0];
@@ -83,24 +84,25 @@ function initializeFireParcel(
   const FFire =
     ((fire.omega * fire.C * fire.spread) / fire.d) * (1 - fire.radiativeLoss);
   const FqFire = (fire.omega * fire.Cq * fire.spread) / fire.d;
-  const FvFire = FFire * (1 + 0.61 * theta * FqFire);
+  const FvFire = FFire * (1 + 0.61 * thetal * FqFire);
 
   // Use cube root as the root may be negative and js will yield NaN for a complex number result
   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
+  // Add excess temperature/humidity
   const dtheta = FFire / (rho * cp * w);
   const dqv = FqFire / (rho * w);
-  theta += dtheta;
+  thetal += dtheta;
   qt += dqv;
 
-  // Thermodynamics
-  const T = saturationAdjustment(theta, qt, p, exner);
+  // Update thetav/qsat accordingly
+  const T = saturationAdjustment(thetal, qt, p, exner);
   const qsat = qsatLiq(p, T);
   const ql = Math.max(qt - qsat, 0);
+  const theta = thetal + (Lv / cp / exner) * ql;
   const thetav = virtualTemperature(theta, qt, ql);
   const rh = ((qt - ql) / qsat) * 100;
   const Td = dewpoint(qt, p / 100);
@@ -118,6 +120,7 @@ function initializeFireParcel(
   return {
     z,
     w,
+    thetal,
     theta,
     qt,
     thetav,
@@ -156,13 +159,15 @@ export function calculatePlume(
     // Mass flux through plume
     const m = parcel.m + (parcel.e - parcel.d) * dz;
     const emz = (parcel.e / parcel.m) * dz;
-    const theta = parcel.theta - emz * (parcel.theta - bg.theta[i - 1]);
+    // FIXME: Line below assumes background state is unsaturated
+    const thetal = parcel.thetal - emz * (parcel.thetal - bg.theta[i - 1]);
     const qt = parcel.qt - emz * (parcel.qt - bg.qt[i - 1]);
 
     // Thermodynamics and buoyancy
-    const T = saturationAdjustment(theta, qt, bg.p[i], bg.exner[i]);
+    const T = saturationAdjustment(thetal, qt, bg.p[i], bg.exner[i]);
     const qsat = qsatLiq(bg.p[i], T);
     const ql = Math.max(qt - qsat, 0);
+    const theta = thetal + (Lv / cp / bg.exner[i]) * ql;
     const thetav = virtualTemperature(theta, qt, ql);
     const rh = ((qt - ql) / qsat) * 100;
     const Td = dewpoint(qt, bg.p[i] / 100);
@@ -191,6 +196,7 @@ export function calculatePlume(
     parcel = {
       z,
       w,
+      thetal,
       theta,
       qt,
       thetav,

packages/class/src/thermodynamics.ts

@@ -8,6 +8,7 @@ const Lv = 2.5e6; // Latent heat of vaporization [J/kg]
 const ep = Rd / Rv;
 
 export function virtualTemperature(t: number, qt: number, ql: number): number {
+  // Note t is theta, not thl as conserved by the plume
   return t * (1.0 - (1.0 - Rv / Rd) * qt - (Rv / Rd) * ql);
 }