The Longevity of Water Ice on Ganymedes and Europas around Migrated Giant Planets

In the blackbody approximation, radiative cooling is given by $\sigma {T}^{4}$ at isothermal temperature T, allowing a straightforward formulation of escape versus radiative cooling. As such, the first order global energy balance for an icy exomoon is between incoming stellar flux versus the energy flux lost to vaporizing the water, lifting molecules out of the gravity well, and radiative cooling, i.e.,

$$\begin{eqnarray}&&\mathop{\underbrace{\displaystyle \frac{1}{4}(1-A){F}_{s}}}\limits_{\mathrm{absorbed}\,\mathrm{stellar}\,\mathrm{flux}}=\mathop{\underbrace{\left(\displaystyle \frac{{GM}}{{r}_{s}}+{L}_{v}\right){\rho }_{s}{u}_{s}}}\limits_{\mathrm{mass}\,\mathrm{loss}\,\mathrm{flux}}+\mathop{\underbrace{\sigma {T}^{4}}}\limits_{\mathrm{radiative}\,\mathrm{cooling}}.\end{eqnarray} \tag{ 2 }$$

Here, F_s is the incoming solar flux available for absorption, A is the Bond albedo, σ is the Stefan–Boltzmann constant, r_s is the surface radius where the atmospheric density is ${\rho }_{s}$, such that ${\text{}}{\rho }_{s}{u}_{s}$ is the mass flux given an outward radial flow velocity at the surface u_s, G is the gravitational constant, and M is the mass of the exomoon.

In addition to the C–C relationship (Equation (1)), three equations are needed to derive the steady state, hydrodynamic atmospheric loss in the isothermal approximation (e.g., Catling & Kasting 2017, Ch. 5). The first is steady state mass continuity, given by

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial r}({r}^{2}\rho u)=0\end{eqnarray} \tag{ 3 }$$

where r is the radial distance from the planet's center, u is the outward radial flow velocity, and ρ is the atmospheric density. Steady state momentum conservation is expressed as

$$\begin{eqnarray}&&u\displaystyle \frac{\partial u}{\partial r}+\displaystyle \frac{1}{\rho }\displaystyle \frac{\partial p}{\partial r}=g\end{eqnarray} \tag{ 4 }$$

with gravity $g=-{GM}/{r}^{2}$ and pressure p. Finally, the equation for energy balance is given by Equation (2). Combining Equations (1)–(4) an analytic expression for the isothermal atmospheric mass loss rate in kg ${{\rm{s}}}^{-1}$ is given by (see Appendix A for the derivation):

$$\begin{eqnarray}&&\dot{M}=\pi {\rho }_{s}\displaystyle \frac{{G}^{2}}{{u}_{0}^{3}}{M}^{2}\exp \left[\displaystyle \frac{3}{2}-\displaystyle \frac{G}{{u}_{0}^{2}}{\left(\displaystyle \frac{3}{4\pi {\rho }_{m}}\right)}^{-\tfrac{1}{3}}{M}^{\tfrac{2}{3}}\right]\end{eqnarray} \tag{ 5 }$$

where u₀ is the isothermal sound speed given by ${u}_{0}^{2}={kT}/m$ with Boltzmann constant k and mean molecular weight m, and ${\rho }_{m}$ is the mean density of the exomoon (assumed $2\,{\rm{g}}\,{\mathrm{cm}}^{-3}$). The atmospheric surface density, ${\rho }_{s}$, is set the by C–C equation for the saturation vapor pressure of water at the prevailing temperature.

For time averaged mass loss rate, $\bar{M}$, the lifetime of the exomoon surface water is given by

$$\begin{eqnarray}&&{\tau }_{\mathrm{Water}}=\displaystyle \frac{{M}_{\mathrm{Water}}}{\bar{M}}.\end{eqnarray} \tag{ 6 }$$

For an upper limit, we assume the total mass of water present on the exomoon surface, M_Water, is 40% of the bulk mass. However, even if 5% water were used (the lower limit for Europa (Schubert et al. 2004)), from Equation (6) we can see that it would translate to a change in ${\tau }_{\mathrm{Water}}$ by a factor of 8, compared to 40% water. From Equation (5) we see that $\dot{M}$, and hence ${\tau }_{\mathrm{Water}}$, has an exponential dependence on mass, so we would anticipate that the difference between 5% and 40% water is not the major factor determining ${\tau }_{\mathrm{Water}}$, which is borne out by our results. In addition, if substantial water vapor is lost the bulk density of the moon, ${\rho }_{m}$ may increase over time. However, from Equation (5) we see that the exponential term scales as $-{\rho }_{m}^{1/3}{M}^{2/3}$ with ${M}^{2/3}$ largely determining the loss rate, so the sensitivity to ${\rho }_{m}$ is small.

