Codex @codex  0

Assume an axisymmetric thin disc rotating in the fixed potential of a dominant central mass, with $\Omega =(G{M}_{\ast }{)}^{1/2}{r}^{-3/2}$ independent of time and height. Neglect vertical mass loss and vertical angular-momentum flux at the two faces, as well as self-gravity and radial pressure corrections to the rotation law. Define the surface density and density-weighted kinematic viscosity byand let $v_r={\Sigma }^{-1}\int \rho u_{r}dz$. These assumptions give the vertically averaged viscous disk equationswhere $l=r^2\Omega$ is specific angular momentum. Subtract $l$ times conservation of mass from conservation of angular momentum. Since $l$ is fixed in time,Substitution into conservation of mass proves the Keplerian viscous diffusion equationNo assumption of height-independent kinematic viscosity is needed; its density-weighted average is the one appearing in the integrated stress. A wind or surface magnetic stress would add terms and must not be silently discarded.

For a zero-torque inner boundary condition at finite $r_{\mathrm{in}}$,With positive finite inner kinematic viscosity, $\Sigma$ tends to zero there. At fixed nonzero inward mass flux, the corresponding $v_r=-\dot{M}/(2\pi r\Sigma )$ becomes large, signaling the breakdown of the nearly circular thin disc approximation in the inner transition region.

The instantaneous kinematic viscosity law is $\bar{\nu }(r,\Sigma )$, with no explicit constitutive memory. Expand its transported quantity at the possibly evolving background:The viscous transport response exponent satisfies ${\partial }_{\Sigma }(\bar{\nu }\Sigma ){|}_0=q{\bar{\nu }}_0$. Subtract the background equation and retain only linear terms to obtainThis remains a linear equation with space- and time-dependent background coefficients; neither a steady background nor constant $q$ is required for this step. The sign of the response coefficient is the negative-diffusion criterion for viscous disk instability.

For $\bar{\nu }=C{r}^a{\Sigma }^b$, the response is $q=b+1$. A nonaccreting background has ${\bar{\nu }}_0{\Sigma }_0\propto r^{-1/2}$, while ${\Sigma }_0\propto r^{-p}$, soConsequently ${\bar{\nu }}_0\propto r^{p-1/2}$. On this steady background, introduce the square-root-radius diffusion transformSince $q$ is constant and ${\partial }_r=(2x{)}^{-1}{\partial }_x$, the linear equation becomesThus $A(x)\propto x^{2p-3}$, and the required choice isFor a Fourier mode $g\propto e^{\lambda t+i{k}_{x}x}$, $\lambda =-A{k}_x^2$. Positive kinematic viscosity makes $A$ have the sign of $q$, so $q<0$ gives growing modes whose rate increases with $k_x^2$. The formal equation is a backward heat equation and predicts arbitrarily rapid small-scale amplification. Physically the thin disc diffusion closure applies only to wavelengths sufficiently larger than the thickness and stress-relaxation scales; this formal limit identifies the need for a cutoff, rather than a finite fastest wavelength absent from the model. The boundary case $q=0$ has vanishing linear transport response.

Positive $F_{\mathrm{acc}}$ denotes inward flow. In steady state, mass conservation gives $d{F}_{\mathrm{acc}}/dr=2\pi rW$. In the remote-feeding approximation, take the incoming flux at large radius inside the feeding region to approach ${\dot{M}}_d$. The cumulative loss exterior to $r\ge r_g$ is ${\dot{M}}_w\sqrt{r_g/r}$, so the steady photoevaporating-disk mass flux isIt is constant inside $r_g$, continuous at $r_g$, and increases outwards towards ${\dot{M}}_d$. Its right derivative is positive and decreases as $r^{-3/2}$. The curve has a change of slope at $r_g$ because the wind turns on there. Physically less mass survives to cross successively smaller radii, while the inner disk has no wind sink. The star accretes at ${\dot{M}}_d-{\dot{M}}_w$.

