12. Boundary and Initial Conditions

This chapter specifies how Pyrolysis.jl closes the coupled energy, species, and gas-pressure equations at the two physical surfaces of the 1D column ($z=0$ substrate / :bottom, $z=L$ exposed surface / :top) and how it seeds the initial state. It derives the surface energy balance and the Newton iteration that solves for the true surface temperature $T_s$ in a cell-centered finite-volume scheme; documents every concrete thermal, mass, pressure, and mesh-motion boundary condition with its exact flux law, temperature derivative, and sign convention as implemented; explains the type-stable BoundaryConditionSet container and the tag$\to$id bridge; and treats the initial-condition closures and their auto-defaults. Throughout, the equations are checked term-by-term against the source so that the math and the code agree exactly — including the boundary-face sign flips, the Kirchhoff emissivity coupling, the neglected $d\varepsilon/dT$ term, and the strict mutual exclusion of the surface transpiration term with the volumetric advection source $S_{\text{conv}}$.

This chapter conforms to the unified nomenclature of §2, Nomenclature, and uses the three binding sign conventions: positive flux is transport in the $+z$ direction (bottom$\to$top); the discrete divergence $(F_R-F_L)A/V$ is positive when the quantity leaves a cell; and for heats of reaction $h>0$ is endothermic with $Q_{\text{rxn}}=-\sum_r h_r r_r$. In this chapter the index $i$ denotes the cell index and $j$ the component index, consistent with the discretization chapters.

12.1 Overview and design philosophy

12.1.1 What a boundary condition must supply

The continuum model of §3, Governing Equations, is a system of conservation laws for the condensed-phase temperature $T(z,t)$ and the per-component mass concentrations $\xi_j(z,t)$ (kg·m⁻³), with the gas-phase pressure $P$ as an algebraic closure in Darcy mode. In flux (divergence) form the spatial operators are

\[\rho c_p^{\text{eff}}\,\frac{\partial T}{\partial t} = -\nabla\!\cdot\mathbf{q}_{\text{cond}} + S_{\text{conv}} + S_{\text{gen}} + Q_{\text{rad}} + Q_{\text{rxn}}, \qquad \frac{\partial \xi_j}{\partial t} = -\nabla\!\cdot\mathbf{J}_j + S_{j,\text{rxn}} .\]

A cell-centered finite-volume discretization (see §13, Finite-Volume Discretization) integrates each equation over a cell $i$ and replaces the divergence by a difference of face fluxes,

\[(\nabla\!\cdot \mathbf{F})_i = \frac{(F_{f_R}-F_{f_L})\,A}{V_i}.\]

For an interior face the flux $F_f$ is a constitutive law evaluated from the two adjacent cells (Fourier, Fick, Darcy–Fick). For a boundary face there is no neighboring cell, so the flux must instead be supplied by a boundary condition. The boundary-condition subsystem therefore has exactly one job: for each boundary face, given the interior state and the time $t$, produce the heat flux $q_{\text{BC}}$, the species fluxes $J_j$, and (in Darcy mode) the gas velocity $v_g$ that enter the divergence operator with the correct sign.

There is one subtlety unique to the thermal boundary: because temperature is stored at cell centers and radiative and convective fluxes depend nonlinearly on the surface temperature $T_s$ (which is not a state variable), the code does not simply evaluate $q_{\text{BC}}$ at the adjacent cell-center temperature. Doing so produces an $\sim$20% energy imbalance for $T^4$ radiation. Instead it solves a small surface energy balance for the true $T_s$ (§12.2). The species and pressure boundaries have no such nonlinearity and are evaluated directly.

12.1.2 Type hierarchy and dispatch

Boundary conditions are organized by physical kind, each rooted in an abstract supertype declared in the Materials module:

• AbstractThermalBC — energy/heat boundary conditions;
• AbstractMassBC — species and pressure-driven (Darcy) boundary conditions;
• AbstractRadiationBC — P1 radiative-transfer boundary conditions (in the Experimental module; see §8, Thermal Radiation).

Concrete types live in the BoundaryConditions module (src/BoundaryConditions/boundary_conditions.jl), flux evaluation lives in src/Discretization/boundary_fluxes.jl, and the surface-temperature Newton solve lives in src/Discretization/boundary_temperature.jl. The whole set is stored in a type-stable, tag-keyed BoundaryConditionSet (§12.7) so that the residual specializes statically with zero dispatch overhead.

Implementation note. Concrete BC types and the BoundaryConditionSet container are defined in src/BoundaryConditions/boundary_conditions.jl; flux laws and their derivatives in src/Discretization/boundary_fluxes.jl and src/Discretization/boundary_temperature.jl; BC application into the flux arrays in src/Discretization/finite_volume.jl (apply_thermal_bc_flux!, apply_mass_bc_flux!); pressure BCs in src/Discretization/darcy_velocity.jl (apply_pressure_bc_to_darcy_velocity!).

12.1.3 The boundary-face sign convention

Two boundaries exist in 1D, with the binding id/tag/position mapping (src/Geometry/topology.jl:77, src/Geometry/mesh_generation.jl:127–128):

tag_to_boundary_id(:top)=1, tag_to_boundary_id(:bottom)=2, boundary_id_to_tag(1)=:top.

Because the discrete divergence is $(F_{f_R}-F_{f_L})A/V_i$ with the convention that a positive face flux points in $+z$, the BC application must orient the computed into-the-domain flux correctly for each boundary. A flux $q_{\text{cond}}$ that carries heat into the slab is in the $+z$ direction at the bottom face (left face of cell 1) but in the $-z$ direction at the top face (right face of cell $n$). Hence the implementation stores

\[q[f_{\text{bottom}}] = +q_{\text{conducted}}, \qquad q[f_{\text{top}}] = -q_{\text{conducted}},\]

and identically for mass fluxes (finite_volume.jl:550–554, 590–596, 712–716). This is not a bug; it is the consistent realization of the signed divergence. We return to it in §12.2.5 and §12.4.4.

12.2 The surface energy balance and the Newton solve

The surface temperature $T_s$ is the keystone of the thermal boundary. It is the temperature of the exposed face itself, distinct from the adjacent cell-center temperature $T_{\text{cell}}$, and it is determined by an instantaneous balance between the externally applied flux and the heat conducted into the slab.

12.2.1 Continuous statement

Consider the exposed face. An energy balance on an infinitesimally thin control volume at the surface (zero thermal mass) equates the net flux applied by the environment to the flux conducted into the material:

\[q_{\text{BC}}(T_s) \;=\; q_{\text{cond}}(T_s, T_{\text{cell}}),\]

with, for the most general thermal BC,

\[q_{\text{BC}}(T_s) = q_{\text{ext}}\;-\;\varepsilon_{\text{eff}}\sigma T_s^4 \;-\;h_{\text{conv}}\,(T_s-T_\infty)\]

(absorbed external flux, re-radiation loss, convective loss; all signed positive into the domain) and the conducted flux written as a one-sided Fourier difference from the surface to the first cell center,

\[q_{\text{cond}}(T_s,T_{\text{cell}}) = k_{\text{eff}}\,\frac{T_s - T_{\text{cell}}}{\Delta z_{1/2}},\]

where $\Delta z_{1/2}$ is the distance from the boundary face to the first cell center and $k_{\text{eff}}$ is the effective conductivity of the boundary cell (§5, Effective Properties & Mixing Rules). Both sides are positive when heat flows into the material. This is Formulation A, the non-conservative convention in which gas sensible enthalpy is transported by the volumetric source $S_{\text{conv}}$ (§7) rather than by a surface blowing term; the gas-storage term is omitted from $\rho c_p^{\text{eff}}$ (see §3, §5, §7). The optional transpiration extension (§12.9) adds a surface enthalpy sink and is off by default.

The conductive distance is exactly half the boundary cell's thickness, computed from state-dependent geometry (AD- and ALE-clean) as

\[\Delta z_{1/2} = \frac{1}{2}\,\frac{V_{\text{cell}}}{A_{\text{face}}} \qquad\text{(`boundary\_temperature.jl:563`)}.\]

Using the half-width (rather than the full cell width) places the surface at the physical face and the datum $T_{\text{cell}}$ at the cell center, which is the correct one-sided gradient for a cell-centered scheme.