It is important to note that in Equation (5) we have assumed the mass loss rate, $\dot{M}$, is sufficiently small that energy balance is dominated by radiative loss. This is indeed the case for exomoons of interest in this paper, where the low temperature water vapor atmospheres last for more than 1 Gyr. The surface pressures are well below ∼500 Pa until the runaway greenhouse limit is reached. For bodies with rapid hydrodynamic escape the numerical approach defined in Appendix A is appropriate.

The saturated case is derived from the same equations as the isothermal case (Equations (1)–(4)) but temperature is allowed to change with altitude. The temperature at the critical point at radius r_c in Figure 1 (where the isothermal sound speed u₀ equals the radial escape speed) is set such that numerically integrating Equations (1)–(4) from the critical point to the surface will result in a surface temperature equivalent to that in equilibrium with incoming solar flux taking into account the evaporative cooling (see Appendix A for details). Once the critical temperature is known the radial outflow velocity is readily calculated, and thus the mass loss rate.

In Case (2.1) we let the isothermal atmospheric temperature be set by just the incoming stellar flux and thus be equal to the effective temperature. However, for a thick water vapor atmosphere the surface will be heated by the greenhouse effect of the overlying atmosphere. In this case, we still used an isothermal atmosphere approximation but increased the atmospheric temperature by the total greenhouse warming of the atmosphere at the surface. The larger isothermal atmospheric temperature under this regime will increase the hydrodynamic loss rate compared to Case (2.1).

To account for the atmospheric greenhouse effect, we used a gray, radiative, plane-parallel approximation where the total gray atmospheric optical depth in the thermal infrared at the surface is given by

$$\begin{eqnarray}&&\tau =\displaystyle \frac{{\kappa }_{\mathrm{ref}}{P}^{2}}{2{{gP}}_{\mathrm{ref}}}\end{eqnarray} \tag{ 7 }$$

for mass absorption coefficient ${\kappa }_{\mathrm{ref}}$, at pressure P_ref and surface pressure P, where pressure broadening causes the ${P}^{2}$ dependency of the optical depth (Catling & Kasting 2017). Here we used ${\kappa }_{\mathrm{ref}}=0.05\,{{\rm{m}}}^{2}\,{\mathrm{kg}}^{-1}$ and ${P}_{\mathrm{ref}}={10}^{4}$ Pa (from Catling & Kasting 2017, Ch. 13). Having $\tau \propto {P}^{2}$ in Equation (7) is appropriate for thick atmospheres, which is the case when the runway greenhouse limit is approached. For thin atmospheres, $\tau \propto P$ is appropriate (Catling & Kasting 2017). In this study, the surface pressures are in the low-pressure regime (less than ∼500 Pa) until the runaway limit is reached. However, the difference between $\tau \propto {P}^{2}$ and $\tau \propto P$ in Equation (7) is small at such low pressures where the total greenhouse warming is less than a few K until the runaway limit is reached. Setting $\tau \propto P$ for such low-pressure moons in Equation (7) has a negligible impact on the calculated mass loss rate, so we approximate the optical depth of all atmospheres in this study with $\tau \propto {P}^{2}$. Once the total optical depth of the atmosphere is known from Equation (7), the first order global energy balance is given by (see Appendix B for derivation)

$$\begin{eqnarray}&&\displaystyle \frac{1}{4}(1-A){F}_{s}\left(1+\displaystyle \frac{\tau }{2}\right)=\left(\displaystyle \frac{{GM}}{{r}_{s}}+{L}_{v}\right){\rho }_{s}{u}_{s}+\sigma {T}_{s}^{4}\end{eqnarray} \tag{ 8 }$$

and from Equations (3) and (4) we derived an expression for ${\rho }_{s}{u}_{s}$ (see Appendix A)