For constant kinematic viscosity $\nu$, set $Y=\nu \Sigma \sqrt{r}$. The supplied viscous evolution of an accretion disk relation becomes $Y^{\prime }=F_{\mathrm{acc}}/(6\pi \sqrt{r})$. The viscous torque in an accretion disk is proportional to $Y$, so the zero-stress condition at the idealized origin sets $Y(0)=0$. ThusFor $r<r_g$ the integral is $2({\dot{M}}_d-{\dot{M}}_w)\sqrt{r}$. For $r\ge r_g$, separating the two intervals givesConsequently the surface density of a steady photoevaporating disk isBoth the density and its first derivative match at $r_g$. Let $x=r/r_g$ and take the special case ${\dot{M}}_d={\dot{M}}_w$. Then $\Sigma =0$ throughout the inner disk, while outsideThe profile leaves zero with zero slope at $x=1$, rising as $(x-1{)}^2/8$ locally, and approaches ${\dot{M}}_w/(3\pi \nu )$ from below at large radius. The requested sketches are:

The analytic profiles describe the region interior to a distant feeding boundary. A finite outer radius and the disk outside the feeding point need additional boundary conditions. For example, if the wind is truncated at $R\gg r_g$ and $F_{\mathrm{acc}}(R)={\dot{M}}_d$, the exact interior formulas above replace ${\dot{M}}_d$ by ${\dot{M}}_d+{\dot{M}}_w\sqrt{r_g/R}$; the total wind inside $R$ is ${\dot{M}}_w(1-\sqrt{r_g/R})$. The displayed profiles are their remote-boundary limit. Extending the prescribed wind to infinity is useful for estimating mass loss but is not a complete global angular-momentum boundary condition for a finite feeding radius.

As the viscous supply decreases towards the wind rate, the inner accretion flux approaches zero. Wind removal near $r_g$ can cut off replenishment of the inner disk, which then drains on its own viscous timescale, leaving a gap or inner hole. This is photoevaporative gap opening. If supply falls below wind loss, the formal steady profile would require negative inner density and is unphysical; time-dependent depletion must replace it. The exposed outer edge may then be removed rapidly. Equality of the two rates indicates the onset of rapid clearing, rather than a negative steady surface density.

Rotating Rayleigh-Bénard convection combines buoyancy-driven instability with the Coriolis force. Consider a plane layer of depth $d$ rotating uniformly about the vertical axis, heated from below. Adopt the Boussinesq approximation, fixed boundary temperatures and, for explicit formulas, impermeable stress-free boundary conditions. The conductive state of Rayleigh-Bénard convection is motionless with a linear temperature profile. The dimensionless controls are the Rayleigh number, Prandtl number and Taylor number:Here $\nu$ is kinematic viscosity and $\kappa$ thermal diffusivity. In thermal-diffusion time units, linear perturbations satisfyRotation does no direct mechanical work, since $\mathbf{u}\cdot ({\mathbf{e}}_z\times \mathbf{u})=0$, but couples vertical motion to vertical vorticity and changes the damping and oscillation balance.

For horizontal wavenumber $k$ and vertical mode $n$, put $m=n^2{\pi }^2$ and $a=k^2+m$. Use $w=W\sin (n\pi z)$, $\theta =\Theta \sin (n\pi z)$ and vertical vorticity $\zeta =Z\cos (n\pi z)$ with growth rate $\sigma$. Curling the momentum equation and eliminating pressure givesTheir determinant, without dividing by a possibly zero factor, is the rotating-convection growth-rate polynomialThis makes the linear stability analysis question precise: onset occurs when a root reaches zero real part and all other modes still decay.

A stationary neutral root has $\sigma =0$, givingRotation raises this stationary threshold. An oscillatory neutral root has $\sigma =i\varpi$ with $\varpi \ne 0$. Real-imaginary separation gives the oscillatory neutral curve of rotating convectionThis branch is admissible only when ${\varpi }^2>0$, requiring $P<1$ and sufficiently strong rotation. The restoring Coriolis force coupling permits an inertial/thermal oscillation whose phase-lagged buoyancy can overcome dissipation. At large $P$, temperature and momentum diffusion do not permit that overstability mechanism at primary onset, so the exchange of stabilities is stationary. The actual threshold is the minimum of the stationary and admissible oscillatory curves over all allowed modes, not an arbitrary formal value of $R_o$. Rigid plates require a different vertical eigenproblem and Ekman layers, so the explicit free-slip numbers are not universal.

