Capacities and latent heat

Air contributes negligible thermal mass: \(f_a \rho_a c_a \approx 1.1 \times 10^3\) J m\(^{-3}\)K\(^{-1}\) against \(f_t \rho_w c_w \approx 1.0 \times 10^6\) J m\(^{-3}\)K\(^{-1}\) even at \(f_t = 0.24\) - about 0.1%. So

\[C_{liq} = f_t \rho_w c_w, \qquad C_{sol} = f_t \rho_i c_i, \qquad (\rho L)_{\text{eff}} = f_t \rho_w L_f\]

The volumetric latent heat scales linearly with tissue fraction. This is exact and model-independent - it is just bookkeeping of how much water is present.

Conductivity: the model that decides everything

Effective conductivity of a two-phase medium is not determined by the volume fractions alone; it depends on the microstructure. We carry three models.

k_air <- 0.026; k_w <- 0.60; k_i <- 2.25

## (1) Bruggeman symmetric effective-medium theory. Random bicontinuous mixture,
##     with a percolation threshold at f = 1/3: below that the conductive phase
##     stops forming an effective backbone and k_eff collapses.
k_emt <- function(ft, ks, ka = k_air) {
  B <- ft*(2*ks - ka) + (1-ft)*(2*ka - ks)
  (B + sqrt(B^2 + 8*ks*ka)) / 4
}
## (2) Maxwell-Eucken with the SOLID as the continuous phase and air as dispersed
##     inclusions. Appropriate when the septal network is structurally continuous.
k_mec <- function(ft, ks, ka = k_air) {
  fa <- 1 - ft
  ks*(2*ks + ka - 2*fa*(ks - ka)) / (2*ks + ka + fa*(ks - ka))
}
## (3) Volume-weighted parallel rule = Wiener upper bound. Physically unattainable
##     for a foam, but it is what much of the tissue-property literature quotes.
k_par <- function(ft, ks, ka = k_air) ft*ks + (1-ft)*ka

mixfun <- function(name) switch(name, emt = k_emt, me = k_mec, par = k_par)

At the FRC tissue fraction of 0.241 the parallel rule gives 0.164 W m\(^{-1}\)K\(^{-1}\) for unfrozen lung - essentially the 0.16-0.20 range quoted in the literature. This is strong circumstantial evidence that those published values are computed, not measured, and it is why we do not treat them as ground truth.

Which model is right for lung?

This matters enormously, so it deserves an argument rather than an assertion.

EMT is the wrong model here. Bruggeman’s symmetric theory describes a random mixture of two grain populations and has a percolation threshold at \(f = 1/3\): below that fraction the minority phase is assumed to lose its connected backbone and conduction collapses. But alveolar septa are not randomly dispersed grains. They form a fully connected, load-bearing, structurally continuous network that carries the pulmonary capillary bed. That network does not percolation-fail at 20% volume fraction - it is continuous by anatomy. Applying EMT to it imposes a disconnection that does not physically occur.

Maxwell-Eucken with the solid continuous is the appropriate baseline, since it encodes exactly that: a continuous conducting skeleton with non-conducting inclusions. It sits, as it should, between the EMT floor and the parallel ceiling.

We nonetheless carry all three, because the reader should see how much rides on it. The consequence is dramatic:

The front speed is governed by \(k_{frozen}/(f_t \rho L)\) - conduction supply divided by latent-heat demand. Under Maxwell-Eucken this ratio is almost flat across the entire physiological range: less water to freeze is almost exactly offset by worse conduction. Under EMT it falls sevenfold below \(f_t \approx 0.35\). So the two models make qualitatively different predictions: one says the ice ball barely cares how aerated the lung is, the other says aeration dominates. We will use the clinical specimen size to decide between them.

Model formulation

Identical to the nodule analysis - a 2-D axisymmetric enthalpy-method Stefan solver - but with composite properties substituted throughout. The probe metal tip is a Dirichlet boundary, the shaft proximal to it is adiabatic, tissue surrounds the tip laterally and distally, and the far field is held at 37 \(^\circ\)C.

