Volatile transport on inhomogeneous surfaces: II. Numerical calculations (VT3D)

Several distant icy worlds have atmospheres that are in vapor-pressure equilibrium with their surface volatiles, including Pluto, Triton, and, probably, several large KBOs near perihelion. Studies of the volatile and thermal evolution of these have been limited by computational speed, especially for models that treat surfaces that vary with both latitude and longitude. In order to expedite such work, I present a new numerical model for the seasonal behavior of Pluto and Triton which (i) uses initial conditions that improve convergence, (ii) uses an expedient method for handling the transition between global and non-global atmospheres, (iii) includes local conservation of energy and global conservation of mass to partition energy between heating, conduction, and sublimation or condensation, (iv) uses time-stepping algorithms that ensure stability while allowing larger timesteps, and (v) can include longitudinal variability. This model, called VT3D, has been used in Young (2012), Young (2013), Olkin et al. (2015), Young and McKinnon (2013), and French et al. (2015).


Introduction
Pluto and Triton have atmospheres whose pressures have been measured by stellar occultations (e.g., Young et al. 2008a, Olkin et al. 1997) and spacecraft (Gurrola 1995, Krasnopolsky et al. 1993, Stern et al. 2015, Gladstone et al. 2016). These measurements reveal atmospheres for Pluto and Triton that are global in extent, almost certainly controlled by vapor-pressure equilibrium of the surface N 2 ice, (Spencer et al. 1997, Yelle et al. 1995, similar to the role of CO 2 on Mars (Leighton and Murray 1996).
Vapor pressure is an exceedingly sensitive function of temperature, and early models predicted that the surface pressures of Pluto and Triton would vary by orders of magnitude over their years (e.g., Paige 1992, 1996;Moore and Spencer 1990;Spencer and Moore 1992). Those early models were based on a single observation of the atmospheric pressure, either the Triton flyby in 1989 or the definitive discovery Pluto occultation in 1988 (e.g., Elliot and Young 1992). Since that time, further occultations have shown a large increase in the atmospheric pressures of both Pluto and Triton since the late 1980's (Elliot et al. 1998;Elliot et al. 2003). Other advances in the past decade include an improved understanding of the surface compositions of Pluto and Triton (Grundy andBuie 2001, Grundy et al. 2010). It is time for new models (Young 2012a, Young 2013, Olkin et al. 2015. This work describes the model used by Young (2012aYoung ( , 2013 and Olkin et al. (2015).
Since the rash of models in the 1990's, the large, volatile-covered ice worlds Pluto and Triton have been joined by other large, volatile-covered bodies in the outer solar system, including the large Kuiper Belt Objects (KBOs) Eris, Sedna, Makemake, Haumea, Quaoar (Schaller and Brown 2007), and 2007 OR10 (Brown et al. 2011). Some of these should have atmospheres at some time in their orbit (Stern & Trafton 2008). In particular, the 98% albedo of Eris argues for a perihelion atmosphere that collapses near aphelion, freshening Eris's surface (Sicardy et al. 2011).  Hansen and Paige 1996). Locally, we balance absorbed insolation, S, emitted thermal energy εσT 4 , and latent heat of sublimation or condensation, L S dm V /dt, where m V is the mass per area of the volatile slab and L S is the latent heat of sublimation. Additionally we balance (i) heat to and from the substrate, k dT/dz, where k is the thermal conductivity and dT/dz is the vertical gradient of temperature, and (ii) the heat capacity of the isothermal ice layer, m V dH V /dt ≈ m V c V dT V /dt, where H V is the enthalpy and c V is the specific heat of the volatile slab (subscript V for volatile). At the lower boundary, there is a heat flow of F. All variables except T V are free to vary with latitude and longitude. Compared with Young (2012a; Paper I), this figure illustrates (i) heating within the substrate for vertically varying k, and (ii) enthalpy of the ice slab, H V , to allow the treatment of solid-phase transitions.
The latent heat of sublimation term of the energy equation depends on the mass flux (dm V /dt in Fig. 2-1). For extremely thin atmospheres, such as on Io or possibly currently on Eris, some atmospheric flow occurs, but is ineffective in changing local surface temperatures (Fig 2-2A). In this case, the volatile slab temperature is controlled by local conditions only. The volatile slab temperature and the local atmospheric pressure are generally higher in areas of high insolation. For thin atmospheres, we assume no atmosphere over the bare areas. This approach allows efficient calculation of surface and subsurface temperatures. Once  I end this section with a few words about the modularity of the techniques described here.
• Anyone interested in this model should read Section 2 (this section) because it is short, and provides an overview.
• If your object has no volatiles, you do not need to read past Section 3. Volatile Transport II (VT3D) • If you want to characterize which processes are important in controlling surface temperatures, you can stop at the calculation of the thermal parameters, or Eq. 3.2-11 for volatile-free bodies, plus Eq. 4.2-7 and 4.2-8 for isolated volatile-covered areas, or Eq. 5.2-2a-d for volatile-covered interacting areas.
• If you want to very quickly approximate a temperature field based on the solar forcing, read Section 3.1-3.2, plus Section 4.1-4.2 if you have isolated volatiles, and Section 5.1-5.2 if you have interacting volatile covered areas. The critical equations are  • If you are calculating temperatures at one volatile-free location at a time, you can stop at Section 3.3. If you are calculating one isolated volatile-covered location at a time, read through Section 3.3, then skip ahead to Sections 4.1-4.3.
• If you are calculating roughly several hundred timesteps per period (e.g., to gain insights at short timescales or to make smooth plots), then the explicit equations will be stable, and the implicit equations will not save much computation time. In that case, you can skip those equations in Section 3.3, 4.3, and 5.3 that are described as implicit (roughly half of them), and all of sections 3.4b and 5.4b.  Table 1 recap the definitions of Areas I and III and their interaction with the atmosphere. For Area I (bare, no mass exchange, no atmosphere), the physics in VT3D is identical to the well-known thermophysical model (TPM) used to interpret thermal emission from airless bodies (e.g., Thomas et al. 2000;Spencer et al. 1989;Harris 1988). Heating in the top-most layer is balanced by thermal emission, insolation, and conduction; heating in interior layers is balanced by conduction only; heating in the lower layer is balanced by conduction and a flux condition at the lower boundary.

VT3D for bare locations (Areas I and III)
Area III (bare, no mass exchange, isobaric atmosphere) represents, for example, the "bedrock" H 2 O on current-day Triton. There is no volatile slab and no sublimation. The difference between the two bare area types are (i) Area III is a potential deposition site, and (ii) an increase in the volatile temperature for Area IV (volatile-covered, isobaric) also increases the pressure over Area III, so the atmosphere above Area III needs to be included in mass balance equation for Area IV. As long as there is no condensation (which will alter the state from bare to volatile-covered), the energy balance for Area III is the same as for Area I. Therefore, both bare areas, I and III, are treated in this section. Volatile Transport II (VT3D) These equations demonstrate several aspects of the numerical power of VT3D. In Section 3.1 and 3.2, I show the analytic expressions for the initial conditions, and show that a simple calculation can approximate the numerical solution. In Section 3.3, I present the explicit and implicit (Crank-Nicholson) numerical solutions for a single bare location, showing that solutions spin up in less than a quarter period. In Section 3.4, I show a compact representation of the linearized, discretized equations. In Section 3.5, I present a worked example of Mimas's diurnal temperatures, with code and output in the supplementary materials.

