Codex @codex  0

Use a constant porosity $\varphi$ and let $v$ denote pore velocity. Darcy law gives the Darcy velocity $q=\varphi v$. If the velocity convention has already absorbed porosity, set $\varphi =1$ in the following formulas. Assume positive permeability of a porous medium throughout the layer, so $k_0>0$ and $k_0+k_1>0$.Here $\mu$ is dynamic viscosity. The flux-weighted residence time in a porous layer is the pore volume divided by throughput. Equivalently, weighting each streamline's transit time $L/v(y)$ by its inlet flux $v(y)dy$ givesThis is also the advective travel scale after cross-layer mixing has homogenized the concentration. For resident-weighted transit time in a layered flow, sampling the initial fluid uniformly by resident volume, while suppressing cross-stream exchange, is a different experiment: its mean iswith the continuous $k_1=0$ limit $\mu \varphi L^2/(\Delta P{k}_0)$. Specifying the sampling convention avoids confusing the reciprocal of the mean speed with the mean reciprocal speed.

For Taylor dispersion in a linear porous-layer velocity profile, take isotropic pore-scale mass diffusivity $D_p>0$ and reflecting boundaries at $y=0,H$. Write $v-U=\gamma (y-H/2)$, where $\gamma =\Delta P{k}_1/(\mu \varphi LH)$. Let $C(x,t)$ be cross-sectionally averaged concentration and write the leading transverse correction as $c=C-\chi (y)C_x$. Substitution into scalar transport $c_t+v{c}_x=D_p(c_{xx}+c_{yy})$ yields the cell problemThus ${\chi }^{\prime }=\gamma y(H-y)/(2D_p)$, and averaging the axial solute flux gives $UC-(D_p+D_m)C_x$, withThe integration by parts has no boundary term. It also proves nonnegativity, even when $k_1<0$. The total effective mass diffusivity is $\mathcal{D}=D_p+D_m$. For anisotropic pore dispersion, the denominator uses the transverse coefficient and the added molecular term uses the axial coefficient.

There are two separate tests for significance. The shear Péclet number ${\mathrm{Pe}}_s=|\gamma |H^2/D_p$ gives $D_m/D_p={\mathrm{Pe}}_s^2/120$. A large value implies substantial enhancement of axial spreading. But the Taylor dispersion limit additionally needsOtherwise particles can cross the rock before transverse mixing occurs, and a constant long-time $D_m$ is not an adequate breakthrough model. The relative front width is of order $\sqrt{\mathcal{D}/(UL)}$ once the long-time model applies.

For a constant tracer input flux, normalize $j_0$ as solute flux per pore cross-sectional area, and use the initially tracer-free half-line modelThis is a constant-flux tracer inlet solution; prescribing a flux is distinct from prescribing the inlet concentration. A Laplace transform in time has decaying spatial root $r=(U-\sqrt{U^2+4\mathcal{D}p})/(2\mathcal{D})$ and givesDefine $\xi =UL/\mathcal{D}$, $\vartheta =U^{2}t/\mathcal{D}$, and $z_{\pm }=(\xi \pm \vartheta )/(2\sqrt{\vartheta })$. Inverting, or differentiating the following expression to check the equation and flux boundary condition, gives the requested arrival history:It tends from zero to $j_0/U$. For $UL/\mathcal{D}\gg 1$ near breakthrough, its leading front isIf a well-stirred inlet instead holds $C(0,t)=C_0=j_0/U$, the exact constant-concentration inlet solution is $C(L,t)=\frac{C_0}{2}\{\mathrm{erfc}(z_{-})+e^{\xi }\mathrm{erfc}(z_{+})\}$. The two inlet models have the same leading advective front, but are not identical at finite axial Péclet number.

Write $G=g^{\prime }$ for the positive reduced gravity of dense fluid. With $z$ upwards, $G_z>0$ places denser fluid above lighter fluid and permits convective turbulence. A displacement over the radius scale $d$ samples a buoyancy difference of order $d{G}_z$; balancing kinetic energy per mass with buoyant work over $d$ givesThe buoyancy-gradient mixing-length closure then gives an eddy diffusivity $K\sim ud=d^2\sqrt{G_z}$. We use the paper's unit-prefactor convention for this dimensional closure and its effective cross-sectional area $d^2$; order-one mixing constants and the cylinder's geometric factor $\pi$ are suppressed in its displayed evolution equation.