P0 <- list(
  rho_w = 1000, c_w = 4180, k_w = 0.60,
  rho_i = 917,  c_i = 2000, k_i = 2.25,
  Lf = 334000, Tf = 0, dTmush = 0.5,
  Tbody = 37, Tprobe_min = -79, tau_cool = 0.75,   # CO2 probe
  w_perf = 0, c_b = 3600, k_ice_Tdep = FALSE,
  ft = 0.241,          # tissue volume fraction (FRC normal lung)
  mix = "me",          # mixing model: "me", "emt" or "par"
  a = 0.85e-3,         # 1.7 mm OD probe
  Ltip = 4.0e-3,       # metallic contact tip length (unpublished; swept 3-5 mm)
  dx = 0.085e-3,       # 10 cells across the probe radius
  Rmax = 8.0e-3, zmin = -5.0e-3, zmax = 12.0e-3,
  tend = 15, safety = 0.15
)

run_cryo <- function(P, snap_times = seq(1, 15, by = 1), n_rec = 60) {

  dr <- dz <- P$dx
  Nr <- round(P$Rmax / dr); Nz <- round((P$zmax - P$zmin) / dz)
  rc <- (seq_len(Nr) - 0.5) * dr
  zc <- P$zmin + (seq_len(Nz) - 0.5) * dz
  Rm <- matrix(rc, Nr, Nz); Zm <- matrix(zc, Nr, Nz, byrow = TRUE)
  Vol <- 2 * pi * Rm * dr * dz

  PROBE <- (Rm < P$a) & (Zm > 0) & (Zm < P$Ltip)
  SHAFT <- (Rm < P$a) & (Zm >= P$Ltip)
  TIS   <- !PROBE & !SHAFT

  ## ---- composite properties ------------------------------------------------
  mf <- mixfun(P$mix)
  ki_solid <- P$k_i
  kF <- mf(P$ft, ki_solid)      # frozen composite conductivity
  kU <- mf(P$ft, P$k_w)         # unfrozen composite conductivity
  Cs <- P$ft * P$rho_i * P$c_i  # frozen volumetric heat capacity
  Cl <- P$ft * P$rho_w * P$c_w  # unfrozen volumetric heat capacity
  rhoL <- P$ft * P$rho_w * P$Lf # volumetric latent heat

  Ts <- P$Tf - P$dTmush; Tl <- P$Tf
  Hs <- Cs * Ts; Hl <- Hs + rhoL

  H2T <- function(H) {
    out <- numeric(length(H))
    s <- H < Hs; m <- !s & (H <= Hl); l <- H > Hl
    out[s] <- H[s] / Cs
    out[m] <- Ts + (Tl - Ts) * (H[m] - Hs) / rhoL
    out[l] <- Tl + (H[l] - Hl) / Cl
    matrix(out, nrow(H), ncol(H))
  }
  H2f <- function(H) matrix(pmin(1, pmax(0, (H - Hs) / rhoL)), nrow(H), ncol(H))

  k_of <- function(f, Tm) {
    kFl <- if (P$k_ice_Tdep)
      mf(P$ft, pmin(4, P$k_i * 273.15 / pmax(120, 273.15 + Tm))) else kF
    km <- kFl + f * (kU - kFl)
    km[PROBE] <- 1e4; km[SHAFT] <- 0
    km
  }
  harm <- function(k1, k2) ifelse(k1 + k2 > 0, 2*k1*k2/(k1+k2), 0)

  Ar_full <- matrix(2*pi*(seq_len(Nr)*dr)*dz, Nr, Nz)
  Ar_m <- Ar_full[-Nr, , drop = FALSE]
  Az_m <- 2*pi*Rm[, -Nz, drop = FALSE]*dr
  Aout <- Ar_full[Nr, ]; Azb <- 2*pi*rc*dr

  mFr   <- PROBE[-Nr, , drop = FALSE] & TIS[-1, , drop = FALSE]
  mFz_a <- PROBE[, -Nz, drop = FALSE] & TIS[, -1, drop = FALSE]
  mFz_b <- TIS[, -Nz, drop = FALSE] & PROBE[, -1, drop = FALSE]

  ## stability from whichever phase diffuses fastest in THIS parameter set
  alpha_max <- max(kF/Cs, kU/Cl, (if (P$k_ice_Tdep) mf(P$ft, 4)/Cs else 0))
  dt <- P$safety * dr^2 / alpha_max
  nst <- ceiling(P$tend / dt); dt <- P$tend / nst

  H <- matrix(Hl + Cl * (P$Tbody - Tl), Nr, Nz)
  Tm <- H2T(H); fm <- H2f(H); km <- k_of(fm, Tm)
  E0 <- sum(H[TIS] * Vol[TIS])

  jmid <- which.min(abs(zc - P$Ltip/2))
  snaps <- vector("list", length(snap_times)); isnap <- 1L
  rec <- vector("list", 0); irec <- 0L

  for (n in seq_len(nst)) {
    tn <- n * dt
    Tm[PROBE] <- P$Tbody + (P$Tprobe_min - P$Tbody)*(1 - exp(-tn/P$tau_cool))

    kf  <- harm(km[-Nr, , drop = FALSE], km[-1, , drop = FALSE])
    Fr  <- kf * Ar_m * (Tm[-Nr, , drop=FALSE] - Tm[-1, , drop=FALSE]) / dr
    kfz <- harm(km[, -Nz, drop = FALSE], km[, -1, drop = FALSE])
    Fz  <- kfz * Az_m * (Tm[, -Nz, drop=FALSE] - Tm[, -1, drop=FALSE]) / dz

    Q <- matrix(0, Nr, Nz)
    Q[-Nr, ] <- Q[-Nr, ] - Fr; Q[-1, ] <- Q[-1, ] + Fr
    Q[, -Nz] <- Q[, -Nz] - Fz; Q[, -1] <- Q[, -1] + Fz
    Q[Nr, ] <- Q[Nr, ] + km[Nr, ]*Aout*(P$Tbody - Tm[Nr, ])/(dr/2)
    Q[, 1]  <- Q[, 1]  + km[, 1] *Azb *(P$Tbody - Tm[, 1]) /(dz/2)
    Q[, Nz] <- Q[, Nz] + km[, Nz]*Azb *(P$Tbody - Tm[, Nz])/(dz/2)
    if (P$w_perf > 0) Q <- Q + P$w_perf*P$c_b*(P$Tbody - Tm)*Vol*(fm > 0.5)

    H[TIS] <- H[TIS] + dt * Q[TIS] / Vol[TIS]
    Tm <- H2T(H); fm <- H2f(H); km <- k_of(fm, Tm)

    if (n %% max(1L, round(nst/n_rec)) == 0 || n == nst) {
      irec <- irec + 1L
      idx <- which(TIS[, jmid]); frow <- fm[idx, jmid]; rrow <- rc[idx]
      i0 <- which(frow >= 0.5)[1]
      rfront <- if (is.na(i0)) max(rrow) else if (i0 <= 1) P$a else
        rrow[i0-1] + (rrow[i0]-rrow[i0-1])*(0.5-frow[i0-1])/(frow[i0]-frow[i0-1])
      Vice <- sum((1 - fm[TIS]) * Vol[TIS])
      dTi <- pmax(0, Ts - Tm); dTi[!TIS] <- 0
      rec[[irec]] <- data.frame(
        t = tn, r_front = rfront, V_ice = Vice,
        power = sum(-Fr[mFr]) + sum(-Fz[mFz_a]) + sum(Fz[mFz_b]),
        E_removed = E0 - sum(H[TIS]*Vol[TIS]),
        E_latent = rhoL * Vice,
        E_sens_ice = sum((1 - fm) * Cs * dTi * Vol))
    }
    while (isnap <= length(snap_times) && tn >= snap_times[isnap] - 1e-12) {
      snaps[[isnap]] <- list(t = snap_times[isnap], T = Tm, f = fm)
      isnap <- isnap + 1L
    }
  }
  list(snaps = snaps, ts = bind_rows(rec), P = P, dt = dt, nst = nst,
       props = c(kF = kF, kU = kU, Cs = Cs, Cl = Cl, rhoL = rhoL),
       grid = list(rc=rc, zc=zc, Rm=Rm, Zm=Zm, Vol=Vol,
                   TIS=TIS, PROBE=PROBE, SHAFT=SHAFT, dr=dr, dz=dz))
}

