Ocean Density (kg/m 3)

Shell Density (kg/m3)

Ocean Thickness (km)a60
40200950

900

850
1100
1080
1060
1040
1020
1000

02040600

Shell Thickness (km)204060 850

Ocean Thickness (km)900

Shell Density (kg/m3)

950
1000

1050

1100

Ocean Density (kg/m 3)

Ocean Density (kg/m 3)

Shell Density (kg/m3)

Ocean Thickness (km)

b60
40
200
950

900

850
1100
1080
1060
1040
1020
1000

0204060

Shell Thickness (km)0204060 850

Ocean Thickness (km)900

Shell Density (kg/m3)

950
1000

1050

1100

Ocean Density (kg/m 3)

Figure S1: The ranges of model parameters that best account for the observed gravitational potential andlibration amplitude, given the shape model of (a) Nimmo et al. (2011) or (b) Tajeddine et al. (2017), andassuming the topography is supported by Airy isostasy (Hemingway and Matsuyama, 2017). Illustrated along
the diagonal are the normalized probability density functions across each of the four model parameters. The
remaining cells in the grid illustrate the joint probability density function for each pair of parameters. The
shaded regions represent 68% (dark), 95% (intermediate), and 99.7% (pale) confidence contours (section 2.6).
The minimum misfits are L = 0.6 and L = 3.2 for the Nimmo et al. (2011) and Tajeddine et al. (2017) shapemodels, respectively.

Mantle Density (kg/m 3)
Mean Shell Density (kg/m3)
Mantle Thickness (km)a6040200950

900

850
1300

1200

1100

1000

2040600204060 850

900950
1000

1100

1200

1300

Mantle Thickness (km) Mean Shell Density (kg/m3) Mantle Density (kg/m 3)

Compensation Depth (km)Mantle Density (kg/m 3)
Mean Shell Density (kg/m3)
Mantle Thickness (km)b6040200950

900850
1300

1200

1100

1000

204060Compensation Depth (km)0204060 850

900950
1000

1100

1200

1300

Mantle Thickness (km) Mean Shell Density (kg/m3) Mantle Density (kg/m 3)

Figure S2: The ranges of model parameters that best account for the observed gravitational potential, given
the shape model of (a) Nimmo et al. (2011) or (b) Tajeddine et al. (2017), and assuming the topography issupported by Pratt isostasy. Illustrated along the diagonal are the normalized probability density functionsacross each of the four model parameters. The remaining cells in the grid illustrate the joint probability
density function for each pair of parameters. The shaded regions represent 68% (dark), 95% (intermediate),and 99.7% (pale) confidence contours (section 2.6). The minimum misfits are L = 6.0 and L = 7.0 for theNimmo et al. (2011) and Tajeddine et al. (2017) shape models, respectively.

aOcean Thickness (km)70605040302010

Elastic Thickness (m)020015010050002040600Shell Thickness (km)b2040600Ocean Thickness (km)50100150200

Elastic Thickness (m)Ocean Thickness (km)70605040302010

Elastic Thickness (m)020015010050002040

Shell Thickness (km)6002040

Ocean Thickness (km)60050

100150200

Elastic Thickness (m)Figure S3: The ranges of model parameters that best account for the observed gravitational potential, given
the shape model of (a) Nimmo et al. (2011) or (b) Tajeddine et al. (2017), and assuming the topography
is supported by Airy isostasy plus elastic flexure, assuming top loading. Illustrated along the diagonal are
the normalized probability density functions across each of the three model parameters. The remaining cells
in the grid illustrate the joint probability density function for each pair of parameters. The shaded regions
represent 68% (dark), 95% (intermediate), and 99.7% (pale) confidence contours (section 2.6). The minimum
misfits are L = 0.6 and L = 3.3 for the Nimmo et al. (2011) and Tajeddine et al. (2017) shape models,
respectively.

aOcean Thickness (km)70605040302010

Elastic Thickness (m)01000800600400200002040600Shell Thickness (km)b2040600

Ocean Thickness (km)200

400600800

1000

Elastic Thickness (m)Ocean Thickness (km)70605040302010

Elastic Thickness (m)

01000
800
600
400
200002040

Shell Thickness (km)6002040

Ocean Thickness (km)

60

0200

400

600

800

1000

Elastic Thickness (m)Figure S4: The ranges of model parameters that best account for the observed gravitational potential, given
the shape model of (a) Nimmo et al. (2011) or (b) Tajeddine et al. (2017), and assuming the topography
is supported by Airy isostasy plus elastic flexure, assuming bottom loading. Illustrated along the diagonal
are the normalized probability density functions across each of the three model parameters. The remaining
cells in the grid illustrate the joint probability density function for each pair of parameters. The shaded
regions represent 68% (dark), 95% (intermediate), and 99.7% (pale) confidence contours (section 2.6). The
minimum misfits are L = 1.0 and L = 3.7 for the Nimmo et al. (2011) and Tajeddine et al. (2017) shapemodels, respectively.

90°N50750015000251250751000

030°N

Core Topography (meters)

00
1060°N

a

500

0°
30°S

0-12-750

025-100-10 50-750-2500-50 0-75 000-10-500

0250
500
750
1000
1250
1500

60°S

-500-250-1000
-1500

90°S
0°
90°N

45°E

b

50
100
150

90°E

135°E

180°E

225°E

270°E

-100-50

0-250

315°E

360°E

-150-200300

0200

35

150
0
10 0
51000000-1-20

0-1

150
100
50
500-10-50-150-1-1050

200°
30°S

0-500-200300 02502030°N-100-50-200

00Core Topography (meters)60°N

60°S
-300

90°S
0°

45°E

90°E

135°E

180°E

225°E

270°E

315°E

360°E

Figure S5: Non-hydrostatic core topography (up to l = 3 only) required to satisfy gravity observations
perfectly given an external shape described by Tajeddine et al. (2017) and assuming (a) no isostatic compensation, such that the ice/ocean interface follows a hydrostatic equilibrium figure; or (b) the topography
is supported by Airy isostasy, as in section 3.2.3, such that most of the compensation occurs through deflection of the ice/ocean interface and that the small remaining gravity anomalies are accommodated by
non-hydrostatic core topography. In this example case, dshell = 21 km, docean = 37 km, ρshell = 925 kg/m3 ,and ρocean = 1020 kg/m3 .

050

51520

10030251040

1560°N

155204035302530°N

110 5
5300°
0520

0030°S90°S
0°10060°S

105
1520
2530
35
4050
80120

Heat Flux (mW/m2 )

90°N

045°E

90°E

135°E

180°E

225°E

270°E

315°E

360°E

Figure S6: Heat anomaly obtained by subtracting the tidal dissipation model results (Figure 11a) from theinferred heat flux (Figure 11b). Integrated over the surface, the total heat anomaly here is approximately
3 GW.