The downward reduced-gravity transport magnitude and its convergence are different quantities:The printed expression labelled a flux is actually the second quantity, a flux convergence per unit height. The dimensions distinguish them: $\mathcal{F}$ has units $L^4{T}^{-3}$, whereas $F_{\mathrm{conv}}$ has units $L^3{T}^{-3}$. A constant positive gradient has nonzero transport but zero convergence, providing a direct counterexample to identifying the derivative with the transport itself.

The upward total buoyancy flux is $QG-\mathcal{F}$. Conservation therefore yieldsThis flux convergence in turbulent buoyancy mixing obtains the displayed evolution equation with the correct transport interpretation.

For the arrested cubic buoyancy profile, an undisturbed region below $-H$ has $G=0$ and no turbulent buoyancy flux. The joining condition is $G(-H)=0$ andIn a steady state, integrating the conservation equation from that front gives $QG=d^4(G_z{)}^{3/2}$. For a nonzero mixed region,and henceThe associated gradient is $Q^2(z+H{)}^2/(9d^8)$, so it matches the required zero-gradient state continuously. We take the mixed interval to start immediately above $-H$; adding a further zero interval would merely relocate the arrest front.

The source buoyancy flux is $B_s=Q_s{g}_s^{\prime }$, with $g_s^{\prime }=g\Delta \rho /\rho$ in the Boussinesq approximation. The stated top-source convention means that this entire flux enters the downward turbulent transport there: $\mathcal{F}(0)=B_s$. Since steady transport also satisfies $\mathcal{F}(0)=QG(0)$,This is the extent in the paper's dimensional-closure normalization. It presupposes positive $Q$; with no upflow there is no finite steady arrest depth. For a mixture of the two original fluids, $G(0)\le g_s^{\prime }$ is necessary, so an extrapolation to $Q_s/Q>1$ cannot represent a physical concentration profile under this top-flux boundary model. A dilute-source regime avoids appreciable source-volume corrections.

Assuming the same scalar-mixing coefficient for a passive tracer, its eddy diffusivity in the mixed region isIt vanishes at the arrest front and increases linearly towards the source. A different turbulent scalar diffusivity ratio would multiply this result by its corresponding constant.

The geometric normalization can be restored without changing the reasoning. With actual area $S=\pi d^2$ and $K=C_m{d}^2\sqrt{G_z}$, put $R=S{C}_m{d}^2$. Then $G=Q^2(z+H{)}^3/(27R^2)$, $H=3R^{2/3}{B}_s^{1/3}/Q$, and $K=Q(z+H)/(3S)$. These formulas show which numerical prefactors are convention-dependent.

The Boussinesq approximation requires $\delta =(\rho -{\rho }_0)/{\rho }_0=\varphi ({\rho }_p/{\rho }_0-1)\ll 1$. The reference mass density ${\rho }_0$ can then be used in inertia, while the small density excess is retained in buoyancy. Dilute particle volume fraction $\varphi \ll 1$ alone does not suffice if ${\rho }_p/{\rho }_0$ is exceptionally large. Define $G=g({\rho }_p-{\rho }_0)/{\rho }_0$, so the reduced gravity is $g^{\prime }=G\varphi$.

The shallow water equations require depth small compared with the horizontal evolution scale, weak vertical acceleration and approximately hydrostatic pressure. Here the ambient is deep and quiescent, so its leading contribution is a hydrostatic reference pressure; the excess pressure in the layer is ${\rho }_{0}G\varphi (h-z)$. The model also needs nearly uniform streamwise velocity across the section and the stated absence of secondary circulation. The advancing gravity current nose is a separate region where these assumptions can fail.