$$\begin{eqnarray}&&{\rho }_{s}{u}_{s}={\rho }_{s}{u}_{0}{\left(\displaystyle \frac{{r}_{c}}{{r}_{s}}\right)}^{2}\exp \left[-\displaystyle \frac{1}{2}+\displaystyle \frac{{GM}}{{u}_{0}^{2}}\left(\displaystyle \frac{1}{{r}_{c}}-\displaystyle \frac{1}{{r}_{s}}\right)\right].\end{eqnarray} \tag{ 9 }$$

Equations (1), (7)–(9) were solved simultaneously to find T and u_s, with ${\rho }_{s}$ being given by the ideal gas law. The mass loss rate is then calculated by

$$\begin{eqnarray}&&\dot{M}=4\pi {\rho }_{s}{u}_{s}{r}_{s}^{2}\end{eqnarray} \tag{ 10 }$$

Using Equation (10), the time averaged loss rate is calculated and surface water lifetime is then obtained via Equation (6).

It is possible that no physically meaningful solution exists to Equations (1), (7)–(9). When the initial surface temperature, and therefore surface pressure, is large (above ∼260 K for this model), the optical depth given by Equation (7) will be significant. This will cause an increase in surface temperature, further increasing the surface pressure and thus the optical depth of the atmosphere. The positive feedback between temperature, pressure, and optical depth will cause Equations (1), (7)–(9) to have no valid solution if the initial surface temperature, set by the incoming stellar flux, is large. The exomoon–star distance where this positive feedback results in no solution is the runaway greenhouse limit, and it is akin the runaway limit found by Ingersoll (1969).

For all three model scenarios, we considered a pure water vapor atmosphere above a surface water reservoir. We looked at icy exomoons with masses ranging from 0.005 to 0.04 Earth masses between 0.9 and 2.0 au from a Sun-like star. This mass range includes bodies slightly smaller than Europa ($0.008\,{M}_{\mathrm{Earth}}$), and slightly larger than Ganymede ($0.025\,{M}_{\mathrm{Earth}}$). We set the Bond albedo to 0.2 for each run. We chose a Bond albedo of 0.2 for two reasons, the first is that it approximately represents the lower bound for icy moon Bond albedos in the solar system (Buratti 1991; Howett et al. 2010). In addition, a Bond albedo of 0.2 approximates the albedo of open ocean with partial cloud cover (Leconte et al. 2013; Goldblatt 2015). Should an icy moon form surface oceans, the 0.2 Bond albedo gives us the best representation when calculating water longevity.

The three key equations for hydrodynamic escape—continuity, momentum, and energy—can be written generally (e.g., multiple species, etc. Koskinen et al. 2013) but we will use a simplified spherically symmetric model with constant mean molecular mass (see Chapter 5 in Catling & Kasting 2017 for a more complete discussion of the topic). We assume the atmospheric density and atmospheric flow velocity only change in the radial direction. As such, the derivatives for mass continuity and momentum conservation are complete. Under these assumptions, the time-dependent and steady-state continuity and mass conservation equations are as follows. Continuity is given by:

$$\begin{eqnarray}&&\displaystyle \frac{\partial \rho }{\partial t}=-\displaystyle \frac{1}{{r}^{2}}\displaystyle \frac{d}{{dr}}({r}^{2}\rho u)\,,\mathrm{steady}\,\mathrm{state}:\,\,\displaystyle \frac{d}{{dr}}({r}^{2}\rho u)=0\end{eqnarray} \tag{ 13 }$$

where ρ is the mass density, r is the radial distance from the planet's center, and u is the atmospheric flow velocity. Momentum conservation is given by:

$$\begin{eqnarray}&&\displaystyle \frac{\partial (\rho u)}{\partial t}=-\rho u\displaystyle \frac{{du}}{{dr}}-\displaystyle \frac{{dp}}{{dr}}+\rho g,\\ &&\quad \mathrm{steady}\,\mathrm{state}:\,\,u\displaystyle \frac{{du}}{{dr}}+\displaystyle \frac{1}{\rho }\displaystyle \frac{{dp}}{{dr}}=g\end{eqnarray} \tag{ 14 }$$

where p is pressure, and g is gravity.

If we assume an isothermal atmosphere, we can relate pressure and density with the isothermal sound speed

$$\begin{eqnarray}&&{u}_{0}^{2}=\displaystyle \frac{{kT}}{m}\end{eqnarray} \tag{ 15 }$$

where k is the Boltzmann constant, T is the isothermal temperature, and m is the mean molecular mass of the atmosphere. From the ideal gas law