Areas I and III : Continuous expressions for bare areas
At the lower boundary, there may be positive (or negative) heat flow, F, which is balanced by upward (or downward) thermal conduction from a negative (or positive) thermal gradient: where k is the thermal conductivity, and T is the temperature. As with Paper I, z is a height coordinate, defined to be zero at the top of the substrate and decreasing downward. Thus, z = 0 at the substrate-volatile interface for locations where there is a volatile slab, or at the surface on volatile-free areas.
Within the substrate, I assume there are no heating sources, so net conductive heat flux is balanced by changes in the temperature of the substrate: where ρ is the density, c is the specific heat at constant pressure for the substrate, t is time, and T is the temperature.
The energy balance at the surface balances net heating with absorbed sunlight, thermal emission, and thermal conduction. There is no latent heat of sublimation or condensation. The total equation is where S is the absorbed solar energy, and ε is the emissivity, and σ is the Stefan-Boltzmann constant.
The first term of Eq. (3.1-3) describes the solar energy absorbed by the volatile slab, in power per area. For Triton, Pluto, Eris and other large KBOs, the fraction of sunlight absorbed by the atmosphere is small, and we do not need to alter S to account for atmospheric absorption. The absorbed solar energy at a particular location and time of day depends on the solar flux at 1 AU, S 1AU , the heliocentric distance, r, the hemispheric albedo, A h , and the cosine of the solar incidence angle, µ 0 (where µ 0 is 0 when the sun is below the horizon).
where S SS is the absorbed insolation at the sub-solar point. µ 0 depends on latitude, λ, subsolar latitude, λ 0 , and the hour angle, h (where h is the difference between the location's longitude and the subsolar longitude, defined to increase with time at any given location).
€ µ 0 = max 0,sin λ sin λ 0 + cos λ cos λ 0 cos h ( ) The hemispheric albedo, A h , is a local quantity, also known as the directional-hemispherical reflectance, hemispherical reflectance, or plane albedo (Hapke 1993). It is defined as the ratio of the total scattered power to the incident collimated power, € (S 1AU /r 2 )µ 0 , and depends on the location on the surface and the incidence angle. It is useful to approximate the hemispheric albedo by its average over all incidence angles, or where A S is known as the spherical reflectance, spherical albedo, or the Bond albedo (note, however, that Bond albedo is strictly defined for an entire surface). For typical phase functions in the outer solar system, substituting A S for A h tends to slightly underestimate solar heating for direct illumination and overestimate solar heating for large incidence angles. Since there is typically large uncertainty in the values of A S or A h due to uncertain phase functions, this distinction is usually ignored. The remainder of the paper uses A for A h , and does not distinguish between A h and A S .
The second term of Eq. (3.1-3) represents thermal energy emitted by the substrate. For a physical surface, this term might include such effects as self-heating from crater sides (Spencer 1990;Rozitis and Green 2011). In VT3D the emissivity, ε, is treated as a parameter Volatile Transport II (VT3D) that defines the power per area lost by thermal emission. Since ε can vary with location and time, it can be used to encompass these more subtle physical effects.
The final term of Eq. (3.1-3) represents thermal conduction from the substrate. If the substrate just below the interface is warmer than the surface temperature (dT/dz < 0), then conduction expressed by this term warms the surface.

Areas I and III : Analytic approximation and initialization bare areas
This section expands on key results of Paper I. The purpose is to introduce variables that will be used later, and to show the equations that will be used to initialize numerical calculations. For more discussion of the derivation, see Paper I.
If the solar insolation, S, at latitude λ and longitude φ, is a known function of time, t, with period P, then it can be approximated as a sum of M+1 sinusoidal terms where ω = 2π / P is the frequency of the diurnal or seasonal forcing, and Ŝ m λ,φ ( ) are the complex sinusoidal coefficients, with the hat indicating complex quantities (note, however, that S 0 is real). The coefficients are derived from the insolation in an expression closely related to the Fourier transform: A common application is diurnal forcing. For areas in permanent darkness, the solution is trivially € ˆ S m (λ,φ) = 0. For others, the diurnally averaged insolation can be expressed analytically (e.g., Levine et al. 1977). One first finds the maximum hour angle of illumination, h max , (h max = π for areas of constant illumination) where λ is the latitude and λ 0 is the sub-solar latitude, as before. The average insolation is a real quantity, so written without the hat, and is given by where the ratio S 0 /S SS is the longitudinal average of µ 0 . The decomposition of the solar forcing can also be written analytically. For a location that has hour angle h 0 at time t = 0, the first term is and, for m > 1, If the latitude of the surface element or the sub-solar latitude are near equatorial, then the solar terms are dominated by the first two terms, then diminish quickly with higher order; at the equator, the magnitudes of the terms are proportional to 1, π/2, 2/3, 0, -2/15, etc., (Paper I). Fig. 3-1 shows an example of the decomposition of insolation for a body with a sub-solar longitude of 2.24° and at a latitude of +30° into a constant plus one term (dashed) or seven terms (dot-dashed).

Insert fig 3-1 here.
1 mu0 = vt3d_solar_mu(lat, lon, lat0, lon0, /lonavg) 2 sol_terms = vt3d_sol_terms_diurnal(dist_sol_au, albedo, lat, h_phase0, lat_sol, n_terms) Fig 3-1. 1 Solid gray: numerical calculation of insolation on a bare spot at 9.5 AU with A = 0.6, at latitude 30°, a sub-solar latitude of 2.24°, and an hour angle at zero phase of -6 hours (-90°). The off-center maximum heating was chosen to force complex coefficients of the sinusoidal expansion. Solid: sinusoidal approximation with M=1, which captures the approximate phase and amplitude of the solar forcing. Dashed: sinusoidal approximation, with M=7, which is hard to see except at the "corners" near dawn and dusk because of the accuracy of this approximation.
As discussed in more detail in Paper I, the temperature can be written in terms of sinusoidal terms as well. If the density, specific heat, and thermal conductivity are constant with depth, then the solution to the diffusion equation (Eq. 3.1-2) with flux specifying the lower boundary condition (Eq. 3.1-1) is the sum of damped waves with wavelength € 2π 2mZ and e-folding distance of € 2mZ ( Fig. 3-2), where Γ = kρc (3.2-5) 2 1 vty16_fig3_1, phase, flux_sol, flux_sol_1, flux_sol_7 2 therminertia = vt3d_thermalinertia(dens, specheat, thermcond) Volatile Transport II (VT3D) is the thermal inertia (in cgs units of erg cm -2 K -1 s -1/2 , or MKS units of tiu = J m -2 K -1 s -1/2 , where tiu, or thermal inertia units, is the SI unit proposed in Putzig 2006), and is the skin depth, as defined by Spencer et al. (1989) and HP96. Other authors use definitions of the skin depths that differ by a constant from Eq. 3.2-6 (e.g., Mellon et al. 2008).
The solution to the conduction equation (3.1-2) can be written as (3.2-7) 2 where € ζ = z /Z is identical to the unitless scaled depth introduced by Spencer et al. (1989).
Temperatures for cases where the thermal-physical properties are variable with depth are treated elsewhere (Fivez & Thoen 1996;Grossel & Depasse 1998;Karam 2000 The goal is to use Eq. (3.2-7) to create initial conditions for numerical calculations in the three dimensions of latitude, longitude, and depth, given the coefficients for the temperature. There are three ways to do this. The simplest is to expand Eq. (3.1-3) to get the Fourier terms of the temperature directly, as described in Paper I and recapped here. The second is to follow this step by an adjustment of the average temperature, to ensure time-averaged energy balance. The third is to expand into Fourier terms of (T 4 ).
The average temperature, € T 0 , is found by substituting the sinusoidal forms of S and T into Eq. (3.1-3) and taking the first-order, time-averaged component, resulting in . This simply states that the mean temperature balances the mean solar insolation and the flux at the lower boundary condition.
The temperature coefficients, € ˆ T m , are found for each m by also substituting S and T into Eq. (3.1-3), taking the appropriate derivatives ( € d / dt →imω , € d /dz → imZ −1 ), and taking only those terms proportional to exp(imωt). The results are most simply expressed by defining the following variables, which represent the derivative of energy flux or heating with respect to temperature (in cgs units of erg cm -2 s -1 K -1 ) for the fundamental frequency: As described in Paper I, a system where Φ S is zero has temperatures that track the solar forcing, while positive Φ S serves to dampen the amplitude of the temperature variation and introduce a lag. The temperature variation (T m ) as a function of the solar variation for bare or volatile-covered area is found from: The temperature is then calculated from Eq. (3.2-8) overestimates mean temperatures, with the discrepancy being worse with larger peak-to-peak temperature variations, because the time average of T 4 is larger than Once an estimate of the peak-to-peak variation is found, the value of € T 0 can be adjusted downward so that the time-average thermal emission equals the sum of the insolation and internal heat flux, iterating over Eq. (3.2-9b) and (3.2-10) until the mean thermal emission converges on the mean absorbed insolation.
As described in Paper I, the time lag and smaller temperature variation can be described by a dimensionless parameter, Θ S .
(3.2-11) 1 phi_s = vt3d_dfluxdtemp_substrate(freq, therminertia) 2 phi_e = vt3d_dfluxdtemp_emit(emis, temp) 3 temp_terms = vt3d_temp_terms_bare(sol_terms,flux_int,emis,freq, therminertia) temp_terms= vt3d_temp_terms_bare_iter(sol_terms,flux_int,emis,freq, therminertia,thermcond) Volatile Transport II (VT3D) The ratio Θ S =4Φ S /Φ E is essentially the thermal parameter of Spencer et al. (1989), but defined at the time-averaged local temperature, rather than at the subsolar temperature. As in Paper I, Θ S can be used to quantify the shift and decrease in amplitude of the response to solar forcing . See Paper I for more discussion on the interpretation of Eq. (3.2-12) in terms of real quantities.   Fig 3-1, with unit emissivity. The result of the numerical integration is shown as solid, thick gray. Initial approximations to the temperature are shown for M = 1, without an adjustment of ! to balance energy fluxes (solid), M = 7, without Volatile Transport II (VT3D) an adjustment of T 0 (triple-dot-dash), and M = 7 with an adjusted value of T 0 (dashed), for a period of 22.6 hours and a thermal inertia of 16 tiu (for a thermal parameter Θ S = 6.0).
In some situations of large variation in the solar forcing and small values of Θ S , the linearization of T 4 is poor, and it is better to expand in the emitted flux instead.
The mean term is found from Eq. (3.2-8): F 0 E = S 0 + F . Before, we expanded the emitted flux in terms of temperature, but now we expand temperature in terms of emitted flux: The conduction is a now a small correction to the thermal emission, so the error in the linearization is confined to the second-order term. Substituting Eq. 3.2-14 into the original equation for energy balance, Eq. (3.1-3), gives: which, with some manipulation, gives an expression similar to Eq. (3.2-10): and one similar to Eq. (3.2-12): From F E , calculate the surface temperature and its Fourier terms from T = (F E / εσ ) 1/4 . In some cases, more Fourier terms (M = 30 to 100) need to be used than when calculating the temperature terms directly, to avoid ringing at sharp transitions in the solar forcing. The analysis in this section can be used for more complex insolation patterns as well. Any insolation pattern, no matter how complex, can be decomposed with a Fast Fourier Transform (FFT) algorithm. If the forcing happens on two different frequencies, such as seasonal and diurnal, then the sums (in e.g., 3.2-7) can be performed over a discrete set of m, not necessarily contiguous. For the specific case of combined seasonal and diurnal variation, we can often decouple the two timescales (Young 2012a). First, calculate the seasonal thermal wave as a function of time, using longitudinally averaged insolation. If the seasonal and diurnal skin depths are sufficiently different, then the diurnal wave is superimposed on the uppermost portion of the seasonal one, and the seasonal wave can be treated as a linear contribution to the diurnal wave. This is mathematically identical to an internal heat flux term, F, already introduced. In other words, the seasonal thermal heat flow to and from deeper layers affects the diurnal temperatures by affecting the energy flux at the lower boundary. This works because the orbital periods in the outer solar system are orders of magnitude longer than the rotational periods. Pluto, for example, has an orbital period of 248 years and a rotational period of only 6.4 days. The seasonal scale height is larger than the diurnal one by (248 years/6.4 days) 1/2 , or a factor of 119. A typical depth for the lower boundary is 6 diurnal skin depths. This is only 0.05 times the seasonal skin depth, or a tenth of a tick in Fig 3-2, clearly in the linear regime of the seasonal wave.

Numerical solution for a single bare location (Area I or III)
For some applications, the results of the analytic calculations may be adequate. For others, higher accuracy is needed. Even for these applications, the analytic solution provides an initial condition that improves convergence.
The continuous equations of Section 3.1 are converted to a form suitable for computation. This is done by discretizing the variables into L locations on the surface (indexed by l), J + 1 layers within the substrate (indexed by j), and choosing time step schemes that take the state from time n to time n + 1 (i.e., no leap-frogging time step schemes). The general approach is to treat the time step as a finite-difference diffusion problem, with flux conditions at both the upper and lower boundaries (Press et al., 2007;Haltiner and Williams 1984). Figure 3-4 represents the discretization of the numerical model. The substrate is divided into J+1 layers, indexed with j = 0 for the top-most layer to j = J for the lowest layer, and defined by a depth z j and a thickness Δ j . Depths (z) are less than or equal to zero, and become more negative with increasing index. Thicknesses of the layers (Δ j ) are positive. Thickness can vary with index j to speed computation ( Table 2). All layers except the top layer extend from € z j − Δ j /2 to € z j + Δ j /2, with temperature T l,j,n defined at the center of the layer. The Volatile Transport II (VT3D) top layer extends from z = 0 to z = -Δ 0 , with the temperature T l,0,n defined at the top of the layer. If the layering is the same across the globe, then The use of layers that are free to vary their thickness 1 with depth improves efficiency, since the computational time is proportional to the number of layers, requiring only a little additional computation at the beginning of a calculation. A common layering approach uses a geometrically increasing thickness, where the thickness of each layer is some factor larger than the layer above (typically a factor of 1.1 to 1.5, e.g., Hansen and Paige 1996;Keiffer 2013). When modeling a diurnal wave, this allows modest computational savings, since geometrically thickening layers can span down to six skin depths with 2-3 times fewer layers than for layers of equal thickness. Unevenly spaced layers is even more important for practical modeling of the diurnal and seasonal wave simultaneously. Because the skin depth is proportional to ω -1/2 , the ratio of diurnal and seasonal skin depths equals the square root of the ratio of their periods, if thermophysical properties are constant with depth. This is important even for Mars, where the orbital period is roughly 669 times the rotation period, so the seasonal skin depth is roughly 25 times the diurnal skin depth (if thermophysical properties are constant with depth). In the outer solar system, the orbital periods can be quite long, so that the equivalent ratio of seasonal to diurnal skin depths is 88 for Enceledus, 118 for Pluto, and 700 for Eris. If the thermal conductivity is greater at depth, these ratios can be even larger. Here the savings for geometrically thickening layers is dramatic, allowing calculation to 100, 1000, or even 10,000 diurnal skin depths with computational savings of ~20, ~100, or ~1000 respectively. For example, layers that begin with a thickness of 0.25 diurnal skin depths can reach 10,000 diurnal skin depths with only 41 layers for a thickening factor of 1.5, or with 87 layers for a thickening factor of 1.2.   The goal is to cast the equations as matrix operations to take advantage of the fast array operations that are available in many modern computer languages. The continuous equations of Section 3.1 can be cast as explicit equations (Fig 3-5), where the new temperature depends explicitly only on the previous temperature (Press et al. 2007;Haltiner and Williams 1984). The explicit expressions for diffusion equations are only accurate to first order in the time step, Δt, and require small time steps for stability. For explicit equations, the timesteps must satisfy (Δt / P) ≤ (Δz / Z ) 2 / 4π , or slightly more than 200 steps per period for a vertical sampling of 4 layers per skin depth.
The explicit linearized problem can be described with a (J +1) × (J +1) tridiagonal matrix ( Fig. 3-5). The new temperatures depend on the current temperatures in the layer above (with matrix element α, mnemonically "a for Above"), the current temperature in that layer (with matrix element η, mnemonically "h for Here"), and the current temperatures in the layer below (with matrix element β, mnemonically "b for Below").  Accuracy and stability can be improved by using implicit (Crank-Nicholson) methods, which solve equations involving both the current and the next temperatures ( Fig. 3-7), at the cost of computational complexity (Press et al. 2007;Haltiner and Williams 1984). The Crank-Nicholson scheme results in an equation that is accurate to second order in the time step, and satisfies von Neumann stability criteria for all sizes of time step. The implicit (Crank-Nicholson) problem uses two (J +1) × (J +1) tridiagonal matrices, with primed elements on the right-hand side of the equation and double-primed elements on the left. The goal of this section is to derive the matrix elements, which are summarized in Tables 3 to 5.

Matrix equation Matrix elements Explicit
To find the energy balance in layer 0, integrate the conduction equation (Eq. 3.1-2) over the top layer, from z = -Δ 0 to z = 0. Add this to the energy balance equation (Eq. 3.1-3) to get: where the overbar indicates the time-averaged value over the time step t n to t n+1 . The subscript for time in the insolation and emission terms is n' to indicate that it varies over the time interval from n to n+1. The change in enthalpy over layer 0 can be approximated as a function of the temperature sampled at the top of the layer: where € Φ l, j H has units of erg cm -2 s -1 K -1 , with the superscript H representing heat or enthalpy.
Eq. 3.3-2 samples the temperature of layer 0 at the top of the layer. If the temperature is integrated over layer 0 instead, then the slope of the temperature through layer 0 needs to be included; this depends on € T l,1,n , and is a second-order effect that I ignore here.
Defining a unitless measure of the time step, τ, (radians per timestep): and a unitless measure of the thickness of layer j expressed as a fraction of the skin depth (c.f., Spencer et al. 1989) (3.2-9a) and only depends on the physical properties of the problem. In contrast, H additionally includes non-dimensional factors that depend on the numerical choices of τ and € δ l, j . In general, I represent the fluxes-per-temperature that depend only on the physical properties with a single subscript for the physical process (e.g., S or H), and the ones that are discretized and depend on τ and € δ l, j with the superscript for the process and a subscript for the indices of location and time.
The average solar insolation between t n and t n+1 , S l,n' , depends on the geometry (heliocentric distance, and subsolar latitude and longitude) and the albedo. If the insolation is evaluated at the start of the timestep ( S l,n' ≈ S l,n ), then the results will be skewed in time by half a timestep, which is acceptable when timesteps are small (e.g., Spencer et al. 1989), but not at the larger timesteps allowed by the Crank-Nicholson method. A simple correction is to average the insolation at the start and end of the timestep Volatile Transport II (VT3D) The average thermal emission at the midpoint of the time interval is found by evaluating the first-order Taylor expansion of € T 4 at the average temperature for the time interval, € (T l,0,n +1 − T l,0,n ) /2, assuming that the emissivity is constant over the time interval.
where € Φ l,n E has units cgs of erg cm -2 s -1 K -1 , with the superscript E representing emission: where Φ E is defined in Eq. (3.2-9b). As with the enthalpy term, € Φ E only depends on the physical properties of the problem, and € Φ i, j E is the value used in the descretized calculation.
Unlike the enthalpy term, The next term in Eq. (3.3-1) is the thermal conduction. For explicit equations, the derivative is evaluated at the start of the time interval: where € Φ l, j K ,B has cgs units of erg cm -2 s -1 K -1 . The superscript K represents thermal conduction, and the superscript B represents conduction from the layer below. The expression for is essentially the distance to the middle of the layer below, modified to ensure continuity of fluxes at layer boundaries: and the unitless distances used for calculating thermal gradients from the layer below is Volatile Transport II (VT3D) Even if z j and Δ j are constant from one location to the next, the dependence on k means that € Δ l, j A and € Δ l, j B may vary with location. Again, € Φ S only depends on the physical properties of the problem, and € Φ i, j K ,B additionally includes non-dimensional factors that depend on the numerical implementation.
The more accurate and more stable Crank-Nicholson scheme (Press et al. 2007) replaces the derivative in  with the average of the derivatives calculated at the start and end of the time step (at time t n and time t n+1 ): The explicit discretized equation for energy balance of layer 0 becomes while the implicit equation is Collecting terms for the explicit equation (only € T l,0,n +1 on the left-hand side) results in: and for the implicit equation ( € T l,0,n +1 and € T l,1,n +1 on the left-hand side) results in: The goal is to now turn Eq. (3.3-16a) and (3.3-16b) into equations that express the matrix multiplication shown in Fig. 3-5 and Fig. 3-6, respectively. For the top layer, the matrix equations are T l,0,n+1 = η l,0,n T l,0,n + β l,0,n T l,1,n +γ l,0,n (3.3-17a) for explicit and for implicit time step schemes. Divide Eq. (3.3-15) by the total flux per temperature with units erg cm -2 s -1 K -1 , where the superscript T represents total, to get the matrix elements for j = 0, Areas I and III (Table 3). The forcing is a function of time, and is subscripted n.
Because the derivative of the thermal emission depends on time, the matrix elements β l,0,n and η l,0,n also depend on time.
For interior layers, the integral form of the diffusion equations (Eq. 3.1-2), averaged over time step n is € ρc ∂T ∂t dz where the overbar indicates the time-averaged value over the time step t n to t n+1 .
In the lowest layer, as in the interior layers, the net change in enthalpy of the layer is balanced by the difference between the flux entering from below and leaving from above ( Fig. 3-4). For layer J, unlike for layers j = 1... J-1, the flux from below is specified as a lower boundary condition. The energy balance equation for the lowest layer is: where F l is the heat flux at the lower boundary for location l.
The change in enthalpy over layer j (j = 1.. J) can be approximated as a function of the temperature sampled at the middle of the layer: Volatile Transport II (VT3D) The expressions for conduction into the layer above for from the layer j are similar to those into layer 0 from layer 1 (Eq. 3.3-9 and 3.3-13). For the explicit scheme, it is: and for the Crank-Nicholson implicit scheme it is: and substituting Eq.
Collect terms and divide by € Φ l, j H , to get the matrix elements (Tables 4 and 5). For the interior layers, the matrix elements α l,j , β l,j and η l,j are independent of time.
Fig. 3-7 compares the sinusoidal, explicit, and implicit calculations (at large and small timesteps) for a bare spot at 9.5 AU with A = 0.6, ε=1, and Γ = 16000 erg cm -2 s -1/2 K -1 = 16 tiu, at latitude of 30°, a sub-solar latitude of 2.24°, P = 22.6 hours, and an hour angle at zero phase of -6 hours (-90°). Calculations were performed on a vertical grid with  For the simplest, the single frequency sine-wave (M = 1) with no adjustment to the mean temperature, the numerical answer agrees well with the converged answer within 60° rotational phase (Fig 3-8, top); the other initial conditions agree with the converged answer even more quickly (within one time step, for the M=7 case with adjusted mean temperature). The calculated temperatures for both M = 1 and M=7 are too warm at the end of one period if the mean temperature for the initial condition was not adjusted (Fig 3-7, top and middle), but reaches the proper temperature with adjustment (Fig 3-7, bottom). All three cases shown have a similar convergence rate. Most of the gain is in the first period, with subsequent periods improving the solution by 12-20% per period.

3.4a. Overview and explicit timesteps
In this section, I present notes on how to solve the matrix equations in Figs. 3-5 and 3-6 in a way that takes advantage of the fact that for many problems substrate properties are often constant with time and location. I show how the implicit and explicit equations can be computed as a single matrix operation for those locations which share common substrate properties. This speeds calculation because it avoids "for-loop" constructions, with a speed savings that depends on the computer language involved. This section also shows how to 1 vty16_fig3_7a vty16_fig3_7b vty16_fig3_7c precompute the matrices associated with the substrate: both the elements for explicit calculations (the light gray elements in Fig 3-5, and the light-gray single-primed elements on the right-hand side of Fig 3-6), and the Lower-Upper (LU) decomposition of the matrix needed for implicit calculation (the light gray double-primed elements on the left-hand side of the equation in Fig 3-6.). Since LU decomposition is the first of the two steps needed in solving a tridiagonal matrix (Press et al. 1997), precomputing the LU decompositon of the substrate portion of the tridiagonal matrix cuts computation time roughly in half.
The key to these efficiencies is to separate the calculations for the uppermost layer (j = 0) from the lower layers (j = 1 to J). In addition to helping with the bare calculations, some of the notions introduced here will be required for implicit calculations of the interacting surfaces.
We separate the temperatures at a given location into a scalar describing the temperature of layer 0, T l,0,n , and a row vector of length J describing the temperatures of interior layers, T l,1.. J,n : With this separation, the timestep for Areas I and III can be written for the explicit timestep (Fig 3-5  € a l is a Jelement column vector with one non-zero element. 1 vt3d_step_expl_1loc, alpha_i, beta_0, beta_i, gamma_0, gamma_J, temp_0, temp_i Volatile Transport II (VT3D) S l is a tridiagonal matrix whose J-1 lower elements are vector with one non-zero element: To simplify the graphic, the time and location subscripts are dropped (e.g., η 0 for η l,0,n ). The temperature array is divided into the uppermost layer, T 0 , the next lower layer, T 1 , and remaining layers for j = 2..J, T j . The elements of the substrate matrix S consist of the three arrays α 2..J , η 1..J , and β 1..J-1 . Darker elements with white lettering correspond to the dark gray elements in Fig. 3-5, and change with each time step. Lighter elements with black lettering correspond to the light gray elements in Fig. 3-5, and are independent of time. White elements are zero.
Computation of Eq. (3.4-2) is displayed graphically in Fig 3-9. The uppermost temperature, € T l,0,n , is calculated by simple scalar arithmetic. The interior temperatures are calculated by matrix multiplication using a matrix that is likely to be time-independent, with additional terms added for T 1 and T J . In many applications, the substrate properties and internal heat flux are assumed to be constant over much of the body. In that case, in Eq. 3.4-2, the substrate arrays, and the temperatures in the interior layers 1 .. J as a J × L matrix with J rows and L columns formed by the concatenation of L temperature arrays of length J: It is admittedly awkward that T {L},0,n is an array, while α {L},1 is a scalar. I hope that context and Appendix A can help.
The surface temperatures are listed as a single 1-D array covering all the locations, rather than as a rectangular matrix of longitude and latitude. This is to simplify the matrix expressions of multiple locations. In addition, this allows for other divisions of the surface rather than simply a rectangular division, which tends to have needlessly small surface elements near the poles. Tiling schemes that maintain similar areas per tile need π/2 fewer tiles than equirectangular tiling schemes.
The new temperatures can be calculated in a way that takes advantage of array arithmetic: The computation represented by Eq.  1 vt3d_step_expl_nloc, alpha_i, beta_0, beta_i, gamma_0, gamma_J, temp_0, temp_i Elements are labeled as in Fig 3-8. . "*" indicates element-by-element multiplication of two arrays (above the dotted line) or the multiplication of each row by a scalar (below the dotted line). "X" indicates matrix multiplication.

3.4b. Implicit timesteps
With the division of temperatures into layer 0 and layer 1 .. J in Eq. (3.4-1), the implicit timestep for a single location for Areas I and III (Fig 3-6) can be written as     We treat Eq. (3.4-9b) as a banded tridiagonal matrix to take advantage of the fact that the terms € " " a l and € " " S l are constant with time. This is a special case of inversion by partitioning, whose solution is presented in Press et al. (2007;section 2.7.4). A similar problem was addressed by Xing-Bo (2009). This allows us to precompute the lower-upper (LU) decomposition of € " " S l . The solution to  can be written by defining two column vectors € y l and z l,n of length J, and two scalars c l,n and d l,n : The solution is shown graphically in Fig 3-12. Note that the only the time-independent substrate matrix needs to be inverted, and this can be done at the start of the computation, rather than for each time step. Furthermore, the array y is also independent of time.  The new temperatures are then Volatile Transport II (VT3D) is an outer product of a J-length column vector and an L-length row vector, yielding a J × L matrix obtained by The graphical schematic is shown in Fig 3- (equivalent to the outer product of two arrays for the multiplication below the lowest dotted line).

Example: Mimas
As a worked example, Fig 3-

VT3D for local volatile-covered locations (Area II)
In this section, I consider locations that have volatiles on their surfaces, but for which the energy balance is essentially local. For worlds where the surface pressure is too low to effectively transport volatiles over the surface, the transport of energy, through latent heat of sublimation and deposition, does not effectively influence on the surface temperatures. This is the case on Io, and almost certainly the case on the large volatile-covered Kuiper-belt objects when far from perihelion. These are the isolated, volatile-covered areas (Area II) in Fig 2-2.
Within the substrate, the physics of thermal conduction and the lower boundary condition for the volatile covered locations (Area II) is identical as for the bare locations (Areas I and III, Section 3), and will not be repeated here. At the surface, on the other hand, the energy equation contains two new terms, one related to the energy needed to heat the volatile slab, and another related to latent heat exchange between the surface and the local gas column via deposition and sublimation. The continuous form is discussed in Section 4.1 and analytic expression for an initial condition is discussed in Section 4.2. Because the energy equations are strictly local, the form of the numerical implementation is very similar to that in Section 3. Only the form of the matrix elements € η 0 and € β 0 change, as discussed in Section 4.3.

Analytic expressions for isolated volatile-covered locations (Area II)
The energy equation at the surface balances net heating or crystalline phase changes with absorbed sunlight, thermal emission, thermal conduction, and latent heat of sublimation/condensation. The total energy equation is where m V is the mass per area of the volatile slab, € ∂H V /∂t is the time derivative of the enthalpy of the volatile slab in energy per mass (equal to € c V ∂T /∂t if there is no phase change, see Eq. 4.1-2, where c V is the specific heat of the volatile slab. Note c V is subscripted V for volatile, not V for constant volume), and L S is the latent heat of sublimation. L S is subscripted with S to distinguish it from the latent heat of crystalline phase change (L C ) and or the number of discrete locations on the surface (L, Section 3.3).
At the surface, a volatile slab is assumed to be isothermal within its vertical extent (See Fig 4.1), with a temperature equal to that at the top of the substrate. As described in Paper I, the isothermal slab was assumed in Hansen and Paige (1992) and Hansen and Paige (1996). This has been justified (David Paige, personal communication) by assuming that if the slab porous, it is in contact with the local atmosphere and the gas can isothermalize the solid; conversely, if the slab is not porous (e.g., from annealing, Eluszkiewicz et al., 1998) then its conductivity will be high, helping to isothermalize a thin enough slab. For very thick deposits, such as the suspected N 2 reservoir seen on Pluto, one approach is to keep track of mass per area of the volatiles available for sublimation as a separate quantity from the mass per area that is isothermalized (Young et al., 2016). Layering within the volatile slab will be treated in a later paper.

Energy Fluxes Thickness or Separation
Temperature temperature of a crystalline phase transition, the derivative of H V with respect to T at constant pressure equals c V , the specific heat of the volatile slab, Eq. (4.1-2a). Adding energy to the slab raises its temperature. At the temperature of a crystalline ice phase transition, the latent heat equals the difference in H V between two phases (L C ); adding energy to the slab converts ice from the low-temperature to the high-temperature phase without changing the temperature. This gives: where T C is the temperature of a crystalline phase transition, L C is the latent heat of crystalline phase change, and X is the mass fraction of the high-temperature phase. If c V is treated as a constant, then we can write € H V = c V T +L C X , which is proportional to the "pseudo temperature" used by John Spencer (personal communication).
Tracking the enthalpy of the slab, rather than its temperature, was introduced because N 2 has a reversible transition between the α and β phase at 35.6 K (e.g., Scott 1976), a relevant temperature for Pluto, Triton, and elsewhere in the outer solar system. Some volatile ices have no solid-state phase transitions at relevant temperatures, which simplifies matters. Others have multiple transitions, or non-reversible transitions. In all cases, the enthalpy is the general quantity that can account for phases as well as temperatures, and Eq. 4.1-2b represents the "special case" of enthalpy change at a phase transition temperature.
Area II satisfies local energy and mass balance. Assuming negligible horizontal transport of mass, any mass lost by the atmosphere either condenses or escapes.
where m A is the mass per area of the atmosphere, and E is the escape rate in mass per area per time. Negative values of E can be used to account for injection into the atmosphere from non-sublimation sources such as geysers (see Paper I). If the atmosphere is in vapor-pressure equilibrium with the surface, then the mass of the atmosphere is a function only of the surface pressure and effective gravity (defined by g = p S /m A , which is smaller than the surface gravity for extended atmospheres by a factor of 1 -2 H/R, where H is the scale height and R is the surface radius, see Paper I): where € p S (T) is the equilibrium vapor pressure at temperature T. The pressure derivative in Eq. (4.1-4) can be evaluated using the Clausius-Clapeyron relation, where m molec is the mass of one molecule and k B is Bolzmann's constant. Substituting Eqs.
(4.1-2a), (4.1-2b), (4.1-3), and (4.1-4) into Eq. (4.1-1) and collecting like terms yields: Eq (4.1-6a) is strikingly similar to the equivalent equation for the bare areas (3.1-3), differing only by the inclusion of the enthalpy and latent heat terms on the left-hand side, and the latent heat of the escaping atmosphere on the right side. The enthalpy and latent heat of sublimation introduce terms proportional to the frequency, ω, in the analytic equations (Section 4.2). They also introduce two additional terms to the total expression for the change in energy flux per temperature for the upper-most layer ( € Φ l,n T ) in the numeric solutions (Section 4.3), but the form of the matrix equations is unchanged. When there is a phase change, (4.1-6b), the analytic and numeric forms are both simpler, as the temperature does not change with time.

Analytic approximation and initialization for isolated volatile-covered areas (Area II)
As in Section 3.2, an analytic form of the continuous equations (Eq. 4.1-6a, b) can be found by decomposing the solar insolation and temperature into a sum of sinusoidal terms of frequency ω (Eqs. 3.2-1, 3.2-7). Additionally, we specify that the temperature of the volatile slab equals the substrate temperature at the substrate-slab interface € T V (λ,φ,t) = T(λ,φ,z = 0,t) (4.2-1) The escape rate is decomposed into a sum of sinusoidal terms in an analogous manner to the solar forcing where is the frequency of the diurnal or seasonal forcing, and € ˆ E m is the complex sinusoidal coefficient (the complexity is indicated by the hat).
As in Section 3.2, the average temperature is found by substituting the sinusoidal forms of S and T into Eqs. (4.1-6a, b) and taking the first-order, time-averaged component.
Φ V is simply related the to specific heat per area of the volatile slab, being the energy per degree per area. Φ A is related to the energy needed for the atmosphere to vary its column mass (atmospheric "breathing"). If the surface temperature rises, the equilibrium pressure rises too. The column mass of the equilibrium atmosphere increases due to sublimation from the surface. This takes energy, through the latent heat of sublimation. The result is that the specific heat of the volatile slab and the atmospheric "breathing" delay and decrease the thermal response (Paper I). The resulting expansion of 4.1-6a is: 1 temp_0 = vt3d_temp_term0_local(sol_0, flux_int, emis, latheat, mflux_esc) 2 phi_v = vt3d_dfluxdtemp_slab(freq, mass_0, specheat) 3 phi_a = vt3d_dfluxdtemp_atm(freq, temp_v, frac_varea, gravacc, name_species) € ω = 2π /P Volatile Transport II (VT3D) If the equilibrium temperature is at a crystalline phase boundary, then the corresponding equation for the change in the slab's state is Latent heat of escaping gas    , T = T C (4.2-6) As described in Paper I, we can define non-dimensional thermal parameters, analogous to the thermal parameter of Spencer et al. (1989), to quantify the importance of heating of the volatile slab and atmospheric breathing. The substrate thermal parameter, Θ S , is defined in Eq. 3.2-11. Two new parameters are: Substituting into Eq. (4.2-5) shows how the amplitude and phase of the thermal response depends on the thermal inertia, the specific heat and depth of the volatile slab, and the extent of the atmospheric "breathing." As for the bare areas (Areas I and III), the expansion can be written in terms of the emitted thermal flux in the case of large temperature variations, giving Volatile Transport II (VT3D)

Numerical solution for isolated volatile-covered areas (Area II)
The discretization for the interior layers (j = 1..J-1) and the lowest layer (j = J) is the same for the isolated, volatile-covered locations (Area II) as it is for the bare locations (Areas I and III). The discretization for the volatile slab and the upper two layers are shown in Fig 4-1. Although the physics is different in the presence of a volatile, the numerics are nearly identical for all calculations on a local level, whether volatiles are present or not.
First consider usual case where the volatile slab is not at a crystalline phase transition temperature. As with Areas I and III, to find the energy balance in layer 0, integrate the conduction equation (Eq. 3.1-2) over the top layer, from z = -Δ 0 to z = 0. Add this to the energy balance equation (Eq. 4.1-6a) to get Eq. (4.3-1), the volatile-covered equivalent to Eq. where the overbar indicates the time-averaged value over the time step t n to t n+1 .
The enthalpy of layer 0, insolation, emission, and conduction are the same as for Areas I and III (Section 3.3).
The second term in Eq (4.3-1) reflects the change in the enthalpy of the volatile slab with temperature. The volatile slab mass, m l,n V , can change over the time interval. However, this change is going to be small unless the slab is about to completely sublime, in which case this term contributes little. Ignoring the change in volatile slab mass during the time interval, this term becomes: Volatile Transport II (VT3D) where c l V is the specific heat of the volatile slab at location l, € Φ l,n V has units of erg cm -2 s -1 K -1 , and the superscript V stands for volatile slab The third term in Eq (4.3-1) is related to the amount of latent heat required sublime the atmospheric mass needed to maintain vapor-pressure equilibrium with a higher surface temperature. Linearizing the change in surface pressure with respect to time gives where € Φ l,n A has units of erg cm -2 s -1 K -1 , and the superscript A stands for atmosphere.
The temperature dependence of pressure is highly non-linear. If this is a dominant source of error, then one either chooses a small τ, or iterates from an initial guess at a temperature € T l,0,n +1 approx to an improved temperature € T l,0,n +1 . In the latter case, by Taylor expansion of p around This can be cast in a form parallel to that of Eq. (4.3-4) by where the derivative in € Φ l,n A is evaluated at the current guess at a temperature € T l,0,n +1 approx . The term € F l,n A has units of erg cm -2 s -1 , and combines mathematically with the solar forcing.
By writing  in terms of the change in temperature relative to the previous time step (i.e., € T l,0,n +1 − T l,0,n ), rather than in terms of the smaller change in temperature relative to the current guess (i.e., € T l,0,n +1 − T l,0,n +1 approx ),  can be simply combined with the other terms in the discretized energy equation. On the first iteration, € T l,0,n +1 approx = T l,0,n , and Eq. (4.3-7) reduces to , so Eq. (4.3-7) can be used with very little added computational complexity.
The escape rate, E, if present, can be calculated at the start or mid time, similarly to the insolation.
Substituting the expressions for the explicit equations gives an equation similar to Eq. 3.3-14: Collecting terms for the explicit equation gives As in Section 3.3, divide by € Φ l,n T , with units erg cm -2 s -1 K -1 , where the superscript T represents total, and the total "flux-per-temperature" now includes terms for enthalpy of the slab and interaction with the atmosphere The explicit equations for Area II can be written in a form that is identical to the explicit equation for the bare areas, Areas I and II (See Fig 3-6), with the resulting matrix elements given in the first row of Table 6. Volatile Transport II (VT3D)
The implicit form of the energy balance equation for Area II away from a crystalline transition temperature is found by substituting the Crank-Nicholson expression for the conduction term into Eq. 4.3-1. The energy balance for the implicit equation is Collecting terms for the implicit equation gives the volatile-covered equivalent to 3.3-16b: Again, divide by € Φ l,n T , with the resulting matrix elements given the second row in Table 5.
The matrix elements for j = 1 to J are identical as for the bare areas, Areas I and III (Section Volatile Transport II (VT3D) Tables 3 and 4). The methods for solving the matrix equations are identical as for the bare areas, Areas I and III (Section 3.4).

Matrix operations for single or multiple isolated volatile-covered locations (Area II)
As with Areas I and III, computation can be sped up considerably by taking advantage of matrix operations to calculate the temperature evolution on multiple locations with a single operation. The form of the matrices for isolated volatile-covered locations (Area II) is the same as for bare locations (Areas I and III). Therefore, once the matrix elements are found, the calculations can proceed identically to Section 3.4.

Example: KBOs with bare areas or locally-supported atmospheres
As an example, consider a point on the equator of a generic KBO with A = 0.7, ε = 0.9, no internal heat flux or mass loss, 50 cm s -2 surface gravity, and an equatorial sub-solar latitude, at a range of heliocentric distances (r) from 30 to 80 AU (Fig 4.2). The thermal parameter for the substrate (Θ S ) ranges from ~4-17 for 5 tiu (similar to those found by Lellouch et al. 2013), and ~1600 to ~7000 for 2100 tiu (pure, compact water ice). The thermal parameter for heating one g cm -2 of a volatile slab (Θ V ) is 7 times larger than Θ S for the 5-tiu case, or 27 to117, so it is not insignificant. Both Θ S and Θ V increase with heliocentric distance, since their numerators stay constant and their denominators (proportional to T 3 ) decrease. Thus, the same object can be a slow rotator at perihelion and a fast rotator at aphelion. The atmospheric thermal parameter (Θ A ), which has equilibrium pressure in the numerator, varies by 5-7 orders of magnitude over the range of r from 5.2 x 10 3 to 8.6 x 10 -2 for N 2 and 1.0 to 1.5 x 10 -7 for CH 4 .
For simplicity, the remainder of Fig 4.2 only contrasts a bare substrate with thermal inertia of 5 tiu, with a surface that is either N 2 -covered or CH 4 -covered. The effect of the decreasing temperature and increasing Θ S with r is clear in the progression for the substrate temperatures in the second panel, which plots the temperature for a bare substrate as a black solid line. The N 2 atmospheric "breathing" (green dashed line) has little effect at 70 AU, modifies the temperatures at 60 AU; by 40 and 30 AU, it nearly flattens out the temperature variation. The atmospheric breathing shifts the maximum by 90° phase, while the thermal conduction into the substrate shifts it by 45°; this is most evident when Θ A is comparable to Θ S , such as for N 2 near 60 AU or CH 4 (red triple-dot-dashed line) near 30 AU. This has the effect of decreasing the peak temperature, and increasing the temperature at both the dawn and dusk limbs.
The third panel of Fig. 4.2 shows the increase in the mean and amplitude of the temperature for a bare substrate (gray fill) with decreasing heliocentric distance. For N 2covered areas (green slanted fill), the temperatures are similar to the bare temperatures beyond ~70 AU. Closer than that, first the maximum temperature decreases while the dusk temperature rises, then the minimum and dawn temperature rise in tandem, until finally the maximum, dusk, dawn, and minimum temperatures all converge inward of 40 AU. For CH 4covered areas (red vertical fill), the temperatures match the bare temperatures beyond ~40 AU; inward of 40 AU, as with the N 2 , the maximum temperature decreases while the dusk temperature rises, with the slight rise in the minimum and dawn temperatures. The corresponding minimum, dawn, dusk, and maximum pressures are shown in the final panel. values of thermal inertia, for 1 g of slab at 1.3e7 erg g -1 K -1 , and for atmospheric "breathing" by N 2 (green dashed) and CH 4 (red triple-dotdashed). Second: Temperature for Γ=5 tiu over a single day for bare (solid, black), N 2 -covered (green dashed, indistinguishable from bare at 70 AU), and CH 4 -covered (red triple-dot-dashed, indistinguishable from bare at 40 AU and farther) at selected distances. Third: Minimum, dawn, dusk, and maximum temperatures over a range of distances for bare (gray), N 2 -covered (green), and CH 4 -covered (red) areas. Fourth: Minimum, dawn, dusk, and maximum pressures over a range of distances for N 2 -covered (green) and CH 4 -covered (red) areas.

VT3D for interacting volatile-covered areas (Area IV)
Currently, Pluto and Triton are expected to have similar surface pressures over the entire globe, independent of local insolation (Trafton & Stern 1983, Trafton 1984, Spencer et al. 1987. N 2 sublimes from areas of high insolation, with latent heat loss balancing the excess insolation. Sublimation winds carry this mass to areas of low insolation, where N 2 is deposited, adding latent heat as well as solid N 2 (Fig 2-2B). As long as the atmosphere is dense enough, transport of mass and latent heat will keep the volatile ice temperatures nearly constant over the globe. Through vapor-pressure equilibrium, the surface pressures will also be nearly constant. If the atmosphere is thin enough so that the sublimation winds are a significant fraction of the sound speed, then the surface pressures will vary over the globe. This case can be handled efficiently by treating the surface as a "splice" between the interaction regions or isobaric regions, which share the same surface pressure, and local regions, for which the surface pressure varies with location (Fig 2-2C In this section, I consider areas that have volatiles on their surfaces and which interact to share the same volatile ice temperature and surface pressure. This includes the entire globe for dense atmospheres, or the interacting portions of the splice for intermediate atmospheres ( See Fig 2-2B, 2-2C). I will discuss the continuous equations in Section 5.1, analytic equations in Section 5.2, the discrete equations in Section 5.3, and efficient solutions to the matrix equations in Section 5.4. In Section 5.5, I present a worked example of Pluto's seasonal activity, with code and output in the supplementary materials.

Continuous expressions for interacting volatile-covered locations (Area IV)
For interacting volatile-covered locations, Area IV, energy is transported between locations through mass transport of volatiles through the atmosphere and the latent heat of sublimation. What ties the multiple locations together is (1) a common volatile-ice temperature, € T V , and (2) conservation of mass over the interacting regions. This latter includes the atmosphere over all areas that share a single surface pressure, whether bare (Area III) or volatile-covered (Area IV), because raising the surface pressure increases the atmospheric mass over all locations that share a common surface pressure. That is, if the surface pressure of the atmosphere increases in the region of effective transport, the mass of the atmosphere will increase above both the volatile-covered areas (Area IV) and the bare areas (Area III). The expression for mass balance in the area of effective transport is found by integrating Eq. 4.1-4 over both Area III and Area IV: where Ω III and Ω IV represent the solid angle of areas III and IV. Both the surface pressure and the temperature of the volatile slab are constant over Areas III and IV; the terms involving gravity, pressure, and temperature can be factored out of the middle integral. Futhermore, the mass flux for Area III is zero, so that the first integral can be evaluated over just Area IV. With these changes, the mass balance equation becomes The areal average of the mass flux over Area IV is: where brackets represent an areal average over Area IV. The atmosphere escapes from above both bare and volatile-covered areas, so the areal average of E is taken over Areas III and IV: where primed brackets represent an areal average over Area III and Area IV.
f V is the fraction of the interacting areas (III and IV) covered with volatiles. In Paper I, which only treated a global atmosphere, this was fraction of the surface covered by volatiles. Here, with the possibility of a spliced atmosphere, the expression is written more generally.
With these definitions, the equation for mass balance over the areas of isobaric surface pressure becomes Eq. 5.1-5 illustrates the significance of the fraction of the surface covered by volatiles, f V . If the volatile ices are confined to a small patch, then that patch has to lose a lot of mass to supply an increase of the entire atmosphere in the isobaric area.
The local energy balance is the same as for localized volatile-covered areas, Eq. (4.1-1). Integrating Eq. 4.1-1 over Area IV, and substituting the equation for conservation of mass over isobaric areas, yields an equation for energy balance over all of Area IV, using the same notation for spatial averages as in Eq. 5.1-3a.
While the temperature of isolated volatile-covered areas depend only on local conditions b), the volatile ice temperature in the interacting areas depends on the spatial average of energy sources and sinks.

Analytic approximation and initialization for interacting volatile-covered locations (Area IV)
The analytic form of the continuous equations (Eq. 5.1-6) is very similar to that for the isolated volatile-covered areas, Area II, described in Section 4. As in Section 4, the solar forcing, the atmospheric escape, and the thermal wave are (1) decomposed into sinusoidal terms (3.2-1 for absorbed insolation, 3.2-7 for temperature, and 4.2-2 for escape), (2) substituted into Eq 5.1-6, and (3) isolated term-by-term. The m=0 term gives the expression for the time-averaged temperature: To find the variation in the temperature (the terms with m = 1 and higher), substitute the expressions for solar forcing, temperature, and escape into Eq. 5.1-6, expand the thermal emission term to first order in T m , and take only those terms proportional to exp(imωt). This expression is simpler with the spatially averaged versions of the "flux-per-temperature" expressions: If the substrate under all of the volatile ices has the same thermophysical properties, then the first two terms reduce to their local equivalents: Eq. 3.2-9a, b. Likewise, if the specific heat of the volatile ices are the same over Area IV, then the third equation (Eq. 5.2-2c) differs from its local equivalent (4.2-4a) simply by replacing the local volatile slab mass with the areal average. If there is no bare ground in the isobaric area (that is, if no locations are Area III), then f V = 1, and the last expression  is identical to its local equivalent (4.2-4b). However, if only part of the isobaric area is volatile-covered, then Φ A (T ) > Φ A (T ) . A change in volatile temperature increases the atmosphere above both bare and volatile-covered locations in the isobaric areas, so more mass is exchanged between the surface and atmosphere, and more latent heat of sublimation is required. This means that the latent heat term is more effective at suppressing the temperature variation when there is a smaller fraction of surface volatiles. For temperatures away from a crystalline phase, with these substitutions, the spatially averaged energy equation is: If the equilibrium temperature is at a crystalline phase transition, then the corresponding equation for the change in the slab's state is

Numerical solution for interacting volatile-covered locations (Area IV)
Fig 5-1 shows the interaction between different locations in Area IV. There is no horizontal heat flow within the substrate. However, the volatile slabs exchange energy through latent heat of sublimation and condensation, and share a single temperature, T V . The temperature of the volatile ice slab therefore depends on the insolation over the entire volatile-covered interacting region, and the conduction from each of the substrate layers (layer 0) that immediately underlie the volatile ice slab. The temperatures of each of the topmost substrate layers depend, in turn, on the single volatile slab temperature, through thermal conduction. Because there is no horizontal heat flow within the substrate, the discretization for layers j = 2 .. J is the same as the other areas, so that much of the matrix form for the explicit equations is tridiagonal (Fig. 5-2). However, because volatile slabs of the areas interact ( Fig  5-1), the explicit discretized equation for the new T V has non-zero coefficients accounting for the conduction upward from each of the j = 1 layers (the upper row of the matrix). Similarly, the explicit discretized equation for each new T 1 has non-zero coefficients accounting for the conduction downward from each of the j = 0 slabs, all assumed to be at T V . The resulting  The elements of the substrate arrays are derived from the discretation of the conductivity equation, Eq. (3.1-2), as before. The matrix elements for the substrate-the light gray elements in Figs 5-2 and 5-3-are unchanged from the previous cases. This holds even for the first layer, j = 1. The dependence of the temperature of the first layer at location l, T l,1 , depends only on the temperature below (T l,2 ) and above (T l,0 ). For Area IV, the assumption is that T l,0 = T V (that is, the upper surface of the substrate equals the volatile slab temperature, Fig 5-1). While this changes the format of the matrices (the line of α's in the left-most column in Fig 5-2 and 5-3), it does not change the value for the α's themselves. To find the elements for the implicit arrays α l, j , η l, j , β l, j (j = 1 .. J) and the lower-boundary element γ l,J , or their explicit counterparts (primed for the right-hand side and double-primed for the left) consult Tables 3 and 4.
The elements for the volatile slab-the dark gray elements on the top row of Figs 5-2 and 5-3-are related to, but different than, the corresponding elements for Area II (volatile- covered, non-interacting). As before, I first solve for temperatures away from the solid phase transition (T ≠ T C ). For Area IV, I integrate the conduction equation (Eq. 3.1-2) over the top layer, average that over Area IV, and add the result to Eq. 5.1-6a to replace the term with conduction at z = 0 (at the slab-substrate interface) with one at −Δ 0 (at the bottom of the first substrate layer). Taking the time average from time n to n+1 (indicated by overbars) yields Eq. 5.3-1, the areal averaged equivalent to Eq. 5.3-1 has areal averages for the thermophysical parameters (density, specific heat, mass of a slab, thermal conduction, emissivity), areal averages for the solar gain and heat lost by escape, and the inclusion of f V , the fraction of the interacting area that is covered by volatiles, in the latent heat and escape terms.
where ρ 0 c 0 is the areal average of the product of density and specific heat in layer 0, with cgs units of erg K -1 cm -3 , and m V c V is the areal average of the product of volatile slab mass and specific heat in the volatile slab, with cgs units of erg K -1 cm -2 .
The treatment of the first term is similar to that in the bare case; see the discussion near Eq. 3.3-2. As before, the temperature of layer 0 is sampled at the top of the layer. Because this is the slab-substrate interface, the temperature of layer 0 equals the volatile slab temperature within Area IV: T l,0,n = T n V . With the assumption that we can sample the temperature at the top of layer 0, the enthalpy term depends only on the change in the volatile slab temperature: , has units of erg cm -2 s -1 K -1 , with the superscript H representing heat or enthalpy. The discrete form for the areal average (cf.  is simply the weighted average of the local values, summed over the locations within Area IV, The weights  are simply the ratio of the solid angle of each location ( Ω l ) to the total solid angle of Area IV: Continuing to treat Eq. 5.3-1 term-by-term, the change enthalpy of the volatile slab also depends on the change in volatile slab temperature; see discussion near .
The latent heat term is the same over all locations, but differs from the local equivalents  by the factor of f V : The insolation terms is simply the areal average of the insolation at each location in Area IV: The thermal emission depends on the areal average of the emissivity: For explicit equations, the expression for the areal average of thermal conduction is found by taking the areal average of Eq. 3.3-9, and making the substitution that T l,0,n = T n V : where Φ l,0 K,B is given by . Similarly, the expression for the implicit (Crank-Nicholson) equations takes the areal average of Eq. 3.3-13: Finally, the escape rate is calculated by the average over all the interacting regions, Area III and Area IV: Substituting the expressions for the explicit equations gives Collecting terms for the explicit equation gives As in Section 3.3 and 4.3, divide by Φ n T , with cgs units erg cm -2 s -1 K -1 , where the superscript T represents total. The total "flux-per-temperature" includes terms for enthalpy of the slab and interaction with the atmosphere The resulting of dividing Eq. 5.3-13 by 5.3-14, and the resulting matrix elemens, are given in Table 8.
The implicit (Crank-Nicholson) equation (5.3-15) differs from equation (5.3-12) only with the substitution of the conduction term: Collecting terms for the implicit equation gives 5.3-16) This equation is used to derive the elements for the matrix elements in Figs. 5-2 and 5-3, given in Table 8.

Matrix equation Matrix elements Explicit
is given by 5.3-14.
The discrete form of the equation for the change in temperature at a crystalline phase is trivial, since the volatile slab temperature does not change from time n to time n+1 in Eq. 5.1-6b.

5.4a. Overview and explicit timesteps
In Section 3.4, I divided the temperature into the upper layer and the interior layers (Eq. 3.4-1), as a means to speeding up calculations in Areas I, II and III. In Area IV, this division is required, as the temperature of each of the upper layers ( The matrix elements € η n V and € γ n V are defined in Table 8 or 9. The b arrays are similar to Eq. (3.4-3), except that the weighting factor is included: The new temperature of the volatile slab depends on the substrate ; this is similar to Eq. 3.4-8a, but slightly simpler. The multi-location matrix equation for the temperatures of the substrate  is also similar to the non-interacting equivalent , differing only in that the topmost substrate temperature equals the temperature of the volatile slab.
Graphically, this is represented by Fig 5-5.  . Elements are labeled as in Fig 5-5. "*" indicates scalar multiplication (above the dotted line) or element-by-element multiplication of two arrays (above and below the dotted line). "X" indicates matrix multiplication.

5.4b. Implicit timesteps
For the implicit case, it is most straight-forward to write the Crank-Nicholson scheme in terms of intermediate temperatures for the volatile slab  T n V and substrate,  T l,1.. J,n .
For the other areas, the banded tridiagonal matrix was a computational convenience. For Area IV, it is the most direct way of solving . The solution to  can be written by defining two column vectors y l and z l of length J (defined as in 3.4-16a, 16b), and two scalars c n and d n : with which the temperatures at the next time step for location l are This solution can be confirmed by direct substitution into . The solution is shown graphically in Fig 5-6. Note that only the time-independent substrate matrix needs to be inverted, and this can be done at the start of the computation, rather than for each time step. Furthermore, the array y is also independent of time.
y 0 Fig 5-6. Graphical schematic of the implementation of an implicit time-step from time n to time n+1 for multiple non-interacting locations .
Elements are labeled as in Fig 3-9. "*" indicates element-by-element multiplication of two arrays (above the dotted line). "X" indicates matrix multiplication (equivalent to the outer product of two arrays for the multiplication below the lowest dotted line).

Example: PNV9 from Young 2013
As a worked example, Fig 5-7 shows the results of the calculations used for case PNV9 (permanent northern volatile #9) from Young 2013. This example was illustrated in Fig 1 of Young (2013) and Fig 3 of Olkin et al. (2015). The format of the figure is a still from the movies that show the seasonal evolution, as shown in various talks (e.g., Young 2012a). The code is included in the supplemental materials as vty16_fig5_7. The code included here is taken from the code actually run for Young (2013), with only superficial changes, to allow myself or others to reproduce the results of Young (2013) and Olkin et al. (2015). ε V ), the Bond albedo and emissivity of the substrate (A S and ε S ), the thermal inertia of the substrate (Γ), and the globally averaged N 2 inventory (N 2 ). Top right: Pluto's temperature and volatile mass for the listed year (2014.6). The subsolar latitude and heliocentric radius are listed (48° and 32.7 AU). The purple line gives an indication of the direction and magnitude of the sublimation winds, running from the North to the South. Blue indicates the volatile mass, where volatiles are present; the thickness of the bars are proportional to the mass, and the maximum mass is indicated (66.6 g/cm 2 ). The thin solid line indicates the surface temperature, which is a uniform 39.0 K for volatile covered areas, and is just above 40 K for bare areas (south of ~20°). Top left: Pluto as seen from the sun. Volatiles and substrate are shaded by their respective albedos. Pluto is tilted by the subsolar latitude. When plotted as a movie, the size of Pluto varies with the inverse of the heliocentric distance. Bottom left: graphical depiction of the seasonal volatile evolution. The shape of the orbit is in scale with Pluto's eccentricity. The 12 light and dark gray "pie pieces" mark out equal durations in the orbit, with the sun at the vertex of the pie pieces. The circles represent Pluto as seen with a zero sub-observer latitude. The pole is a squat bar running behind the circles. The circles and the pole bar are oriented so that the pole is perpendicular to the Pluto-sun line at the two equinoxes and so the summer hemisphere is oriented toward the sun. The red line and the circle outlined in red represent Pluto's position and state at the listed year (2014.6). Within the circles, lighter gray shows the location of volatiles, and darker gray shows substrate. Lower right: 1 vty16_fig5_7, which calls res = pluto_mssearch_func(run, av, ev, as, es, ti, mvtot, n_off, res_all) vty16_plutostill_mssearch_func, run, res, yr_still vty16_pluto_mssearch_resub_mat, flag_frostslab,time_delta,n_loc,n_z,emis,temp_surf,eflux_sol,mass_slab,specheat_frost,z_delta,z_delta_bot,dens,specheat,thermcond,eflux_int,beta,alpha_bot,alpha_mid vty16_pluto_mssearch_resub_timestep,flag_atm,freq,time_delta,gravacc,name_species,flag_stepscheme,n_loc,lat,n_z,z_delta,z_delta_bot,is_xport,angarea_delta,emis,temp_surf,eflux_sol,mass_slab,specheat_volatile,dens,specheat,thermcond,alpha_top,alpha_mid,alpha_bot,beta,denom,eflux_net,temp,temp_volatile,temp_next,temp_volatile_next,angarea_atm,mflux vty16_plutowrite_mssearch_func,run,res_all Volatile Transport II (VT3D) Surface pressure (log scale) and geometric albedo (linear) as a function of year, with the listed year marked by a small red circle. The pressure at the listed year is indicated (33 µbar), and the years of the equinoxes and solstices are marked.
In order to relate skin depth to depth with physical units, the substrate is assumed to have a density, ρ, of 0.93 g cm -3 and the skin depth, Z, is assumed to be 15 m; from the specified thermal inertia, Γ, Eq. 3.2-5 and 3.2-6 define the specific heat, c, and the thermal conductivity, k. The specific heat of the volatile, c V , (I remind the reader that the V is for volatile, not volume) is assumed to be that of N 2 (β), or 1.3e7 erg g -1 K -1 (Spencer & Moore 1992).
The run is initialized with the entire surface of Pluto covered with N 2 at aphelion, and the initial surface and subsurface temperatures are calculated assuming that the entire surface was volatile-covered and interacting over the previous Pluto year. The solar forcing is calculated assuming orbital elements of eccentricity of 0.254, inclination of 23.439°, Longitude of Ascending Node of 43.960°, argument of perifocus of 183.994°, last periapsis at Julian date 2447899.597, mean motion of 0.00392581°/day, a semi-major axis of 39.79700 AU, and a pole with right ascension 132.993° and declination -6.163° (see code for full precision). The diurnally averaged absorbed solar flux was calculated at 240 time steps over Pluto's orbit, at each of 60 latitude bands, and expanded to M = 2 (constant and two sinusoidal terms). The initial temperature field is calculated from the sinusoidal expansion of the absorbed solar flux, assuming a flux from the interior of 6 erg cm -2 s -1 . This follows the prescription of Section 5.2, except that the atmospheric "breathing" term is ignored (it is small on the seasonal timescales, Young 2013). The substrate uses a "medium" grid, with 19 layers of width 0.4 Z, where Z is the skin depth. The top layer is half that, or 0.2 Z.

Conclusions
A variety of mathematical techniques for speeding up thermophysical models or volatile transport models have been presented. They include an improved initial condition, an implicit time-step step scheme, and a matrix formulation that allows for the calculation of several locations at once. These can be used separately or in combination.
This formulation described here has been previously applied to Pluto's diurnal cycle with volatile distributions and albedos that vary with both latitude and longitude (Young 2012a). The speed gains allowed me to perform a wide parameter-space search of Pluto's seasonal cycle in anticipation of New Horizons (Young 2013). This work has also been used to study KBO seasons (Young and McKinnon 2013) and the first Pluto volatile transport models to include an N 2 reservoir .
Other routines are taken from elsewhere in the layoung IDL library, and from the astron library.