1. Introduction {#se0010}
===============

Jupiter\'s moon, Europa, and Saturn\'s moon, Enceladus, are two of only a few moons in the solar system that bear the possibility of extraterrestrial life (e.g., [@br0090], [@br0270]). Europa has a deep (∼100 km) ocean that underlies an icy shell, more than several kilometers deep (e.g., [@br0070], [@br0060], [@br0280], [@br0170], [@br0160], [@br0250], [@br0390], [@br0320], [@br0420]), where chemical interactions at the rocky bottom of the ocean may enable the existence of a habitable environment (e.g., [@br0090], [@br0100], [@br0130], [@br0200]). The Voyager and Galileo (and to a lesser extent, the Cassini-Huygens and New Horizons) spacecrafts/missions discovered many interesting features of Europa including chaos terrains ([@br0330], [@br0410]) and craters ([@br0190], [@br0230], [@br0120], [@br0340]). More recently, based on the Hubble telescope observations, scientists raised the possibility of water vapor plumes at Europa\'s south pole ([@br0310], [@br0350]). Europa is one of the youngest, largest, and brightest moons in the solar system ([@br0290]). A basic property of Europa is its surface temperature; surface temperature is needed to calculate the properties of its icy shell and ocean dynamics. Thus, an accurate estimation of Europa\'s surface temperature is required.

The annual mean surface temperature of Europa was previously estimated by [@br0260]. These authors took into account the obliquity of Europa with respect to the plane of rotation of Jupiter around the Sun. Using a surface albedo of 0.5, they found that the annual mean temperature varies from ∼110 K at the equator to ∼52 K at the poles; the global (and annual) mean surface temperature was found to be ∼100 K (see the gray curve in [Fig. 2](#fg0020){ref-type="fig"}d below). However, the following factors were not taken into account when calculating the annual mean surface temperature of Europa: the eccentricity of Jupiter\'s orbit around the Sun, the effect of the emissivity of Europa\'s ice, the heat capacity of the surface ice, the effect of Europa\'s eclipse (i.e., the time that Europa is within the shadow of Jupiter), the longwave radiation of Jupiter that is absorbed by Europa, and Europa\'s internal heating. The eccentricity \[see Eqs. [(18)](#fm0180){ref-type="disp-formula"} and [(21)](#fm0210){ref-type="disp-formula"} below\], emissivity, Jupiter\'s radiation, and internal heating factors increase the incoming radiation to Europa\'s surface from above and below, while the eclipse factor reduces the absorbed incoming solar radiation. In addition, an updated surface Bond albedo of $0.68 \pm 0.05$ ([@br0140]) that is based on New Horizons measurements should be considered; this by itself reduces the absorbed incoming solar radiation by more than 20%. Moreover, to our knowledge, the seasonal cycle of the surface temperature of Europa has not been previously estimated. The goal of this study is to develop a more accurate estimation of the diurnal, seasonal and annual mean surface temperature of Europa, taking all the above factors into account; here, we systematically investigate the role of the different parameters on the surface temperature of Europa. In addition, the analytic approximation developed here may be used to estimate the surface temperature of other moons in the solar system.

The brightness temperature of Europa was measured by [@br0370] and [@br0300], based on measurements performed by the Galileo spacecraft. They measured the diurnal temperature cycle at the low latitudes to be between 86 K and 132 K and also provided spatial snapshots of surface temperatures up to, roughly, 70^∘^ latitude. Still, the diurnal mean temperature for latitudes poleward of 15^∘^ and the temperatures at the high latitudes were not measured. The polar region temperatures are important as in these regions, the internal heating significantly affects the surface temperatures, and measurements of surface temperatures in these regions may help to estimate the internal heating rate and, consequently, the thickness of the ice (see, e.g., [@br0020]). The dependence of surface temperature on the internal heating is studied here.

Below we first discuss the diurnal, seasonal, and annual mean incoming solar radiation to Europa (Section [2.1](#se0020){ref-type="sec"}). Some of Jupiter\'s longwave radiation is absorbed by Europa, and this effect is quantified in Section [2.4](#se0070){ref-type="sec"}. We then quantify the effect of Europa\'s eclipse, i.e., the relative time that Europa is within the shadow of Jupiter (Section [2.5](#se0080){ref-type="sec"}). Next we calculate the surface temperature of Europa (Section [3.1](#se0100){ref-type="sec"}). A summary and discussion close the paper (Section [4](#se0120){ref-type="sec"}). A very rough estimation of the mean thickness of Europa\'s icy shell as a function of the internal heating rate is detailed in Section [3.2](#se0110){ref-type="sec"}. The parameters that are used in this study are listed in [Table 1](#tbl0010){ref-type="table"}.Table 1List of parameters.Table 1ParameterDescriptionValue*S*~0~Jupiter solar constant51 W m^−2^*e*eccentricity of Jupiter0.048*ε*obliquity of Europa3^∘^*α*~*p*~bolometric Bond albedo of Europa0.68 ± 0.05*p*Europa eclipse relative time0.033*σ*Stefan-Boltzmann constant5.67 × 10^−8^ W m^−2^ K^−4^*ϵ*emissivity of Europa0.94*ρ*~*I*~density of ice917 kg m^−3^*c*~*p*,*I*~heat capacity of ice2000 J kg^−1^ K^−1^*κ*deep ice heat diffusion constant1.54 × 10^−6^ m^2^ s^−1^*κ*~*s*~surface ice heat diffusion constant7.7 × 10^−10^ m^2^ s^−1^*g*surface gravity of Europa1.314 m s^−2^*r*~*j*~mean Jupiter-Sun distance7.785 × 10^11^ m*r*~*e*~mean Europa-Jupiter distance6.71 × 10^8^ m*a*~*s*~radius of Sun6.96 × 10^8^ m*a*~*j*~radius of Jupiter6.99 × 10^7^ m*a*~*e*~radius of Europa1.561 × 10^6^ mΩ~*e*~rotational frequency of Europa2.05 × 10^−5^ s^−1^*ω*Jupiter\'s longitude of the perihelion14.7285^∘^*J*~0~Jupiter\'s radiation constant0.176 W m^−2^

2. Theory {#se0180}
=========

2.1. The incoming solar radiation {#se0020}
---------------------------------

Below, we present results regarding the diurnal cycle, seasonal cycle, and annual mean solar radiation reaching the surface of Europa. Details regarding the derivations of the daily mean insolation can be found in Section [2.2](#se0030){ref-type="sec"} and regarding the seasonal and annual mean insolation in Section [2.3](#se0040){ref-type="sec"}.

[Fig. 1](#fg0010){ref-type="fig"}a depicts the diurnal cycle of the incoming solar radiation at the equator during the northern hemisphere (NH) vernal (spring) equinox and NH summer solstice and at 89^∘^S during the southern hemisphere (SH) summer solstice. The incoming solar radiation peaks at the equator during the NH summer solstice and does not exceed the value of the solar constant of Europa/Jupiter, $S_{0} = 51\ \text{W}\ \text{m}^{- 2}$; the maximal diurnal solar radiation is $\sim S_{0}$ when the eccentricity is small \[see Eq. [(1)](#fm0010){ref-type="disp-formula"}\], as for Europa. The NH summer solstice radiation is higher than the NH vernal equinox radiation as Jupiter (and hence Europa) is closer to the Sun at the NH summer solstice such that the eccentricity effect overcomes the obliquity effect. In other words, at the present time, the solstices are close to the aphelion and perihelion (see [Fig. 2](#fg0020){ref-type="fig"}b below); this situation may be different at different times due to the long (many thousands of years) precession cycle of Jupiter. Close to the pole, even during the NH summer solstice (blue curve), the radiation is much smaller, as is the radiation range.Figure 1(a) Diurnal cycle of the incoming solar radiation at the equator during the equinox (*λ* = 0, NH vernal equinox, green) and solstice (*λ* = *π*/2, NH summer solstice, red), and at 89^∘^S during the solstice (*λ* = 3*π*/2, NH winter solstice, blue). (b) Diurnal cycle of the surface temperature at the equator during the equinox (*λ* = 0, NH vernal equinox, green) and the solstice (*λ* = *π*/2, NH summer solstice, red). The dotted and dashed horizontal lines indicate the daily mean surface temperature calculated using the depth-dependent model \[Eq. [(29)](#fm0290){ref-type="disp-formula"}\] (dotted) and the daily mean radiation \[Eq. [(32)](#fm0320){ref-type="disp-formula"}\] (dashed). (c) The temperature versus depth at the equator at the equinox (*λ* = 0, NH vernal equinox) at several times during the diurnal cycle. Note the fast decline in temperature oscillations with depth. (d) Same as b but for 89^∘^S, at the solstice (*λ* = 3*π*/2, NH winter solstice). The temperature was obtained by integrating the depth-dependent model \[Eq. [(29)](#fm0290){ref-type="disp-formula"}\] at a specific latitude.Figure 1Figure 2(a) Annual mean (blue and red) and NH summer (SH winter) solstice (yellow) incoming solar radiation as a function of latitude. Both the numerically integrated (blue) and the analytic (red) approximation \[Eqs. [(18)](#fm0180){ref-type="disp-formula"}, [(21)](#fm0210){ref-type="disp-formula"}\] are shown although the two are almost indistinguishable (the maximal difference between the two is 0.01 W m^−2^). (b) Seasonal variations of the surface temperature versus time (in Jupiter\'s years) for the equator (blue) and for the NH pole (red). The internal heating rate is *Q* = 0.05 W m^−2^ and the curves are the solution of Eq. [(29)](#fm0290){ref-type="disp-formula"}. The dashed horizontal lines indicate the annual mean temperature based on the annual mean insolation \[Eqs. [(18)](#fm0180){ref-type="disp-formula"}, [(21)](#fm0210){ref-type="disp-formula"}\] while the dotted lines indicate the annual mean temperature of the solid lines. *t* = 0 is set at the spring (vernal) equinox and the dotted vertical lines indicate the time of the aphelion and perihelion. (c) Surface temperature at the NH summer solstice for different internal heating rates. (d) Annual mean surface temperature for different internal heating rates. The gray line indicates the surface temperature as estimated by [@br0260] based on [@br0240].Figure 2

In [Fig. 2](#fg0020){ref-type="fig"}a, we plot the annual mean and NH summer solstice incoming solar radiation as a function of latitude. Both the numerically integrated and the analytically approximated annual mean curves are presented, and the two are almost indistinguishable. The numerical integration was based on the daily mean insolation \[Eqs. [(6)](#fm0060){ref-type="disp-formula"}, [(7)](#fm0070){ref-type="disp-formula"}\] using Eq. [(12)](#fm0120){ref-type="disp-formula"} and $\langle W\rangle = \frac{1}{\tau}\int_{0}^{2\pi}W\frac{dt}{d\lambda}d\lambda$ with $d\lambda = 0.01\ \text{rad}$. Note the very small level of incoming solar radiation ($\sim 1\ \text{W}\ \text{m}^{- 2}$) that reaches the poles, suggesting that the internal heating cannot be ignored at the high latitudes. No solar radiation reaches the winter pole during the corresponding winter solstice; thus, we expect the internal heating to have a larger influence on polar region surface temperature during this time.

2.2. Derivation of the daily mean insolation {#se0030}
--------------------------------------------

Below we provide the mathematical details regarding the derivations of the daily mean insolation. The reader that is not interested in the mathematical details may skip this section.

The incoming insolation *W* at a certain latitude *ϕ* and a certain time of the day on Europa ([@br0150]) is$$W = S_{0}\left( \frac{\overline{d}}{d} \right)^{2}\cos\theta_{s},$$ where $\theta_{s}$ is the solar zenith angle, $\overline{d}$ is the mean distance of Jupiter (Europa) from the Sun, *d* is the distance from the Sun, and $S_{0}$ is the solar constant of Jupiter. $\theta_{s}$ is the solar zenith angle that depends on latitude, season, and time of day and$$\frac{\overline{d}}{d} = \frac{1 + e\cos(\lambda - \omega - \pi)}{1 - e^{2}},$$ where *e* is the eccentricity of Jupiter, *ω* is the precession (longitude of the perihelion), and *λ* is the longitude of Jupiter with respect to its orbit around the Sun. Seasons are expressed through the declination angle *δ*, which is the latitude of the line connecting the center of Europa and the Sun during noontime. The hour angle *H* indicates the longitude of the subsolar point relative to its position at noon. Then,$$\cos\theta_{s} = \sin\phi\sin\delta + \cos\phi\cos\delta\cos H.$$ At sunrise and sunset, the zenith angle is $90^{\circ}$ such that $\cos\theta_{s} = 0$ and$$\cos H_{0} = - \tan\phi\tan\delta.$$ The declination angle *δ* is$$\sin\delta = \sin\varepsilon\sin\lambda,$$ where *ε* is the obliquity of Europa with respect to the plane of rotation of Jupiter around the Sun; the maximal obliquity can be calculated as the sum of the axial tilt of Europa, 0.1^∘^, the inclination angle of Europa, 0.47^∘^, and the axial tilt of Jupiter, 3.13^∘^, yielding an angle of 3.7^∘^. Yet, since the axes of rotation of both Europa and Jupiter exhibit precession, the above angles can either add up or subtract. We thus estimated the mean present day obliquity as the subsolar latitude at the solstices based on the JPL Horizons web-interface (<https://ssd.jpl.nasa.gov/horizons.cgi#top>); it is $\varepsilon \approx 3^{\circ}$. We performed the calculations presented in this paper also using $\varepsilon = 3.7^{\circ}$ and obtained almost identical results for the low latitudes and global mean values, while at the polar regions, the temperature was higher by a few degrees for $\varepsilon = 3.7^{\circ}$; see [Table 2](#tbl0020){ref-type="table"}. At the solstice, $\sin\lambda = \pm 1$ such that $\delta = \pm \varepsilon$, while at the equinox, $\delta = 0$.Table 2Sensitivity of surface temperature (K) of Europa to the different parameters when using the numerical model \[Eq. [(29)](#fm0290){ref-type="disp-formula"}\]. Control run parameters: eccentricity (*e* = 0.048), obliquity (*ε* = 3^∘^), emissivity (*ϵ* = 0.94), albedo (*α*~*p*~ = 0.68), Europa eclipse relative time (*p* = 0.033), surface ice heat diffusion constant (*κ*~*s*~ = 7.7 × 10^−10^ m^2^ s^−1^), Jupiter\'s longwave radiation constant of Europa (*J*~0~ = 0.176 W m^−2^), and internal heating (*Q* = 0.05 W m^−2^).Table 2ParameterEquatorPoleGlobal mean15^∘^S--15^∘^N meancontrol96.146.190.195.9*e* = 095.845.989.995.6*ε* = 096.231.290.195.9*ε* = 3.7^∘^96.147.990.195.9*ϵ* = 194.645.489.495.1*α*~*p*~ = 0.5105.849.794.8100.9*p* = 096.846.490.596.2*κ*~*s*~ = 7.7 × 10^−8^ m^2^ s^−1^99.349.091.497.5*J*~0~ = 095.846.190.095.7*Q* = 095.842.889.995.7*Q* = 0.2 W m^−2^97.052.890.796.3

2.3. Derivation of the annual mean insolation {#se0040}
---------------------------------------------

Below we provide the mathematical details regarding the derivations of the annual mean insolation. The reader that is not interested in the mathematical details may skip this section.

### 2.3.1. Small obliquity {#se0050}

The daily mean insolation ([@br0220], [@br0040], [@br0050], [@br0150], [@br0010]) can be found by integrating Eq. [(1)](#fm0010){ref-type="disp-formula"} from $- H_{0}$ to $H_{0}$ and dividing by 2*π* and is given by$$W = \frac{S_{0}}{\pi}\frac{{\lbrack 1 + e\cos(\lambda - \omega - \pi)\rbrack}^{2}}{{(1 - e^{2})}^{2}}(H_{0}\sin\phi\sin\delta + \cos\phi\cos\delta\sin H_{0}),$$ where $H_{0}$ is the hour angle at sunrise and sunset \[Eq. [(4)](#fm0040){ref-type="disp-formula"}\]; see [Table 1](#tbl0010){ref-type="table"}. Equation [(6)](#fm0060){ref-type="disp-formula"} is relevant for latitudes $|\phi| < |\pi/2 - |\delta||$, while outside this latitude range, the insolation during the polar night is zero and during the polar day is equal to$$W = S_{0}\frac{{\lbrack 1 + e\cos(\lambda - \omega - \pi)\rbrack}^{2}}{{(1 - e^{2})}^{2}}\sin\phi\sin\delta.$$

Below we will use the following mathematical relations:$$\tan\delta = \frac{\sin\delta}{\sqrt{1 - \sin^{2}\delta}} = \frac{\sin\varepsilon\sin\lambda}{\sqrt{1 - \sin^{2}\varepsilon\sin^{2}\lambda}}.$$ Since the obliquity is small, one can perform the following approximations$$\tan\delta \approx \sin\delta \approx \delta \approx \varepsilon\sin\lambda$$$$\sin H_{0} = \sqrt{1 - \tan^{2}\phi\tan^{2}\delta} \approx \sqrt{1 - \varepsilon^{2}\tan^{2}\phi\sin^{2}\lambda}$$ and, based on Eqs. [(4)](#fm0040){ref-type="disp-formula"} and [(9)](#fm0090){ref-type="disp-formula"}$$H_{0} = \arccos( - \varepsilon\tan\phi\sin\lambda).$$

To calculate the annual mean insolation, one must take into account the derivative of the time, *t*, with respect to the longitude, *λ* ([@br0220], [@br0150], [@br0010])$$\frac{dt}{d\lambda} = \frac{{(1 - e^{2})}^{3/2}}{{\lbrack 1 + e\cos(\lambda - \omega - \pi)\rbrack}^{2}}.$$ Then, the relative time, *τ*, of one cycle of Jupiter around the Sun may be approximated as$$\tau = \int\limits_{0}^{2\pi}\frac{dt}{d\lambda}d\lambda \approx 2\pi{(1 - e^{2})}^{3/2},$$ as the eccentricity is much smaller than one. The relative time, *τ*, is the ratio of the time passed since the NH summer solstice ($\lambda = 0$) to the time it takes for Jupiter to complete one cycle around the Sun.

The total insolation during the year is then$$\int\limits_{0}^{2\pi}W\frac{dt}{d\lambda}d\lambda = \frac{S_{0}}{\pi\sqrt{1 - e^{2}}}\int\limits_{0}^{2\pi}(H_{0}\sin\phi\sin\delta + \cos\phi\cos\delta\sin H_{0})d\lambda.$$ We solve each part of the integral separately. The first part is$$\sin\phi\int\limits_{0}^{2\pi}H_{0}\sin\delta d\lambda = \varepsilon\sin\phi\int\limits_{0}^{2\pi}\arccos( - \varepsilon\tan\phi\sin\lambda)\sin\lambda d\lambda\quad \approx \varepsilon\sin\phi\int\limits_{0}^{2\pi}\left( \frac{\pi}{2} + \varepsilon\tan\phi\sin\lambda \right)\sin\lambda d\lambda = \varepsilon^{2}\tan\phi\sin\phi\int\limits_{0}^{2\pi}\sin^{2}\lambda d\lambda\quad = \pi\varepsilon^{2}\sin\phi\tan\phi.$$

The second part of the integral of Eq. [(14)](#fm0140){ref-type="disp-formula"} is$$\cos\phi\int\limits_{0}^{2\pi}\cos\delta\sin H_{0}d\lambda = \cos\phi\int\limits_{0}^{2\pi}\sqrt{1 - \sin^{2}\varepsilon\sin^{2}\lambda}\sqrt{1 - \varepsilon^{2}\tan^{2}\phi\sin^{2}\lambda}d\lambda \approx \cos\phi\int\limits_{0}^{2\pi}\left( 1 - \frac{1}{2}\sin^{2}\varepsilon\sin^{2}\lambda \right)\left( 1 - \frac{1}{2}\varepsilon^{2}\tan^{2}\phi\sin^{2}\lambda \right)d\lambda \approx \cos\phi\int\limits_{0}^{2\pi}\left( 1 - \frac{1}{2}\sin^{2}\lambda(\sin^{2}\varepsilon + \varepsilon^{2}\tan^{2}\phi) \right)d\lambda = \frac{1}{2}\pi\cos\phi(4 - \sin^{2}\varepsilon - \varepsilon^{2}\tan^{2}\phi).$$

Following the above, the total insolation during the year can be approximated as:$$\int\limits_{0}^{2\pi}W\frac{dt}{d\lambda}d\lambda \approx \frac{S_{0}}{2\sqrt{1 - e^{2}}}(4\cos\phi + \varepsilon^{2}\sin\phi\tan\phi - \cos\phi\sin^{2}\varepsilon).$$ Finally, the annual mean insolation outside the polar regions is$$\langle W\rangle = \frac{1}{\tau}\int\limits_{0}^{2\pi}W\frac{dt}{d\lambda}d\lambda \approx \frac{S_{0}}{4\pi{(1 - e^{2})}^{2}}\left( 4\cos\phi + \varepsilon^{2}\sin\phi\tan\phi - \cos\phi\sin^{2}\varepsilon \right).$$

To find the annual mean insolation in the polar region, we first find the total annual insolation at the pole ($\phi = \pi/2$),$$\int\limits_{0}^{2\pi}W\frac{dt}{d\lambda}d\lambda = \frac{S_{0}\sin\varepsilon}{\sqrt{1 - e^{2}}}\int\limits_{0}^{\pi}\sin\lambda d\lambda = \frac{2S_{0}\sin\varepsilon}{\sqrt{1 - e^{2}}}.$$ Then, using Eq. [(13)](#fm0130){ref-type="disp-formula"}, we find the annual mean isolation at the poles:$$\langle W(\phi = \pm \pi/2)\rangle = \frac{S_{0}\sin\varepsilon}{\pi{(1 - e^{2})}^{2}}.$$ Then we match the above polar annual mean insolation to the annual mean insolation at the edge of the polar region (i.e., $|\phi| = \frac{\pi}{2} - \varepsilon$) to obtain the annual mean insolation at the polar regions$${\langle W\rangle}_{p} \approx \frac{S_{0}}{\pi{(1 - e^{2})}^{2}}\left( \sin\varepsilon + \frac{\cos^{2}\phi}{4\sin\varepsilon} \right).$$

### 2.3.2. Large obliquity -- the case of Enceladus {#se0060}

The approximation we developed above is based on the assumption that the obliquity, *ε*, is small. This situation, however, is not the case, for example, for Enceladus, the moon of Saturn, whose obliquity with respect to the orbit around the Sun is 27^∘^. Enceladus is also a moon with an underlying ocean ([@br0380]), and, like Europa, is one of the most probable places in the solar system to find extraterrestrial life. Surprisingly, our approximation for the low and mid-latitudes \[Eq. [(18)](#fm0180){ref-type="disp-formula"}\] holds also for the relatively large obliquity of Enceladus, while the solar radiation in the polar regions is better approximated using$${\langle W\rangle}_{p} \approx \frac{S_{0}}{\pi{(1 - e^{2})}^{2}}\left( \sin\varepsilon + \frac{\cos^{2}\phi}{2\pi\sin\varepsilon} \right).$$ [Fig. 3](#fg0030){ref-type="fig"} depicts the numerically integrated and analytically approximated annual mean \[Eq. [(22)](#fm0220){ref-type="disp-formula"}\] solar radiation for Enceladus; for reference, we also include the solar radiation during the solstice.Figure 3The incoming solar radiation as a function of latitude for Enceladus, the moon of Saturn. Both the numerically integrated and the analytic approximation are shown, and the two are almost indistinguishable except for the region around latitudes −*π*/2 + *ε* and *π*/2 − *ε* (63^∘^S and 63^∘^N), at the transition to the polar regions. The parameter values we used to calculate the incoming solar radiation to Enceladus are: *S*~0~ = 15.14 W m^−2^, *e* = 0.055, *ε* = 27^∘^, and *α*~*p*~ = 0.8. For simplicity, to estimate the incoming solar radiation during the solstice, we assume that *e* = 0.Figure 3

2.4. The effect of Jupiter\'s longwave radiation {#se0070}
------------------------------------------------

Jupiter emits longwave radiation that may affect the surface temperature of Europa. Europa is phased-locked to Jupiter such that one side of Europa always absorbs the longwave radiation of Jupiter while the other side never absorbs this radiation. Below, we show that Jupiter\'s longwave radiation only slightly affects the surface temperature of Europa, mainly at the low latitudes -- the effect of this radiation is very small at the high latitudes of Europa. For simplicity, we assume that: (i) the axes of rotation of Jupiter and Europa are parallel (in reality, the difference between the two is about $0.5^{\circ}$), (ii) the eccentricity of Europa around Jupiter is 0 (where in reality is 0.009), (iii) Europa absorbs all longwave radiation of Jupiter reaching its surface without any reflection, (iv) the longwave radiation of Jupiter is distributed evenly during Europa\'s diurnal cycle (i.e., the Jupiter-facing side of Europa receives as much radiation as the side that never faces Jupiter), and (v) ignore the spatial and temporal variations in the emission temperature (and hence longwave radiation) of Jupiter. Given the assumptions above (especially assumption iv), the analysis below is not exact. We note that the justification for assumption (iv) is our aim to calculate the zonal mean temperature.

Consider [Fig. 4](#fg0040){ref-type="fig"} below. One can view the total radiation emitted from Jupiter as emitted from a point source concentrated at the center of Jupiter; this assumption is valid if one assumes that the radiation emitted from Jupiter\'s surface is perpendicular to the surface. This "source" emits radiation that propagates spherically. The radiation decreases as $1/r^{2}$ where *r* is the distance from the center of Jupiter. This radiation at Europa is $J_{0} = (4\pi a_{j}^{2}\sigma T_{j}^{4})/(4\pi r_{e}^{2}) = 0.176\ \text{W}\ \text{m}^{- 2}$. Here, $a_{j}$ is the radius of Jupiter, *σ* is the Stephan-Boltzmann constant, $T_{j} = 130\ \text{K}$ is the emission temperature of Jupiter ([@br0210]), and $r_{e}$ is the distance between Jupiter and Europa. The maximum of Jupiter\'s longwave radiation that is absorbed by Europa at a certain latitude $\phi_{e}$ is $J_{0}\cos(\phi_{e} + \beta)$, where *β* is indicated in [Fig. 4](#fg0040){ref-type="fig"}. Then, by calculating the average Jovian radiation over all longitudes at a given latitude, the mean Jovian radiation over a diurnal cycle is $J_{0}\cos(\phi_{e} + \beta)/\pi$ where *β* can be approximated as $\tan\beta = \frac{a_{e}\sin\phi_{e}}{r_{e} - a_{e}\cos\phi_{e}} \approx \frac{a_{e}}{r_{e}}\sin\phi_{e} \approx 0$ since $r_{e} \approx 430a_{e}$. Thus, the daily mean radiation absorbed by Europa at a certain latitude is, to a good approximation,$$W_{j} = \frac{J_{0}}{\pi}\cos\phi_{e}.$$ It follows that this radiation is zero at the poles and maximal at the equator where it is comparable ($0.056\ \text{W}\ \text{m}^{- 2}$) to the contribution of the internal heating *Q*. Since the daily mean absorbed solar radiation at the equator is $W(1 - \alpha_{p})(1 - p) \approx 5\ \text{W}\ \text{m}^{- 2}$, the effect of Jupiter\'s longwave radiation on the surface temperature of Europa is small.Figure 4A drawing showing the longwave radiation emitted from a point source of Jupiter, reaching Europa. The relative dimensions are not realistic as the distance between Jupiter and Europa is 10 times the radius of Jupiter and the radius of Europa is 45 times smaller than the radius of Jupiter. The inset shows a scaled configuration of Jupiter (black circle) and Europa (red dot that is indicated by the arrow).Figure 4

2.5. The effect of Europa\'s eclipse {#se0080}
------------------------------------

Since the obliquity of both Jupiter and Europa is very small and since Europa is close to Jupiter (the distance between Europa and Jupiter is only ∼10 times the radius of Jupiter) but much smaller than Jupiter (the radius of Europa is about 45 times smaller than that of Jupiter), Europa passes in the shadow of Jupiter during each of Europa\'s days. Since Europa is phased-locked to Jupiter, only the side of Europa that faces Jupiter will experience the eclipse -- below, when simulating the diurnal cycle in temperature, we ignore this fact and assume that the eclipse reduces the overall absorbed daily solar radiation by some factor. The relative time that Europa is within the shadow of Jupiter (eclipse conditions) can be approximated as the ratio between the diameter of Jupiter and the perimeter of Europa\'s orbit around Jupiter; it is $p = 2a_{j}/2\pi r_{e} \approx 0.033$, where $a_{j}$ is the radius of Jupiter and $r_{e}$ is the distance of Europa from Jupiter. Here we consider only the umbra effect and ignore the penumbra effect since (i) Jupiter is relatively very far from the Sun and since (ii) Europa is relatively close to Jupiter. See Section [2.6](#se0090){ref-type="sec"} for a more accurate estimation of the eclipse effect.

It is possible to approximate the decrease in surface temperature at the end of the eclipse by considering only the internal heating and outgoing longwave radiation, using a time-dependent energy balance equation. In this case, the decrease in temperature is primarily controlled by the heat capacity of the surface ice and much less by the internal heating. The decrease in temperature is then less than 0.3 K.

2.6. More accurate estimation for the effect of the eclipse {#se0090}
-----------------------------------------------------------

Below we provide the details regarding the more accurate estimation for the effect of the eclipse. The reader that is not interested in the mathematical details may skip this section.

A more accurate estimation is based on [Fig. 5](#fg0050){ref-type="fig"}. The term $a_{s,j,e}$ indicates the radius of the Sun, Jupiter, and Europa, $r_{p}$ is the distance between the Sun and point *p*, $r_{j}$ is the distance of Jupiter from the Sun, $r_{e}$ is the distance between Europa and Jupiter, *d* is the segment of the perimeter of Europa around Jupiter that is in Jupiter\'s shadow, and 2*α* is the angle between two rays from the Sun that intercept at point *p*. Thus,$$\frac{d}{2} = a_{s}\left( 1 - \frac{r_{j}}{r_{p}} - \frac{r_{e}}{r_{p}} \right)$$$$\tan\alpha = \frac{a_{j}}{r_{p} - r_{j}}$$$$r_{p} = \frac{r_{j}a_{s}}{a_{s} - a_{j}},$$ which implies that$$d = 2a_{j} - 2r_{e}(a_{s} - a_{j})/r_{j}.$$ The relative time that Europa is in the shadow of Jupiter is approximately$$p = \frac{d - a_{e}}{2\pi r_{e}} = \frac{a_{j} - r_{e}(a_{s} - a_{j})/r_{j} - a_{e}/2}{\pi r_{e}}.$$ For the parameter values listed in [Table 1](#tbl0010){ref-type="table"}, $d \approx 2a_{j}$ and $p \approx a_{j}/(\pi r_{e}) = 0.033$, as indicated in Section [2.5](#se0080){ref-type="sec"}.Figure 5A drawing depicting the different measures that are used to calculate the time that Europa spends in Jupiter\'s shadow. The different measures do not reflect the real relative measures.Figure 5

3. Model {#se0190}
========

3.1. Calculation of the surface temperature of Europa {#se0100}
-----------------------------------------------------

We first focus on the diurnal cycle of the surface temperature. Here one has to take into account the heat capacity of the surface of Europa that leads to absorption of solar radiation only during the day and emission of longwave radiation during both day and night. For this purpose, it is possible to construct a simple diffusion equation for ice temperature, *T*, that is forced from below by the internal heat of Europa and from above by the incoming solar (and Jupiter) radiation and outgoing longwave radiation:$$\frac{\partial T}{\partial t} = \kappa_{s}\frac{\partial^{2}T}{\partial z^{2}}$$ where $\kappa_{s}$ is the surface ice temperature diffusion coefficient. We note that the surface ice temperature diffusion constant $\kappa_{s}$ is 2--3 orders of magnitude smaller than the value of Earth; this follows previous studies ([@br0370], [@br0300]) that estimated the thermal inertia of the surface ice of Europa to be tens of times smaller than the value of Earth; the thermal inertia is $\Gamma = \sqrt{k_{s}\rho_{I}c_{p,I}}$ where $k_{s}$ is the thermal conductivity constant which is related to the ice surface diffusion coefficient $\kappa_{s}$ through $k_{s} = \rho_{I}c_{p,I}\kappa_{s}$. The value of $\kappa_{s}$ was chosen to fit the diurnal variations of surface temperature reported in [@br0300]. The surface boundary condition is$$\rho_{I}c_{p,I}\kappa_{s}\frac{\partial T_{s}}{\partial z} = W(1 - \alpha_{p})(1 - p) + W_{j} - \epsilon\sigma T_{s}^{4},$$ where $\rho_{I}$ is the density of the ice, $c_{p,I}$ is the heat capacity of the ice, $T_{s}(t) = T(z = 0,t)$ is Europa\'s surface temperature, *W* is the diurnal cycle of the insolation \[Eq. [(1)](#fm0010){ref-type="disp-formula"}\], $W_{j}$ is the daily mean absorbed Jupiter\'s longwave radiation \[Eq. [(23)](#fm0230){ref-type="disp-formula"}\], $\alpha_{p}$ is the planetary albedo of Europa, *p* is the relative time that Europa passes through Jupiter\'s shadow, and *ϵ* is the emissivity of Europa ([@br0360]). The boundary condition at the bottom of the surface ice layer is$$\rho_{I}c_{p,I}\kappa_{s}\frac{\partial T}{\partial z} = - Q,$$ meaning that the internal heating rate *Q* is proportional to the temperature gradient deep enough within the ice.

Under the assumptions that Europa has no atmosphere and that there is no heat capacity at Europa\'s surface (i.e., $c_{p,I} = 0$), it is possible to calculate the daily mean surface temperature of Europa based on the energy balance between the incoming heat fluxes (i.e., shortwave solar radiation, incoming longwave radiation of Jupiter, and internal heating) and the outgoing longwave radiation. More specifically,$$W(1 - \alpha_{p})(1 - p) + W_{j} + Q = \epsilon\sigma T_{s}^{4},$$ where *W* is the daily mean insolation \[Eqs. [(6)](#fm0060){ref-type="disp-formula"}, [(7)](#fm0070){ref-type="disp-formula"}\], and the other parameters are explained above. Consequently, Europa\'s surface temperature can be expressed as$$T_{s}(\phi) = \left\lbrack \frac{W(1 - \alpha_{p})(1 - p) + W_{j} + Q}{\epsilon\sigma} \right\rbrack^{\frac{1}{4}}.$$ The annual mean surface temperature, $\langle T_{s}\rangle$, can be calculated using Eqs. [(33)](#fm0330){ref-type="disp-formula"}, [(12)](#fm0120){ref-type="disp-formula"}, [(13)](#fm0130){ref-type="disp-formula"} and is$$\langle T_{s}\rangle = \frac{1}{\tau}\int\limits_{0}^{2\pi}T_{s}\frac{dt}{d\lambda}d\lambda.$$ All the parameters in Eq. [(33)](#fm0330){ref-type="disp-formula"} are well constrained except the internal heating *Q*. The internal heating, most probably, has spatial dependence, due to, for example, tidal heating within the ice. However, due to a lack of knowledge and large uncertainties, we assume that it is spatially constant. In addition, [@br0390] indicated that the heat flux at Europa\'s surface is almost constant. \[We note that it is easily possible to use spatially variable internal heating and surface albedo in Eq. [(33)](#fm0330){ref-type="disp-formula"}.\]

A more common way to estimate the surface temperature is using the annual mean insolation, i.e.,$${\langle T_{s}^{4}\rangle}^{1/4} = \left\lbrack \frac{\langle W\rangle(1 - \alpha_{p})(1 - p) + W_{j} + Q}{\epsilon\sigma} \right\rbrack^{\frac{1}{4}},$$ where $\langle W\rangle$ is the annual mean insolation given in Eq. [(18)](#fm0180){ref-type="disp-formula"}. Such an approach was taken by [@br0260]. Still, this is a non-trivial assumption as $\langle T_{s}\rangle \neq {\langle T_{s}^{4}\rangle}^{1/4}$. We discuss the validity of Eq. [(35)](#fm0350){ref-type="disp-formula"} below.

To find the surface temperature, it is necessary to numerically integrate Eq. [(29)](#fm0290){ref-type="disp-formula"} with respect to depth and time using the surface and bottom boundary conditions \[Eqs. [(30)](#fm0300){ref-type="disp-formula"}, [(31)](#fm0310){ref-type="disp-formula"}\]. The surface layer is subjected to diurnally and seasonally varying incoming solar radiation where a Jupiter year is equivalent to ≈1220 Europa days. It is thus necessary to use a sufficiently fine temporal and vertical resolution with a long enough integration time.

We used 3500 vertical levels with a resolution of 0.87 mm, covering a depth of 3.04 m. The effect of the diurnal surface oscillations (due to the diurnal variations of solar radiation) decay within the upper 5 cm while the seasonal surface variations decay within the upper 1 m. The integration time step (in terms of hour angle) is $10^{- 3}\ \text{rad}$. The initial conditions were uniform temperature (60 K) with depth, and the model \[Eq. [(29)](#fm0290){ref-type="disp-formula"}\] was integrated for 150 Jupiter years -- this integration time was long enough to achieve convergence.

We first discuss the diurnal cycle of the surface temperature of Europa. The results are summarized in [Fig. 1](#fg0010){ref-type="fig"}. In [Fig. 1](#fg0010){ref-type="fig"}b, we present the diurnal surface temperature versus time at the equator, at the NH vernal equinox and NH summer solstice. Consistent with the incoming solar radiation presented in [Fig. 1](#fg0010){ref-type="fig"}a, the temperature during the NH summer solstice is higher than the NH vernal equinox temperature by a few degrees. The temperature peaks a few hours after the maximum in radiation due to the heat capacity of the ice. The difference between the minimal and maximal temperatures spans about 40 K, consistent with previous studies ([@br0370], [@br0300]). The ice surface temperature diffusion constant, $\kappa_{s}$, plays an important role in the range of variations -- a smaller constant would permit larger variations while a larger value that is similar to Earth\'s ice diffusion constant would yield much smaller variations. The relatively fast warming followed by slower cooling is also apparent and similar to that of [@br0370] and [@br0300]. \[Yet, we note that these studies did not use a temperature diffusion equation to estimate the surface, temperature, and did not include in their calculations all the factors considered here; they only concentrated on the diurnal cycle of the surface temperature and did not calculate the seasonal cycle and annual mean temperature.\] The daily mean values are also plotted in [Fig. 1](#fg0010){ref-type="fig"}b where the dotted horizontal lines indicate the mean of the solid (numerically calculated) curves while the dashed horizontal lines depict the daily mean temperature calculated based on the daily mean incoming solar radiation \[Eqs. [(6)](#fm0060){ref-type="disp-formula"}, [(33)](#fm0330){ref-type="disp-formula"}\]. The difference between these two daily mean temperature estimations is about 3 K, small compared to the diurnal variations in temperature.

The temperature profile within the surface layer of the ice is presented in [Fig. 1](#fg0010){ref-type="fig"}c. The temperature fluctuations decay very fast with depth, and the different temporal temperature profiles converge to a single profile at a depth of about 5 cm. This is due to the small surface ice temperature diffusion constant, $\kappa_{s}$, where higher values of $\kappa_{s}$ would allow deeper convergence together with smaller diurnal fluctuations. The temperature gradient below the fluctuation level is influenced by both the seasonal cycle and the endogenic heating.

We repeated the surface temperature estimation at the SH summer solstice close to the pole (89^∘^S) in [Fig. 1](#fg0010){ref-type="fig"}d. Here the temperature variations are much smaller than the equatorial ones presented in [Fig. 1](#fg0010){ref-type="fig"}b, consistent with the much smaller incoming solar radiation presented in [Fig. 1](#fg0010){ref-type="fig"}a. The difference between the daily mean temperature calculated based on the numerically calculated temperature (dotted line) and the daily mean radiation energy-balance-based temperature (dashed line) is about 1.5 K.

We now switch to the seasonal cycle of the surface temperature. In [Fig. 2](#fg0020){ref-type="fig"}b, we plot the equatorial (blue) and polar (red) daily mean surface temperature as a function of time. These curves are based on the numerically calculated daily mean temperatures (dotted line in [Fig. 1](#fg0010){ref-type="fig"}b,d). At the equator, the difference between the maximal and minimal temperatures is about 4 K, much lower than the maximum-minimum difference at the pole (∼30 K). This is expected due to the relatively large variations in the incoming solar radiation at the high latitudes across the year. The minimum and maximal equatorial temperatures occur at the aphelion and perihelion, close to the equinox, and consistent with [1](#fg0010){ref-type="fig"}b, as at the present day, Jupiter is farther from (or closer to) the Sun during the NH vernal (NH autumnal) equinox, and this eccentricity effect overcomes the effect of obliquity, which is relatively small. Following the above, the equatorial solstice temperature is close to the annual mean equatorial temperature; see also [Fig. 2](#fg0020){ref-type="fig"}c,d. At the poles, as expected, the temperature peaks close to the summer solstice and is minimal at the vernal equinox. There is a change in the rate of cooling starting at the autumnal equinox, when sunlight does not reach the winter pole.

The NH summer solstice temperature as a function of latitude for several internal heating rates *Q* is plotted in [Fig. 2](#fg0020){ref-type="fig"}c. It is clear that the internal heating rate more drastically affects the surface temperature at the poles, and especially at the winter pole where the incoming solar radiation is minimal; the difference between the temperature of the different heating rates exceeds 2 K. Naturally, as the internal heating is the only source of radiation at the winter pole, the temperature difference due to the internal heating is larger than that of the summer pole. The solstice temperatures at the winter/summer poles vary by $\sim 3\ \text{K}/ \sim 1\ \text{K}$ for a internal heat difference of $0.05\ \text{W}\ \text{m}^{- 2}$. We note that even in the absence of internal heating ($Q = 0$) and the absence of solar radiation, the solstice winter polar temperature is much greater than zero, due to the heat capacity of the ice. Yet, if the heat capacity of the ice is neglected (as in [@br0260]), the temperature should drop to zero at the winter pole region; still, in this case, the annual mean polar temperature is not zero as the polar regions absorb solar radiation during the summer.

The numerically calculated annual mean temperature for several internal heating rates *Q* is plotted in [Fig. 2](#fg0020){ref-type="fig"}d. The equator to pole temperature difference is about 50 K, depending on the internal heat rate. As expected, the poles are more drastically affected by the internal heating rate, about 2 K for a difference in internal heating of $0.05\ \text{W}\ \text{m}^{- 2}$.

The internal heating in Europa is not well constrained. In addition, it may vary spatially due to tidal heating within the ice, in both the meridional and longitudinal directions (e.g., [@br0260], [@br0390], [@br0240]). There are several sources that contribute to this heating: rocky mantle (metallic core and silicate mantle) radiogenic heating ($6\text{–}8\ \text{mW}\ \text{m}^{- 2}$, [@br0030]), tidal heating of Europa\'s core (e.g., $30\text{–}230\ \text{mW}\ \text{m}^{- 2}$, [@br0080]), and tidal heating of the icy shell ([@br0390]). The tidal heating of Europa\'s ocean is negligible ([@br0080]). In addition, when the ice is sufficiently shallow (typically, shallower than 10 km), it has only one conductive layer without a bottom convective layer, and the tidal heating of the icy shell is relatively small ([@br0240], [@br0030]). Here we assume that the icy shell is only conductive such that the internal heating originates in Europa\'s rocky mantle. Yet, it is easy to calculate the surface temperature at least when the latitudinal dependence of the tidal heating is specified.

In [Fig. 6](#fg0060){ref-type="fig"} (solid lines), we plot the global mean (weighted by the cosine of the latitude), minimum (NH polar), and maximum (equatorial) annual mean surface temperatures as a function of the internal heating, *Q*. The mean and maximal (equator) surface temperatures are not drastically affected by the internal heating as the main heating source is the incoming solar radiation; the mean and maximal surface temperatures increase by ∼1 K for an increase of internal heating from 0 to $0.2\ \text{W}\ \text{m}^{- 2}$. However, the polar temperature is more drastically affected by the internal heating as the solar radiation is very weak in these locations. As shown in [Fig. 2](#fg0020){ref-type="fig"}c, this increase is much larger at the poles during the (winter) solstice. [Fig. 6](#fg0060){ref-type="fig"} also depicts the temperatures that are based on the annual mean radiation \[Eqs. [(18)](#fm0180){ref-type="disp-formula"}, [(21)](#fm0210){ref-type="disp-formula"}, [(35)](#fm0350){ref-type="disp-formula"}\] (dashed lines). These annual mean temperature estimations that are based on an energy balance consideration are higher by a few degrees than the numerically calculated temperatures, probably since the heat capacity of the ice moderates the seasonal variations in temperatures and since the outgoing radiation is proportional to $T_{s}^{4}$, giving more weight to higher temperatures.Figure 6(a) The annual mean global mean (blue) and maximum (equator, yellow) of Europa\'s surface temperature (in K) as a function of the internal heating rate, *Q* (in W m^−2^). The solid lines indicate temperatures that are based on the numerical solution of Eq. [(29)](#fm0290){ref-type="disp-formula"}, while the dashed lines indicate those that are based on the annual mean radiation and energy balance equation \[Eqs. [(18)](#fm0180){ref-type="disp-formula"}, [(21)](#fm0210){ref-type="disp-formula"}, [(35)](#fm0350){ref-type="disp-formula"}\]. (b) Same as (a) for the minimum (NH polar) surface temperature.Figure 6

The sensitivity of the annual mean surface temperature with regards to the different parameters is summarized in [Table 2](#tbl0020){ref-type="table"}. The equator, NH pole, global mean, and low latitude (15^∘^S--15^∘^N) annual mean temperatures are given for control values and for specific parameter values that are different from the control values. The parameters include the eccentricity, *e*, obliquity, *ε*, emissivity, *ϵ*, surface albedo, $\alpha_{p}$, relative time of Europa\'s eclipse, *p* (as estimated in Section [2.5](#se0080){ref-type="sec"}), surface ice heat diffusion, $\kappa_{s}$, longwave radiation from Jupiter $J_{0}$, and internal heating, *Q*.

As shown in [Figs. 2](#fg0020){ref-type="fig"} and [6](#fg0060){ref-type="fig"}, the internal heating more drastically affects the polar regions, as there the incoming solar radiation heating is relatively small compared to the low latitudes. The effect of the eccentricity is small. This also follows from Eqs. [(18)](#fm0180){ref-type="disp-formula"}, [(21)](#fm0210){ref-type="disp-formula"} from which one can see that the annual mean incoming solar radiation is approximately proportional to $1 + 2e^{2}$ since $e \ll 1$. Thus, the eccentricity increases the annual mean incoming solar radiation, but since for Europa, $e \approx 0.05$, this increase does not exceed 0.5%. The obliquity, *ε*, has a significant effect on the high latitudes and thus must be taken into account. We also study the case of maximum obliquity ($\varepsilon = 3.7^{\circ}$) in comparison to present day obliquity of the control case ($\varepsilon = 3^{\circ}$), and it is apparent that only the polar region is affected by this parameter, by about $3^{\circ}$K. When ignoring the effect of ice emissivity, i.e., taking $\epsilon = 1$, the surface temperature decreases by more than 1 K. The effect of Europa\'s eclipse on surface temperatures reduces the incoming solar radiation by up to 3.3%, and the associated drop in surface temperature does not exceed 1 K.

The effect of the surface ice heat diffusion, $\kappa_{s}$, on the annual mean temperatures is large as a much larger value of $\kappa_{s}$ yielded surface temperatures that are higher by ∼3 K compared to the control temperature values. Moreover, the diurnal cycle and, at the high latitudes, also the seasonal cycle in the surface temperatures are drastically affected by this parameter, exhibiting much smaller variations. More specifically, at the equator, the range of the diurnal variation reduces to about 7 K (compared to 40 K for the control case) while the range of seasonal variations remains 4 K. At the 89^∘^, the diurnal cycle variation is less than 0.2 K (compared to 2 K for the control case) and the seasonal cycle variation is 12 K (compare to 30 K for the control case); see [Figs. 1](#fg0010){ref-type="fig"}, [2](#fg0020){ref-type="fig"}.

We also examine the effect of the longwave radiation of Jupiter that is absorbed by Europa by setting $J_{0} = 0$. It is clear that the contribution of Jupiter\'s longwave radiation on Europa\'s surface temperature is small as the drop in temperature compared to the control case is less than 0.3 K. Thus, this effect can be neglected. [Table 2](#tbl0020){ref-type="table"} also includes cases of no ($Q = 0$) and increased internal heating ($Q = 0.2\ \text{W}\ \text{m}^{- 2}$). Consistent with the results presented above, the internal heating mainly affects the high latitudes.

In [Table 3](#tbl0030){ref-type="table"}, we present the annual mean surface temperatures under the assumption of an energy balance between the incoming solar radiation, internal heating, incoming longwave radiation of Jupiter, and emitted longwave radiation of Europa \[Eq. [(35)](#fm0350){ref-type="disp-formula"}\]. The setup of this table is similar to [Table 2](#tbl0020){ref-type="table"}. Generally speaking, the surface temperatures that are based on the energy balance assumption are warmer than the numerically calculated ones by several degrees, especially at the poles. It is interesting to compare the case of the increased heat diffusion in ice coefficient (increased $\kappa_{s}$ in [Table 2](#tbl0020){ref-type="table"}) to the control case of the energy balance calculation ([Table 3](#tbl0030){ref-type="table"}), as very high $\kappa_{s}$ leads to much smaller diurnal and seasonal variations, which leads to a temporally constant temperature that is equivalent to the temperature calculated based on energy balance assumptions. Indeed the surface temperature estimations of these two cases are very similar. Thus, when the heat diffusion in ice coefficient (or the thermal inertia) becomes larger, the energy-balance-based temperature estimation becomes better.Table 3Sensitivity of surface temperature (K) of Europa to the different parameters using the energy balance assumption \[Eq. [(35)](#fm0350){ref-type="disp-formula"}\]. Control run parameters are as in [Table 2](#tbl0020){ref-type="table"}.Table 3ParameterEquatorPoleGlobal mean15^∘^S--15^∘^N meancontrol98.949.392.298.6*e* = 098.849.292.198.5*ε* = 098.931.192.298.6*ε* = 3.7^∘^99.151.592.498.8*ϵ* = 197.448.590.897.1*α*~*p*~ = 0.5110.454.3102.9110.1*p* = 099.749.692.999.4*J*~0~ = 098.949.392.298.6*Q* = 098.647.291.998.3*Q* = 0.2 W m^−2^99.654.393.199.3

Given the above, we conclude that an accurate albedo map is essential in estimating the low and mid-latitude surface temperatures of Europa while the obliquity and internal heating rate are essential in accurately determining the temperature at the high latitudes. The surface ice heat diffusion has a profound effect on the range of the variability of the diurnal and seasonal cycles of the surface temperatures. To a first approximation, the effects of the eccentricity, Europa\'s eclipse, and the effect of longwave radiation of Jupiter may be neglected. It is plausible that the polar regions\' temperatures will be estimated/measured in the future, similar to the measurement of the low and mid-latitude temperatures that followed the Galileo spacecraft\'s measurements ([@br0370], [@br0300]). If this occurs, for example during future missions (like Europa Clipper of NASA) that will examine the water plumes over the south pole of Europa ([@br0310], [@br0350]), it will be possible to roughly estimate the internal heating rate based on the surface temperature, especially over the winter pole.

It is possible to provide a very rough estimate of Europa\'s ice thickness based on the above; see Section [3.2](#se0110){ref-type="sec"}. This estimate may be a reference thickness for studying the effect of both the vertical and horizontal ice flows due to both tidal heating within the ice and ice flow due to pressure gradients that are associated with variations in ice thickness (see, e.g., [@br0020]).

3.2. Ice thickness {#se0110}
------------------

Below we provide the details regarding the rough estimation of the ice thickness of Europa. The reader that is not interested in the details of the ice thickness may skip this section.

It is possible to obtain a lower limit estimate for Europa\'s ice thickness based on the energy balance between the incoming solar radiation, the outgoing longwave radiation, and the internal heating. This rough estimate may be a reference thickness for studying the effect of both vertical and horizontal ice flow due to both tidal heating within the ice and the ice flow due to pressure gradients that are associated with variations in ice thickness. Here we assume that (i) the internal heating, *Q*, is uniform in space and time, (ii) the icy shell is conductive but not convective such that the temperature within the ice varies linearly with depth, and (iii) the tidal heating within the ice is negligible. These assumptions are probably valid for relatively shallow ice (shallower than 10 km) when there is only one conductive layer ([@br0390]); when there are two layers, an upper conductive layer and a lower convective layer, the temperature is approximately uniform within the convective layer (where the temperature is close to the melting temperature), and tidal heating is not negligible within this layer.

Under the above assumptions, it is possible to calculate the mean thickness of the icy shell as follows:$$h = \rho_{I}c_{p,I}\kappa\frac{T_{f} - T_{s}}{Q},$$ where $\rho_{I}$ is the ice density, $c_{p,I}$ is the heat capacity of the ice, *κ* is the deep ice temperature diffusion constant, $T_{f}$ is the freezing (or melting) temperature of the ice, and *h* is the ice thickness. $T_{s}$ (annual mean energy-balance-based) and *h* depend on latitude and the freezing temperature, and $T_{f}$, depends on the thickness of the ice (through the pressure at the bottom of the ice) and on the salinity of the water as follows ([@br0110], [@br0180]):$$T_{f} = 273.16 + 0.0901 - 0.0575 \times S - 7.61 \times 10^{- 8} \times P_{b},$$ where *S* is the salinity of the ocean water, $P_{b}$ is the pressure (in Pa) at the bottom of the ice (i.e., $P_{b} = g\rho_{I}h$), and $T_{f}$ is given in K. Thus, $T_{f}$ depends linearly on ice thickness.

The dependence of the thickness of the ice on latitude is not trivial as ice may flow due to gradients in ice thickness (see, for example, [@br0400], [@br0020]). To bypass this complexity, we only estimate the mean ice thickness by performing a spatial mean on Eq. [(36)](#fm0360){ref-type="disp-formula"}. In addition, since $T_{f}$ depends on *h*, the global mean thickness of the ice can be estimated iteratively as follows:$${\langle h\rangle}_{j + 1} = \rho_{I}c_{p,I}\kappa\frac{{\langle T_{f}\rangle}_{j} - \langle T_{s}\rangle}{Q},$$ where *j* is a counter that indicates the number of the iteration. We start the process from a typical mean freezing temperature (e.g., $T_{f} \approx 270\ \text{K}$), find the mean thickness, and then use it to estimate the new mean freezing temperature. Convergence is achieved after a few iterations.

The global mean ice thickness and freezing temperature as a function of internal heating are shown in [Fig. 7](#fg0070){ref-type="fig"}. First, as predicted by Eq. [(38)](#fm0380){ref-type="disp-formula"}, the mean ice thickness, $\langle h\rangle$, inversely depends on the internal heating, *Q*. Second, the thickness hardly depends on the salinity of the water. Third, the mean freezing temperature converges to a constant value for realistic internal heating values (i.e., $Q > 40\ \text{mW}\ \text{m}^{- 2}$). Fourth, an increase of salinity by 100 ppt decreases the freezing temperature by 5.75 K; this is a direct consequence of Eq. [(37)](#fm0370){ref-type="disp-formula"}.Figure 7(a) Global mean ice thickness as a function of the internal heating rate, *Q*, when the underlying ocean is fresh (blue) and salty (100 and 200 ppt, green and red). The ice thickness is hardly affected by the salinity of the ocean water. The salinity unit is "parts per thousand" (ppt), or grams of salt per kilogram of seawater. (b) Global mean freezing temperature of seawater as a function of the internal heating rate, *Q*, for fresh (blue) and salty (100 and 200 ppt, green and red) water. (c) Ice thickness (in km), *h*, and freezing temperature (in ^∘^C), *T*~*f*~, as a function of latitude for internal heating of *Q* = 0.05,0.1 W m^−2^. As mentioned in the text, this is a very rough estimate of the ice thickness as our underlying assumption is that the ice is stagnant and that the tidal heating is negligible.Figure 7

A mean ice thickness of $\langle h\rangle = 10\ \text{km}$ was obtained for an internal heating of $Q \approx 0.05\ \text{W}\ \text{m}^{- 2}$. Below this value, the estimated mean ice thickness should be regarded as a lower limit, as an almost uniform temperature bottom convective layer may be formed, violating our assumption of a linear increase of temperature within the ice with depth. For a mean ice thickness shallower than $\sim 10\ \text{km}$, we expect one conductive layer, and the estimated mean ice thickness is more accurate.

It is possible to calculate the ice thickness and the freezing temperature as a function of latitude under the assumption that the ice does not flow and is stagnant. In this case, there is a simple energy balance between the internal heating, the incoming solar radiation and the outgoing longwave radiation at each latitude. [Fig. 7](#fg0070){ref-type="fig"}c depicts the ice thickness as a function of latitude for different internal heating rates, using Eq. [(36)](#fm0360){ref-type="disp-formula"}. As expected, the ice thickness increases poleward and the equator to pole gradient becomes smaller as the internal heating increases. The gradient in thickness is on the order of a few kilometers. The freezing temperature as a function of latitude is shown in [Fig. 7](#fg0070){ref-type="fig"}c. As expected, the freezing temperature is higher for shallower ice, and the thickness increases toward the poles. The typical equator to pole freezing temperature gradient is about 0.2^∘^C and is larger for thicker ice and a smaller internal heating rate.

4. Discussion and conclusions {#se0120}
=============================

The Galileo mission triggered many studies regarding the moon Europa (see, [@br0290]). One of the observations made by the Galileo spacecraft (Photopolarimeter-Radiometer) was used to measure Europa\'s surface temperature ([@br0370], [@br0300]). [@br0370] concentrated on low latitude temperatures, which were also relatively high, and consequently suggested either a low local albedo (0.5 compared with $0.68 \pm 0.14$ of [@br0140]) or a very high local endogenic heating (of $1\ \text{W}\ \text{m}^{- 2}$). The mean surface temperature measurement of [@br0370] was limited to latitudes equatorward of $\sim \pm 70^{\circ}$, and these approximately correspond to an internal heating of $\sim 0.05\ \text{W}\ \text{m}^{- 2}$. The temperature measurements of [@br0370] and [@br0300] are diurnal temperatures -- these helped us to tune the range of diurnal variations through the heat diffusion in ice coefficient. We performed the sensitivity tests of the heat diffusion in ice parameter since this parameter is not well constrained in the upper part of the ice and since these tests helped us to understand the role of the surface heat capacity of the ice. The other parameters are fairly constrained except the internal heating. We further approximated the seasonal and annual mean temperature for a range of internal heating.

We used two ways to estimate the surface temperature of Europa. The first method is by using a diffusion equation for the top (several meters) layer of the ice. This diffusion equation is subject to top and bottom boundary conditions. The second method is based on an energy balance between the incoming and outgoing radiation. The first method is more accurate as it is able to simulate, fairly reasonably, the diurnal cycle while the second method may be suitable to estimate the annual mean temperature. We have also examined an alternative way, similar to the approach used in previous studies ([@br0370], [@br0300]), where we use a one layer simple ordinary differential equation to estimate the surface temperature of Europa; i.e., we solved the following differential equation: $\rho_{I}c_{p,I}d\frac{\partial T_{s}}{\partial t} = W(1 - \alpha_{p})(1 - p) + W_{j} - \epsilon\sigma T_{s}{(t)}^{4}$ where $T_{s}$ is the surface temperature and the scale depth is $d = \sqrt{2\kappa_{s}/\Omega_{e}}$. However, this approach yielded much colder polar temperatures in comparison to the more fundamental diffusion equation model. We thus refrained from presenting the results of the single layer model.

Future measurements of Europa\'s surface temperature may include the polar regions. Exact polar temperatures may help to better estimate the rate of the internal heating as these regions are only partly affected by the solar radiation, and the contribution of the internal heating in these regions is large, especially during the winter solstice. Polar temperatures have hardly any diurnal variations, thus simplifying the estimation of the diurnal mean temperature. Large uncertainties are associated with the internal heating and the ice depth of Europa, and an exact estimation of the polar region temperatures may significantly reduce these uncertainties ([@br0020]).

In summary, we discussed the diurnal and seasonal variations of the incoming solar radiation to Europa and developed a mathematical approximation for the global mean incoming solar radiation at Europa\'s surface (Section [2.3](#se0040){ref-type="sec"}). Based on the incoming solar radiation, we estimated the diurnal, seasonal and annual mean surface temperature in two ways: 1) based on the numerical integration of a temperature diffusion equation and 2) based on the energy balance between the incoming solar radiation, internal heating, Jupiter\'s radiation, and outgoing longwave radiation. Our estimates take into account the eccentricity of Jupiter, as well as Europa\'s obliquity, emissivity, eclipse, surface ice heat diffusion, Jupiter\'s longwave radiation, and internal heating. We showed that the temperature varies moderately at the low latitudes and much more drastically at the high latitudes and that the high latitudes are more drastically affected by the internal heating, especially during the winter solstice. The diurnal and seasonal variations are controlled by the heat diffusion in ice coefficient. For a typical internal heating rate of $0.05\ \text{W}\ \text{m}^{- 2}$ (e.g., [@br0390], [@br0240]), the equator, pole, and global mean annual surface temperatures are 96 K, 46 K, and 90 K, respectively. These values are not far from the, much simpler, energy-balance-based temperatures (99 K, 49 K, and 92 K, respectively) such that energy balance based temperatures can be used when studying the ice flow and ocean dynamics of Europa. Based on the internal heating rate, we provide a very rough estimate for the mean thickness of Europa\'s icy shell (Section [3.2](#se0110){ref-type="sec"}). We also estimate the incoming solar radiation to Enceladus, the moon of Saturn (Section [2.3.2](#se0060){ref-type="sec"}). The approach we developed here may be applicable to other moons in the solar system.

Declarations {#se0130}
============

Author contribution statement {#se0140}
-----------------------------

Yosef Ashkenazy: Conceived and designed the experiments; Performed the experiments; Analyzed and interpreted the data; Contributed reagents, materials, analysis tools or data; Wrote the paper.

Funding statement {#se0150}
-----------------

This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.

Competing interest statement {#se0160}
----------------------------

The authors declare no conflict of interest.

Additional information {#se0170}
----------------------

No additional information is available for this paper.

We thank Roiy Sayag and Eli Tziperman for helpful discussions.