$$\begin{eqnarray}&&p=\rho {u}_{0}^{2}.\end{eqnarray} \tag{ 16 }$$

Integrating Equation (13) in the steady state, we get the mass escape rate per steradian of ${r}^{2}\rho u$, which, when combined with Equations (13) and (14), gives an expression for the isothermal planetary wind from a body with mass M:

$$\begin{eqnarray}({u}^{2}-{u}_{0}^{2})\displaystyle \frac{1}{u}\displaystyle \frac{{du}}{{dr}} & = & \displaystyle \frac{2{u}_{0}^{2}}{r}+g,\\ \mathrm{or}\,({u}^{2}-{u}_{0}^{2})\displaystyle \frac{1}{u}\displaystyle \frac{{du}}{{dr}} & = & \displaystyle \frac{2{u}_{0}^{2}}{r}-\displaystyle \frac{{GM}}{{r}^{2}}\end{eqnarray} \tag{ 17 }$$

Equation (17) is analogous to Parker's solar wind equation. For a strongly bound atmosphere at some critical distance from the planet's surface, the right-hand side of Equation (17) reaches zero, indicating that either the flow reaches the speed of sound or ${({du}/{dr})}_{c}=0$. The subsonic solution, ${({du}/{dr})}_{c}=0$, requires a finite background pressure that inhibits escape, so we will focus on the transonic solution where ${u}^{2}={u}_{0}^{2}$. The transonic solution has ${du}/{dr}\gt 0$ at all times and is consistent with a strongly bound atmosphere at the surface and zero pressure at infinity.

The critical distance r_c occurs in Equation (17) when ${u}^{2}={u}_{0}^{2}$ which gives us:

$$\begin{eqnarray}&&0=\displaystyle \frac{2{u}_{0}^{2}}{{r}_{c}}-\displaystyle \frac{{GM}}{{r}_{c}^{2}}.\end{eqnarray} \tag{ 18 }$$

Solving for r_c in Equation (18) we find

$$\begin{eqnarray}&&{r}_{c}=\displaystyle \frac{{GM}}{2{u}_{0}^{2}}.\end{eqnarray} \tag{ 19 }$$

If we integrate Equation (17) from the surface radius r_s to r_c and ignore the ${u}^{2}$ term near the surface, where it is negligible for bodies of interest in this study, we get the equation:

$$\begin{eqnarray}&&{\rho }_{s}{u}_{s}={\rho }_{s}{u}_{0}{\left(\displaystyle \frac{{r}_{c}}{{r}_{s}}\right)}^{2}\exp \left[-\displaystyle \frac{1}{2}+\displaystyle \frac{{GM}}{{u}_{0}^{2}}\left(\displaystyle \frac{1}{{r}_{c}}-\displaystyle \frac{1}{{r}_{s}}\right)\right].\end{eqnarray} \tag{ 20 }$$

As the radial distance from the moon increases the mass flux, $\rho u$ (in kg ${{\rm{m}}}^{-2}$ ${{\rm{s}}}^{-1}$), decreases. The steady state continuity given by Equation (13), when integrated gives $4\pi {r}^{2}\rho u=C$ where the constant of integration C is just the total rate of mass loss (in kg ${{\rm{s}}}^{-1}$) through a spherical surface. As r goes to infinity, $\rho u$ goes to 0 since $\rho u\propto 1/{r}^{2}$. Therefore, the outflowing wind loses kinetic energy as $r\to \infty$. Thus, the energy flux required to drive the escaping mass flux is given by the energy required to remove the mass flux from the gravity well of the moon, ${\rho }_{s}{u}_{s}{GM}/{r}_{s}$.

A first-order global energy balance between insolation and cooling via mass loss is then given by:

$$\begin{eqnarray}&&\mathop{\underbrace{\displaystyle \frac{1}{4}(1-A){F}_{s}}}\limits_{\mathrm{absorbed}\,\mathrm{stellar}\,\mathrm{flux}}=\mathop{\underbrace{\left(\displaystyle \frac{{GM}}{{r}_{s}}+{L}_{v}\right){\rho }_{s}{u}_{s}}}\limits_{\mathrm{mass}\,\mathrm{loss}\,\mathrm{flux}}+\mathop{\underbrace{\sigma {T}^{4}}}\limits_{\mathrm{radiative}\,\mathrm{cooling}}\end{eqnarray} \tag{ 21 }$$