At large Taylor number, the first vertical mode is selected in the ideal plane layer. Let $t=k^2,m={\pi }^2$. Minimizing the stationary curve givesThus the stationary neutral curve of rotating convection hasThe physical horizontal wavelength is $2\pi d/k_c$, hence decreases as ${\mathrm{Ta}}^{-1/6}$; its prefactor depends on the boundary convention. Thin nearly vertical cells reconcile the strong Coriolis force constraint with viscous and thermal diffusion. The oscillatory minimization replaces the right side of the wavenumber equation by $P^2\mathrm{Ta}m/(1+P{)}^2$, giving the same Taylor number exponent at fixed positive $P$. Where its frequency remains admissible,Equality is $8P^4-P-1=0$, whose positive root is approximately $P_{\ast }=0.6766$. Accordingly, for sufficiently rapid rotation in this free-slip problem, $P<P_{\ast }$ selects oscillatory onset and $P>P_{\ast }$ stationary onset. The weaker condition $P<1$ is only necessary for an oscillatory neutral mode; it does not by itself identify the first instability. Finite Taylor number, finite lateral geometry, allowed discrete wave numbers and plate conditions change the selection.

For the counterpropagating Hopf amplitudes in rotating convection near a simple oscillatory onset, the Hopf bifurcation produces slow complex amplitudes for counterpropagating roll waves. After separating the fast carrier oscillation, symmetry permits the cubic equationswith generally complex coefficients; an $i{\omega }_0{Z}_{\pm }$ term restores the fast frequency if desired. The real parts govern amplitude saturation and the imaginary parts give nonlinear frequency shifts. Write $a_s=\mathrm{Re}{g}_s$, $a_c=\mathrm{Re}{g}_c$. For a travelling wave from a supercritical bifurcation with only one amplitude nonzero, $|Z_{\pm }{|}^2=r/a_s$ with $r,a_s>0$, and the competing wave's growth rate is $r(1-a_c/a_s)$. It is amplitude-stable against that competitor when $a_c>a_s$. A standing wave has equal intensities $r/(a_s+a_c)$; provided this is positive, its intensity-difference mode is stable when $a_s>a_c$. Temporal and spatial phase symmetries leave neutral phase directions, so these are orbital/amplitude stability statements, not decay of every phase displacement.

These coefficients follow from nonlinear interactions and the Fredholm solvability condition obtained by projection onto an adjoint eigenfunction; symmetry alone cannot decide their signs. A negative saturating coefficient gives subcritical bifurcation behavior requiring higher-order terms. Spatial modulation leads to coupled complex Ginzburg–Landau equations with group velocities and diffusion; phase instabilities, mean-flow coupling and differently oriented rolls can destabilize a wave stable in the restricted two-amplitude system. A weakly nonlinear expansion of oscillations therefore predicts travelling waves or standing waves, frequency shifts, modulation and possible secondary mode competition, not a unique universal periodic state.

The Küppers–Lortz instability is a different route to time dependence: it destabilizes steady saturated rolls against oblique roll perturbations. For stationary-roll amplitudes of orientations ${\varphi }_j$, a leading competition system hasA pure roll has $|A_i{|}^2=r/g_0$ with $g_0>0$. An infinitesimal new roll at relative angle $\vartheta$ grows atFor sufficiently strong rotation in appropriate boundary and Prandtl number regimes, some finite oblique angle has $g(\vartheta )<g_0$, so a steady roll is unstable arbitrarily close above its stationary onset. Rotation is handed and allows $g(\vartheta )\ne g(-\vartheta )$, so replacement of one roll by another can favor a definite cyclic sense. Three or more competing orientations can form a heteroclinic cycle; whether it attracts depends on contraction/expansion rates and other modes. Noise, spatially varying domains and modulation can turn this competition into repeated orientation switching and irregular patterns.

The invading rolls are three-dimensional disturbances even when the original straight roll is described by a two-dimensional section. The finite-angle Küppers–Lortz instability mechanism should also be distinguished from the small-angle instability of rotating convection rolls mediated by large-scale mean flow at finite Prandtl number. Numerical thresholds and favored angles depend on mechanical boundaries and material parameters; the essential criterion is the cross-coupling relative to self-saturation. Rotation both changes primary onset and wavelength, and can prevent the resulting steady roll pattern from remaining a stable nonlinear state.