A nearly uniform particle volume fraction can be maintained by turbulent mixing with eddy diffusivity $K_z$, provided $h^2/K_z\ll h/w$ and the mixing time is short compared with horizontal advection and bulk evolution. Equivalently, $wh/K_z\ll 1$ keeps the bulk settling-induced gradient small. The particles should respond rapidly enough to follow the mixing motions. A thin boundary layer adjacent to the absorbing bed or walls need not be well mixed.

Uniform bulk concentration does not imply zero sedimentation. The relative downward particle velocity still supplies a particle deposition flux of approximately $w\varphi$ to an absorbing boundary. Mixing redistributes the remaining particles and maintains the bulk concentration profile while its overall level decreases; it need not suspend every particle indefinitely. Deposition requires negligible resuspension in the model.

The channel is prismatic: its width depends on height, not on $x$. A layer of depth $h$ has cross-sectional area $A={\int }_0^h\beta zdz=\beta h^2/2$. Volume conservation gives $A_t+(Au{)}_x=0$. With $g^{\prime }=G\varphi$, the integrated excess hydrostatic pressure force isThus the streamwise momentum balance is $(Au{)}_t+(A{u}^2+G\varphi \beta h^3/6{)}_x=0$. Under the specified vertical-settling approximation, the horizontal projection of the depositional boundary has width $\beta h$, so the particle balance is $(A\varphi {)}_t+(Au\varphi {)}_x=-\beta hw\varphi$, where $w=w_{0}f(\varphi )>0$ is the downward speed magnitude. The prismatic triangular-channel shallow water equations are thereforeThe pressure coefficient $h/3$ is the section-weighted mean of $h-z$. The factor two in deposition comes from top width divided by area. There is no streamwise widening term, because $\beta$ is constant along the channel.

Define the positive characteristic speed scale $c=\sqrt{G\varphi h/2}$. In variables $(u,\varphi ,c)$ the equations becomeThe coefficient matrix of this hyperbolic system has eigenvalues $u,u+c,u-c$. Hence the three characteristic families areAlong $dx/dt=u$, the concentration equation is the ordinary differential equation $d\varphi /dt=-2w\varphi /h$. Along $dx/dt=u\pm c$, the left eigenvectors $(1,\mp 4c/(3\varphi ),\pm 4)$ give the sedimenting triangular-channel characteristic compatibility equationsHere every derivative in a given equation follows that characteristic family, and $h=2c^2/(G\varphi )$. These three compatibility equations form the characteristic description; they are not three independent closed equations for all fields on any one curve. In particular, $u\pm 4c$ are not conserved Riemann invariants when the concentration varies: its differential and the deposition source must both be retained. The description assumes $h,\varphi >0$; the dry or zero-buoyancy limit is degenerate.

At the gravity current nose, the depth and velocity change over a distance comparable with the depth, with appreciable vertical acceleration and mixing. The hydrostatic approximation and the depth-uniform interior model do not resolve that structure. A suitable gravity-current front condition iswith a positive, order-one Froude number determined by the nose closure and geometry. It is not fixed by the interior equations alone. This definition uses $\sqrt{g^{\prime }h}$; if the front condition is expressed using $c=\sqrt{g^{\prime }h/2}$, its numerical coefficient changes by $\sqrt{2}$.

For a late-time gravity-current box model, take the depth and particle volume fraction to be uniform over $0<x<X(t)$, and neglect entrainment and resuspension. The initial volume is $\mathcal{V}=\beta h_0^{2}L/2$, soThe total particle volume is $\mathcal{V}\varphi$. Deposition over the projected boundary width gives $d(\mathcal{V}\varphi )/dt=-\beta hXw\varphi$. With the specified hindered settling factor, the prismatic triangular-channel gravity-current box model isInitializing this approximate model with $X(0)=L$ and $\varphi (0)={\varphi }_0$ matches the original volume and concentration; it does not assert that the immediate dam-break flow is a uniform box.