where A is the Bond albedo, F_s is the incident stellar flux, and σ is the Stefan–Boltzmann constant. The escape flux is given by ${\rho }_{s}{u}_{s}$ and is multiplied by the energy required for that flux to escape the planet. The energy includes a gravitational potential energy term, and the latent heat of vaporization L_v (for this model ${L}_{v}=2.5\times {10}^{6}\,{\rm{J}}\,{\mathrm{kg}}^{-1}$). In Equation (21) we assume the atmosphere is transparent to both shortwave and infrared radiation.

Equations (20) and (21) can be solved simultaneously for the two unknowns u_s and T. Once solved, we can calculate the total escaping mass rate by:

$$\begin{eqnarray}&&\dot{M}=4\pi {\rho }_{s}{u}_{s}{r}_{s}^{2}\end{eqnarray} \tag{ 22 }$$

with ${\rho }_{s}$ being calculated from ${\rho }_{s}={P}_{s}/{u}_{0}^{2}$ with surface pressure P_s. We calculate surface pressure with Equation (1) for a C–C atmosphere given the surface temperature of water, where reference parameters are at the triple point: ${P}_{0}=611.73\,\mathrm{Pa}$, ${T}_{0}=273.16\,{\rm{K}}$. For our model, we only consider water worlds with pure ${{\rm{H}}}_{2}{\rm{O}}$ atmospheres, so estimating the surface density from the saturation vapor pressure is valid (Adams et al. 2008). We refer to this approach, where Equations (20) and (21) are solved numerically, as the Numerical Model.

For the slowly evaporating moons of interest in this study, those with surface water lasting more than 1 Gyr, the escape is so slow that it does not appreciably cool the moon. Thus, an analytic model can be derived by neglecting the mass loss flux cooling term in Equation (21). With the simplified Equation (21) our isothermal temperature is simply calculated from incoming stellar flux. And, from our assumption that the exomoons have an average bulk density of $2\,{\rm{g}}\,{\mathrm{cm}}^{-3}$, we can calculate the surface radius ${r}_{s}={(3M/(4\pi {\rho }_{\mathrm{exomoon}}))}^{1/3}$. Substituting these two equations into Equation (20) we find that

$$\begin{eqnarray}&&{u}_{s}={u}_{0}\displaystyle \frac{{r}_{c}^{2}}{{r}_{s}^{2}}\exp \left[\displaystyle \frac{3}{2}-\displaystyle \frac{G}{{u}_{0}^{2}}{\left(\displaystyle \frac{3}{4\pi {\rho }_{\mathrm{exomoon}}}\right)}^{\tfrac{-1}{3}}{M}^{2/3}\right].\end{eqnarray} \tag{ 23 }$$

By plugging Equation (23) into Equation (22), we get the analytic expression for the mass loss due to hydrodynamic escape given in Equation (5).

Calculating u_s in this manner assumes the temperature in our energy balance equation is a constant and set solely by the incoming stellar flux and the emitted thermal flux from the surface. For the low temperature bodies (<273 K) we are interested in for this study, Equation (23) gives identical results as the previously defined numerical model until the runaway limit is approached. See Figure 3 for a comparison.

Zoom In Zoom Out Reset image size

We also modeled hydrodynamic escape from a non-isothermal atmosphere, the saturated case. To model escape in the saturated case we start with Equations (13)–(16). Instead of using Equation (1) to relate temperature and pressure, we will approximate the Clausius–Clapeyron relation with an expression similar to Tentens' formula, given by

$$\begin{eqnarray}&&p={p}_{w}\exp (-{T}_{w}/T)\end{eqnarray} \tag{ 24 }$$

for reference temperature T_w and pressure p_w. A reasonable approximation for $250\lt T\lt 400\,{\rm{K}}$ over water takes ${T}_{w}\,=5200\,{\rm{K}}$ and ${p}_{w}=1.13\times {10}^{6}$ bar. A very good approximation for $150\lt T\lt 273\,{\rm{K}}$ over ice takes ${T}_{w}=6140\,{\rm{K}}$ and ${p}_{w}=3.53\times {10}^{7}$ bar. From Wexler (1977), whose expression we have approximated, the simple exponential fit is likely good to within a few percent for the temperatures in our model. This simplified expression is desirable because we want to work with an analytic expression for dT/dr.

We can eliminate p from Equation (14) using Equations (15) and (16), giving us