ice_metrics <- function(sn, g, P) {
  ice <- g$TIS & (sn$f < 0.5)
  if (!any(ice)) return(NULL)
  rmax <- max(g$Rm[ice]) + g$dr/2
  zlo  <- min(g$Zm[ice]) - g$dz/2
  zhi  <- max(g$Zm[ice]) + g$dz/2
  Vice <- sum((1 - sn$f[g$TIS]) * g$Vol[g$TIS])
  data.frame(t = sn$t,
             D_max = 2*rmax*1e3,
             wall_ice = (rmax - P$a)*1e3,
             distal = -zlo*1e3,
             proximal = (zhi - P$Ltip)*1e3,
             L_axial = (zhi - zlo)*1e3,
             V_ice = Vice*1e9,
             tissue_mg = Vice * P$ft * P$rho_w * 1e6,
             A_xsect = pi*(rmax*1e3)^2)
}

t_run <- system.time(sim <- run_cryo(P0))
cat(sprintf("Baseline: f_t=%.3f (%s), tip=%.1f mm, %d x %d cells, dt=%.2e s, %d steps, %.0f s\n",
            P0$ft, P0$mix, P0$Ltip*1e3, length(sim$grid$rc), length(sim$grid$zc),
            sim$dt, sim$nst, t_run[["elapsed"]]))