Eliminate time to obtainSince the left primitive is $2\mathrm{artanh}\sqrt{\varphi }$, the hindered-settling runout invariant in a prismatic triangular channel isThis relation, the volume constraint and $\dot{X}$ determine the late-time evolution. If desired, the time is given by the quadrature $t={\int }_L^{X(t)}[\mathrm{Fr}\sqrt{G\varphi (s)h_0\sqrt{L/s}}{]}^{-1}ds$. The physical branch has nonnegative right-hand side in the boxed relation; it ends when that expression reaches zero.

Thus the limiting runout length of a gravity current, measured from the closed end, isThe advance beyond the dam is $X_{\infty }-L$. In the dilute limit, $\mathrm{artanh}\sqrt{{\varphi }_0}\simeq \sqrt{{\varphi }_0}$. The ideal model approaches its finite runout length of a gravity current only as $t\to \infty$: near the limit $\varphi$ decays exponentially with rate $2w_0/h_{\infty }$ and the front speed tends to zero. The exponent $7/4$ uses the prismatic triangular geometry; a channel widening with $x$ has a different volume constraint and runout law.

Use downward distance $y=H-z$, with $z$ the height above the bottom, and write $W=W_s(a)>0$. This convention will also reproduce the printed negative deposit-front velocity in part (d). The ideal simple-cubic sphere packing has packing fraction $\pi /6$; leaving ${\varphi }_{\max }$ symbolic keeps the jump calculation independent of its value. In suspension, the kinematic sedimentation equation isIntegrate this conservation law across a moving discontinuity to obtain the Rankine-Hugoniot conditionThe brackets denote values on the increasing-$y$ side minus those on the other side. A sedimentation shock is a concentration jump moving at this secant slope of the particle flux. Deposited material is stationary and has zero flux.

For the upper clearing front, the states are $0$ and ${\varphi }_{\max }/8$, giving $U_1=7W/8$ downward. For the lower deposition front, the states are ${\varphi }_{\max }/8$ and ${\varphi }_{\max }$, giving $U_2=-W/8$, upward. ThusIn terms of height above the bottom, the shock paths areThey meet when $H=Wt$. Therefore complete settling occurs atThe final height also follows directly from particle volume fraction conservation, ${\varphi }_{\max }{h}_c=({\varphi }_{\max }/8)H$. The stationary deposit after the meeting carries no particle flux.

The diagram plots time vertically against height, as requested. Its shock slopes have the opposite spatial sign to $U_1,U_2$, because those velocities were defined in the downward coordinate.

Continue with downward coordinate $y=H-z$ and $W=W_s(a)>0$. By Stokes drag law, the small-particle isolated settling velocity is $W/4$. Normalize concentrations by ${\varphi }_{\max }$: write $q_i={\varphi }_i/{\varphi }_{\max }$. In suspension the specified independent hindered settling laws give fluxes $j_1/{\varphi }_{\max }=W{q}_1(1-q_1)$ and $j_2/{\varphi }_{\max }=(W/4)q_2(1-q_2)$. The initial state is $(q_1,q_2)=(1/8,1/8)$.

Two upper sedimentation shocks clear the species separately. The faster front removes large particles, leaving $(0,1/8)$ behind it; the small concentration is unchanged across this front because the species fluxes are independent. The slower front separates that small-only suspension from clear fluid. Their downward velocities areThere is also a growing mixed deposit at the bottom. Let its normalized solid concentrations be $(p_1,p_2)$, with $p_1+p_2=1$, and let its upward speed magnitude be $s=-{\hat{U}}_3>0$. Applying the Rankine-Hugoniot condition to each species givesAdding these relations provesThe deposited flux is zero even though each species occupies only part of the packing: the settling law applies to the suspension, while the deposit is a separate stationary state. These jump relations determine the deposit composition from sedimentation jump conditions.

The large-particle clearing front has height $H-7Wt/8$ and meets the mixed-deposit front at height $35Wt/192$. ThusAt this merger all large particles have entered the deposit. Only the small particles remain suspended, at normalized concentration $1/8$. The new material deposited above the mixed layer is pure small-particle solid, at total concentration ${\varphi }_{\max }$. Its sedimentation shock therefore has downward speedFor $t>t_1$, the deposit top and the remaining clearing front have heightsTheir meeting, or the total solid-volume balance, givesThis is two-stage bidisperse batch sedimentation: the fast species finishes first, followed by deposition of the remaining slow species.