$$\begin{eqnarray}&&u\displaystyle \frac{{du}}{{dr}}+\displaystyle \frac{{u}_{0}^{2}}{\rho }\displaystyle \frac{d\rho }{{dr}}+\displaystyle \frac{{u}_{0}^{2}}{T}\displaystyle \frac{{dT}}{{dr}}=-\displaystyle \frac{{GM}}{{r}^{2}}.\end{eqnarray} \tag{ 25 }$$

We use Equation (24) to express dT/dr in terms of $d\rho /{dr}$

$$\begin{eqnarray}&&\displaystyle \frac{{u}_{0}^{2}}{\rho }\displaystyle \frac{d\rho }{{dr}}=\left(\displaystyle \frac{{T}_{w}-T}{T}\right)\displaystyle \frac{{u}_{0}^{2}}{T}\displaystyle \frac{{dT}}{{dr}}\end{eqnarray} \tag{ 26 }$$

and Equation (13) eliminates $d\rho /{dr}$ in terms of du/dr giving us our saturated wind equation

$$\begin{eqnarray}&&\left(u-\displaystyle \frac{{u}_{0}^{2}}{u}\left(\displaystyle \frac{{T}_{w}}{{T}_{w}-T}\right)\right)\displaystyle \frac{{du}}{{dr}}=\displaystyle \frac{2{u}_{0}^{2}}{r}\left(\displaystyle \frac{{T}_{w}}{{T}_{w}-T}\right)-\displaystyle \frac{{GM}}{{r}^{2}}\end{eqnarray} \tag{ 27 }$$

or equivalently as an expression for du/dr

$$\begin{eqnarray}&&\displaystyle \frac{1}{u}\displaystyle \frac{{du}}{{dr}}=\displaystyle \frac{N}{D}=\displaystyle \frac{(2{u}_{0}^{2}/r)({T}_{w}/({T}_{w}-T))-{GM}/{r}^{2}}{{u}^{2}-{u}_{0}^{2}({T}_{w}/({T}_{w}-T))}\end{eqnarray} \tag{ 28 }$$

where the numerator $N(r,T)$ is

$$\begin{eqnarray}&&N=\displaystyle \frac{2{u}_{0}^{2}}{r}\left(\displaystyle \frac{{T}_{w}}{{T}_{w}-T}\right)-\displaystyle \frac{{GM}}{{r}^{2}}\end{eqnarray} \tag{ 29 }$$

and the denominator $D(r,T,u)$ is

$$\begin{eqnarray}&&D={u}^{2}-{u}_{0}^{2}\left(\displaystyle \frac{{T}_{w}}{{T}_{w}-T}\right).\end{eqnarray} \tag{ 30 }$$

Equation (28) is the form we will use to numerically integrate u(r). Equation (27) can be written equivalently as

$$\begin{eqnarray}&&D\displaystyle \frac{1}{u}\displaystyle \frac{{du}}{{dr}}=N.\end{eqnarray} \tag{ 31 }$$

Recall from the isothermal case that, for hydrodynamic escape from a strongly bound atmosphere, ${du}/{dr}\gt 0$. Near the surface of the moon the numerator, $N(r,T)$, will be negative as the gravity term will dominate given that our atmosphere is strongly bound. At some distance r_c the $2{u}_{0}^{2}/r$ term will equal the force of gravity, so $N(r,T)=0$ at r_c. Since $N(r,T)=0$ and ${du}/{dr}\gt 0$, from Equation (31), $D(r,T,u)=0$ at r_c as well. At the critical point, ${N}_{c}=0$ provides a simple relation between T_c and r_c

$$\begin{eqnarray}&&{r}_{c}=\displaystyle \frac{{GMm}({T}_{w}-{T}_{c})}{2{{kT}}_{c}{T}_{w}}.\end{eqnarray} \tag{ 32 }$$

Similarly, ${D}_{c}=0$ relates u_c and T_c by

$$\begin{eqnarray}&&{u}_{c}^{2}={({u}_{0}^{2})}_{c}\left(\displaystyle \frac{{T}_{w}}{{T}_{w}-{T}_{c}}\right)=\displaystyle \frac{{GM}}{2{r}_{c}^{3}}.\end{eqnarray} \tag{ 33 }$$

The transonic solution is obtained by numerically integrating Equation (28) from the critical point to the surface. The first step is to solve for ${({du}/{dr})}_{c}$ at the critical point. This is obtained from Equation (28) by using L'Hopital's rule.