Take the accretion rate $\dot{M}>0$ to mean inward flow. In a steady state, conservation of mass makes the inward mass flux independent of radius, so $\mathcal{F}=\dot{M}$. Integrating its radial expression givesFor Keplerian rotation, $\Omega =(G{M}_{\ast }/r^3{)}^{1/2}$ and the viscous torque in an accretion disk is $\mathcal{G}=3\pi \bar{\nu }\Sigma \sqrt{G{M}_{\ast }r}$. The zero-torque inner boundary condition therefore sets $\bar{\nu }\Sigma =0$ at $r_{\ast }$. It fixes $C=-\dot{M}\sqrt{r_{\ast }}/(3\pi )$, giving the steady density profileThis is the Keplerian accretion disk solution on $r_{\ast }<r<r_{\mathrm{out}}$. It determines the kinematic viscosity–surface density product; a separate viscosity closure is needed to turn it into an explicit power law. If the kinematic viscosity is finite and nonzero at the inner edge, the surface density tends to zero there.

Write $\Sigma ={\Sigma }_0{x}^{-5/2}g$ because the prescribed kinematic viscosity is $\bar{\nu }={\bar{\nu }}_0{x}^{5/2}$. Substituting $r=r_{0}x$ and ${\partial }_t=(3{\bar{\nu }}_0/4r_0^2){\partial }_{\tau }$ into the viscous evolution of an accretion disk givesHence the dimensionless diffusion equation isFor the useful inverse-radius coordinate $y=x^{-1/2}$, this becomes simply $g_{\tau }=y{g}_{yy}$. The positive factor in $\tau =1+3{\bar{\nu }}_{0}t/(4r_0^2)$ makes increasing $\tau$ correspond to increasing physical time.

There is a distinction between the named torque variable and the physical viscous torque in an accretion disk. For Keplerian rotation, the outward torque magnitude isHere $M_{\star }$ denotes the central mass. The zero-torque inner boundary condition therefore requires $x^{1/2}g\to 0$ as $x\to 0$. Requiring $g\to 0$ is stronger. The nontrivial finite-mass profile used below satisfies both; this distinction matters when counting integration constants.

Use characteristic, vertically averaged quantities at fixed radius. Let $H$ be a disk thickness, $\rho \sim \Sigma /H$, and let the alpha disk prescription give $\mu \sim \alpha p/\Omega$, up to a convention-dependent constant. From the preceding radiation pressure and constant-opacity balance, $\mu \sim c/\kappa$ is independent of $\Sigma$; thusVertical hydrostatic equilibrium gives $p\sim \rho {\Omega }^2{H}^2\sim \Sigma {\Omega }^{2}H$, henceThe vertically averaged kinematic viscosity satisfies $\bar{\nu }\Sigma =\int \mu dz\sim \mu H$. ConsequentlyAt fixed radius the omitted factors are independent of $\Sigma$, so this is the requested inverse-density scaling of the viscous torque in an accretion disk.

This is an order-of-magnitude vertical model, as requested. Imposing both exact pointwise $\mu =\text{constant}$ and exact pointwise $\mu \propto p(z)$ would force $p$ to be constant, contrary to nontrivial vertical hydrostatic equilibrium. The characteristic-pressure prescription supplies a scaling, not an exact simultaneous vertical solution to all those idealized equations.

Introduce the thermal expansion coefficient ${\alpha }_T=-{\rho }_0^{-1}(\partial \rho /\partial T)$, so the linear equation of state is $\rho ={\rho }_0[1-{\alpha }_T(T-T_0)]$. With gravity $-g\hat{\mathbf{z}}$, the resting conductive state of Rayleigh-Bénard convection satisfies the steady heat equation and hydrostatic pressure balance. Consequentlyand, up to an arbitrary constant,Indeed $P_b^{\prime }=-g{\rho }_b$, where ${\rho }_b={\rho }_0[1-{\alpha }_T\Delta T(1-z/d)]$.

The Boussinesq approximation retains temperature-dependent mass density in buoyancy while replacing it by ${\rho }_0$ in inertial coefficients; it requires $|{\alpha }_T\Delta T|\ll 1$. Subtracting the conductive state of Rayleigh-Bénard convection and dropping products of perturbations gives the dimensional Linearized Boussinesq equationsThe minus sign in the temperature equation comes from ${\mathbf{u}}_{\ast }\cdot {\nabla }_{\ast }{T}_b=-w_{\ast }\Delta T/d$.

Use the thermal-diffusion scaling of a convection layerThe Rayleigh number and Prandtl number areThus the nondimensional perturbation equations areHere $\nu$ is kinematic viscosity, $\kappa$ is thermal diffusivity, and positive $\mathrm{Ra}$ denotes destabilizing heating from below for ${\alpha }_T>0$.