## Baseline: f_t=0.241 (me), tip=4.0 mm, 94 x 200 cells, dt=1.15e-03 s, 13039 steps, 24 s

cat(sprintf("Composite properties: k_frozen=%.3f  k_unfrozen=%.3f W/m/K   rho*L=%.2e J/m3\n",
            sim$props["kF"], sim$props["kU"], sim$props["rhoL"]))

## Composite properties: k_frozen=0.416  k_unfrozen=0.128 W/m/K   rho*L=8.05e+07 J/m3

Shape

front_df <- function(sn, g, P, lab = NULL) {
  fc <- sn$f; fc[g$PROBE] <- 0
  cl <- contourLines(x = g$rc, y = g$zc, z = fc, levels = 0.5)
  if (!length(cl)) return(NULL)
  bind_rows(lapply(seq_along(cl), function(i) {
    d  <- data.frame(r = cl[[i]]$x, z = cl[[i]]$y)
    dd <- d[rev(seq_len(nrow(d))), , drop = FALSE]
    rbind(transmute(d,  x = -z*1e3, y =  r*1e3),
          transmute(dd, x = -z*1e3, y = -r*1e3)) |>
      mutate(grp = paste0(sn$t, "_", i))
  })) |> mutate(t = sn$t, lab = lab %||% "")
}
`%||%` <- function(a, b) if (is.null(a)) b else a

probe_layer <- function(P) {
  a <- P$a*1e3; L <- P$Ltip*1e3
  list(annotate("rect", xmin = -L-4, xmax = -L, ymin = -a, ymax = a,
                fill = "grey70", colour = "grey40", linewidth = 0.3),
       annotate("rect", xmin = -L, xmax = 0, ymin = -a, ymax = a,
                fill = "grey25", colour = "grey15", linewidth = 0.3))
}

