Adding the ordinary differential equations for all enzyme states makes the net rates telescope:With initially unbound enzyme, this sum is $E_0$ for all time. HenceThis is enzyme conservation, not conservation of each complex individually. The reaction changes site occupancy and makes reaction product, but no event creates or destroys an enzyme molecule.

Use $S,C_i,P$ for concentrations and define net complex-formation rates $R_i=k_i{C}_{i-1}S-(k_{-i}+k_i^{+})C_i$, with $R_{n+1}=0$. The law of mass action givesEach binding event removes one free substrate molecule; each unbinding returns one; each catalytic event releases one reaction product and reduces the number of occupied sites by one. Suitable initially unbound enzyme and substrate data are $S(0)=S_0$, $P(0)=0$, $C_0(0)=E_0$, $C_i(0)=0$ for $i\ge 1$. All rates are nonnegative on the boundary of the nonnegative concentration region in the direction required to preserve it. These equations also conserve $S+P+{\sum }_{i=1}^{n}i{C}_i=S_0$, accounting for bound substrate as well as free substrate.

Write $A(x)=x-x^2$, $q=1/(n+2)$ and $p=(n+1)/(n+2)$. Integrating the steady mass conservation equation from the head givesFor the time-dependent first-order equation $H_t+H^{n+1}{H}_x=1-2x$, method of characteristics gives $\dot{x}=H^{n+1}$ and $\dot{H}=1-2x$. Consequently $d[H^{n+2}/(n+2)-A(x)]/dt=0$. On a characteristic starting at $(\xi ,0)$ with zero thickness,On a characteristic emitted by the boundary $x=0$, the invariant is zero and $H=H_0(x)$. The travel time of the one emitted at time zero isThe incomplete beta function is finite at both endpoints because $0<p<1$ for the physical exponents $n>0$. Before characteristic crossings, the initial-data branch and boundary-data branch have the displayed forms in the source, with the appropriate emission time $t-f(x)$ for the latter.

Glacier characteristic admissibility gives a substantive domain restriction. Factoring the invariant gives $A(x)-A(\xi )=(x-\xi )(1-x-\xi )$. A positive forward branch requiresEquivalently it is the image of this parameter domain under $(\xi ,x)\mapsto (x,t(\xi ,x))$, restricted further to regions where that map is single-valued. Characteristics cannot start with positive thickness growth at $\xi >1/2$, where the prescribed initial rate $H_t=1-2\xi$ is negative. This is why the printed piecewise expression is a formal characteristic construction, not a global nonnegative solution from a bare bed. For example for $x>1/2$, putting $L_{\xi }=1-2\xi$ givesBoth factors increase with $\xi$ there, so no positive initial-data branch exists at that $x$ for $0<t<f(x)$, contrary to a literal global reading of the first case. Compression can also make characteristics intersect; entropy solutions then require shock waves.

One can nevertheless give the exact zero-thickness endpoint locus of these formal characteristics. Let $K=f(1)=(n+2{)}^{-p}B(q,q)$. At $x=1-\xi$ the change of variable $z=\xi +(1-2\xi )y$ givesThis is the formal snout locus inferred from the given characteristic family; it retreats from $1$ towards $1/2$. It is not justified as the physical terminus of an accumulating glacier formed from zero initial ice. In particular it says nothing about a positive advancing snout at $t<K$. A physical free boundary problem must specify no ablation on a bare bed, maintain $H\ge 0$, and select an entropy solution or retain the small slope regularization. Those additional conditions are absent from the printed first-order construction. Thus a unique global physical snout trajectory cannot be claimed from that construction alone. The steady endpoint remains $x=1$. The degenerate $n=0$ limit would instead have the same endpoint time $K$ for every initial characteristic, so the inverted formula presumes $n>0$.

A typical net ice accumulation curve decreases downstream: snow gain at high elevation, zero at the equilibrium line, and negative values in the ice ablation zone. In a steady glacier with zero volume flux at its head, $Q(x)={\int }_0^{x}a(X)dX$. If its terminus also has zero volume flux, a natural glacier length $L$ is the positive root of ${\int }_0^{L}a(X)dX=0$, with the integral positive between the endpoints. Such a root requires a suitable given accumulation profile; it is not guaranteed by dimensional analysis alone.