The final solid composition by volume, and also by mass because the particle densities are identical, isThe upper pure-small layer has thickness $9H/116$. As checks, $(29/40)h_1=H/8$ accounts for all large-particle solid volume, and $(11/40)h_1+(h_2-h_1)=H/8$ accounts for all small-particle solid volume, in units of ${\varphi }_{\max }$. If “percentage of particles” means number rather than volume, a large sphere occupies eight times the volume of a small sphere: the mixed bottom layer contains $29/117\simeq 24.8\%$ large particles by number and $88/117\simeq 75.2\%$ small. The upper layer remains entirely small under either convention.

Use the Gaussian units convention implicit in the $4\pi$ coefficient. For a uniform dielectric and a counterion-only electrolyte, the Boltzmann distribution and Poisson-Boltzmann equation areThe reference electric potential is $\varphi (a)=0$, so $n_0=n(a)$; it is not a prescribed nonzero bulk concentration at infinity. A salt reservoir with additional ion species would require their charge densities as well. With $u=\log (r/a)$, the cylindrical Laplacian is $r^{-2}{d}^2/d{u}^2$. Since $\beta q\varphi =\Phi +2u$ and $r^2=a^2{e}^{2u}$, the factors of $e^{2u}$ cancel:Here primes mean differentiation with respect to $u$, and the surface condition is $\Phi (0)=0$ in this coordinate.

Multiply by ${\Phi }^{\prime }$ to obtain a first integral:The two specified far-field limits set $K=0$. For a solution tending upwards to infinity, take the positive square root. Writing $b=\sqrt{C/2}$ gives $(e^{\Phi /2}{)}^{\prime }=b$, hence the cylindrical counterion-only Poisson-Boltzmann profileFor $b>0$, $\Phi \to \infty$ and ${\Phi }^{\prime }=2b/(1+bu)\to 0$. The corresponding number density is

Apply Gauss's law to a coaxial surface per unit cylinder length. Define the positive accumulated cylindrical screening charge byThe signed outward electric field is $E_r=2[-\lambda +Q(r)]/(ϵr)$, so its magnitude is $2|\lambda -Q(r)|/(ϵr)$. In the branch here, $Q(r)<\lambda$, the field is inward andIn this display, ${\varphi }^{\prime }$ is a radial derivative. The quantities $Q$ are electric charges per unit axial length, not total charges of an infinite cylinder.

At the surface, $a{\varphi }^{\prime }(a)=2\lambda /ϵ$. Define the charge-$q$ Bjerrum length ${\ell }_B=\beta q^2/ϵ$ and the Manning parameter $\eta =\lambda {\ell }_B/q$. Matching the derivative of the solution givesThus the condensed branch exists only for $\eta >1$, and its normalization isThe sign condition $b>0$ is essential: simply squaring $\eta -1$ would produce a spurious subcritical density. Direct integration using $rn(r)dr=n_0{a}^{2}du/(1+bu{)}^2$ yieldsEquivalently, $E(r)=2{\lambda }_c[1+b/(1+b\log (r/a))]/(ϵr)$: the remaining far-field line-charge magnitude is ${\lambda }_c$. The Manning condensed fraction is $1-1/\eta$.

Below threshold, no positive-density profile can satisfy those same far-field conditions. Indeed ${\Phi }^{\prime \prime }\le 0$, so a negative initial slope $2(\eta -1)$ cannot increase to zero. The physical interpretation is the infinite-dilution limit of cylindrical counterions: start with a neutral finite cylindrical cell and send its outer radius to infinity. For $\eta \le 1$, counterions escape to arbitrarily large radii, leaving $n_0\to 0$ and no finite condensed charge. The bare-cylinder Boltzmann distribution $n\propto r^{-2\eta }$ also shows why its radial normalization $\int r^{1-2\eta }dr$ diverges in this regime. In the local zero-density limit, $\Phi =2(\eta -1)u$; this does not obey the condensed branch's imposed limits. At $\eta =1$ the branch has zero amplitude and $\Phi =0$, rather than $\Phi \to \infty$.