show_t <- c(1,2,3,4,5,6,8,10,12,15)
fr <- bind_rows(lapply(sim$snaps[sapply(sim$snaps, `[[`, "t") %in% show_t],
                       front_df, g = sim$grid, P = P0))

ggplot(fr, aes(x, y, group = grp)) +
  geom_polygon(fill = "#9EC9EC", colour = "#12508F", linewidth = 0.55) +
  probe_layer(P0) +
  facet_wrap(~ factor(paste0(t, " s"), levels = paste0(show_t, " s")), ncol = 5) +
  coord_fixed(xlim = c(-8, 4), ylim = c(-3.6, 3.6)) +
  labs(title = "Frozen lung attached to a 1.7 mm cryoprobe (4 mm tip, normal aeration)",
       subtitle = "Longitudinal section. Distal tip apex at x = 0; probe extends to the left.",
       x = "Axial distance from distal tip apex (mm)", y = "Radius (mm)") +
  theme(strip.text = element_text(face = "bold"))

sn6 <- sim$snaps[[which(sapply(sim$snaps, `[[`, "t") == 6)]]
g <- sim$grid
fld <- expand.grid(r = g$rc, z = g$zc) |>
  mutate(Temp = as.vector(sn6$T), inp = as.vector(g$PROBE | g$SHAFT),
         x = -z*1e3, y = r*1e3) |> filter(!inp)
fld <- bind_rows(fld, mutate(fld, y = -y))

ggplot(fld, aes(x, y, fill = pmax(pmin(Temp, 37), -79))) +
  geom_raster() +
  scale_fill_gradientn(
    colours = c("#08306B","#2171B5","#6BAED6","#C6DBEF","#FFFFFF",
                "#FEE0D2","#FC9272","#CB181D"),
    values = scales::rescale(c(-79,-60,-40,-20,0,10,25,37)),
    name = "T (C)", limits = c(-79, 37)) +
  geom_path(data = filter(fr, t == 6), aes(x, y, group = grp),
            inherit.aes = FALSE, colour = "black", linewidth = 0.7) +
  probe_layer(P0) +
  coord_fixed(xlim = c(-8, 4.5), ylim = c(-4, 4)) +
  labs(title = "Temperature field at t = 6 s",
       subtitle = "Black line = 0 C isotherm (specimen boundary)",
       x = "Axial distance from distal tip apex (mm)", y = "Radius (mm)")

Using the clinical data to discriminate

Reported transbronchial cryobiopsy with a 1.7 mm probe in DPLD gives a median specimen diameter of 6.7 mm. That is an external, independent constraint - but applying it requires care about which dimension it refers to.

Our specimen is a capsule, markedly longer than it is wide. Pathology reports of cryobiopsy specimens conventionally quote the greatest dimension, which for such a shape is the axial length, not the transverse diameter. We therefore test both readings rather than assume one.

val_runs <- bind_rows(lapply(names(mixes), function(mm) {
  Pi <- modifyList(P0, list(mix = mm, dx = 0.17e-3))
  s <- run_cryo(Pi, snap_times = 1:15, n_rec = 6)
  bind_rows(lapply(s$snaps, ice_metrics, g = s$grid, P = Pi)) |>
    mutate(model = mixes[[mm]])
}))

## first time a monotone series reaches `target`; NA if it never does
cross_t <- function(tt, yy, target) {
  ok <- !is.na(yy); tt <- tt[ok]; yy <- yy[ok]
  if (max(yy) < target) return(NA_real_)
  i <- which(yy >= target)[1]
  if (i == 1) return(tt[1])
  tt[i-1] + (tt[i]-tt[i-1])*(target-yy[i-1])/(yy[i]-yy[i-1])
}
target <- 6.7