Let $a_0$ be a representative accumulation rate and $K=A(\rho g\sin \alpha {)}^n$. ChooseUsing $x=L\hat{x}$, $H=H_{\ast }\hat{H}$, $Q=Q_{\ast }\hat{Q}$, $t=t_{\ast }\hat{t}$, $a=a_0\hat{a}$, both terms in mass conservation have scale $a_0$. The dimensionless flux isDropping hats gives the required partial differential equation. The parameter measures the surface-thickness slope correction relative to the bed slope. For the linear example $a=1-2x$, accumulation changes sign at $x=1/2$ and its integrated balance vanishes at $x=1$.

Set $G=\rho g(\sin \alpha -\cos \alpha H_x)$. For downstream shear $G\ge 0$, Glen's flow law and the no-slip boundary condition giveHere $A$ is the ice rate factor, distinct from the cross-sectional area in the river problem. Integrating this velocity profile through the thickness gives the volume flux per unit valley width:The factor $H^{n+2}$ combines the depth-dependent shear stress, integration to obtain velocity, and integration to obtain flux. If reverse shear is admitted for noninteger $n$, the signed version of Glen's flow law is $u_z=A|\tau {|}^{n-1}\tau$; the displayed unsigned power presumes the downstream regime.

Measure pressure relative to atmospheric pressure. At leading order the stress-free upper surface gives $p(H)=0$, $\tau (H)=0$. Integrating the normal balance first yieldsAt fixed $z$, $p_x=\rho g\cos \alpha H_x$. Therefore ${\tau }_z=\rho g(\cos \alpha H_x-\sin \alpha )$, independent of $z$. Integrating from $z$ to $H$ givesThis is linear in $z$, is zero at the surface, and reaches its basal value at $z=0$. In particular the surface slope modifies the shear driving; a flat upper surface relative to the horizontal eliminates it in this approximation.

Resolve gravity into components $g\sin \alpha$ down the bed and $-g\cos \alpha$ normal to it. For slow incompressible flow, inertia is negligible and Cauchy momentum equation reduces to the divergence of the stress tensor plus the gravitational force density being zero. Write the dominant shear component as $\tau ={\sigma }_{xz}$ and the leading normal stresses as $-p$.

The shallow-ice approximation assumes a small thickness-to-length ratio. Vertical variation of the shear is retained in the along-bed equation, while longitudinal deviatoric-stress gradients and the shear contribution to the normal equation are higher order. The leading balances are consequentlyThe second is a hydrostatic pressure balance; the first balances gravitational driving and the surface-induced pressure gradient against vertical shear. This approximation requires shallow geometry and sufficiently slow flow, not a Newtonian constitutive equation for ice.

The water mass conservation equation inherited from the earlier model is $ϵh_t+(hu{)}_x=0$. It is essential: without it the printed displayed equations are not closed. Take the slow bed branch $ϵ=0$, and perturb the uniform state by $h=1+\hat{h}{e}^{\lambda t+ikx}$, $u=1+U{e}^{\lambda t+ikx}$, $c=1+C{e}^{\lambda t+ikx}$, $s=S{e}^{\lambda t+ikx}$ and $\eta =1+\hat{\eta }{e}^{\lambda t+ikx}$. For $k\ne 0$, mass conservation gives $\hat{h}=-U$ and $\hat{\eta }=S-U$.

Put $m=E^{\ast \prime }(1)+1$; the extra $1$ comes from linearizing $c/u$. The linearized momentum and suspended-sediment equations areEliminating $U,C$ gives the dispersion relationAway from $F=1$, for fixed nonzero $k$ and negligible $\delta$,Thus growth requires $m/(F^2-1)>0$. For $F>1$, any $E^{\ast \prime }(1)>-1$ gives an unstable branch; for $F<1$ growth would require $E^{\ast \prime }(1)<-1$. A physically increasing erosion function therefore gives the supercritical instability in this simplified model, rather than proving growth of ordinary subcritical dunes. Retaining drag giveswhich shows its stabilizing correction and the need to reconsider the limit near critical flow.

The negative bed-to-surface ratio for subcritical flow means bed crests lie beneath surface troughs: the dune phase relation. Its positive sign for supercritical flow gives aligned crests: an antidune. In the unstable supercritical case the imaginary part gives upstream migration, since $-\mathrm{Im}\lambda /k=-m/[(F^2-1)(1+k^2)]$. The figure shows these two phase relations. They identify the bedforms but do not by themselves establish growth. The qualitative limitation of the suspended-sediment mechanism is also discussed in the official ^[ref].