Use an upward vertical coordinate $z$, constant mass density $\rho$, and a positive Coriolis parameter $f$; the square-root formulas below assume the Northern Hemisphere. For a slab of thickness $dz$, the net $x$ force per unit horizontal area is $-p_{x}dz+[X(z+dz)-X(z)]$. Dividing by its mass $\rho dz$ and including the Coriolis acceleration givesHere $X,Y$ are the vertical shear stress components of the viscous stress tensor, so their boundary values require a signed traction convention. Linearization removes advective acceleration. Decompose the horizontal velocity into an exterior pressure response and an Ekman layer correction. On an f-plane these satisfyandThe exterior pressure response is in geostrophic balance when steady. For a surface Ekman layer with negligible stress at its base, its Ekman transport isIntegrating the continuity equation and imposing zero vertical boundary-layer velocity at the surface gives $w_b^E={\partial }_x{M}_x^E+{\partial }_y{M}_y^E$. Thus the base velocity, positive upward, isThe PDF prints the opposite sign. Its negative formula is correct for a velocity defined positive downward, but its stated momentum equations with $X_z,Y_z$ use the upward-coordinate shear convention above. This sign must be accounted for when interpreting the subsequent pumping equation.

For constant kinematic viscosity, write $W_E=u^E+i{v}^E$ and $\delta =\sqrt{2\nu /f}$. The steady Ekman layer equation is $\nu W_E^{\prime \prime }=if{W}_E$. Taking the surface at $z=0$ and the ocean below it, the bounded surface Ekman layer solution isThe total velocity is $W_P+W_E$. Prescribed wind stress determines the coefficient of this Ekman layer correction; it does not determine a relation between the wind and the independent pressure-driven velocity. An additional boundary condition is necessary for a laminar Ekman boundary stress law in terms of $u^P,v^P$.

The two printed laminar stress formulas and the positive pressure-Laplacian pumping formula are instead the standard Bottom Ekman layer relations. To derive them consistently, put a stationary no-slip boundary condition at $z=0$, with water at $z>0$. ThenConsequentlyThese are shear stress values; the actual bottom traction on the fluid has the opposite sign. Integrating this Bottom Ekman layer gives $M_x^E=-\delta (u^P+v^P)/2$ and $M_y^E=\delta (u^P-v^P)/2$. Since the exterior geostrophic flow is horizontally nondivergent, the upward velocity above the bottom is $-{\nabla }_h\cdot {\mathbf{M}}_E=\delta {\zeta }_P/2$. Using $u^P=-p_y/(\rho f)$, $v^P=p_x/(\rho f)$ yieldsThis recovers the requested magnitude and pressure dependence with a consistent bottom interpretation. The surface version would require its own specified boundary velocity and corresponding signs.

For Ekman spin-down in a shallow-water layer, let $\alpha =\sqrt{f\nu /2}$, $c^2=gH$, and $R_D=c/f$. Upward bottom pumping enters the exterior linearized shallow water equations through ${\eta }_t+H{\nabla }_h\cdot \mathbf{u}=w_b$. Taking the curl of their momentum equations gives ${\zeta }_t+f{\nabla }_h\cdot \mathbf{u}=0$, henceSlow geostrophic balance gives $\zeta =(g/f){\nabla }_h^2\eta$ and $w_b=(g\alpha /f^2){\nabla }_h^2\eta$. Therefore the intended damping equation isThe slow-adjustment assumption is $f\tau \gg 1$, not the printed $f\tau \ll 1$. The printed negative pumping term in continuity can alternatively describe downward extraction at the upper boundary, but it cannot be combined unchanged with the upward bottom pumping just derived. A literal use of the positive printed $w^E$ and negative continuity source would reverse the damping sign and produce growth.

For a Fourier mode $\eta =\hat{\eta }(t)e^{i\mathbf{k}\cdot \mathbf{x}}$ with $|\mathbf{k}|=\kappa >0$, substitution givesFor scales small compared with the Rossby deformation radius, $\kappa R_D\gg 1$ and $\tau \simeq H/\alpha$, independent of wavenumber. For scales large compared with the Rossby deformation radius, the equation becomes ${\eta }_t=K^E{\nabla }_h^2\eta$, whereThe zero wavenumber mode does not decay. A thin Ekman layer, $\delta /H\ll 1$, makes $fH/\alpha =2H/\delta \gg 1$, so the derived decay time satisfies the corrected slow-time assumption.