val_runs |> group_by(model) |>
  summarise(`t to reach 6.7 mm transverse (s)` = round(cross_t(t, D_max, target), 1),
            `t to reach 6.7 mm greatest dim (s)` = round(cross_t(t, L_axial, target), 1),
            `Transverse D at 15 s (mm)` = round(max(D_max), 1),
            `Greatest dim at 15 s (mm)` = round(max(L_axial), 1),
            .groups = "drop") |>
  kable(caption = paste("Freeze time needed to reach the reported 6.7 mm specimen size,",
                        "under each reading of 'diameter'. NA = never reached within 15 s.",
                        "f_t = 0.24."))

val_long <- val_runs |>
  select(model, t, `Transverse diameter` = D_max,
         `Greatest dimension (axial)` = L_axial) |>
  pivot_longer(c(`Transverse diameter`, `Greatest dimension (axial)`),
               names_to = "dim", values_to = "mm")

ggplot(val_long, aes(t, mm, colour = model)) +
  annotate("rect", xmin = 5, xmax = 8, ymin = -Inf, ymax = Inf,
           fill = "grey70", alpha = 0.20) +
  geom_hline(yintercept = target, colour = "black", linewidth = 0.7) +
  geom_line(linewidth = 0.9) +
  facet_wrap(~ dim) +
  scale_colour_manual(values = c("Bruggeman EMT" = "#B2182B",
                                 "Maxwell-Eucken (baseline)" = "#12508F",
                                 "Parallel (Wiener bound)" = "#F0AD4E"), name = NULL) +
  labs(title = "Model predictions against the reported 6.7 mm clinical specimen size",
       subtitle = "1.7 mm probe, 4 mm tip. Black line = reported median; grey band = typical clinical freeze window",
       x = "Freeze time (s)", y = "Specimen dimension (mm)") +
  theme(legend.position = "right")

The comparison must be read with care: a frozen ice ball and a fixed, thawed pathology specimen are not the same object, and the ice ball is an upper bound on what detaches. With the tip length now fixed at its true 4 mm, the result is sharper than before - and less comfortable:

• Read as a transverse diameter, no model reaches 6.7 mm within 15 s. The best case is the parallel bound at 5.8 mm, and the baseline Maxwell-Eucken model reaches only 5.4 mm.
• Read as the greatest dimension, the model still falls short at clinical freeze times. Maxwell-Eucken does not reach 6.7 mm until 14.4 s, and the parallel bound not until 10.4 s. At a clinically typical 5-8 s, the predicted greatest dimension is 6.0 mm - roughly 11% short of the reported median.
• Bruggeman EMT remains firmly falsified, reaching only 4.8 mm greatest dimension at 6 s and stalling at 3.4 mm transverse. Since the ice ball is an upper bound on the specimen, this is a genuine falsification.

This is a real and instructive discrepancy, and fixing the tip length is what exposed it. In the earlier version of this analysis the unknown tip length could absorb the gap - a 5 mm tip reconciled the model with the clinical figure at about 5 s. With the tip known to be 4 mm that escape route is closed, and the model under-predicts the reported specimen size at realistic freeze times under every conductivity model tested.

Three explanations are compatible with the physics, and they are not mutually exclusive:

1. The specimen is not just the ice ball. Retraction tears tissue, and the fracture need not follow the 0 \(^\circ\)C isotherm. Frozen lung is brittle and welded to surrounding parenchyma; some tissue beyond the frozen margin plausibly comes with it. The model computes the frozen zone, which is a lower bound on what tears even while being an upper bound on what freezes.
2. The tissue at the probe is denser than the whole-lung average. Wedging, atelectasis, blood and oedema all raise the local tissue fraction. Note, however, that this does not rescue the size: the ice ball is nearly independent of aeration (see the tissue-fraction sensitivity above), so a denser target yields more tissue at much the same dimensions. Density explains yield, not size.
3. The reported dimension may not be measured on a comparable object. A 6.7 mm median is quoted for fixed specimens whose relationship to the frozen mass - collapsed foam, irregular tags, greatest-dimension convention - is not specified.

Explanation 2 is the one that can be excluded on the model’s own evidence, which is a genuinely useful narrowing. The remaining candidates are mechanical and methodological rather than thermal, which is to say: the thermal model is probably not the thing that is wrong here.