Let the original area be $A_0>0$ and the increased upstream area be $A_1>A_0$. At $ϵ=0$, the conservation law flux $q(A)=A^{5/4}$ is strictly convex. Characteristics behind the increase travel faster than those ahead, forming an admissible shock wave. Rankine-Hugoniot conditions give its speedAt $x=1$ the idealized hydrograph jumps from $Q_0=A_0^{5/4}$ to $Q_1=A_1^{5/4}$ at $t=1/c_s$. With positive $ϵ$, the diffusion term spreads the rise over a finite interval. The hydrograph is rounded, with an earlier small response and a later approach to the new discharge, rather than an instantaneous jump.

The physical cause is the longitudinal pressure gradient: a local increase in depth influences neighbouring sections and redistributes water before pure downstream transport would deliver a discontinuity. For a small increase, linearization about $A_0$ has speed $c=\frac{5}{4}{A}_0^{1/4}$ and diffusivity $D=ϵA_0^{3/4}/4$. The smooth curve in the figure is the exact half-line step response for this linearized advection-diffusion equation; it illustrates the shape without claiming to solve a large nonlinear jump exactly.

First neglect inertia relative to the $O(ϵ)$ pressure correction, requiring $F\ll 1$ for smooth fields on these scales. The positive velocity isMultiplying by area gives the discharge $Q=A^{5/4}-(ϵ/4)A^{3/4}{A}_x+O({ϵ}^2)$. HenceThe leading equation is a nonlinear kinematic wave equation. The pressure correction has positive diffusion coefficient for $A>0$ and smooths steep gradients. Dropped terms are $O({ϵ}^2)$ and inertial terms $O(ϵF)$ on the assumed scales; sufficiently thin fronts require a separate balance rather than an unrestricted Taylor expansion.

Write $R=k_R\sqrt{A}$, $\bar{h}=k_H\sqrt{A}$, where the given geometry givesTake a Chézy drag law $\tau =\rho c_{d}u|u|$ and restrict to downstream flow $u>0$. Choose an area scale $A_0$, length scale $L$, and speed scale $U_0$ satisfyingThis makes gravitational driving and drag equal in the reference uniform steady flow. Dividing the printed momentum equation by $gS$ givesMass conservation retains $A_t+(Au{)}_x=0$. In this notation $F$ is the square of the velocity-to-gravity-wave-speed ratio based on the chosen depth scale; authors sometimes instead denote its square root by Froude number. If the actual triangular pressure head from the preceding solution is used, replace $k_H$ by $2k_H$ in both parameter definitions. The displayed dimensionless form remains the same, with consistently redefined parameters.

Assume constant mass density, a prismatic channel, approximately uniform cross-sectional velocity $u$, a hydrostatic pressure distribution, small bed slope $S$, and no lateral inflow. Integrating mass conservation over a slice of length $dx$ gives $\rho A_{t}dx=\rho Au{|}_x-\rho Au{|}_{x+dx}$, hence $A_t+(Au{)}_x=0$.

Let $I(A)={\int }_A(H-z)dA$ be the hydrostatic first moment, where $H$ is the actual free-surface depth. The downstream momentum balance isUsing $A_t+(Au{)}_x=0$ reduces its left side to $A(u_t+u{u}_x)+g{I}_x$. For a fixed prismatic cross-section, $dI/dH=A$, so $I_x/A=H_x$. Therefore the exact depth-averaged result isHere hydraulic radius converts perimeter shear stress into force per unit fluid volume. The printed notation $\bar{h}$ yields the requested equation if it denotes this hydraulic pressure head, whose derivative equals $H_x$.

The geometrical definition printed in the next part needs care. Triangular-channel hydrostatic pressure head gives the factor explicitly. For the stated symmetric triangular notch, the actual depth is $H=(\ell /2)\sin \theta$, whereas its printed $\bar{h}=(\ell /4)\sin \theta$ is $H/2$. Indeed $A=H^2\cot \theta$ and $I=AH/3$, giving $I_x/A=H_x=2{\bar{h}}_x$. Thus the equation with precisely the printed geometrical $\bar{h}$ has pressure term $-2g{\bar{h}}_x$. A first-principles derivation cannot remove that factor. The subsequent intended scaling can be carried out algebraically using the printed convention, as below, or with the pressure scale doubled for the actual triangular channel.