12.2.2 Residual and Newton iteration

Define the surface residual

\[f(T_s) \;=\; q_{\text{BC}}(T_s)\;-\;q_{\text{cond}}(T_s)\;-\;q_{\text{transpire}}(T_s),\]

where the transpiration term $q_{\text{transpire}}=\sum_g \dot m_g\,\Delta h_g(T_s,T_\infty)$ is zero in the default configuration (§12.9). The root $f(T_s)=0$ is found by damped Newton iteration

\[T_s^{(k+1)} = T_s^{(k)} - \omega\,\frac{f(T_s^{(k)})}{f'(T_s^{(k)})}, \qquad f'(T_s) = \frac{dq_{\text{BC}}}{dT_s} - \frac{k_{\text{eff}}}{\Delta z_{1/2}} - \frac{dq_{\text{transpire}}}{dT_s},\]

with under-relaxation $\omega=1$ by default (boundary_temperature.jl:168–205). The initial guess is a warm start from the previously solved surface temperature ($T_s^{(0)} =$ state.T_surface, the persistent per-boundary BoundaryState; the very first call starts from the 300 K mesh-construction seed) — on smooth time-stepping this saves Newton iterations relative to restarting from $T_{\text{cell}}$. (The non-mutating convenience wrapper solve_boundary_temperature accepts a T_guess keyword defaulting to $T_{\text{cell}}$.) The conduction derivative $dq_{\text{cond}}/dT_s = k_{\text{eff}}/\Delta z_{1/2}$ is exact and always present; the BC derivatives are given per-type in §12.2.4.

Convergence. The residual is compared against a relative tolerance scaled by the dominant physical flux,

\[|f(T_s)| < \texttt{tol}\cdot \max\!\bigl(|q_{\text{BC}}|,\,|q_{\text{cond}}|,\,|q_{\text{transpire}}|,\,1\bigr), \qquad \texttt{tol}=10^{-10},\]

with a maximum of 50 iterations (boundary_temperature.jl:184–193). Scaling by the dominant flux makes the criterion meaningful whether the surface flux is $\sim 10^4$ W·m⁻² (cone) or near zero (near-adiabatic). Nothing aborts on non-convergence: on max-iterations or a vanishing Jacobian ($|f'|<10^{-30}$ — e.g. an adiabatic boundary with no temperature dependence anywhere) the loop emits a throttled @warn and falls through to the final flux evaluation at the unconverged iterate; the ODE integrator's step-error control is the backstop that rejects a step built on a bad surface state.

Clamping. Each Newton iterate is clamped to the physical window $[200,\,3000]$ K (boundary_temperature.jl:202–204) to prevent unphysical excursions during the transient. The clamp constructs the clamped value with zero ForwardDiff partials (ForwardDiff convert semantics) — derivatives do not ride along through a clamped iterate. This self-heals in the normal path because a subsequent un-clamped Newton update restores the implicit- function derivatives; only when the final iterate is pinned at 200/3000 K (root outside the window) are the returned $T_s$/$q$ sensitivities zero — which is then the mathematically correct clamp derivative, and the solver warns (throttled) with the AD consequence spelled out.

For the physical BCs, convergence is fast and robust: convection contributes a constant negative slope $-h_{\text{conv}}$ (one Newton step is exact), and radiation contributes a monotone cubic slope $-4F\varepsilon_{\text{eff}}\sigma T_s^3$ that Newton handles in 2–5 iterations. Combined BCs sum monotone contributions, so the residual is monotone decreasing in $T_s$ and the root is unique.

12.2.3 Automatic differentiation through the solve

The solver is written to carry whatever scalar type its inputs carry. It promotes to the widest type among $(T_{\text{cell}}, k_{\text{eff}}, \Delta z_{1/2}, T_\infty, t, \varepsilon_{\text{eff}}, \dot m_g)$ (boundary_temperature.jl:154–162); under sensitivity analysis (§17, Sensitivity Analysis) this is a ForwardDiff.Dual. The Newton convergence test uses ForwardDiff.value(f) (the primal residual), so the iteration count is identical with or without AD, while the returned tuple $(T_s,\,q_{\text{applied}},\, q_{\text{conducted}})$ is fully Dual-typed and carries sensitivities into the Jacobian (boundary_temperature.jl:207–227). One guard is needed: on the very first iteration under AD the perturbed residual may pass the primal tolerance while the dual part has not yet propagated, so convergence is suppressed for dual_mode && i==1 (boundary_temperature.jl:190). The diagnostic BoundaryState stores plain Float64 (it is bookkeeping, never differentiated).

12.2.4 Per-type flux and derivative

The Newton kernel calls evaluate_bc_flux_and_derivative(bc, T_s, t, ε_eff, radiation_model), which returns $(q_{\text{BC}}, dq_{\text{BC}}/dT_s)$ for each concrete thermal BC. The exact pairs (verified against boundary_temperature.jl:307–429) are collected below; the $\varepsilon_{\text{eff}}$ argument is the composition-weighted effective surface emissivity from the property cache (§12.2.6).

Each BC is derived and motivated in §12.3.

Neglected emissivity slope. For a temperature-dependent emissivity $\varepsilon=\varepsilon(T_s)$ the exact radiative derivative would be

\[\frac{dq_{\text{BC}}}{dT_s} = -4F\varepsilon(T_s)\sigma T_s^3 \;+\; F\,\frac{d\varepsilon}{dT_s}\,\sigma\,(T_\infty^4 - T_s^4),\]

but the code keeps only the first term (boundary_temperature.jl:373–377). The omitted term is part of the Newton Jacobian of the iteration, not the residual itself — the residual still uses the true $\varepsilon(T_s)$ — so the converged $T_s$ is unbiased; only the convergence rate degrades from quadratic to super-linear when $d\varepsilon/dT$ is large. This is documented as a known limitation (§12.13).

12.2.5 Application into the flux array and the sign flip

After the Newton solve, _apply_thermal_bc_static! writes the conducted flux $q_{\text{conducted}}$ (not the applied flux) into the energy face-flux array, with the boundary-orientation sign

\[q[f] = \begin{cases} +q_{\text{conducted}}, & \text{bottom (}z=0\text{)},\\[2pt] -q_{\text{conducted}}, & \text{top (}z=L\text{)}, \end{cases}\]

(finite_volume.jl:590–596). Writing the conducted flux (rather than the applied flux) is what makes the scheme energy-conserving: the divergence operator then sees the flux that actually enters the first interior cell, and the surface energy balance $q_{\text{BC}}=q_{\text{cond}}$ guarantees that the energy delivered at the surface equals the energy conducted inward. The diagnostic check check_boundary_energy_balance returns the relative imbalance $|q_{\text{applied}}-q_{\text{conducted}}|/|q_{\text{applied}}|$ (the absolute difference when $|q_{\text{applied}}|\le10^{-10}$) per boundary, which is $\approx$ machine precision when Newton has converged (boundary_temperature.jl).

Recall (§7, §13) that the interior conductive fluxes are explicitly zeroed at boundary faces in compute_conductive_fluxes! before BC application, so the overwrite here is the sole source of the boundary face flux; there is no double-counting.

12.2.6 Effective surface emissivity (composition-weighted, Kirchhoff)

When a RadiativeBC or RadiativeFluxBC is constructed with $\varepsilon=\texttt{nothing}$ (resp. $\alpha=\texttt{nothing}$), the surface radiative property is taken from the boundary cell's effective emissivity $\varepsilon_{\text{eff}}$, computed in the property cache as the volume-weighted mean over the solid/liquid components,

\[\varepsilon_{\text{eff}} = \sum_j v_j\,\varepsilon_j,\]

clamped to $[0,1]$ (see §5.8). Two physical principles are encoded here. First, Kirchhoff's law $\alpha=\varepsilon$ at thermal equilibrium lets the same quantity serve as both emissivity (re-radiation) and absorptivity (absorbed incident flux). Second, surface emissivity is a geometric property of the exposed mixture, so a volume-weighted (not mass-weighted) average is appropriate and it automatically tracks composition changes during pyrolysis (e.g. polymer $\to$ char as $v_j$ shifts). The boundary solver reads this as cache.ε[cell_id] (boundary_temperature.jl).

12.2.7 BoundaryState storage

The per-boundary diagnostic record (src/Geometry/mesh.jl:29–39) holds the solved surface temperature, the applied flux, the conducted flux, the transpiration flux, the inputs $(T_{\text{cell}}, k_{\text{eff}}, \Delta z_{1/2})$, and convergence metadata (converged, iterations). All fields are Float64; they are diagnostics and are never differentiated. The live (possibly Dual-typed) quantities are returned from the solve and flow into the residual.

12.3 Thermal boundary conditions

This section derives each concrete thermal BC, its physical meaning, and how it plugs into §12.2. All flux laws are signed positive into the domain.

12.3.1 AdiabaticBC

Perfect insulation: $q_{\text{BC}}=0$, $dq_{\text{BC}}/dT_s=0$. The surface energy balance reduces to $q_{\text{cond}}=0$, hence $T_s=T_{\text{cell}}$ — the surface floats to the adjacent cell temperature with zero conduction across the face. Typical for the substrate face in a cone-calorimeter idealization with a well-insulated back. (boundary_fluxes.jl:43, boundary_temperature.jl:307–314.)

12.3.2 HeatFluxBC

A prescribed flux $q_{\text{BC}}=q(t)$, constant (HeatFluxBC{<:Real}) or a function of time (HeatFluxBC{<:Function}), independent of $T_s$ so $dq_{\text{BC}}/dT_s=0$. The Newton solve still runs to find the $T_s$ that balances this against conduction. Sign: positive $q$ heats the material. (boundary_fluxes.jl:46–47, boundary_temperature.jl:319–335.)

12.3.3 ConvectiveBC (Newton's law of cooling)

\[q_{\text{BC}} = h_{\text{conv}}(t)\,\bigl(T_\infty(t) - T_s\bigr), \qquad \frac{dq_{\text{BC}}}{dT_s} = -h_{\text{conv}}(t).\]

Both $h_{\text{conv}}$ and $T_\infty$ may be constants or functions of $t$ (evaluated through evaluate_param). The flux is positive (heating) when $T_\infty>T_s$ and negative (cooling) otherwise, so for a hot surface in a cool ambient this is a convective loss. Because the derivative is the constant $-h_{\text{conv}}$, the contribution to the Newton residual is linear and exactly resolved in a single step. (boundary_fluxes.jl:50–54, boundary_temperature.jl:340–352.)

Blowing (transpiration) correction — opt-in. With use_blowing_correction = true (a solve/to_problem_def switch, default off) every ConvectiveBC in the surface balance is rescaled by the classical Couette-film (Spalding) factor

\[h \;\to\; h\,\frac{B}{e^{B}-1}, \qquad B \;=\; \frac{\dot m''\,c_{p,g}}{h}, \qquad \dot m''\,c_{p,g} \;=\; \sum_{g}\dot m_g\,c_{p,g}(T_s),\]

where $\dot m_g$ is the signed outward gas mass flux at the boundary face, read from the same J_gas[:, f] the mass-flux pass populates (film-law diffusion plus, in Darcy–Fick mode, the pressure-BC venting flux). This is the momentum-side companion of the transpiration enthalpy term: the gas efflux thickens the boundary layer and suppresses convective heat exchange. Unlike the (rejected) surface transpiration term it modifies only the convective coefficient, not the gas enthalpy accounting, so it composes exactly with the $S_{\text{conv}}$-only convention (§7.6) — no double-count is possible. Limits: $B = 0$ recovers the bare $h$ exactly; $B = 4$ (peak-PMMA outgassing against canonical cone $h$) gives $h_{\text{eff}}/h = 0.075$, i.e. convective cooling essentially extinguished; $B < 0$ (net inflow) is the smooth suction enhancement. The factor is evaluated by the numerically stable blowing_factor (Taylor branch at $|B|<10^{-5}$, overflow-safe $B\,e^{-B}/(1-e^{-B})$ form for large $B$), and the blowing conductance is rebuilt at every Newton iterate with $c_{p,g}(T_s)$. The Newton derivative uses the frozen-$B$ form $-h\,f(B)$ — first-order consistent, the same convention as the transpiration derivative (§12.9), so convergence is marginally sub-quadratic and the converged root is exact. Each convective member of a CombinedThermalBC forms its own $B$ against its own $h$; all other BC types are untouched, and the back face is unaffected automatically ($\dot m''\approx 0$ there). The converged conductance is stored as BoundaryState.blowing_mcp and consumed by the flux-decomposition diagnostics. The structured Jacobian's BCThermalOperator and BCConvectionOperator rebuild $\dot m''$ from the stencil state via the shared boundary-face kernel, so the $\partial(h\chi)/\partial(T_b,\xi_b)$ chain is exact (support-complete strict verification passes with blowing enabled in Fickian, Darcy–Fick, and ALE modes). Enabling blowing alone on a case calibrated without it will raise $T_s$ and MLR; introduce it together with flame heat-flux feedback before any recalibration pass (the two partially cancel — see the physics roadmap, docs/Physics_Audit_2026-06-09.md). (boundary_temperature.jl blowing_factor/_blowing_conductance, adapter.jl use_blowing_correction.)

12.3.4 RadiativeBC (gray surface exchange with surroundings)

\[q_{\text{BC}} = F\,\varepsilon\,\sigma\,\bigl(T_\infty^4 - T_s^4\bigr), \qquad \frac{dq_{\text{BC}}}{dT_s} = -4\,F\,\varepsilon\,\sigma\,T_s^3,\]

with the Stefan–Boltzmann constant $\sigma=5.670374419\times10^{-8}$ W·m⁻²·K⁻⁴ (Materials/constants.jl, exported as STEFAN_BOLTZMANN). The fields are: the view factor $F\in[0,1]$ (fraction of the surface's hemisphere that exchanges with surroundings at $T_\infty$; default $1$), the emissivity $\varepsilon$ (either a supplied constant/function or, when nothing, the composition-weighted $\varepsilon_{\text{eff}}$ of §12.2.6), and the surroundings temperature $T_\infty$. This term is the surface's own re-radiation to and absorption from a gray environment; it is the dominant cooling mechanism for a hot pyrolyzing surface. The cubic temperature derivative is the source of the $\sim$20% error that would result from naively evaluating radiation at the cell-center temperature — and the reason the Newton solve exists. (boundary_fluxes.jl:58–68, boundary_temperature.jl:357–380.)

The neglected $d\varepsilon/dT$ term (§12.2.4) applies here when $\varepsilon$ is temperature-dependent.

12.3.5 RadiativeFluxBC (external incident radiation; model-dependent)

This BC represents incident radiation from an external source (cone heater, furnace, radiant panel) — not the surface's own emission. The absorbed fraction is $\alpha\,q_{\text{ext}}$, where $q_{\text{ext}}=$ flux (constant or $q(t)$) and the absorptivity $\alpha$ is either supplied or, when nothing, defaults to $\varepsilon_{\text{eff}}$ via Kirchhoff's law (§12.2.6).

Its behavior is dispatched on the global radiation_model passed to solve() (§8, Thermal Radiation). Note the defaults diverge: solve() defaults to NO_RADIATION, while several internal helpers (apply_thermal_bc_flux!, lower-level RHS entry points) default to SURFACE_ABSORPTION — a trap for direct callers of the internals; always pass radiation_model explicitly when bypassing solve(). The crucial design point is that incident radiation must be deposited exactly once:

• NO_RADIATION and SURFACE_ABSORPTION: the radiation is treated as a surface heat flux. The Newton solve receives

  \[q_{\text{BC}} = \alpha\,q_{\text{ext}}(t), \qquad dq_{\text{BC}}/dT_s = 0,\]

  so the absorbed flux enters the surface energy balance and correctly raises $T_s$ above $T_{\text{cell}}$ when there is net heating (boundary_temperature.jl:414–423).
• BEER_LAMBERT and P1_QUASI_STEADY: the radiation is deposited volumetrically by compute_radiation_source! (Beer–Lambert) or by the P1 residual path. To avoid double-counting, the boundary flux is forced to zero,

  \[q_{\text{BC}} = 0, \qquad dq_{\text{BC}}/dT_s = 0\]

  (boundary_temperature.jl:425–429, boundary_fluxes.jl:84–99).

A note on consistency. The Newton-solve evaluator (evaluate_bc_flux_and_derivative) applies $\alpha q_{\text{ext}}$ as a boundary flux for SURFACE_ABSORPTION, while compute_radiation_source! deposits zero volumetric radiation for SURFACE_ABSORPTION. These are consistent: SURFACE_ABSORPTION absorbs at the boundary face only, so the energy enters via the boundary flux and the volumetric source is zero. Only BEER_LAMBERT (and the experimental P1 path) routes energy through the volume, and only then is the boundary flux suppressed. The diagnostic decomposition compute_boundary_heat_fluxes follows the same dispatch: its net is the flux the Newton balance actually applies — $\alpha q_{\text{ext}}$ is inside net under NO_RADIATION/SURFACE_ABSORPTION and absent under BEER_LAMBERT/P1, where the volumetric deposit is reported separately as extras.volumetric_absorption (§12.10.1).

Implementation note. The incident intensity is extracted at compile time from the (possibly combined) BC by extract_radiative_flux_from_bcs(bc_set, Val(tag), t), which returns (I, has_explicit_absorptivity) — the flux with an explicit absorptivity already folded in, or the bare incident flux whose $\varepsilon_{\text{eff}}$ scaling happens at injection (boundary_conditions.jl, §8.3.4). The intensity is passed as I_radiation to compute_rhs! together with the radiation_kirchhoff flag.

12.3.6 DirichletThermalBC (prescribed surface temperature)

A prescribed surface temperature $T=T_{\text{bc}}(t)$ is enforced by one of two mechanisms, selected by how the BC appears, rather than a hard algebraic constraint:

1. In the Newton solve (boundary_temperature.jl): a penalty "flux" with a very large conductance,

   \[q_{\text{BC}} = h_{\text{large}}\,(T_{\text{bc}} - T_s),\qquad dq_{\text{BC}}/dT_s = -h_{\text{large}}, \qquad h_{\text{large}}=10^{8} \;\text{W·m}^{-2}\text{·K}^{-1}.\]

   The large slope drives $T_s\to T_{\text{bc}}$ to high accuracy in one step. The value $10^8$ is large enough to pin $T_s$ to $T_{\text{bc}}$ within Newton tolerance while keeping the iteration well conditioned (a still-larger penalty would degrade conditioning with no accuracy benefit).
2. Standalone, in the RHS assembly (finite_volume.jl): a bare DirichletThermalBC is special-cased and bypasses the Newton path entirely. The surface temperature is the prescribed value, $T_s = T_{\text{bc}}(t)$, and the conducted flux is written directly,

   \[q_{\text{bc}} = k_{\text{cell}}\,\frac{T_{\text{bc}} - T_{\text{cell}}}{\Delta z_{1/2}}, \qquad q[f]= \begin{cases}+q_{\text{bc}},&\text{bottom},\\ -q_{\text{bc}},&\text{top},\end{cases}\]

   with $\Delta z_{1/2}=\tfrac12\,V_{\text{cell}}/A_{\text{face}}$ — the boundary face's own area, identical to the Newton path's distance, so the two mechanisms agree under a variable cross-section. The bypass also updates the diagnostic BoundaryState ($T_{\text{surface}}=T_{\text{bc}}$, $q_{\text{applied}}=q_{\text{conducted}}=q_{\text{bc}}$, converged = true), so surface-temperature and flux diagnostics report the prescribed state. The penalty mechanism of item 1 is used only when the Dirichlet BC is nested inside a CombinedThermalBC (the composite goes through the Newton path; the two mechanisms agree to $\sim10^{-4}$ K).

The penalty approach is acknowledged as less clean than a direct constraint (a genuine algebraic constraint would require a DAE index reduction); it is robust in practice because the prescribed temperature pins the boundary tightly.

12.3.7 CombinedThermalBC (additive superposition)

Real surfaces experience several mechanisms at once. The + operator builds an additive composite, e.g. the canonical cone-calorimeter top BC

bc = RadiativeFluxBC(flux = 50e3)            # absorbed incident cone flux
   + RadiativeBC(ε = nothing, T_inf = 300)   # surface re-radiation loss
   + ConvectiveBC(h = 10.0, T_inf = 300)     # convective loss

Fluxes and derivatives sum termwise,

\[q_{\text{BC}} = \sum_m q_{\text{BC},m}(T_s,t),\qquad \frac{dq_{\text{BC}}}{dT_s} = \sum_m \frac{dq_{\text{BC},m}}{dT_s}.\]

The sum is built by a @generated function that unrolls the tuple at compile time for type stability (boundary_temperature.jl:460–514, boundary_fluxes.jl:109–120). Because the radiation model can affect a RadiativeFluxBC member, the combined evaluator threads radiation_model through to every child. The + operator overloads handle every nesting case (BC+BC, combined+BC, BC+combined, combined+combined) so associativity holds (boundary_conditions.jl:231–236).

12.3.8 Ambient (reference) temperature resolution

Convective, radiative, and (when enabled) transpiration terms need a reference ambient $T_\infty$. This is resolved per-boundary by ambient_temperature(bc, t, T_{\text{default}}) (boundary_conditions.jl:331–357):

• ConvectiveBC and RadiativeBC return their own T_inf (evaluated at $t$);
• all other thermal BCs (HeatFluxBC, AdiabaticBC, RadiativeFluxBC, DirichletThermalBC) carry no ambient: the two-argument form ambient_temperature(bc, t) returns nothing for them, and the three-argument form substitutes the supplied default (typically the package constant $T_{\text{REF}}=298.15$ K);
• CombinedThermalBC walks its children in order and returns the first ambient-carrying child's value, which preferentially surfaces the ConvectiveBC/RadiativeBC ambient over the flux-only members.

The nothing sentinel (rather than comparison against a default value) makes a user-set ambient of exactly 298.15 K distinguishable from "this BC has no ambient"; resolution never silently falls through just because the value equals the default.

This lets each boundary carry its own effective ambient without a per-cell datum. The same resolver, applied to the top boundary first then the bottom, also selects the single scalar $T_{\text{ref}}$ used by the optional volumetric generation sink (finite_volume.jl, _resolve_volumetric_T_ref); the limitation of that single scalar under asymmetric ambients is discussed in §12.13.

12.4 Mass (species) boundary conditions

The species boundaries supply the gas mass flux $J_j$ at each boundary face. Unlike the thermal boundary, there is no surface-temperature nonlinearity, so the flux is evaluated directly from the adjacent cell state. The internal diffusive flux model is the volume-fraction (partial-pressure) gradient form of §9, Gas-Phase Mass Transport & Pressure; here it is closed at the surface by an external resistance.

12.4.1 ImpermeableBC (default)

$J_j = 0$ for all components — no gas crosses the boundary. This is the default for both faces (boundary_fluxes.jl) and is the physically correct choice for the substrate (sealed back face) and any non-venting surface. It is also the implicit closure for solid/liquid components everywhere, since only gas components are transported.

12.4.2 MassFluxBC

A prescribed species flux $J_j = J(t)$ (constant or function), positive into the domain (boundary_fluxes.jl). Used to impose a known gas influx/efflux, typically zero or negative (gas escaping).

12.4.3 ConvectiveMassBC (external film transfer)

\[J_j = h_m(t)\,\bigl(\xi_{j,\infty} - \xi_{j,s}\bigr),\]

the mass-transfer analog of Newton's cooling law, with the film coefficient $h_m$ (m·s⁻¹, velocity-based) and the ambient concentration $\xi_{j,\infty}$ (typically $0$ for pyrolysis gas escaping to atmosphere) (boundary_fluxes.jl). The coefficient may be estimated by the heat–mass analogy,

\[h_m \approx \frac{h_{\text{conv}}}{\rho_{\text{air}}\,c_{p,\text{air}}},\]

e.g. $h_{\text{conv}}\approx10$ W·m⁻²·K⁻¹ gives $h_m\approx0.008$ m·s⁻¹. With $\xi_{j,\infty}=0$ and $\xi_{j,s}>0$, $J_j<0$ (escaping), consistent with the flux-into-domain sign.

The ambient $\xi_{j,\infty}$ accepts three forms, resolved per gas species $j$ by ambient_concentration(bc, t, j): a Real (the same ambient for every gas species — typically 0, the atmosphere carrying no pyrolyzate), a Function of $t$ (time-varying scalar, broadcast to every gas species), or an NTuple{NC} of per-component ambient concentrations indexed by the material's component order (condensed-phase entries are ignored; entries are validated non-negative at construction and the tuple length against $N_C$ at problem validation). The per-species form is what represents ambient air ingress without injecting pyrolyzate: a scalar $\xi_\infty>0$ would feed that concentration of every gas species into the domain.

12.4.4 DiffusiveMassBC (combined internal/external resistance)

This BC couples the internal Fickian diffusion resistance from the surface cell to the face with an external film resistance, in series. Working in the volume-fraction formulation $J=-\lambda\,\partial\xi/\partial x$ with the surface gradient approximated as $(\xi_s-0)/d$ (zero exterior concentration), the two resistances are

\[R_{\text{internal}} = \frac{d}{\lambda}, \qquad R_{\text{external}} = \frac{1}{h_m},\]

where $d=\Delta z_{1/2}$ is the cell-center-to-face distance and $\lambda$ is the effective gas-transfer coefficient. The combined flux is

\[J_j = -\frac{\xi_{j,s}}{R_{\text{internal}}+R_{\text{external}}} = -\frac{\lambda\,h_m\,\xi_{j,s}}{\lambda + h_m\,d}.\]

The minus sign places escaping gas at $J_j<0$ (flux-into-domain convention). Two limits follow:

• No external resistance ($h_m\to\infty$, the default DiffusiveMassBC() with h_m = Inf): pure internal diffusion,

  \[J_j = -\frac{\lambda\,\xi_{j,s}}{d};\]

• Impermeable ($h_m\to0$): $J_j\to0$.

Implementation note. The flux law is in mass_bc_flux(bc::DiffusiveMassBC, ξ_surface, λρ, ρ, d, t, j) (boundary_fluxes.jl) — every mass-BC flux method carries the trailing species index j so per-species ambients resolve correctly. The discretization passes the combined transport coefficient $\lambda\rho$ (consistent with the rigorous gas-transport closure) and the density $\rho$; for $\rho \le 0$ the flux is zero, and otherwise the function recovers $\lambda = \lambda\rho/\rho$ by exact division (no floor — the callers construct $\lambda\rho$ as the product, so division is exact for any positive density), because $h_m$ is defined per concentration, not per mass fraction. The surface-cell quantities $(\xi_{j,s},\,\lambda,\,\rho,\,d)$ are read in _apply_mass_bc_static! (finite_volume.jl).

Sign reconciliation. Three sign conventions coexist and must be reconciled (boundary_fluxes.jl):

The discretization bridges these by orienting per boundary (_apply_mass_bc_static!, finite_volume.jl): at the bottom $J[j,f]\mathrel{+}=J_{\text{bc}}$, at the top $J[j,f]\mathrel{-}=J_{\text{bc}}$. With $J_{\text{bc}}<0$ for escaping gas, the top-face stored flux is positive — gas leaving in $+z$ at the exposed surface — exactly what the divergence operator needs. The application accumulates rather than overwrites: in Darcy–Fick mode the boundary face already carries the pressure-BC advective venting flux (§9.5.3), and the film-law flux adds to it; in Fickian mode the boundary faces arrive zero so the accumulation is equivalent to assignment.

12.4.5 CombinedMassBC

Like the thermal composite, $J_j=\sum_m J_{j,m}$, built by a @generated unrolled sum (boundary_fluxes.jl), e.g. DiffusiveMassBC(Inf) + ConvectiveMassBC(h_m=0.01, ξ_inf=0.0). Pressure BCs (§12.5) participate in the composite by contributing zero diffusive flux (they act on the Darcy velocity, not the diffusive flux), so they can be combined freely with DiffusiveMassBC (boundary_fluxes.jl).

12.5 Pressure boundary conditions (Darcy flow)

In Darcy–Fick transport mode (§9), gas advection is driven by the Darcy velocity $v_g = -(\kappa/\mu)\nabla P$, and the pressure boundaries set $v_g$ at the boundary faces. These BCs are subtypes of AbstractMassBC so they live in the mass slot of the BoundaryConditionSet; they do not contribute a diffusive flux (returning $0$ from mass_bc_flux, §12.4.5) — they act through apply_pressure_bc_to_darcy_velocity! (darcy_velocity.jl), which runs after the interior Darcy velocities are computed and sets the boundary face velocity $v_b$. That velocity then enters the species equations as the boundary advective venting flux $J_j=v_b\,\xi_j^{\text{int}}/\varphi_{\text{int}}$ (outflow only; §9.5.3), written by compute_gas_fluxes_darcy_fick! before the diffusive mass BC is added.

12.5.1 DirichletPressureBC (prescribed pressure)

Prescribes the boundary pressure $P_{\text{bc}}$ (default $101325$ Pa, atmospheric). The boundary face velocity follows from Darcy's law with a one-sided pressure gradient between the boundary and the interior cell:

\[v_g = -\frac{\kappa_{\text{int}}}{\mu_{\text{int}}}\,\nabla P, \qquad \nabla P = \begin{cases} (P_{\text{bc}} - P_{\text{int}})/d, & \text{top},\\[2pt] (P_{\text{int}} - P_{\text{bc}})/d, & \text{bottom}, \end{cases}\]

with $d$ the cell-center-to-face distance, $\mu_{\text{int}}$ the Sutherland viscosity at the cell temperature, and the velocity clamped to zero if $\kappa_{\text{int}}\le0$ or $\mu_{\text{int}}<10^{-20}$ Pa·s (darcy_velocity.jl:369–410). The orientation of the gradient ensures the correct outflow direction at each face: at the top a higher interior pressure ($P_{\text{int}}>P_{\text{bc}}$) yields $v_g>0$ (outflow in $+z$); at the bottom the same drives $v_g<0$ (outflow in $-z$).

12.5.2 NeumannPressureBC (prescribed gradient)

Prescribes the pressure gradient $\partial P/\partial n$ along the outward normal (default $0$ = impermeable). For zero gradient the boundary velocity is set to zero; for a nonzero prescribed gradient, $v_g=-(\kappa_{\text{int}}/\mu_{\text{int}})\,\partial P/\partial n$ at the top and $v_g=+(\kappa_{\text{int}}/\mu_{\text{int}})\,\partial P/\partial n$ at the bottom (where $\hat n=-\hat z$, so $\partial P/\partial z=-\partial P/\partial n$), with the same $\kappa,\mu$ guards (boundary_darcy_velocity, darcy_velocity.jl). The zero-gradient default makes a NeumannPressureBC() equivalent to an impermeable wall for Darcy flow.

12.5.3 ConvectivePressureBC (Robin / surface resistance)

A Robin-type law for a surface with finite resistance to gas escape (char layer, protective coating, blowing resistance):

\[v_g\cdot\hat n = h_P\,\bigl(P_{\text{int}} - P_{\text{ext}}\bigr),\]

with $h_P$ the pressure-transfer coefficient (m·Pa⁻¹·s⁻¹) and $P_{\text{ext}}$ the exterior pressure, gated on $\kappa_{\text{int}}>0$ — an impermeable surface cell cannot vent regardless of $h_P$. The implementation orients the outward normal per boundary (boundary_darcy_velocity, darcy_velocity.jl):

\[v_g[f] = \begin{cases} +h_P\,(P_{\text{int}}-P_{\text{ext}}), & \text{top},\\[2pt] -h_P\,(P_{\text{int}}-P_{\text{ext}}), & \text{bottom}. \end{cases}\]

The limits bracket the other two pressure BCs: $h_P\to\infty$ pins the boundary pressure (Dirichlet) and $h_P\to0$ seals it (zero-gradient Neumann / impermeable).

Implementation note. All three pressure BCs are scalar methods of the shared kernel boundary_darcy_velocity(bc, P_int, κ_int, μ_int, d, Val(tag)) (darcy_velocity.jl); a @generated overload sums the contributions of a CombinedMassBC, and the AbstractMassBC fallback contributes $0$ (e.g. a DiffusiveMassBC in the mass slot does not set a Darcy velocity). The face loop _apply_pressure_bc_faces! writes the kernel value per boundary face; the orchestration entry apply_pressure_bc_to_darcy_velocity! applies top then bottom. The same kernel is evaluated (with Dual arguments) inside the BCMassOperator Jacobian and (with the saved state) in the MLR diagnostics, so the three paths cannot drift apart.

12.6 Mesh-motion boundary conditions (ALE)

In ALE mode (§11, The ALE Framework) the node positions $z_i$ are state variables and the mesh-velocity boundary treatment is fixed by the formulation, not user-configurable: the bottom node is pinned ($w_1 = 0$, hard-coded in the residual's mesh-velocity integration) and the velocity at every other node follows cumulatively from the species dilation $\theta_i$ and the column-global radial strain $\theta_A$, $w_{i+1} = w_i + (\theta_i - \theta_A)\,\Delta z_i$ (§11.3, §15.3). The surface node is therefore automatically material-following (Lagrangian) — it recedes or expands with the integrated dilation.

There are no user-facing mesh-motion BC types: BoundaryConditionSet carries only thermal and mass slots, and the ALE bottom anchor ($w_1 = 0$) is hard-coded in the residual's mesh-velocity assembly. See §11 for the full ALE mesh-velocity derivation and the geometric conservation law.

12.7 The BoundaryConditionSet container

12.7.1 Structure and tag keying

All boundaries are bundled into a single type-stable container (boundary_conditions.jl:674–677):

struct BoundaryConditionSet{TT<:NamedTuple, MT<:NamedTuple}
    thermal::TT   # NamedTuple{tags, Tuple{<:AbstractThermalBC,...}}
    mass::MT      # NamedTuple{tags, Tuple{<:AbstractMassBC,...}}
end

Each boundary tag maps to exactly one thermal BC and one mass BC. The 1D convention uses :top ($z=L$) and :bottom ($z=0$); the field shape is designed so that a future 2D/3D constructor can add tags without changing the runtime contract. Storing the BCs in a typed NamedTuple (rather than, say, a Dict) means a symbol lookup compiles to a direct field read and the residual specializes statically.

12.7.2 Constructor, defaults, and validation

The keyword constructor (boundary_conditions.jl:702–722) merges the user-supplied partial NamedTuples with the 1D defaults

\[\texttt{thermal} = (\,\texttt{top}=\texttt{AdiabaticBC()},\ \texttt{bottom}=\texttt{AdiabaticBC()}\,), \qquad \texttt{mass} = (\,\texttt{top}=\texttt{ImpermeableBC()},\ \texttt{bottom}=\texttt{ImpermeableBC()}\,),\]

so a caller may specify only the slots it cares about. After merging it validates that every thermal payload is an AbstractThermalBC and every mass payload is an AbstractMassBC, throwing an ArgumentError if a BC was placed in the wrong slot (boundary_conditions.jl:712–720). The default — adiabatic and impermeable on both faces — is the inert closure: a column with no defined BCs neither gains energy nor loses mass.

12.7.3 Compile-time accessors

Access is through Val{tag} dispatch:

@inline get_thermal_bc(bc_set, ::Val{tag}) where {tag} = bc_set.thermal[tag]
@inline get_mass_bc(bc_set,    ::Val{tag}) where {tag} = bc_set.mass[tag]

At a call site such as get_thermal_bc(bc_set, Val(:top)), the Val(:top) resolves at compile time to a direct NamedTuple field read with zero overhead (boundary_conditions.jl:732–741). The RHS assembly uses these accessors so that each boundary's BC type is known statically and the per-type flux law is inlined. The tag$\to$id bridge tag_to_boundary_id (§12.1.3) connects the user-facing symbol API to the integer-keyed Mesh1D.boundary_faces / boundary_states dictionaries.

12.7.4 Constructor parameter validation

Individual BC constructors optionally validate their physical parameters at build time (validate_bc_params, boundary_conditions.jl:28–61): $h\ge0$, $T>0$ K, $\varepsilon\in[0,1]$, $F\in[0,1]$, $\alpha\in[0,1]$, $h_m\ge0$, $P>0$. Function-valued parameters cannot be checked at construction and are skipped (_get_scalar_value returns nothing for functions and property functions; boundary_conditions.jl:70–73).

12.8 Initial conditions

12.8.1 The state seed

The initial state is specified by a temperature closure $T_{\text{initial}}(z)$ and per-component concentration closures $\xi_{j,\text{initial}}(z)$, evaluated at each cell center $z_i^c$ (initial_conditions.jl:28–47):

\[T_i(0) = T_{\text{initial}}(z_i^c), \qquad \xi_{j,i}(0) = \xi_{j,\text{initial}}(z_i^c), \quad j=1,\dots,N_C .\]

The closures are evaluated and stored in the mutable mesh's cell states, then packed into the flat state vector (§15, Temporal Integration & Structured Jacobian). Cell volumes are not altered by initialization.

12.8.2 Auto-defaults

If the user supplies no initial conditions, the PyrolysisProblem constructor synthesizes them (pyrolysis_problem.jl:106, 114–134):

• $T_{\text{initial}}(z)=300$ K uniformly;
• $\xi_{1,\text{initial}}(z)=\rho_1(300\,\text{K})$ if the first component is a solid (the virgin solid at its reference density), and $\xi_{j,\text{initial}}=0$ for all other components.

This "first solid at reference density, everything else zero" default models a fresh, uniform virgin slab at room temperature — the standard cone/TGA starting point. The uniform-IC variant apply_initial_conditions!(mesh, material; T, ξ) follows the same rule when $\xi=\texttt{nothing}$ (initial_conditions.jl:103–112). Profile helpers (step_profile, linear_profile, gaussian_profile, tanh_profile) are provided for non-uniform seeds (e.g. a pre-heated layer).

12.8.3 The dry-solid mass denominator $m_{\text{dry},i}$

A second, BC-adjacent quantity is fixed at initialization: the per-cell initial dry-solid mass, which is the constant denominator of the pyrolysis-progress state $\chi_i = m_{\text{released,dry},i}/m_{\text{dry},i}$ used by the variable cross-section law (§10, Volume Change & Lateral Shrinkage). It is populated by populate_m_dry_solid_initial! (initial_conditions.jl:256–272):

\[m_{\text{dry},i} = V_i(0)\,\sum_{j\in\text{solid}} \xi_{j,i}(0),\]

summing only the solid components (the is_solid predicate — not liquids, not gases) at the initial state. This is computed once and stored on the mutable mesh as Float64; under AD the $\chi$ numerator is Dual-typed but this denominator stays Float64 (see §15, §17).

12.8.4 Verification

verify_initial_conditions checks that every $T_i>0$ and finite, every $\xi_{j,i}\ge0$ and finite, and that the total concentration $\sum_j\xi_{j,i}>0$ in each cell (initial_conditions.jl), catching an empty or negative seed before the solve begins.

12.9 Surface transpiration (optional; off by default)

Do not confuse this with the blowing correction. This section is the enthalpy-side transpiration term (rejected in production, §7.6). The momentum-side effect of the same gas efflux — suppression of the convective coefficient, $h \to h\,B/(e^B-1)$ — is the supported, opt-in use_blowing_correction documented in §12.3.3; it is independent of this term and composes with the default $S_{\text{conv}}$-only convention.

12.9.1 The transpiration term and its derivative

The surface energy balance of §12.2 can optionally include the sensible enthalpy carried out of the slab by escaping gas (the "blowing" or transpiration term):

\[q_{\text{BC}}(T_s) = k_{\text{eff}}\,\frac{T_s-T_{\text{cell}}}{\Delta z_{1/2}} \;+\; \sum_g \dot m_g\,\Delta h_g(T_s, T_\infty),\]

where $\dot m_g$ is the per-component outward gas mass flux at the surface (positive when leaving) and

\[\Delta h_g(T_s,T_\infty) = \int_{T_\infty}^{T_s} c_{p,g}(\tau)\,d\tau\]

is the sensible enthalpy difference, computed by the midpoint rule in delta_enthalpy (the same quadrature used by the volumetric generation sink $S_{\text{gen}}$, §7.5). The analytic Newton derivative of the transpiration term, used in $f'(T_s)$, is

\[\frac{dq_{\text{transpire}}}{dT_s} = \sum_{g:\ \text{is\_gas}} \dot m_g\,c_{p,g}(T_s) \qquad (\dot m_g>0,\ \text{outflow}).\]

This is the derivative of the true integral, which differs from the exact derivative of the midpoint-rule residual by $O\bigl((T_s-T_\infty)\,c_{p,g}'\bigr)$; the two agree to first order, so Newton still converges to the residual's root (unbiased $T_s$, marginally sub-quadratic rate) (boundary_temperature.jl:125–132).

12.9.2 Inflow/outflow upwinding

The transpiration term is upwinded on the sign of $\dot m_g$ (boundary_temperature.jl:233–263):

• Outflow ($\dot m_g>0$): the donor is the surface at $T_s$; the term carries $\dot m_g\,\Delta h_g(T_s,T_\infty)$ with derivative $\dot m_g\,c_{p,g}(T_s)$.
• Inflow ($\dot m_g\le0$): the donor is the ambient at $T_\infty=T_{\text{ref}}$; the surface enthalpy flux is $\dot m_g\,\Delta h_g(T_{\text{ref}},T_{\text{ref}})=0$ with zero derivative. Entering gas is heated to $T_s$ volumetrically (by $S_{\text{conv}}$) once inside, not at the face.

Upwinding to the donor temperature is what prevents an unphysical artifact in which inflow enthalpy would spuriously heat or cool the surface. The outward flux $\dot m_g$ itself is extracted from the gas-flux matrix with the boundary orientation: $\dot m_g = +J_{g,f}$ at the top, $\dot m_g=-J_{g,f}$ at the bottom (so that "outward" is positive at both faces) (finite_volume.jl). Note this requires the gas fluxes to be computed before the thermal BC solve — hence Step 3 (gas transport) precedes Step 4 (thermal BCs) in the RHS sequence (§7.8).

12.9.3 Mutual exclusivity with $S_{\text{conv}}$ — and why it is forbidden

The transpiration term is mutually exclusive with the volumetric gas-advection source $S_{\text{conv}}=-\sum_g c_{p,g}(T_i)\,\bar J_{g,i}\,(\partial T/\partial z)_i$, which is always active in the default :standard energy form (§7.4). The algebraic identity (with constant $c_p$)

\[\int \bigl(S_{\text{conv}} + S_{\text{gen}}\bigr)\,dz = -\sum_g \dot m_g\,\Delta h_g(T_s, T_{\text{ref}})\]

shows that the volumetric sources already account for the gas-enthalpy transport: $S_{\text{conv}}$ alone integrates to $-\sum_g\dot m_g\Delta h_g + \int S_{\text{gen}}$, i.e. it carries almost all of the surface outflow enthalpy. Adding the surface transpiration term on top would apply the full outflow a second time, double-counting the advective enthalpy and producing a spurious sink of order $5$–$10$ kW·m⁻² at typical pyrolysis rates.

For this reason use_transpiration_bc = true is rejected at problem construction in to_problem_def (adapter.jl), and the unsupported combination is guarded by a regression test (test/discretization/test_transpiration_energy_balance.jl). In the default configuration the caller zeroes $J_{\text{boundary,outward}}$ (finite_volume.jl), so the transpiration term vanishes and the surface balance reduces to $q_{\text{BC}}(T_s)=q_{\text{cond}}$. The supported way to move toward a conservative formulation is the :with_generation_sink energy form (§7.5), which adds $S_{\text{gen}}$; only if a future solver disables $S_{\text{conv}}$ would the surface transpiration term become consistent (see §16, Conservation Diagnostics, for the exact telescoping ledger). The supported convention is S_conv-only; parameter sets calibrated against solvers that additionally apply a surface transpiration term embed a $\sim$5–10 kW·m⁻² spurious sink and should be re-validated when imported (§18).

12.10 Boundary-flux diagnostics

12.10.1 Flux decomposition

For auditing and post-processing, compute_boundary_heat_fluxes(bc, T_s, t, ε_eff, radiation_model[, mcp]) decomposes the boundary flux into named components (boundary_fluxes.jl). The optional trailing mcp is the blowing conductance $\dot m'' c_{p,g}$ at the converged $T_s$ (BoundaryState.blowing_mcp); when supplied, the convective component is reported with the same $h\,B/(e^B-1)$ rescaling the surface Newton applied (§12.3.3), keeping the decomposition consistent with net. Components:

• absorbed — incident radiation absorbed at the boundary face [W·m⁻²];
• reradiation — re-radiation loss (positive = leaving);
• convection — convective loss (positive = leaving);
• prescribed — prescribed flux;
• net — the flux the surface energy balance actually applies at the face (positive into the domain).

For ConvectiveBC and RadiativeBC, the loss components are reported with the opposite sign to the into-domain convention (so a hot surface re-radiating shows a positive reradiation). Per component, net matches evaluate_bc_flux_and_derivative — the Newton-applied flux — for every radiation model, so the ledger identity $\texttt{net} = q_{\text{cond}} + q_{\text{transp}}$ holds at convergence (§12.10.2, §16.7). For combined BCs the components and net are summed over children; net is never recomputed from the components. Consequences of the model dispatch:

• NO_RADIATION / SURFACE_ABSORPTION: a RadiativeFluxBC contributes absorbed = net = α·q_ext (α the explicit BC absorptivity or ε_eff) — the flux is applied at the face.
• BEER_LAMBERT / P1: a RadiativeFluxBC contributes zero absorbed and zero net — the radiation enters volumetrically (§8.3) and is reported separately as extras.volumetric_absorption $= \alpha_{\text{surf}}\,q_{\text{ext}}\, (1-e^{-\tau_{\text{tot}}})$, the exact telescoped sum of the Beer–Lambert deposit on the saved state.

12.10.2 Energy-balance check

check_boundary_energy_balance(mesh) returns, per boundary, the relative imbalance $|q_{\text{applied}}-q_{\text{conducted}}|/|q_{\text{applied}}|$ (the absolute difference when $|q_{\text{applied}}|\le10^{-10}$), which should be $\approx$ machine precision ($\sim10^{-15}$) when the Newton solve has converged (boundary_temperature.jl). This is a direct continuous self-check that the surface energy balance $q_{\text{BC}}=q_{\text{cond}}$ holds. The accessors get_surface_temperature(mesh) and get_back_face_temperature(mesh) return the solved $T_s$ at the top (id 1) and bottom (id 2) boundaries (boundary_temperature.jl).

12.11 Worked example: the cone-calorimeter top boundary

Assembling the canonical top boundary makes the pieces concrete. A 50 kW·m⁻² cone with surface re-radiation and natural convection to a 300 K ambient, on a sealed-but-venting slab, is

bc_set = BoundaryConditionSet(
    thermal = (
        top    = RadiativeFluxBC(flux = 50e3)            # absorbed incident, α = ε_eff
               + RadiativeBC(ε = nothing, T_inf = 300)   # re-radiation loss
               + ConvectiveBC(h = 10.0, T_inf = 300),    # convective loss
        bottom = AdiabaticBC(),                          # insulated back
    ),
    mass = (
        top    = DiffusiveMassBC(),                      # gas vents, internal diffusion limit
        bottom = ImpermeableBC(),                        # sealed back
    ),
)

With radiation_model = SURFACE_ABSORPTION, the surface energy balance solved by Newton at the top is

\[\underbrace{\alpha\,(50{,}000)}_{\text{absorbed}} \;-\;\underbrace{\varepsilon_{\text{eff}}\,\sigma\,(T_s^4 - 300^4)}_{\text{re-radiation}} \;-\;\underbrace{10\,(T_s - 300)}_{\text{convection}} \;=\; k_{\text{eff}}\,\frac{T_s - T_{\text{cell}}}{\Delta z_{1/2}},\]

with $\alpha=\varepsilon_{\text{eff}}$ (Kirchhoff, both nothing). As the surface chars, $\varepsilon_{\text{eff}}$ tracks the composition automatically, raising both the absorbed and re-radiated terms. The Newton derivative is $f'(T_s) = -4\varepsilon_{\text{eff}}\sigma T_s^3 - 10 - k_{\text{eff}}/\Delta z_{1/2}$ (the absorbed term and convection-ambient are $T_s$-independent), which is strictly negative — guaranteeing a unique root reached in a few iterations. The conducted flux $q_{\text{conducted}}$ is written to the top face as $q[f_{\text{top}}]=-q_{\text{conducted}}$, and the vented gas leaves the top with $\dot m_g$ accumulated by DiffusiveMassBC and oriented to $+z$ via the top-face sign flip. If radiation_model = BEER_LAMBERT is selected instead, the RadiativeFluxBC contributes zero to this balance and the 50 kW·m⁻² enters the energy equation as a volumetric Beer–Lambert source $q'''=\alpha I$ (§8).

12.12 Comparison to ThermaKin2Ds, Gpyro, and FDS

The three reference solvers express the same surface physics in their own notation; the subsections below map their symbols to ours on first use and contrast the modeling choices. (ThermaKin/Gpyro publish $h>0=$ exothermic; our internal convention is $h>0=$ endothermic with $Q_{\text{rxn}}=-\sum h_r r_r$, but this sign issue is confined to the volumetric reaction heat, not the surface BCs.)

12.12.1 ThermaKin2Ds

ThermaKin2Ds applies piecewise external/flame boundary conditions and supports an external radiant source with surface re-radiation. The exposed-face energy balance, in ThermaKin's form, is

\[-k\frac{\partial T}{\partial x}\Big|_{\text{surf}} = \alpha\,\dot q''_{\text{ex}} - \varepsilon\sigma\bigl(T_s^4 - T_\infty^4\bigr) - h_{\text{conv}}\bigl(T_s - T_\infty\bigr),\]

with $\dot q''_{\text{ex}}$ the incident external flux, $\alpha$ the absorptivity and $\varepsilon$ the emissivity (Kirchhoff $\alpha=\varepsilon$). Pyrolysis.jl's default cone BC RadiativeFluxBC + RadiativeBC + ConvectiveBC is exactly this balance, with our $q_{\text{ext}}=\dot q''_{\text{ex}}$ and our $\varepsilon_{\text{eff}}$ playing the role of ThermaKin's surface $\varepsilon$. ThermaKin also tracks a reflected component $I_{\text{ex}}^0$; our composition-weighted $\varepsilon_{\text{eff}}$ supplies $\alpha=\varepsilon_{\text{eff}}$ and the reflected fraction $(1-\alpha)$ is implicitly discarded (opaque-surface assumption), consistent with the SURFACEABSORPTION model. The shared idea — that both the absorbed external flux and the surface re-radiation depend on the true surface temperature — is realized by ThermaKin through its Crank–Nicolson surface node and by Pyrolysis.jl through the explicit Newton surface-balance solve of §12.2. ThermaKin carries gas enthalpy via a convective term ``cg(Jg\cdot\partial T/\partial x)$(its Eq. 17), which is precisely our$S{\text{conv}}`` and the reason transpiration is omitted by default (§12.9.3).

12.12.2 Gpyro

Gpyro (Lautenberger & Fernandez-Pello) writes the exposed-surface balance as a radiative–convective balance with an impermeable back face by default. In Gpyro notation,

\[-\bar k\frac{\partial T}{\partial z}\Big|_{\text{surf}} = \dot q''_e - \varepsilon\sigma\bigl(T_s^4 - T_\infty^4\bigr) - h_c\bigl(T_s - T_\infty\bigr),\]

and gas escape at the surface uses a heat–mass transfer analogy for the film coefficient, $h_m \sim h_c/(\rho c_p)$ — the same estimate we recommend for ConvectiveMassBC (§12.4.3). Gpyro's control-volume formulation places the surface temperature at the face via its harmonic interface interpolation; our Newton solve achieves the equivalent by solving the one-sided Fourier balance for $T_s$. Gpyro permits a non-zero back-face pressure/mass boundary in its pressure-evolution model; our DirichletPressureBC / NeumannPressureBC / ConvectivePressureBC (§12.5) are the analogous Darcy boundary closures, with Gpyro's permeability $K$ mapping to our $\kappa$. Gpyro's default impermeable back face corresponds to our default ImpermeableBC + NeumannPressureBC() (zero gradient).

12.12.3 FDS

The FDS solid-phase model imposes the 1D surface boundary condition

\[-k\frac{\partial T_s}{\partial x}\Big|_{x=0} = \dot q''_c + \dot q''_r,\]

i.e. the conducted flux equals the sum of the convective and net radiative surface fluxes, with $\dot q''_r$ including absorbed incident and emitted re-radiation. This is the same content as our $q_{\text{cond}}=q_{\text{BC}}$ with $q_{\text{BC}}=q_{\text{ext}}-\varepsilon\sigma T_s^4-h_{\text{conv}}(T_s- T_\infty)$. FDS additionally performs a surface oxygen Newton iteration for oxidative pyrolysis (coupling surface $O_2$ to the gas phase); Pyrolysis.jl is a condensed-phase-only solver with no gas-phase $O_2$ feedback, so our surface Newton iterates on $T_s$ alone. FDS recovers Arrhenius parameters from TGA reference temperature/rate (a kinetics, not a BC, concern; see §6, §18). The FDS in-depth radiation (a two-flux Schuster–Schwarzschild model) corresponds to our BEER_LAMBERT / P1 volumetric handling, under which the RadiativeFluxBC boundary flux is suppressed to avoid double-counting (§12.3.5), just as FDS deposits the in-depth absorbed fraction volumetrically rather than at the face.

12.13 Limitations and open issues

1. Asymmetric ambients with a permeable back face. The volumetric generation-sink reference $T_{\text{ref}}$ is resolved as a single scalar (top ambient, else bottom, else $T_{\text{REF}}$) and applied to every cell (finite_volume.jl:601–620). If the top and bottom ambients differ and the back face passes gas, the :with_generation_sink energy ledger fails to close at the bottom by exactly $\sum_g \dot m_{g,\text{bottom}}\,c_{p,g} (T_{\infty,\text{top}}-T_{\infty,\text{bottom}})$. This is harmless for the usual impermeable back face (the bottom gas flux vanishes) and for equal ambients; supporting it would require a per-cell datum rather than a scalar.

2. Temperature-dependent emissivity derivative. The radiative Newton derivative omits the $F\,(d\varepsilon/dT)\,\sigma\,(T_\infty^4-T_s^4)$ term (§12.2.4). The converged $T_s$ is unbiased, but convergence can become sub-quadratic for strongly temperature-sensitive $\varepsilon(T)$.

3. Transpiration BC integration. The surface transpiration term is implemented and tested but is opt-in and currently forbidden in production because it double-counts with $S_{\text{conv}}$ (§12.9.3). A consistent conservative-form path would require disabling $S_{\text{conv}}$ when transpiration is on; that combination is not yet supported by the solver.

4. Dirichlet thermal BC via penalty. DirichletThermalBC uses a large conductance $h_{\text{large}}=10^8$ in the Newton solve and a special-cased direct flux in the RHS (§12.3.6). A genuine algebraic constraint (DAE) would be cleaner and would remove the residual conditioning sensitivity of the penalty.

5. View-factor meaning in 1D. The RadiativeBC view factor $F$ is a scalar surrogate for a 3D geometric quantity; in a strictly 1D model it is best read as a tuning factor for partial enclosure rather than a rigorously computed exchange factor.

6. Mass-transfer coefficient estimation. The film coefficient $h_m$ for ConvectiveMassBC / DiffusiveMassBC is supplied by the user, typically via the heat–mass analogy $h_m\approx h_{\text{conv}}/(\rho_{\text{air}} c_{p,\text{air}})$; there is no built-in correlation, so its accuracy depends on the user's estimate of $h_{\text{conv}}$ and the flow regime.

7. Pressure-BC integration maturity. The three pressure BCs are defined and wired into apply_pressure_bc_to_darcy_velocity!, but their interaction with the full Darcy–Fick velocity solve and the structured Jacobian is still maturing relative to the thermal and diffusive-mass paths; in Eulerian/Fickian mode they are inert.