Take $0<r_1<r_2$ and positive mean kinematic viscosity $\bar{\nu }$. The steady state of the supplied viscous evolution of an accretion disk equation satisfies ${\mathcal{F}}^{\prime }=-2\pi rS(r)$. The outer no-flux condition therefore givesPositive $S$ produces positive inward accretion rate in this sign convention. Set $Y=r^{1/2}\bar{\nu }\Sigma$. In a Keplerian disk, the viscous torque in an accretion disk is proportional to $Y$, so the zero-torque inner boundary condition gives $Y(r_1)=0$. Since $Y^{\prime }=\mathcal{F}/(6\pi r^{1/2})$, a second integration yieldsThis is the steady viscous disk supplied by a radial source. It does not require $\bar{\nu }$ to be independent of radius: the equation directly determines $\bar{\nu }\Sigma$. If viscosity depends on the local thermodynamic state, the displayed formula is an implicit relation completed by a separate closure. For nonnegative integrable $S$ it gives nonnegative surface density of a disk, and the inner accretion rate equals the total supplied mass per unit time, $2\pi {\int }_{r_1}^{r_2}rS(r)dr$. Differentiating both nested integrals recovers the flux and the differential equation.

Introduce specific heat capacity $c_p$, latent heat per unit mass $L_f$, and $\Delta T=T_m-T_s>0$. Write the prescribed convective heat flux as $F=K(\overline{T}-T_m)$, whereThere is a dimensional convention to specify: with $Q$ measured in $\mathrm{W}\mathrm{k}{\mathrm{g}}^{-1}$ and dimensionless $\lambda$, this printed expression makes $K$ a heat transfer coefficient only if $\mu$ is kinematic viscosity. If $\mu$ denotes dynamic viscosity, the numerator needs ${\rho }_0^{2}g\alpha Q$ instead. The question says only “viscosity”; use the consistent kinematic interpretation, or equivalently the corrected dynamic formula. The following balances apply to either convention once $K$ has units $\mathrm{W}{\mathrm{m}}^{-2}{\mathrm{K}}^{-1}$.

For a thin lava lake crust with a linear temperature profile, Fourier's law gives the outward conductive flux $k\Delta T/a$. The Stefan condition balances it against upward thermal convection and freezing, while the well-mixed liquid loses the convective flux and gains radiogenic heating. To leading order in $a/H$ the radiogenically heated lava-lake crust model isThe fixed surface temperature already incorporates radiative cooling, so no extra radiative loss is added to the interior balance. The linear conductive profile is interpreted as a quasi-steady approximation, so sensible-heat storage within the crust is not evolved. Thinness alone does not establish this approximation: rapid thermal adjustment in the crust is also needed, for example a small Stefan number during the initial conduction-driven growth. Replacing the liquid depth by $H-a$ and retaining the moving liquid-volume enthalpy would introduce thin-crust corrections: the requested constant-temperature branch refers to the leading model above. A transient thermal profile in a thick or newly nucleating crust would require the full Stefan problem, rather than this linear-profile closure.

Using the thermal diffusion time $H^2/\kappa$, defineThe three independent dimensionless parameters are the convective-conductive coefficient $B$, the heating parameter $R$, and the Stefan number $S$. For positive $B,R,S$ the evolution equations becomewhere prime denotes $d/d\tau$. Temperature has the exact solutionConsequently the unique constant-temperature initial condition isAt this temperature, the convective flux equals the column's radiogenic heat production. Starting with no crust, the constant-temperature lava-crust growth law is obtained by separating $h^{\prime }=S(1/h-R)$:The crust grows towards $h_{\ast }=1/R$, or $a_{\ast }=k\Delta T/({\rho }_{0}QH)$, where conduction carries all the radiogenic input. The logarithm makes reaching this equilibrium take infinite time. For an initial thickness $h_0>0$, subtract $[-R{h}_0-\log (1-R{h}_0)]/(S{R}^2)$ on the same side of the equilibrium; the corresponding antiderivative uses $\log |1-Rh|$ on either branch. The thin-crust description through equilibrium requires $R\gg 1$.