Initially $a=1/2$ and $T=T_{+}/T_0>1$, so $\epsilon T_t=1-T^4<0$. On the fast time $s=t/\epsilon$, albedo changes by only $O(\epsilon )$ andThus temperature rapidly falls towards $1$. At the same time $a_t=a^{\ast }(T)-1/2<0$ while $T>1$, which seeds selection of the warm branch. Afterwards the approximate slow manifold is $T^4=2(1-a)$: decreasing albedo makes temperature rise slowly back towards $T_{+}/T_0$. The preceding figure separates these two time scales; a single time axis would compress the initial fall into a nearly vertical segment.

In the sharp-transition regime, after selection of the warm branch $a^{\ast }(T)\simeq 1/4$. Solving the linear ordinary differential equation with the leading matched value $a(0^{+})=1/2$ givesThese are outer solutions: their initial temperature is $1$, so the fast cooling layer is required to satisfy the actual initial temperature.

There is a useful qualification to the two-limit argument. For fixed finite $\mu >2$, put $x=1-2a$ on the slow manifold. Then $x_t=-x+\frac{1}{2}\tanh (\mu x)$, and $x=0$ is an unstable equilibrium with growth rate $\mu /2-1$. The small positive $O(\epsilon )$ seed from the fast layer can therefore produce a delay of order $\log (1/\epsilon )/(\mu /2-1)$. The simple exponential approximation assumes that the sharp albedo transition is sufficiently fast to make this delay negligible on the displayed slow scale. It is not a uniform approximation for every fixed $\mu >2$ as $\epsilon \to 0$.

Use dimensional quantities $T_d,t_d$ and set $T_d=T_{0}T$, $t_d=t_{0}t$. Dividing the energy balance by $\sigma \gamma T_0^4$ and the albedo equation by its time scale givesThus $\epsilon$ is the ratio of the atmospheric thermal response time to the much longer albedo response time. Converting kilometres to metres and years to seconds gives $C=10^7\mathrm{J}{\mathrm{m}}^{-2}{\mathrm{K}}^{-1}$, $T_0\simeq 266.3\mathrm{K}$, $t_0\simeq 3.156\times 10^{10}\mathrm{s}$, andThe thermal time is about half a year, compared with a millennium for albedo adjustment. This is a singular perturbation: setting $\epsilon =0$ removes the independent initial condition for temperature.

Set $x=T^4/T_0^4-1$. Since $\sigma \gamma T_0^4=Q/8$, the steady state energy balance is $2(1-a_0)=1+x$. Substitution givesThe two sides have slopes $\mu$ and $2$ at the origin. For $x>0$, $\tanh (\mu x)$ is strictly concave, its initial slope is $\mu$, and it remains below $1$. Consequently no positive root exists for $\mu \le 2$, while for $\mu >2$ there is exactly one positive root $x_{+}\in (0,1/2)$. Oddness supplies the corresponding negative root. At $\mu =2$ the only root is still zero. Hence multiple steady temperatures occur precisely when $\mu >{\mu }_c=2$.

The largest root corresponds monotonically to the largest positive absolute temperature:For large $\mu$, the positive root tends to $1/2$, so $a_{+}\simeq 1/4$ and $T_{+}/T_0\simeq (3/2{)}^{1/4}$. The negative branch has larger albedo and lower temperature.

The heat capacity per unit horizontal area is $C=\rho c_{p}d$: mass density times specific heat capacity times atmospheric depth. Thus $C\dot{T}$ is the rate of stored thermal energy per unit area. The solar constant $Q$ is the incident radiative flux on a surface normal to sunlight; the ratio of the Earth's intercepting disk area to its spherical area is $1/4$. The absorbed fraction is $1-a$, where albedo $a$ is the reflected fraction. The outgoing radiative flux is $\sigma \gamma T^4$, a Stefan–Boltzmann law with effective radiation factor $\gamma$ incorporating the simplified greenhouse effect. Every term therefore has units of power per area.

The albedo curve decreases smoothly, passes through $(T_0,1/2)$, and approaches $1/4$ at high temperature. At $T=0$ it is $(2+\tanh \mu )/4$, close to $3/4$ when $\mu$ is large. Its derivative isAn ice-covered planet reflects more sunlight than exposed land or ocean. Warming removes reflective ice, lowers albedo, and increases absorbed radiative flux: this is ice-albedo feedback. The parameter $\mu$ controls the sharpness of this transition; the formula is an idealized global response, rather than a local melting law.