Consequently the intended counterion condensation result, interpreted as a local infinite-dilution density, isEvery finite neutral cell still contains enough counterions to neutralize the cylinder; the missing residual charge in the local limiting profile is carried to infinity. Taking the infinite-cell limit before the charge integral differs from integrating the entire finite cell first. In SI units, the same reference-charge Bjerrum length is $q^2/(4\pi {ϵ}_{\mathrm{abs}}{k}_{B}T)$; the electrostatic prefactors must be changed consistently.

Take $z$ positive downwards from the roof. Write $T_0=T_m-m{C}_0$, $C_B=(T_m-T_B)/m$, $\Delta T=T_0-T_B=m(C_B-C_0)>0$ and $\Delta C=C_B-C_0$. The cooled magma follows the decreasing liquidus from $(C_0,T_0)$ towards $(C_B,T_B)$; the initially superheated liquid starts at $(C_0,T_{\infty })$ above that liquidus. The solid has zero solute concentration. At each temperature in the mushy layer, a horizontal tie line joins the pure solid to the residual liquid on the liquidus. The lever rule gives the liquid fraction from the bulk concentration. Since $T_B>T_E$, this path stops before the eutectic temperature and does not cross the eutectic composition.

The porosity $\varphi$ denotes the liquid fraction, so $1-\varphi$ is the solid fraction. Equal phase mass densities remove volume-change flow; the stable density stratification suppresses convection. In the stagnant mushy-layer model, define bulk specific enthalpy and bulk solute concentration byThe arbitrary reference of specific enthalpy has been omitted. The two conservation laws in $0<z<h(t)$ areHere $\kappa =k_T/(\rho c_p)$ is the thermal diffusivity, with $k_T$ the thermal conductivity. The positive $L{\varphi }_t$ term records increased specific enthalpy on melting; during freezing it supplies latent heat. In the liquid $z>h(t)$, the heat equation and solute diffusion equation arewhere $D$ is the solutal diffusivity. Neglecting solute diffusion inside the mushy layer does not permit dropping $D{C}_{zz}$ in the liquid at this stage.

Use $T(0,t)=T_B$, liquid $T\to T_{\infty }$ and $C\to C_0$ as $z\to \infty$, and initially $h=0$, $T=T_{\infty }$, $C=C_0$ for $z>0$. These are the half-space boundary conditions appropriate before cooling reaches the bottom of the chamber. For a closed chamber of finite depth $H$, replace the far-field conditions by the specified bottom conditions, for example $T_z=C_z=0$ at an insulated impermeable floor; the half-space similarity solution then ceases to apply when the thermal penetration depth approaches $H$.

At $z=h(t)$, let $C_i$ and $T_i=T_m-m{C}_i$ be the common liquid concentration and temperature, and let ${\varphi }_i$ be the limiting porosity on the mush side. Continuity of $T,C$, the Rankine-Hugoniot condition for solute, and the interfacial conservation of energy giveThe first condition has the sign of solute rejection into the liquid. For the sharp-edge stagnant mushy-layer model with nonzero liquid solutal diffusivity, close the free-boundary problem by marginal equilibrium at a mush-liquid boundary:This prevents constitutional supercooling immediately ahead of the mushy layer. A small solid fraction jump at its edge is allowed: imposing ${\varphi }_i=1$ as well as these finite-$D$ conditions would generally overdetermine the reduced model. At the roof, $C=C_B$ follows from the liquidus, and no additional solute-flux boundary condition is needed in the diffusion-free mush.