$$\begin{eqnarray}&&\displaystyle \frac{1}{{u}_{c}}{\left(\displaystyle \frac{{du}}{{dr}}\right)}_{c}=\displaystyle \frac{{N}_{c}}{{D}_{c}}=\displaystyle \frac{0}{0}=\displaystyle \frac{{({dN}/{dr})}_{c}}{{({dD}/{dr})}_{c}}.\end{eqnarray} \tag{ 34 }$$

The numerator becomes

$$\begin{eqnarray}{\left(\displaystyle \frac{{dN}}{{dr}}\right)}_{c} & = & \displaystyle \frac{{GM}}{{r}_{c}^{3}}-\displaystyle \frac{4{({u}_{0}^{2})}_{c}}{{r}_{c}^{2}}\displaystyle \frac{{T}_{w}^{2}{T}_{c}}{{({T}_{w}-{T}_{c})}^{3}}\\ & & -\displaystyle \frac{2{({u}_{0}^{2})}_{c}}{{r}_{c}}\displaystyle \frac{{T}_{w}^{2}{T}_{c}}{{({T}_{w}-{T}_{c})}^{3}}\displaystyle \frac{1}{{u}_{c}}{\left(\displaystyle \frac{{du}}{{dr}}\right)}_{c}\end{eqnarray} \tag{ 35 }$$

and the denominator becomes

$$\begin{eqnarray}&&{\left(\displaystyle \frac{{dD}}{{dr}}\right)}_{c}=\left(2+\displaystyle \frac{{T}_{w}{T}_{c}}{{({T}_{w}-{T}_{c})}^{2}}\right){u}_{c}\left(\displaystyle \frac{{du}}{{dr}}\right)+\displaystyle \frac{2{u}_{c}^{2}}{{r}_{c}^{2}}\displaystyle \frac{{T}_{w}{T}_{c}}{{({T}_{w}-{T}_{c})}^{2}}.\end{eqnarray} \tag{ 36 }$$

If we let $x\equiv {u}_{c}^{-1}{({du}/{dr})}_{c}$, simplify Equation (35) to replace ${GM}/{r}_{c}$ with Equation (33), and divide out the common factor ${u}_{c}^{2}$, then Equation (34) can be written as the quadratic equation

$$\begin{eqnarray}&&\left(2+\displaystyle \frac{{T}_{w}{T}_{c}}{{({T}_{w}-{T}_{c})}^{2}}\right){x}^{2}+\displaystyle \frac{4}{{r}_{c}}\displaystyle \frac{{T}_{w}{T}_{c}}{{({T}_{w}-{T}_{c})}^{2}}x\\ &&\quad +\displaystyle \frac{4}{{r}_{c}^{2}}\displaystyle \frac{{T}_{w}{T}_{c}}{{({T}_{w}-{T}_{c})}^{2}}-\displaystyle \frac{2}{{r}_{c}^{2}}=0.\end{eqnarray} \tag{ 37 }$$

The positive root of this equation corresponds to the accelerating flow at the critical point.

To find the mass flux loss, the first step is to guess an initial temperature at the critical point, T_c. Given T_c, we know ${u}_{c}^{2}$, p_c is given from Equation (24), and ${\rho }_{c}$ is given by the ideal gas law. From Equations (32) and (33), we get r_c and u_c respectively, which allows us to solve Equation (37) for the critical slope ${({du}/{dr})}_{c}$. Density can then be found at the new point from continuity, $\rho {{ur}}^{2}={\rho }_{c}{u}_{c}{r}_{c}^{2}$. Given ρ, we can solve for T and p from Equation (24) with the help of the ideal gas equation. This integration proceeds to the surface. The guess for T_c is adjusted numerically until the desired surface temperature (in balance with incoming stellar flux and mass loss given by Equation (21)) is achieved. Once the correct values are found, Equation (22) will give the mass loss rate.

For both isothermal and non-isothermal models, the surface temperature is assumed to be set by the incident solar flux averaged over time and hemisphere, which is given by Equation (2) for a rapidly rotating body. The isothermal case represents the warmest possible atmosphere neglecting greenhouse effects under the case of hydrodynamic escape. The non-isothermal case represents a minimum possible temperature for a water vapor atmosphere at r_c since it is saturated at all points based on the surface temperature set from the solar flux. These two models represent the extremes of atmospheric temperature profiles for a water vapor atmosphere, with the real solution likely somewhere between them.