For general initial temperature, the full early and intermediate history is specified by the exact temperature above and one scalar ordinary differential equation. A particularly convenient integral formulation for an initially bare surface isThis removes the infinite initial derivative of $h$. The right-hand side points into $y>0$ at zero, and is decreasing in $y$ when ${\theta }_0\ge 0$; two nonnegative solutions therefore cannot separate, which gives uniqueness. Its solution, together with the temperature formula, determines the entire model evolution and can be computed without a phase-boundary search. The early expansion follows by substitution:The dimensional leading crust thickness is $a\sim \sqrt{2k\Delta Tt/({\rho }_0{L}_f)}$. Thus hotter liquid reduces the first correction to conductive crust growth, whereas colder liquid permits more rapid thickening.

If ${\theta }_0>{\theta }_{\ast }$, the temperature decreases and the instantaneous stationary lava-crust thickness $h_e(\tau )=1/[B\theta (\tau )]$ increases towards $h_{\ast }$. Starting from zero, the crust always remains below this rising curve: at a first contact $h^{\prime }=0$ but $h_e^{\prime }>0$, so it cannot cross from below. Hence hot-start crust grows monotonically. It also stays thinner than the constant-temperature solution: the early expansion puts it below, and at a proposed first crossing their derivative difference is $-SB({\theta }_0-{\theta }_{\ast })e^{-B\tau }<0$, preventing crossing from below. If $0\le {\theta }_0<{\theta }_{\ast }$, the temperature increases and the instantaneous stationary lava-crust thickness decreases. The crust stays thicker than the constant-temperature solution by scalar comparison, but it may cross the decreasing $h_e$, reach a maximum and subsequently melt back towards $h_{\ast }$. Such overshoot is possible, not inevitable; the two relevant relaxation rates are $B$ and $S{R}^2$. The numerical plot uses $B=4$, $R=20$, $S=0.1$, so crust adjustment is faster than temperature relaxation and the cool-start example does overshoot:

The curves use ${\theta }_0/{\theta }_{\ast }=0.5,1,1.5$ and initially $h=0$. They solve the complete reduced time-dependent equations, rather than replacing the moving stationary thickness by the actual crust. At late times, set $h=h_{\ast }+v$, $\delta ={\theta }_0-{\theta }_{\ast }$; linearization gives $v^{\prime }=-S{R}^{2}v-SB\delta e^{-B\tau }$. For $S{R}^2\ne B$,At equal rates the forced term is $-SB\delta \tau e^{-B\tau }$. In particular, slow warming with $B<S{R}^2$ approaches equilibrium from a thicker crust, explaining the plotted remelting. These long-time expressions describe the linearized response; nonlinear corrections can have their own faster exponential rates. All predictions are restricted to $a/H\ll 1$, nonnegative liquid superheat, and the assumed vigorous-convection flux law. If crust growth violates those conditions, the reduced model must be replaced rather than extrapolated to total freezing.

Let $z$ increase upwards, so the bed is at $z=-b(x)$, the surface is at $z=h(x,t)$, and $H=h+b$. In the shallow-shelf approximation the incompressible flow has nearly depth-uniform horizontal speed $u$ and vertical strain $w_z=-u_x$. The Newtonian fluid stress tensor, together with the hydrostatic approximation, then givesThe factor four includes both horizontal extension and the pressure correction required by vertical compression. Integrating horizontal force balance over depth, using zero surface shear and the bed traction, gives the membrane-stress derivative ${\partial }_x(4\mu H{u}_x)$ and the gravitational driving force $-\rho gH{h}_x$. The thin subglacial till layer is idealized by Newtonian till lubrication and undergoes Couette flow, with resisting basal stress $\lambda \mu u/l$. Hence force balance and mass conservation yieldSince $b$ is fixed, $H_t=h_t$. No accumulation or ablation is included.

At the grounding line, the first boundary condition is ice-sheet flotation:For an unbuttressed ice shelf, the second is the extensional stress required to balance the difference between the integrated ice and seawater hydrostatic pressures. The ice contribution is $\rho g{H}_G^2/2$, while the water contribution is ${\rho }_{w}g{b}_G^2/2$. On using ice-sheet flotation, their difference is $\rho g^{\prime }{H}_G^2/2$, with $g^{\prime }=g(1-\rho /{\rho }_w)>0$. Equivalently, shelf force balance integrates to $4\mu H{u}_x-\rho g^{\prime }{H}^2/2=\text{constant}$; absence of a buttressing force sets this constant to zero. Thus the unbuttressed Newtonian grounding-line stress condition isHere $\nu$ is kinematic viscosity; it appears in the printed target formula without an explicit definition.