To obtain the large-concentration similarity solution, putThe positive superheat specifies the physically intended liquid initial state. On a similarity solution, $C_i,{\varphi }_i$ are constant, so every newly incorporated level leaves the same constant $B={\varphi }_i{C}_i$. Write $C_i=C_0+\Delta C\delta$ and $B=\Delta C(\mathcal{C}+a)$. Then throughout the mushy layerFor $a,\theta =O(1)$, $\mathcal{C}\gg 1$, and $\mathcal{S}/\mathcal{C}=O(1)$, the effective heat capacity is constant to relative error $O({\mathcal{C}}^{-1})$. Thus setIntegration of the similarity solution ordinary differential equation gives the leading analytical profilesThe temperature is $T=T_0+\Delta T\theta$. In the mushy layer, the concentration and porosity are the displayed $C,\varphi$; in the liquid,These retain the liquid solutal boundary layer. For completeness, finite-$r$ values of $\delta ,a,\lambda$ are determined analytically, to the same large-$\mathcal{C}$ order, by settingand solvingThe solute jump gives ${\varphi }_i=1-b/(\mathcal{C}+\delta )$ exactly; replacing $\mathcal{C}$ by $\mathcal{C}+\delta$ on the right of the last equation retains that jump's next-order factor, but does not make the linearized thermal profiles exact at finite $\mathcal{C}$. The branch has $0\le \delta <1$, $0<\varphi \le 1$. All these profiles are leading asymptotic expressions, not exact solutions of the nonlinear finite-$\mathcal{C}$ heat equation.

For $r\to 0$ at fixed positive $\lambda$, the complementary error function asymptotic gives $K(\lambda /\sqrt{r})\sim \sqrt{\pi }\lambda /\sqrt{r}$, hence $\delta =O(r)$, $b=O(r)$, $a\to 0$, and ${\varphi }_i\to 1$. The result simplifies toIn particular, $\varphi (0,t)=C_0/C_B=\mathcal{C}/(\mathcal{C}+1)$ and $\varphi (h^{-},t)=1$. There is no finite interfacial latent heat jump in this limit: the latent heat is released continuously throughout the mushy layer. Matching the two temperature gradients givesThe printed final equation is missing the square on $\Omega$ in the exponential. The corrected exponent follows directly from differentiating $\mathrm{erf}(\Omega \eta )$ and is consistent with ${\Omega }^2=1+\mathcal{S}/\mathcal{C}$. The positive root is unique: the mush-side gradient decreases from infinity to zero with $\lambda$, whereas ${\theta }_{\infty }K(\lambda )$ increases. The zero-superheat limit is singular and is not described by a finite positive $\lambda$ in this two-region construction.

For freezing, the liquid concentration decreases away from the interface, so the liquidus $T_L=-mC$ increases there. Constitutional supercooling means $T_l<T_L$ somewhere ahead of the interface. Since these temperatures agree at the interface, its local onset isFor the similarity solution and salt rejection condition, this is exactlyIt is also the global criterion for this freezing solution when $D<\kappa$ and $T_{\infty }>T_i$: the ratio $T_L^{\prime }(x)/T_l^{\prime }(x)$ is a positive constant times $\exp [-x^2(1/D-1/\kappa )/(4t)]$. For $x\ge a\ge 0$ this ratio decreases with $x$. If the ratio is at most one at the interface, $T_l-T_L$ cannot become negative; if it is greater than one there, a supercooled interval appears immediately ahead. Equality marks onset, not a finite supercooled interval.

For fixed positive $\Delta \theta$ and $\lambda =O(1)>0$, the left side of the inequality is $O(ϵ)$ and the right side is positive of order one, so freezing is constitutionally supercooled for sufficiently small $ϵ$. To locate the boundary between the two regimes, resolve the smaller scale $\Delta \theta =O(ϵ)$, $\lambda =O(ϵ)$. The salt equation and liquidus then giveHere $ϵ\lambda =O({ϵ}^2)$, so the thermal equation improves to $T_i=(T_{\infty }+T_{-\infty })/2+O({ϵ}^2)$; this justifies the second expansion including latent heat. Also $T_{\infty }-T_i={\theta }_{+}+O(ϵ)$ and the left side of the exact criterion is $ϵ{\theta }_{+}/\sqrt{\pi }+O({ϵ}^2)$. Therefore the small-diffusivity constitutional-supercooling threshold isThis printed condition is an asymptotic onset condition, rather than an exact finite-$ϵ$ inequality. In its $O({ϵ}^2)$ transition window the complete algebraic system and exact gradient criterion must be used.