Headline answer

tbl |> filter(t %in% c(1,2,3,4,5,6,8,10,12,15)) |>
  transmute(`Freeze (s)` = t,
            `Specimen D x L (mm)` = sprintf("%.1f x %.1f", D_max, L_axial),
            `Volume (mm3)` = round(V_ice,1),
            `Tissue mass (mg)` = round(tissue_mg,1),
            `Front velocity (mm/s)` =
              round(c(NA, diff(approx(ts$t, ts$r_mm, t, rule=2)$y))/c(NA, diff(t)), 2)) |>
  kable(caption = "Baseline: 1.7 mm probe, 4 mm tip, f_t = 0.24, Maxwell-Eucken, -79 C.")

For a 1.7 mm Erbe flexible cryoprobe in normally aerated lung with a 4 mm tip:

• The frozen specimen reaches roughly 3.9 mm diameter at 5 s and 4.8 mm at 10 s, as a capsule 5.7-6.4 mm long.
• The 4 mm tip sets the specimen’s length; the probe diameter and freeze time set its width. With the tip length known rather than assumed, every dimension above is a prediction rather than a family of predictions.
• The model under-predicts the reported 6.7 mm clinical specimen at realistic freeze times (see the clinical validation above). Because the ice ball is nearly aeration-independent, a denser target cannot account for the gap; the likely explanations are mechanical tearing beyond the frozen margin and the measurement convention, not the thermal model.
• Aeration barely changes the ice ball but strongly changes the yield. Going from fully inflated to densely fibrotic lung multiplies retrieved tissue mass roughly 6.6-fold at nearly constant frozen size.
• Diminishing returns are as severe as in the nodule case. Volume roughly 1.7\(\times\) from 5 to 10 s, then only 34% more from 10 to 15 s.

What would change the answer

• The mixing model. Everything above assumes a structurally continuous septal network. If the relevant microstructure were closer to a disconnected one - severe emphysema, say, where septa are actually destroyed - conduction would drop toward the EMT curve and the ice ball would be far smaller. That is a testable, disease-specific prediction: cryobiopsy in emphysematous lung should freeze markedly less tissue than in fibrotic lung at the same settings.
• Local aeration at the probe. Wedging, atelectasis, blood and oedema all raise the local tissue fraction well above the CT-derived whole-lung average, and all push toward larger specimens.
• The tip temperature is prescribed, not solved. A real Joule-Thomson probe has finite cooling power; nothing here models gas flow inside the probe. With CO\(_2\) this matters less than it might, because the achievable surface temperature is clamped near the sublimation point rather than being free to vary - the uncertainty is whether the tip reaches \(-79\,^\circ\)C, not how far past it the probe could go.
• The frozen zone is not the fracture surface. The model stops at the 0 \(^\circ\)C isotherm. What actually detaches is decided by where brittle frozen lung fails in tension during retraction, which this model does not attempt.

Limitations specific to the lung case

Beyond the limitations of the nodule study, three are new and material:

• Homogenisation below the alveolar scale is not physical. Alveoli are 200-300 \(\mu\)m across and the grid cell is 85 \(\mu\)m, so the model resolves a continuum that the real tissue does not possess at that scale. The freezing front in lung does not advance as a smooth surface; it propagates septum by septum, with air spaces freezing only at their walls. The continuum treatment is defensible for the millimetre-scale ice ball as a whole, but the computed boundary should be read as the average position of a ragged, microscopically irregular front.
• Air is treated as static. Real alveolar gas cools, contracts and can move, and water vapour redistributes. These carry negligible energy but mean the frozen foam is not quite a sealed composite.
• The specimen is a foam, and foams collapse. The frozen mass is up to 90% air by volume. On thawing and fixation it partially collapses, so a pathologist’s measured specimen diameter is not directly comparable to the frozen diameter computed here - a caveat that applies to the 6.7 mm validation above and cuts in the direction of the model being conservative.