Membrane Worlds

Preamble

Author: Mikael Beuthe
Purpose: computation of tidal/load Love numbers of ‘membrane worlds’, or icy satellites having a subsurface ocean under a thin crust
Version: 1.0
Date: 16 June 2015

The notebook is based on the following papers:
MP (for ‘Membrane Paradigm’) = Beuthe (2015a), Tides on Europa: The membrane paradigm, Icarus 248, 109-134; arXiv.org:1410.4735v1
MW (for ‘Membrane Worlds’) = Beuthe (2015b), Tidal Love numbers of membrane worlds: Europa, Titan and Co., Icarus (in press, available online) doi:10.1016/j.icarus.2015.06.008; arXiv.org:1504.04574v2 (revised version)
If you use this notebook to obtain results for a publication, or include part of it in your own code, please cite the papers above.

Copyright (C) 2015, Mikael Beuthe
All rights reserved.

Table of contents

Part I: Pedestrian example

      1. Love numbers
          First thing to do
           Outline of the methos
           Physical model: general structure
           Choose a fluid-crust model
           Fluid-crust model: density and gravity
           Fluid-crust model: mantle rheology
           Fluid-crust model Love numbers
           Adding the crust
           Physical Love numbers
           Compact formulation
       
       2. Dissipation
           Total dissipation
            Dissipation within the crust
            Dissipation within the core and mantle
            Heat budget

Part II: Definitions

      1. Notation
      
      2. Formulas
          Maxwell rheology
           Two-layer crust with Maxwell rheology
           Fluid-crust Love numbers
           Relations between Love numbers
           Static Love numbers (gravitational and radial)
           Dynamical Love numbers (gravitational and radial)
           A simple thick shell model: the homogeneous crust model

Part III: Other examples

      0. Constants
      
      1. Effective viscoelastic parameters
          Example 1.1: viscoelastic parameters of a uniform layer with Maxwell rheology (MP-Fig 1)
           Example 1.2: effective viscoelastic parameters of a two-layer crust (MP-Fig 2)
           Example 1.3: effective viscoelastic parameters of a two-layer crust (MW-Fig 3)
       
       2. (Subscript[l, 2])-(Subscript[h, 2]) relation
           Example 2.1: (Subscript[l, 2])-(Subscript[h, 2]) ratio as a function of the inverse frequency/inverse viscosity of the crust (MP-Fig 3)
             Example 2.2: (Subscript[l, 2])-(Subscript[h, 2]) ratio as a function of crust thickness (MP-Fig 4)
             Example 3.1: normalized tilt factor in function of the inverse frequency/inverse viscosity of the crust (MP-Fig 5)
       
       3. (Subscript[k, 2])-(Subscript[h, 2]) relation
           Example 3.1: normalized tilt factor in function of the inverse frequency/inverse viscosity of the crust (MP-Fig 5)
             Example 3.2: normalized tilt factor as a function of crust thickness (MP-Fig 6)
             Example 3.3: slope of tilt factor at zero crust thickness (MW-Fig 7)
       
       4. Static Love numbers
           Example 4.1: Radial Love number of Europa as a function of the inverse frequency/inverse viscosity of the crust, for diurnal tides (MP-Fig 7A)
             Example 4.2: Radial Love number of Europa as a function of the inverse frequency/inverse viscosity of the crust, for NSR tides (MP-Fig 7B)
             Example 4.3: Tidal Love numbers of Titan as functions of crust thickness for different crust/ocean density models (MW-Fig 5)
             Example 4.4: Tidal Love numbers of Titan in function of crust thickness: choose another fluid-crust density (not in MW)
             Example 4.5: Load Love numbers of Titan in function of crust thickness for different crust/ocean density models (not in MW)
             Example 4.6: Tilt factor of Europa as a function of crust thickness (MW-Fig 8)
       
       5. Dynamical Love numbers
           Example 5.1: Dynamical tilt factor of Europa as a function of ocean thickness (MW-Fig 14)
            Example 5.2: Dynamical resonance (MW-Fig 15)
       
       6. A simple thick shell model
           Example 6.1: Dependence of geometrical factors on crust thickness (MP-Fig 13; MW-Table 11)

Part I: Pedestrian example

Let’s say you want to check the effect of crustal rheology on Europa’s tidal deformation and dissipation. For this, you need to compute Europa’s tidal Love numbers.
As an application, we will evaluate the heat dissipated in the whole satellite and determine the specific contributions of the crust and mantle.

1. Love numbers

First thing to do

Evaluate the section II.2. Formulas.
You will be asked to do it anyway, because the formulas section is an Initialization Cell that should be evaluated before any other.
Also, clear the variables used in the example:

Clear[y,ξ0,μmHat]
Clear[ω,ecc,R,Rm,ρb,g,ρbar,ρ,δρ,d,μEm,μE,νE]
Clear[ϵ,δm,η,μ,ν,χ,f,δ1,δ2]
Clear[obli,Gconstant]

Outline of the method

The  membrane formulas for the gravitational Love number of degree n is

kn[μhat,ν,χ,δρhat,ϵ,ξ0,hn0,n]

The formula for the radial Love number of degree n is similar:

hn[μhat,ν,χ,δρhat,ϵ,ξ0,hn0,n]

The parameters (μhat, ν, χ, δρhat, ϵ) are relevant to the shallow interior (crust and crust-ocean boundary).
The parameters (ξ0, hn0) are relevant to the deep interior (ocean and below).

In the membrane approach, computing Love numbers is a two-step procedure:
1) evaluate the fluid-crust Love number hn0, i.e. the radial Love number of a model in which the crust is fluid and may differ in density from the physical crust.
2) add the crust modeled as a membrane, which means that its viscoelasticity has been integrated over crust thickness.
    Any arbitrarily complicated depth-dependent rheology is then summed up into three complex numbers: the effective viscoelastic parameters.

Physical model: general structure

For Europa, we know that the surface radius is R=1560.8 km, the bulk density is ρb=3013 kg/m3, and the surface gravity is g=1.315 m/s2.
The bottom of the ocean must be be at a depth of less than 170 km; assume here that this bound corresponds to the radius of the mantle.
The density of the crust (ρbar) and ocean (ρ) are unknown: assume that ρbar=950 kg/m3 and ρ=1050 kg/m3.
The crust-ocean density contrast is defined as δρ= ρbar-ρ=-100kg/m3. For the crust thickness, assume that d=25 km (no attempt is made here to fit the observed moment of inertia).
For later use, note that the tidal angular frequency is ω=2.048 x (10^(-5)) rad/s and that the orbital eccentricity is ecc=0.0094.
These values define the following replacement rule:

rulePedestrian={ω->2.048*10^-5,ecc->0.0094,R->1560.8*10^3,Rm->1560.8*10^3-170.*10^3,ρb->3013.,g->1.315,ρbar->950.,ρ->1050.,δρ->-100.,d->25.*10^3}

Choose a fluid-crust model

The first step consists in choosing a model for the deep interior. The subsurface ocean partially decouples the mantle from the crust, so that the core and mantle mainly affect the Love numbers through their bulk mass. For this reason, we consider here a very simple deep interior model: core and mantle are treated as one incompressible viscoelastic layer, below a homogeneous and incompressible ocean (density stratification of the ocean is not so important for Europa because the ocean is not very deep). However, nothing forbids you to use a more complicated interior model, with any number of layers or even a continuous radial variation of the physical parameters.
Our choice is thus the VL model, where ‘VL’ stands for ‘Viscoelastic/Liquid’: this is a two-layer incompressible model made of viscoelastic mantle and a surface ocean. This choice implies that the fluid-crust density is equal to the ocean density.
The degree-two radial Love number of the VL model is given by:

hVL0[y,ξ0,μmHat,2]

In this formula, the variables are:
y=Rm/R: relative radius of the mantle with respect to the surface,
ξ0 = ρ/ρb0 : ocean-to-bulk density ratio,
μmHat=μm/(ρb0*g0*R) : nondimensional shear modulus of the mantle.

Fluid-crust model: density and gravity

Now, the bulk density ρb0 and the surface gravity g0 of the fluid-crust model are not the same as those of the physical model.
This is because the fluid crust has the density of the ocean, instead of having the density of the physical crust.
This changes the bulk density and surface gravity as follows:

ρb0[ρb,δρ,ϵ]
g0[g,ρb,δρ,ϵ]

where ϵ is the relative crust thickness: ϵ=d/R. These relations are correct up to first order in ϵ=d/R.
For the parameters chosen above, the bulk density and surface gravity of the fluid-crust model are

{ρb0[ρb,δρ,(d/R)],g0[g,ρb,δρ,(d/R)]}/.rulePedestrian

Fluid-crust model: mantle rheology

Before computing the fluid-crust Love numbers, we need to consider mantle rheology.
For the elastic shear modulus of the mantle, assume that μEm=40 GPa.
If the mantle is viscoelastic, its shear modulus is not 40 GPa but a complex number depending on the forcing frequency and on the rheological model.
You could choose any rheological model (Burgers, Andrade etc), but assume here for simplicity that the rheology is Maxwell:

μV[μE,δm]

where (Subscript[δ, m]) is a dimensionless parameter inversely proportional to tidal frequency ω and viscosity (Subscript[η, m]): (Subscript[δ, m])=(Subscript[μ, E])/(ω*(Subscript[η, m])).
If (Subscript[δ, m])=0, the body is purely elastic;
if (Subscript[δ, m])=1, the body is in the critical state;
if (Subscript[δ, m])>>1, the body is fluid-like.
The forcing angular frequency is ω=2.048 x (10^(-5)) rad/s;
the viscosity of the core-mantle system could be (10^21)Pa.s, as the Earth’s mantle. In that case, the (Subscript[δ, m]) parameter is

((μEm)/(ω*ηm))/.rulePedestrian/.{μEm->40.*10^9,ηm->1.*10^21}

which is so small (δ<<1)  that the core-mantle system can be treated in practice as being purely elastic at tidal frequencies (i.e. you can set δ = 0).
Just for the fun, let’s assume that the viscosity of the core-mantle system is such that its rheology is not far from critical: δ = 0.1:

ruleMantleRheology={μEm->40.*10^9,δm->0.1};

Fluid-crust Love numbers

The nondimensionless parameters are expressed in terms of physical parameters with the rule

ruleMantleParam={y->(Rm/R),ξ0->(ρ/(ρb0[ρb,δρ,(d/R)])),μmHat->((μV[μEm,δm])/(ρb0[ρb,δρ,(d/R)]*g0[g,ρb,δρ,(d/R)]*R))};

Now evaluate the radial fluid-crust Love number:

h2FluidCrust=hVL0[y,ξ0,μmHat,2]/.ruleMantleParam/.ruleMantleRheology/.rulePedestrian

In comparison with the elastic case, the real part of the fluid-crust Love number only changes at the 5th decimal, but the imaginary part (representing dissipation) is non-zero at the 3rd decimal. Note that mantle dissipation is still rather low for δ=0.1: it increases until δ=54 before decreasing again.

As the surface is hydrostatic, the gravitational fluid-crust Love number is related to the radial one by

k2FluidCrust=h2FluidCrust-1

Adding the crust

The second step consists in adding the crust.
Assume that the crust is made of an elastic top layer (or lithosphere) and of a bottom layer that may be convecting.
The relative thickness of the top layer is f=(Subscript[d, top])/d.
The rheology of each layer is Maxwell, the top and bottom layers being parameterized by (Subscript[δ, 1]) and (Subscript[δ, 2]), respectively.
Thus the effective viscoelastic parameters of the crust are given by the following replacement rule

ruleEffective={μ->μbar[μE,f,δ1,δ2],ν->νbar[νE,f,δ1,δ2],χ->χbar[νE,f,δ1,δ2]};

For example, the two-layer effective shear modulus is given by

μbar[μE,f,δ1,δ2]

The elastic shear modulus of the icy crust is 3.5 GPa. The elastic Poisson’s ratio is (Subscript[ν, E])=0.33.
For the boundary between top and bottom layers, assume that f = 0.3.
Regarding the rheology, the top layer has (Subscript[δ, 1])=0 (because it is elastic) and the bottom layer has (Subscript[δ, 2])=1 (critical state).

ruleCrustRheology={μE->3.5*10^9,νE->0.33,f->0.3,δ1->0,δ2->1};
μbar[μE,f,δ1,δ2]/.ruleCrustRheology
χbar[νE,f,δ1,δ2]/.ruleCrustRheology
νbar[νE,f,δ1,δ2]/.ruleCrustRheology

Physical Love numbers

It is now possible to compute the gravitational Love number (Subscript[k, 2]) with the formula kn[μhat, ν, χ, δρhat, ϵ, ξ0, hn0, 2]

k2=kn[(μ/(ρ*g*R)),ν,χ,((δρ)/ρ),(d/R),(ρ/(ρb0[ρb,δρ,(d/R)])),h2FluidCrust,2]/.rulePedestrian/.ruleEffective/.ruleCrustRheology

as well as the radial Love number (Subscript[h, 2]) with the formula hn[μhat, ν, χ, δρhat, ϵ, ξ0, hn0, 2]

h2=hn[(μ/(ρ*g*R)),ν,χ,((δρ)/ρ),(d/R),(ρ/(ρb0[ρb,δρ,(d/R)])),h2FluidCrust,2]/.rulePedestrian/.ruleEffective/.ruleCrustRheology

The tangential Love number (Subscript[l, 2]) is related to (Subscript[h, 2]) by the (Subscript[l, 2])-(Subscript[h, 2]) relation:

l2=lnhn[ν,2]*h2/.ruleEffective/.ruleCrustRheology

The (Subscript[l, 2])-(Subscript[h, 2]) relation is only valid at order (ϵ^0), i.e. in the limit of zero crust thickness, whereas the formulas above for (Subscript[k, 2])and (Subscript[h, 2]) are valid up to order ϵ.
Thus the imaginary part of (Subscript[l, 2]) (which is of order ϵ) cannot be fully trusted.

Compact formulation

If you want to evaluate the above model for other sets of parameters, define a (Subscript[k, 2]) function that depends directly on the physical parameters:

k2VLVV[R_,d_,Rm_,ρb_,ρbar_,ρ_,g_,μE_,νE_,f_,δ1_,δ2_,μEm_,δm_]:=kn[((μbar[μE,f,δ1,δ2])/(ρ*g*R)),νbar[νE,f,δ1,δ2],χbar[νE,f,δ1,δ2],((ρbar-ρ)/ρ),(d/R),(ρ/(ρb0[ρb,ρbar-ρ,(d/R)])),hVL0[(Rm/R),(ρ/(ρb0[ρb,ρbar-ρ,(d/R)])),((μV[μEm,δm])/(ρb0[ρb,ρbar-ρ,(d/R)]*g0[g,ρb,ρbar-ρ,(d/R)]*R)),2],2]

In the function name, ‘VLVV’ denotes Viscoelastic/Liquid/Viscoelastic/Viscoelastic for the 4-layer structure: viscoelastic core-mantle, (liquid) ocean, and two viscoelastic crust layers.
Do the same for (Subscript[h, 2]):

h2VLVV[R_,d_,Rm_,ρb_,ρbar_,ρ_,g_,μE_,νE_,f_,δ1_,δ2_,μEm_,δm_]:=hn[((μbar[μE,f,δ1,δ2])/(ρ*g*R)),νbar[νE,f,δ1,δ2],χbar[νE,f,δ1,δ2],((ρbar-ρ)/ρ),(d/R),(ρ/(ρb0[ρb,ρbar-ρ,(d/R)])),hVL0[(Rm/R),(ρ/(ρb0[ρb,ρbar-ρ,(d/R)])),((μV[μEm,δm])/(ρb0[ρb,ρbar-ρ,(d/R)]*g0[g,ρb,ρbar-ρ,(d/R)]*R)),2],2]

As an example, plot (Subscript[k, 2])and (Subscript[h, 2]) in terms of crust thickness.
Assume that the crust top layer and the core-mantle system are purely elastic ((Subscript[δ, 1])=0 and (Subscript[δ, m])=0).
Consider three rheology models for the crust bottom layer: (Subscript[δ, 2])=0 (elastic, black), (Subscript[δ, 2])=1 (critical, red), (Subscript[δ, 2])=100 (fluid-like, blue).
Plot the real parts:

{dMin,dMax}={0,170};
ruleCommon=Join[rulePedestrian,ruleCrustRheology,ruleMantleRheology];
Fig1=Show[
Plot[Re@k2VLVV[R,10^3*dcrust,Rm,ρb,ρbar,ρ,g,μE,νE,f,δ1,0,μEm,0]/.ruleCommon,{dcrust,dMin,dMax},PlotStyle->{Black,Thick}],
Plot[Re@k2VLVV[R,10^3*dcrust,Rm,ρb,ρbar,ρ,g,μE,νE,f,δ1,1,μEm,0]/.ruleCommon,{dcrust,dMin,dMax},PlotStyle->{Red,Thick}],
Plot[Re@k2VLVV[R,10^3*dcrust,Rm,ρb,ρbar,ρ,g,μE,νE,f,δ1,100,μEm,0]/.ruleCommon,{dcrust,dMin,dMax},PlotStyle->{Blue,Thick}],
FrameLabel->{"Crust thickness [km]","Re ((Subscript[k, 2]))"},PlotRange->{0.155,0.29}];
Fig2=Show[
Plot[Re@h2VLVV[R,10^3*dcrust,Rm,ρb,ρbar,ρ,g,μE,νE,f,δ1,0,μEm,0]/.ruleCommon,{dcrust,dMin,dMax},PlotStyle->{Black,Thick}],
Plot[Re@h2VLVV[R,10^3*dcrust,Rm,ρb,ρbar,ρ,g,μE,νE,f,δ1,1,μEm,0]/.ruleCommon,{dcrust,dMin,dMax},PlotStyle->{Red,Thick}],
Plot[Re@h2VLVV[R,10^3*dcrust,Rm,ρb,ρbar,ρ,g,μE,νE,f,δ1,100,μEm,0]/.ruleCommon,{dcrust,dMin,dMax},PlotStyle->{Blue,Thick}],
FrameLabel->{"Crust thickness [km]","Re ((Subscript[h, 2]))"},PlotRange->{0.78,1.32}];
Show[GraphicsRow[{Fig1,Fig2}]]

Plot the imaginary parts:

{dMin,dMax}={0,170};
ruleCommon=Join[rulePedestrian,ruleCrustRheology,ruleMantleRheology];
Fig1=Show[
Plot[Im@k2VLVV[R,10^3*dcrust,Rm,ρb,ρbar,ρ,g,μE,νE,f,δ1,0,μEm,0]/.ruleCommon,{dcrust,dMin,dMax},PlotStyle->{Black,Thick}],
Plot[Im@k2VLVV[R,10^3*dcrust,Rm,ρb,ρbar,ρ,g,μE,νE,f,δ1,1,μEm,0]/.ruleCommon,{dcrust,dMin,dMax},PlotStyle->{Red,Thick}],
Plot[Im@k2VLVV[R,10^3*dcrust,Rm,ρb,ρbar,ρ,g,μE,νE,f,δ1,100,μEm,0]/.ruleCommon,{dcrust,dMin,dMax},PlotStyle->{Blue,Thick}],
FrameLabel->{"Crust thickness [km]","Im ((Subscript[k, 2]))"},PlotRange->{-0.026,0.001}];
Fig2=Show[
Plot[Im@h2VLVV[R,10^3*dcrust,Rm,ρb,ρbar,ρ,g,μE,νE,f,δ1,0,μEm,0]/.ruleCommon,{dcrust,dMin,dMax},PlotStyle->{Black,Thick}],
Plot[Im@h2VLVV[R,10^3*dcrust,Rm,ρb,ρbar,ρ,g,μE,νE,f,δ1,1,μEm,0]/.ruleCommon,{dcrust,dMin,dMax},PlotStyle->{Red,Thick}],
Plot[Im@h2VLVV[R,10^3*dcrust,Rm,ρb,ρbar,ρ,g,μE,νE,f,δ1,100,μEm,0]/.ruleCommon,{dcrust,dMin,dMax},PlotStyle->{Blue,Thick}],
FrameLabel->{"Crust thickness [km]","Im ((Subscript[h, 2]))"},PlotRange->{-0.11,0.005}];
Show[GraphicsRow[{Fig1,Fig2}]]

2. Dissipation

Total dissipation

As an application of the Love numbers computed above, we will now evaluate the heat dissipated in the whole satellite and determine the specific contributions of the crust and mantle. The satellite is synchronously rotating with spin rate ω, orbital eccentricity ecc and obliquity obli.

The global heat flow due to tidal dissipation in the whole satellite is equal to the product of Im((Subscript[k, 2])) and of a global factor independent of the internal structure (see MP-99; MW-100):

Psi0=((21/5)(ecc^2)+(3/5)(((Sin[obli]))^2))/.{obli->0}/.rulePedestrian
GlobalFac=-(5/2)*((((ω*R))^5)/Gconstant)*Psi0/.{Gconstant->6.67*10^-11}/.rulePedestrian

Here the obliquity has been set to zero.
Psi0 represents the nondimensional surface average of the squared norm of the tidal potential.

The global heat flow (in W) due to tidal dissipation in the whole satellite is

TotalDissip=GlobalFac*Im[k2]/.rulePedestrian

Dissipation within the crust

The membrane approach also provides a formula for the dissipation within the crust (MP-98; MW-101), which is proportional to |(Subscript[h, 2])((|)^2) and to Im(Λ), where Λ is the membrane spring constant given by

Λ2=Λ[(μ/(ρ*g*R)),ν,(d/R),2]/.rulePedestrian/.ruleEffective/.ruleCrustRheology

The global heat flow (in W) due to dissipation within the crust is given by (MP-98 and 138; MW-101)

CrustDissip=GlobalFac*(-(3/5))*(ρ/(ρb))*(((Abs[h2]))^2)*Im[Λ2]/.rulePedestrian

Dissipation within the core and mantle

Finally, the membrane approach predicts that the global heat flow due to dissipation in the core and mantle is proportional to the imaginary part of the gravitational fluid-crust Love number.
The exact relation involves the global factor given above and the squared norm of the scaling factor ζ. The latter represents the relative reduction in radial displacement at the mantle-ocean boundary due to the presence of the crust (MW-95 and 96):

ζ=(1-3*(d/R)*((δρ)/(ρb)))((k2+1)/(k2FluidCrust+1))/.rulePedestrian

The global heat flow (in W) due to dissipation within the mantle is given by (MW-103):

MantleDissip=GlobalFac*(((Abs[ζ]))^2)*Im[k2FluidCrust]/.rulePedestrian

Heat budget

Now check that the crust and mantle contributions add up to the total heat flow:

CrustDissip+MantleDissip
TotalDissip

The heat budget is correct up to an error of 0.2%:

100*((TotalDissip-(CrustDissip+MantleDissip))/TotalDissip)

As noted above, dissipation within the mantle is in reality much smaller because the mantle is nearly elastic at tidal frequency.

Part II: Definitions

This part includes
1) a short description of the various symbols and
2) a section listing all the membrane formulas you could need to compute Love numbers

1. Notation

Beware:
- There is one important notation difference between the papers MP and MW.
  (OverBar[ρ]) denotes the bulk density in the former paper but the mean crust density in the latter.
  This notebook follows the latter convention.
- Equation numbers in MW refer to the published version or, equivalently, to the revised version of arXiv

A list of the parameters relevant to the shallow interior is given in MW-Table 4

Size (dimensional)
R: surface radius
Rm: mantle radius
d: crust thickness
dO: ocean thickness

Size (nondimensional)
ϵ=d/R: relative crust thickness (MW-28)
ϵO=dO/R: relative ocean thickness
fϵ=(Subscript[d, top])/d if the crust is made of two uniform layers of thicknesses (Subscript[d, top])and (Subscript[d, bot])
x=(R-d)/R: relative radius of the crust-ocean boundary with respect to the surface
y=Rm/R: relative radius of the mantle with respect to the surface
z=Rm/Rϵ: relative radius of the mantle with respect to the crust-ocean boundary (Rϵ=R-d)

Density
ρb: bulk density
ρ: density of the top layer of the ocean
ρbar: mean crust density (MW-31)
δρ=ρbar-ρ: crust-ocean density contrast (MW-32)
δρhat=δρ/ρ: crust-ocean density contrast, nondimensionalized by the density of the top layer of the ocean.

Gravity
g: surface gravity (MW-34)

Elasticity
μEmantle: elastic shear modulus of the mantle
μE: elastic shear modulus of the icy crust
νE: elastic Poisson’s ratio of the icy crust
χE: elastic compressibility factor of the icy crust (MW-12)
μEhat=μE/(ρ*g*R) : elastic shear modulus of the crust, μbar, nondimensionalized by ρgR (ρ is the ocean density)

Viscoelasticity
δ=(Subscript[μ, E])/(ω*η) : dimensionless parameter inversely proportional to tidal frequency ω and viscosity η (as in Wahr et al. 2009).
μ (or μbar if considered as a function): effective shear modulus of the icy crust (MW-42)
ν (or νbar if considered as a function): effective Poisson’s ratio of the icy crust (MW-43)
χ (or χbar if considered as a function): effective compressibility factor of the icy crust (MW-38 and 45)
μhat=μbar/(ρ*g*R) : the effective shear modulus of the crust, μbar, nondimensionalized by ρgR (MW-46)
μmHat=μm/(ρb0*g0*R) : nondimensional shear modulus of the mantle

Varia
n: harmonic degree
qω=(ω^2)R/g : dynamical parameter

Notation specific to fluid-crust models
In general, attaching ‘0’ to a variable means that we are dealing with a fluid-crust model instead of the physical model.
ρ0 : density of fluid crust.
ρb0 : bulk density of fluid-crust model.
ξ0 = ρ/ρb0 : ocean-to-bulk density ratio for the fluid-crust model (MW-33)
ξbar0 = ρ0/ρb0 : fluid-crust-to-bulk density ratio (MW-33)

2. Formulas

Initialization

Clear function names:

Clear[μV,νV,χV,pV,qV,average,μbar,χbar,νbar]
Clear[hRL0,kRL0,fn,pA,pB,An,Bn,hVL0,kVL0,pC,hRLL0,kRLL0]
Clear[xn,lnhn,fμ,fχ,fρ,Λ,Λχ,Λρ,ΛT,tilt,fρL,fχρ,psi]
Clear[ρb0,g0,Kρ,Hρ,kn,hn,knLoad,hnLoad]
Clear[Λωn,Kρr,Hρr,knrDyn,hnrDyn,Dresonant]
Clear[poly,Xa,Xb,Xc,Xd,Xe,k2HC,h2HC,k2DynHC,h2DynHC]

Maxwell rheology

Remark: ‘E’ and ‘V’ refer to ‘Elastic’ and ‘Viscoelastic’, respectively.

Viscoelastic shear modulus μ in the Fourier domain (MP-123; MW-A.2):

μV[μE_,δ_]:=((μE)/(1-I δ))

Viscoelastic Poisson’s ratio ν in the Fourier domain (MP-125; MW-A.2):

νV[νE_,δ_]:=((νE-I δ (1+νE)/3)/(1-2I δ (1+νE)/3))

Viscoelastic compressibility factor χ in the Fourier domain (MW-12):

χV[νE_,δ_]:=((1-2νV[νE,δ] )/(1-νV[νE,δ]))

Viscoelastic parameters (p, q) appearing in viscoelastic-gravitational equations (MW-36):

pV[μE_,νE_,δ_]:=((μV[μE,δ])/(1-νV[νE,δ]))
qV[μE_,νE_,δ_]:=pV[μE,νE,δ]*νV[νE,δ]

Two-layer crust with Maxwell rheology

Aim: rough model for conductive/convective crust in stagnant lid regime.
Elastic parameters (μE,νE) are uniform throughout the whole crust (thickness d).
Top layer with relative thickness f=(Subscript[d, top])/d and Maxwell viscoelasticity parameterized by (Subscript[δ, 1]).
Bottom layer with relative thickness (1-f)=(Subscript[d, bot])/d and Maxwell viscoelasticity parameterized by (Subscript[δ, 2]).

Weighted average of x (value in top layer) and y (value in bottom layer):

average[x_,y_,f_]:=f*x+(1-f)*y

Effective shear modulus (OverBar[μ]) of the membrane (MP-10; MW-42):

μbar[μE_,f_,δ1_,δ2_]:=average[μV[μE,δ1],μV[μE,δ2],f]

Effective Poisson’s ratio (OverBar[ν]) of the membrane (MP-7; MW-43):
Remark: I set μE=1 in the right-hand side because the ratio does not depend on μE

νbar[νE_,f_,δ1_,δ2_]:=((average[qV[1,νE,δ1],qV[1,νE,δ2],f])/(average[pV[1,νE,δ1],pV[1,νE,δ2],f]))

Effective compressibility factor (OverBar[χ]) of the membrane (MW-38):

χbar[νE_,f_,δ1_,δ2_]:=average[χV[νE,δ1],χV[νE,δ2],f]

Effective ‘load’ compressibility factor (Subscript[(OverBar[χ]), μ]) of the membrane (MW-48):

χμbar[μhat_,νE_,f_,δ1_,δ2_]:=average[((χV[νE,δ1])/(μV[μhat,δ1])),((χV[νE,δ2])/(μV[μhat,δ2])),f]

Fluid-crust Love numbers

Reminder:
The fluid-crust model is defined as having the same internal structure as the body to be modeled (called the physical model), except that the crust is fluid and of density ρ0.
There are two obvious choices for the fluid-crust density ρ0: it is equal either to the original crust density or to the density of the top layer of the ocean.
The latter choice makes it simpler to compute the Love numbers of the fluid-crust model, but it changes the bulk density of the body (except if crust density = ocean density).

Method used to obtain the formulas:
Incompressible propagator matrix method.

Relation between Love numbers:
If the surface is fluid, static Love numbers satisfy the hydrostatic relation (MP-35; MW-25):
h = k+1

Remark: the ' 0' attached to the functions serves to remind us that we are dealing with the fluid-crust model instead of the physical model.

RL (for Rigid/Liquid): Rigid mantle/Surface ocean (MP-2; MW-C.1):

hRL0[ξ0_,n_]:=(((1-(3/(2n+1))ξ0))^(-1))  
kRL0[ξ0_,n_]:=hRL0[ξ0,n]-1

VL (for Viscoelastic/Liquid): Viscoelastic mantle/Surface ocean (MP-118 to 120; MW-C.3 to C.7):

(*Auxiliary functions*)
fn[n_]:=(n/(2(n-1)(3+4n+2(n^2))))
pA[y_,ξ0_,n_]:=(2(n-1)+3(y^(2n+1)))(1-ξ0)+(2n+1)(y^3) ξ0
pB[y_,ξ0_,n_]:=(2n+1-3ξ0)(2(n-1)(1-ξ0)+(2n+1)(y^3)ξ0)-9(1-ξ0)ξ0 (y^(2n+1))
An[y_,ξ0_,n_]:=fn[n]*(2n+1)*(1-ξ0)*pA[y,ξ0,n]
Bn[y_,ξ0_,n_]:=fn[n]*(1-ξ0)*pB[y,ξ0,n]

(*Love numbers*)
hVL0[y_,ξ0_,μmHat_,n_]:=((An[y,ξ0,n]+(2n+1)(y^4)μmHat)/(Bn[y,ξ0,n]+(2n+1-3ξ0)(y^4)μmHat))
kVL0[y_,ξ0_,μmHat_,n_]:=hVL0[y,ξ0,μmHat,n]-1

RLL (for Rigid/Liquid/Liquid): Rigid mantle/Ocean/Fluid crust (MW-C.9 to C.11):

(*Auxiliary function*)
pC[x_,ξ0_,ξbar0_,n_]:=(2n+1)(1-(1-(x^3))ξbar0)+3(ξbar0-ξ0)(x^3)

(*Love numbers*)
hRLL0[x_,ξ0_,ξbar0_,n_]:=(2n+1)*(pC[x,ξ0,ξbar0,n]-3(ξbar0-ξ0)*(x^(2n+4)))/((2n+1-3ξbar0)*pC[x,ξ0,ξbar0,n]+9ξbar0(ξbar0-ξ0)*(x^(2n+4)))
kRLL0[x_,ξ0_,ξbar0_,n_]:=hRLL0[x,ξ0,ξbar0,n]-1

Relations between Love numbers

Eigenvalue of the operator -(Δ+2) (MW-19):

xn[n_]:=(n-1)(n+2)

(Subscript[l, n])-(Subscript[h, n]) relation (MP-25; MW-65):

lnhn[ν_,n_]:=((1+ν)/(xn[n]+1+ν))

Dimensionless coefficients appearing in the (Subscript[k, n])-(Subscript[h, n]) and (Subscript[k, n])'-(Subscript[h, n])' relations (MW-Table 5):

fμ[ν_,n_]:=((2xn[n](1+ν))/(xn[n]+1+ν))
fχ[ν_,n_]:=-((xn[n](1-ν))/(xn[n]+1+ν))
fρ[ν_,n_]:=(((n-1)((n^2)-3+(n+3)ν))/(xn[n]+1+ν))
fρL[ν_,n_]:=((((n^2)-1)(n+1-ν))/(xn[n]+1+ν))
fχρ[ν_,n_]:=((n(n+1)(1-ν))/(2(xn[n]+1+ν)))

Membrane spring constant (MP-28; MW-71):

Λ[μhat_,ν_,ϵ_,n_]:=fμ[ν,n]*μhat*ϵ

Compressibility correction (MW-72):

Λχ[ν_,χ_,ϵ_,n_]:=fχ[ν,n]*χ*ϵ

Minor density correction (MW-73):

Λρ[δρhat_,ν_,ϵ_,n_]:=fρ[ν,n]*δρhat*ϵ

Membrane spring constant plus compressibility, density, and dynamical corrections (MW-70):

ΛT[μhat_,ν_,χ_,δρhat_,Λω_,ϵ_,n_]:=Λ[μhat,ν,ϵ,n]+Λχ[ν,χ,ϵ,n]+Λρ[δρhat,ν,ϵ,n]+Λω

Tilt factor  (MW-69):

tilt[μhat_,ν_,χ_,δρhat_,Λω_,ϵ_,hnValue_,n_]:=ΛT[μhat,ν,χ,δρhat,Λω,ϵ,n]*hnValue-(2n+1)*δρhat*ϵ

Density and compressibility corrections for the (Subscript[k, n])'-(Subscript[h, n])' relation (MW-79):

psi[δρhat_,μhat_,ν_,χ_,χμ_,ϵ_,n_]:=(-fχ[ν,n]*χ+fρL[ν,n]*δρhat+((χμ)/2)+fχρ[ν,n]*((((δρhat+χ))^2)/(μhat)))*ϵ

Static Love numbers (gravitational and radial)

Bulk density and surface gravity of the fluid-crust model (MW-85):

ρb0[ρb_,δρ0_,ϵ_]:=ρb-3*δρ0*ϵ
g0[g_,ρb_,δρ0_,ϵ_]:=g*(1-3*((δρ0)/(ρb))*ϵ)

Density corrections in the Love number formulas (MW-99):

Kρ[δρhat0_,ϵ_,hn0_,n_]:=2*((n-1)hn0-(2n+1))*δρhat0*ϵ
Hρ[δρhat_,δρhat0_,ν_,ϵ_,ξ0_,hn0_,n_]:=Λρ[δρhat,ν,ϵ,n]-((2n+1)/hn0)*δρhat*ϵ+((3ξ0)/(2n+1))*Kρ[δρhat0,ϵ,hn0,n]

Formulas for gravitational and radial Love numbers (MW-98):
Assumption: the fluid-crust density is equal to the density of the top layer of the ocean: ρ0=ρ

kn[μhat_,ν_,χ_,δρhat_,ϵ_,ξ0_,hn0_,n_]:=hn0*(((1+((3ξ0)/(2n+1))*(hn0*((Λ[μhat,ν,ϵ,n])/(1+Λ[μhat,ν,ϵ,n]))+Kρ[δρhat,ϵ,hn0,n])))^(-1))-1
hn[μhat_,ν_,χ_,δρhat_,ϵ_,ξ0_,hn0_,n_]:=hn0*(((1+(1+((3ξ0)/(2n+1))*hn0)*Λ[μhat,ν,ϵ,n]+Λχ[ν,χ,ϵ,n]+Hρ[δρhat,δρhat,ν,ϵ,ξ0,hn0,n]))^(-1))

Load Love numbers (MW-80 and 81):

knLoad[μhat_,ν_,χ_,δρhat_,ϵ_,ξ0_,hn0_,n_]:=-1+ΛT[μhat,ν,χ,δρhat,0,ϵ,n]*hn[μhat,ν,χ,δρhat,ϵ,ξ0,hn0,n]-(2n+1)*δρhat*ϵ
hnLoad[μhat_,ν_,χ_,χμ_,δρhat_,ϵ_,ξ_,ξ0_,hn0_,n_]:=(ΛT[μhat,ν,χ,δρhat,0,ϵ,n]*hn[μhat,ν,χ,δρhat,ϵ,ξ0,hn0,n]-((2n+1)/(3ξ))*(1+psi[δρhat,μhat,ν,χ,χμ,ϵ,n]))/(1+ΛT[μhat,ν,χ,δρhat,0,ϵ,n])

Dynamical Love numbers (gravitational and radial)

Dynamical correction if rigid mantle model (MW-116):

Λωn[qω_,z_,n_]:=-((qω)/(n(n+1)))*((n+1+n*z^(2n+1))/(1-z^(2n+1)))

Density corrections in the dynamical Love number formulas, assuming the rigid mantle model (MW-120):

Kρr[δρhat_,ϵ_,ξ0_,n_]:=2*((n-1)*(1/(1-3ξ0/(2n+1)))-(2n+1))*δρhat*ϵ
Hρr[δρhat_,ν_,ϵ_,ξ0_,n_]:=Kρr[δρhat,ϵ,ξ0,n]+Λρ[δρhat,ν,ϵ,n]+3(1-ξ0)δρhat*ϵ

Formulas for the dynamical gravitational and radial Love numbers, assuming the rigid mantle model (MW-119):

knrDyn[μhat_,ν_,χ_,δρhat_,ϵ_,qω_,z_,ξ0_,n_]:=(1/(1-3ξ0/(2n+1)))*(((1+((3ξ0)/(2n+1))*(1/(1+Λ[μhat,ν,ϵ,n]+Λωn[qω,z,n]))*(((Λ[μhat,ν,ϵ,n]+Λωn[qω,z,n])/(1-3ξ0/(2n+1)))+Kρr[δρhat,ϵ,ξ0,n])))^(-1))-1
hnrDyn[μhat_,ν_,χ_,δρhat_,ϵ_,qω_,z_,ξ0_,n_]:=(1/(1-3ξ0/(2n+1)))*(((1+((Λ[μhat,ν,ϵ,n]+Λωn[qω,z,n])/(1-3ξ0/(2n+1)))+Λχ[ν,χ,ϵ,n]+Hρr[δρhat,ν,ϵ,ξ0,n]))^(-1))

Approximate resonant ocean thickness (MW-123):
Remark: neglect crust-ocean density contrast and compressibility factor, but keep the dependence of Λ on ν because it can be significant if Λ is large

Dresonant[R_,qω_,ξ_,μhat_,ν_,ϵ_,n_]:=R*((qω)/(n(n+1)))*(((1-((3ξ)/(2n+1))+Re[Λ[μhat,ν,ϵ,n]]))^(-1))

A simple thick shell model: the homogeneous crust model

Reference: MP-Appendix A and MW-Appendix E
Characteristics: homogeneous and incompressible crust, no density contrast at the crust-ocean boundary.
In the static case, one can solve the model in terms of fluid-crust Love numbers. which means that the structure below the crust can be freely chosen (MP-Appendix A).
In the dynamical case, an explicit solution can be obtained for the rigid mantle model: infinitely rigid mantle + homogeneous and incompressible ocean (MW-Appendix E).

Geometrical factors appearing in the Love number relations (MW-Table 11):
Recall that x=((R-d)/R) is the relative radius of the crust-ocean boundary

poly[x_]:=24+40 (x^3)-45 (x^7)-19 (x^10)
Xa[x_]:=(24/5)(19-75 (x^3)+112 (x^5)-75 (x^7)+19 (x^10))/poly[x]
Xb[x_]:=(x^2)(19+45 (x^3)-40 (x^7)-24 (x^10))/poly[x]
Xc[x_]:=(1/5)(36-100 (x^3)+308 (x^5)-225 (x^7)-19 (x^10))/poly[x]
Xd[x_]:=(1/2)*(x^2)(3+5 (x^3)-10 (x^7)+2 (x^10))/poly[x]
Xe[x_]:=5(x^2) (1-(x^3)-(x^7)+(x^10))/poly[x]

Factors appearing in the (Subscript[k, 2])-(Subscript[h, 2]) and (Subscript[l, 2])-(Subscript[h, 2]) formulas (MW-E.3):

zh[μhat_,Λω_,x_]:=((Xa[x]*μhat+Xb[x]*Λω)/(μhat+Xe[x]*Λω))
zl[μhat_,Λω_,x_]:=((Xc[x]*μhat+Xd[x]*Λω)/(μhat+Xe[x]*Λω))

Static Love numbers of degree 2 in terms of fluid-crust Love number h20 (MP-112):

k2HC[ξ_,μhat_,x_,h20_]:=h20*(((1+(3/5)*ξ*h20*((Xa[x]*μhat)/(1+Xa[x]*μhat))))^(-1))-1
h2HC[ξ_,μhat_,x_,h20_]:=h20*(((1+(1+(3/5)*ξ*h20)*Xa[x]*μhat))^(-1))

Dynamical Love numbers of degree 2 if rigid mantle model (MW-E.7):

k2DynHC[ξ_,μhat_,Λω_,x_]:=(3/5)ξ*(1/(1-(3/5)*ξ+zh[μhat,Λω,x]*μhat))
h2DynHC[ξ_,μhat_,Λω_,x_]:=(1/(1-(3/5)*ξ+zh[μhat,Λω,x]*μhat))

Part III: Other examples

Most examples are taken from papers MP and MW, to which the reader is referred for more explanations.

0. Constants

Set plotting options:

SetOptions[{Plot,LogLinearPlot},PlotStyle->{Thick,Black},Frame->True,FrameStyle->Directive[Thickness[0.003]],ImageSize->300,AspectRatio->1,Axes->None,BaseStyle->{14,FontFamily->"Arial"},PlotRange->All];

Clear variables that are left unspecified in this notebook:

Clear[μ,ν,χ,μE,νE,χE,fϵ,δ1,δ2,μEm,μEhat]
Clear[ρb,ρ,ρbar,δρ,δρhat]
Clear[ω,Rm,R,g,n]

Units:
radius: m
density: kg/(m^3)
surface gravity: m/(s^2)
shear modulus: Pa
angular frequency: (s^(-1))

Bulk and orbital parameters of Europa (MW-Table 2):

ωEuropa=2.048*10^-5;
RadiusEuropa=1560.8*10^3;
BulkDensityEuropa=3013.;
SurfaceGravityEuropa=1.315;
MantleRadiusEuropa=RadiusEuropa-170.*10^3(*MW-caption of Table 6*);
ruleEuropa={ω->ωEuropa,Rm->MantleRadiusEuropa,R->RadiusEuropa,ρb->BulkDensityEuropa,g->SurfaceGravityEuropa,μEm->μEmantle,n->2};
qωEuropa=(((ω^2)*R)/g)/.ruleEuropa

Bulk and orbital parameters of Titan (MW-Table 2):

ωTitan=0.456*10^-5;
RadiusTitan=2574.76*10^3;
BulkDensityTitan=1881.5;
SurfaceGravityTitan=1.354;
MantleRadiusTitan=RadiusTitan-350.*10^3(*MW-caption of Table 6*);
ruleTitan={ω->ωTitan,Rm->MantleRadiusTitan,R->RadiusTitan,ρb->BulkDensityTitan,g->SurfaceGravityTitan,μEm->μEmantle,n->2};
qωTitan=(((ω^2)*R)/g)/.ruleTitan

Elastic parameters (MP-Table 3 and MW-Table 7):

μEmantle=40.*10^9;
μEmantleSoft=2.*10^8;
μEice=3.5*10^9; 
νEice=0.33;
χEice=((1-2νEice)/(1-νEice));
ruleElastic={μ->μEice,ν->νEice,χ->χEice}

Density of pure water:

ρWater=1000.;

MP: Replacement rules and rheology models (MP-Table 4)

μEicehat=((μEice)/(ρWater*RadiusEuropa*SurfaceGravityEuropa));
ruleMP1={μEhat->μEicehat,νE->νEice,fϵ->0.4,n->2};
ruleMP2={{δ1->10^-7,δ2->0.01},{δ1->10^-7,δ2->1},{δ1->10^-7,δ2->100}};
μMP=μbar[μEhat,fϵ,δ1,δ2]/.ruleMP1/.ruleMP2
νMP=νbar[νE,fϵ,δ1,δ2]/.ruleMP1/.ruleMP2
χMP=χbar[νE,fϵ,δ1,δ2]/.ruleMP1/.ruleMP2

MW: Density models for the crust of Europa and Titan (MW-Table 6):

{ρI,ρO}={930.,1020.};
ruleL={ρbar->ρI,ρ->ρO,δρ->ρI-ρO,δρhat->((ρI-ρO)/(ρO))}
{ρI,ρO}={930.,1280.};
ruleM={ρbar->ρI,ρ->ρO,δρ->ρI-ρO,δρhat->((ρI-ρO)/(ρO))}
{ρI,ρO}={1167.,1280.};
ruleD={ρbar->ρI,ρ->ρO,δρ->ρI-ρO,δρhat->((ρI-ρO)/(ρO))}
ruleLMD={ruleL,ruleM,ruleD};

MW: Rheology models for Europa’s crust (MW-Tables 7 and 8):

{δ1E,δ1C,δ1F}={1,1,1}*10^-7;
{δ2E,δ2C,δ2F}={0.1,1.,10.};
{fE,fC,fF}={0.4,0.4,0.1};
ruleE={μE->μEice,νE->νEice,fϵ->fE,δ1->δ1E,δ2->δ2E};
ruleC={μE->μEice,νE->νEice,fϵ->fC,δ1->δ1C,δ2->δ2C};
ruleF={μE->μEice,νE->νEice,fϵ->fF,δ1->δ1F,δ2->δ2F};
ruleECF={ruleE,ruleC,ruleF};

The corresponding effective viscoelastic parameters are those of MW-Table 8:

μECF=μbar[μE,fϵ,δ1,δ2]/.ruleECF;
νECF=νbar[νE,fϵ,δ1,δ2]/.ruleECF;
χECF=χbar[νE,fϵ,δ1,δ2]/.ruleECF;
TableForm[Transpose@{μECF/μEice,νECF,χECF},TableHeadings->{{"Elastic-like","Critical","Fluid-like"},{"(OverBar[μ])/μE","(OverBar[ν])","(OverBar[χ])"}}]

The associated replacement rules are:

ruleEbar={μ->μECF[[1]],ν->νECF[[1]],χ->χECF[[1]]};
ruleCbar={μ->μECF[[2]],ν->νECF[[2]],χ->χECF[[2]]};
ruleFbar={μ->μECF[[3]],ν->νECF[[3]],χ->χECF[[3]]};

1. Effective viscoelastic parameters

Example 1.1: viscoelastic parameters of a uniform layer with Maxwell rheology (MP-Fig 1)

Recall that δ=(Subscript[μ, E])/(ω*η) is a dimensionless parameter inversely proportional to tidal frequency ω and viscosity η;
δ  quantifies viscous effects in Maxwell rheology.

Here are the formulas for the viscoelastic shear modulus μ, Poisson’s ratio ν, and compressibility factor χ in the Fourier domain (‘E’ refers to ‘Elastic’ and ‘V’ to ‘Viscoelastic’):

{μV[μE,δ],νV[νE,δ],χV[νE,δ]}//Simplify

Plot the viscoelastic parameters (μ, χ, ν) as functions of δ (MP-Fig 1).
Remark: the shear modulus is divided its elastic value.

{δMin,δMax}={10^-2,10^2};
Fig1=LogLinearPlot[{Re@μV[1, δ], Im@μV[1, δ]}, {δ, δMin,δMax}, PlotStyle -> {Blue, Red},FrameLabel->{"δ","μ/(Subscript[μ, E])","(blue,red)=(Re,Im)",None}];
Fig2=LogLinearPlot[{Re@χV[νEice, δ], Im@χV[νEice, δ]}, {δ,δMin,δMax}, PlotStyle -> {Blue, Red},FrameLabel->{"δ","χ","(blue,red)=(Re,Im)",None}];
Fig3=LogLinearPlot[{Re@νV[νEice, δ], Im@νV[νEice, δ]}, {δ, δMin,δMax}, PlotStyle -> {Blue, Red},FrameLabel->{"δ","ν","(blue,red)=(Re,Im)",None}];
Show[GraphicsRow[ {Fig1,Fig2,Fig3}]]

Example 1.2: effective viscoelastic parameters of a two-layer crust (MP-Fig 2)

In this example, elastic parameters (μ,ν) are uniform throughout the whole crust (thickness d);
Top layer with relative thickness f=(Subscript[d, top])/d and Maxwell viscoelasticity parameterized by (Subscript[δ, 1]).
Bottom layer with relative thickness (1-f)=(Subscript[d, bot])/d and Maxwell viscoelasticity parameterized by (Subscript[δ, 2]).

The top layer thickness is 0.4 times the total thickness:

fd=0.4;

The viscosity of the bottom layer is (10^6) smaller than the viscosity of the top layer, so that (Subscript[δ, 2])=mult*(Subscript[δ, 1]) with:

mult=10^6;

Now, plot effective viscoelastic parameters ((OverBar[μ]), (OverBar[χ]), (OverBar[ν])) as functions of (Subscript[δ, 1]) (MP-Fig 2):

{δ2Min,δ2Max}={10^-8,10^2};
Fig1=LogLinearPlot[{Re@μbar[1,fd,δ1,mult*δ1], Im@μbar[1,fd, δ1,mult*δ1]}, {δ1, δ2Min,δ2Max}, PlotStyle -> {Blue, Red},FrameLabel->{"(Subscript[δ, 1])","\!\(\*OverscriptBox[\(\(μ\)\(\\ \)\), \(_\)]\)/ (Subscript[μ, E])","(blue,red)=(Re,Im)",None}];
Fig2=LogLinearPlot[{Re@χbar[νEice,fd, δ1,mult*δ1], Im@χbar[νEice,fd, δ1,mult*δ1]}, {δ1,δ2Min,δ2Max}, PlotStyle -> {Blue, Red},FrameLabel->{"(Subscript[δ, 1])","(OverBar[χ])","(blue,red)=(Re,Im)",None}];
Fig3=LogLinearPlot[{Re@νbar[νEice,fd, δ1,mult*δ1], Im@νbar[νEice,fd, δ1,mult*δ1]}, {δ1, δ2Min,δ2Max}, PlotStyle -> {Blue, Red},FrameLabel->{"(Subscript[δ, 1])","(OverBar[ν])","(blue,red)=(Re,Im)",None}];
Show[GraphicsRow[ {Fig1,Fig2,Fig3}]]

Example 1.3: effective viscoelastic parameters of a two-layer crust (MW-Fig 3)

As in Example 1.2, except that the top layer remains purely elastic while the (Subscript[δ, 2]) parameter varies.

The top layer thickness is 0.4 times the total thickness:

fd=0.4;

The top layer is elastic so (Subscript[δ, 1])=0:

δ1elastic=0;

Critical transitions for χ and ν (MW-Section 3.4):

{δp,δpp}={3*((1-ν)/(1+ν)),(Sqrt[(3/fd)*((1-ν)/(1+ν))])}/.ruleElastic

Now, plot effective viscoelastic parameters ((OverBar[μ]), (OverBar[χ]), (OverBar[ν])) as functions of (Subscript[δ, 2]) (MW-Fig 3):

{δ2Min,δ2Max}={10^-2,10^3};
Fig1=LogLinearPlot[{Re@μbar[1,fd, δ1elastic,δ2], Im@μbar[1,fd, δ1elastic,δ2]}, {δ2, δ2Min,δ2Max}, PlotStyle -> {Blue, Red},FrameLabel->{"(Subscript[δ, 2])","\!\(\*OverscriptBox[\(\(μ\)\(\\ \)\), \(_\)]\)/ (Subscript[μ, E])","(blue,red)=(Re,Im)",None},Axes->True,AxesOrigin->{1,0}];
Fig2=LogLinearPlot[{Re@χbar[νEice,fd, δ1elastic,δ2], Im@χbar[νEice,fd, δ1elastic,δ2]}, {δ2,δ2Min,δ2Max}, PlotStyle -> {Blue, Red},FrameLabel->{"(Subscript[δ, 2])","(OverBar[χ])","(blue,red)=(Re,Im)",None},Axes->True,AxesOrigin->{δp,0}];
Fig3=LogLinearPlot[{Re@νbar[νEice,fd, δ1elastic,δ2], Im@νbar[νEice,fd, δ1elastic,δ2]}, {δ2, δ2Min,δ2Max}, PlotStyle -> {Blue, Red},FrameLabel->{"(Subscript[δ, 2])","(OverBar[ν])","(blue,red)=(Re,Im)",None},Axes->True,AxesOrigin->{δpp,-0.05}];
Show[GraphicsRow[ {Fig1,Fig2,Fig3}]]

For completeness, plot the effective ‘load’ compressibility parameter (Subscript[(OverBar[χ]), μ]) as a function of (Subscript[δ, 2]) (this figure is not in the paper MW);
assume that (Subscript[μ, E])/(ρgR)=1:

LogLinearPlot[{Re@χμbar[1,νEice,fd, δ1elastic,δ2], Im@χμbar[1,νEice,fd, δ1elastic,δ2]}, {δ2, δ2Min,δ2Max}, PlotStyle -> {Blue, Red},FrameLabel->{"(Subscript[δ, 2])","(Subscript[(OverBar[χ]), μ])","(blue,red)=(Re,Im)",None},Axes->True,AxesOrigin->{δp,-0.05}]

2. (Subscript[l, 2])-(Subscript[h, 2]) relation

Example 2.1: (Subscript[l, 2])-(Subscript[h, 2]) ratio as a function of the inverse frequency/inverse viscosity of the crust (MP-Fig 3)

The crust is divided into two layers, the top layer thickness being 0.4 times the total thickness:

fd=0.4;

The viscosity of the bottom layer is (10^6) smaller than the viscosity of the top layer, so that (Subscript[δ, 2])=mult*(Subscript[δ, 1]) with:

mult=10*^6;

Now, plot the (Subscript[l, 2])/(Subscript[h, 2]) ratio as a function of (Subscript[δ, 1]) (MP-Fig 3):

{δ1Min,δ1Max}={10^-8,10^2};
Fig1=LogLinearPlot[Re@lnhn[νbar[νEice,fd, δ1,mult*δ1],2], {δ1, δ1Min,δ1Max}, PlotStyle ->Blue,FrameLabel->{"(Subscript[δ, 1])","Re ((Subscript[l, 2])/(Subscript[h, 2]))"},AspectRatio->0.75,ImageSize->370];
Fig2=LogLinearPlot[Im@lnhn[νbar[νEice,fd, δ1,mult*δ1],2], {δ1, δ1Min,δ1Max}, PlotStyle -> Red,FrameLabel->{"(Subscript[δ, 1])","Im ((Subscript[l, 2])/(Subscript[h, 2]))"},AspectRatio->0.75,ImageSize->370];
Show[GraphicsRow[{Fig1,Fig2}]]

Example 2.2: (Subscript[l, 2])-(Subscript[h, 2]) ratio as a function of crust thickness (MP-Fig 4)

The membrane formula for the (Subscript[l, 2])/(Subscript[h, 2]) ratio is:

lnhn[ν,2]

First, compute the membrane approximation of the  (Subscript[l, 2])/(Subscript[h, 2]) ratio for the rheology models (H,E,C,F) of the paper MP.
In MP and MW, the membrane approximation of (Subscript[l, 2])/(Subscript[h, 2]) is only computed at order 1 in ϵ=d/R.

zlH=lnhn[1/2,2];
zlE=lnhn[νMP[[1]],2];
zlC=lnhn[νMP[[2]],2];
zlF=lnhn[νMP[[3]],2];
{zlH,zlE,zlC,zlF}

For comparison, consider (Subscript[l, 2])/(Subscript[h, 2]) ratio of the homogeneous crust model.
This model assumes that the crust is homogeneous and incompressible, and is not appropriate for the 2-layer rheology.
The static (Subscript[l, 2])/(Subscript[h, 2]) ratio for this model is given by zl[1,0,x]=Xc[x]

zl[1,0,x]

Finally, plot the (Subscript[l, 2])/(Subscript[h, 2]) ratio as a function of crust thickness for a two-layer crust (MP-Fig 4).

R0=RadiusEuropa/1000.;
{dMin,dMax}={0.01,170.};
Show[
Plot[Re@zl[1,0,1-d/R0],{d,dMin,dMax},PlotStyle->{Darker@Green,Dashed,Thick}],
Graphics[{PointSize[0.025],Darker@Green,Point[{0.,Re@zlH}]}],
Graphics[{PointSize[0.025],Black,Point[{0,Re@zlE}]}],
Graphics[{PointSize[0.025],Red,Point[{0.,Re@zlC}]}],
Graphics[{PointSize[0.025],Blue,Point[{5.,Re@zlF}]}],
FrameLabel->{"Crust thickness [km]","Re ((Subscript[l, 2])/(Subscript[h, 2]))","Elastic/Critical/Fluid = black/red/blue
Homogeneous and incompressible = green",None},
PlotRange->{{-5,185},{0.24,0.28}},AspectRatio->0.75,ImageSize->350]

3. (Subscript[k, 2])-(Subscript[h, 2]) relation

Example 3.1: normalized tilt factor in function of the inverse frequency/inverse viscosity of the crust (MP-Fig 5)

For safety, clear variables to be replaced:

Clear[ν,χ,δρhat,Λω,ϵ,n,μEhat,νE,fϵ,δ1,δ2];

Assume static tides (Λω=0) of degree 2 and crust density=ocean density=1000 kg/m3 (MP-Table 3);
Moreover there is no compressibility factor in the paper MP: χ=0.
Under these assumptions, the tilt factor divided by (Subscript[h, 2]) (called normalized tilt factor, see MP-37) is given by

ΛT[(OverHat[μ]),ν,χ,δρhat,Λω,ϵ,n]/.{χ->0,δρhat->0,Λω->0,n->2}

which is equal to Λ:

Λ[(OverHat[μ]),ν,ϵ,n]/.{n->2}

First, define replacement rules for rheology model of Europa’s crust (MP-Table 3)

μEicehat=((μEice)/(ρWater*RadiusEuropa*SurfaceGravityEuropa));
mult=10^6;(*the viscosities of the two crust layers are in a constant ratio*)
ruleMPfig5={μEhat->μEicehat,νE->νEice,fϵ->fE,δ2->mult*δ1,δρhat->0,Λω->0,ϵ->1,n->2}(*ϵ=1 because Λ is in units of d/R*);

Next, define Λ as a one-variable function just for this plot:

Clear[ΛMP];
ΛMP[δ1_]=Λ[μbar[μEhat,fϵ,δ1,δ2],νbar[νE,fϵ,δ1,δ2],ϵ,n]/.ruleMPfig5;

For comparison, plot also the curves for (Subscript[Λ, T]) , i.e. with the correct compressibility correction (χ≠0).
For this purpose, define (Subscript[Λ, T]) as a one-variable function just for this plot:

Clear[ΛMPcomp];
ΛMPcomp[δ1_]=ΛT[μbar[μEhat,fϵ,δ1,δ2],νbar[νE,fϵ,δ1,δ2],χbar[νE,fϵ,δ1,δ2],δρhat,Λω,ϵ,n]/.ruleMPfig5;

Finally, plot the normalized tilt factor as a function of (Subscript[δ, 1]) (MP-Figure 5):

{δ1Min,δ1Max}={10^-8,10^2};
Fig1=Show[
LogLinearPlot[Re@ΛMP[δ1],{δ1, δ1Min,δ1Max},PlotStyle->{Blue,Thick}],
LogLinearPlot[Re@ΛMPcomp[δ1],{δ1, δ1Min,δ1Max},PlotStyle->{Darker@Green,Thick}],
FrameLabel->{"(Subscript[δ, 1])","Re ((Subscript[γ, 2])/(Subscript[h, 2]))  [units of d/R]","χ=0 (blue), χ≠0 (green)",None},PlotRange->{-0.3,3.7},AspectRatio->0.75,ImageSize->370
];
Fig2=Show[
LogLinearPlot[Im@ΛMP[δ1],{δ1,δ1Min,δ1Max},PlotStyle->{Red,Thick}],
LogLinearPlot[Im@ΛMPcomp[δ1],{δ1,δ1Min,δ1Max},PlotStyle->{Orange,Thick}],
FrameLabel->{"(Subscript[δ, 1])","Im ((Subscript[γ, 2])/(Subscript[h, 2])) [units of d/R]","χ=0 (red), χ≠0 (orange)",None},PlotRange->{-0.1,1.1},AspectRatio->0.75,ImageSize->370
];
Show[GraphicsRow[{Fig1,Fig2}]]

Example 3.2: normalized tilt factor as a function of crust thickness (MP-Fig 6)

Same assumptions as for Example 3.1.
Plot normalized tilt factor as a function of crust thickness for different crustal rheology models.
For comparison, plot also the curves for (Subscript[Λ, T]) , i.e. with the correct compressibility correction (χ≠0):

{δρhat0,Λω0,n0}={0,0,2};
{dMin,dMax}={0.01,170.};
R0=RadiusEuropa/1000;
Show[
Plot[Re@Λ[μEicehat,0.5,d/R0,n0],{d,dMin,dMax},PlotStyle->Darker@Green],
Plot[Re@Λ[μMP[[1]],νMP[[1]],d/R0,n0],{d,dMin,dMax},PlotStyle->Black],
Plot[Re@Λ[μMP[[2]],νMP[[2]],d/R0,n0],{d,dMin,dMax},PlotStyle->Red],
Plot[Re@Λ[μMP[[3]],νMP[[3]],d/R0,n0],{d,dMin,dMax},PlotStyle->Blue],
Plot[Re@ΛT[μMP[[1]],νMP[[1]],χMP[[1]],δρhat0,Λω0,d/R0,n0],{d,dMin,dMax},PlotStyle->{Black,Thick,Dotted}],
Plot[Re@ΛT[μMP[[2]],νMP[[2]],χMP[[1]],δρhat0,Λω0,d/R0,n0],{d,dMin,dMax},PlotStyle->{Red,Thick,Dotted}],
Plot[Re@ΛT[μMP[[3]],νMP[[3]],χMP[[3]],δρhat0,Λω0,d/R0,n0],{d,dMin,dMax},PlotStyle->{Blue,Thick,Dotted}],
FrameLabel->{"Crust thickness [km]","Re ((Subscript[γ, 2])/(Subscript[h, 2]))","χ=0 (solid), χ≠0 (dotted)"},
PlotRange->{-0.02,0.46},AspectRatio->0.8,ImageSize->350
]

Example 3.3: slope of tilt factor at zero crust thickness (MW-Fig 7)

For safety, clear variables to be replaced:

Clear[ν,χ,Λω,ϵ,n];
Clear[Rm,R,ρ,ρb,g,μEm];

Here are the formulas for the degree-two dimensionless coefficients appearing in the tilt factor formula (MW-Table 5):

{fμ[ν,n],fχ[ν,n],fρ[ν,n]}/.{n->2}

Here is the formula for the degree-two membrane spring constant modified by corrections (MW-70) :

ΛT[(OverHat[μ]),ν,χ,(OverHat[δρ]),Λω,ϵ,n]/.{n->2}

First, compute the fluid-crust Love numbers at zero crust thickness with the VL model (MW-Figure 7, see caption):
Remark: since the crust is of zero thickness, the fluid-crust model and the physical model have the same bulk density.

h20=hVL0[(Rm/R),(ρ/(ρb)),((μEm)/(ρb*g*R)),n]/.ruleEuropa;
{h20L,h20M,h20D}=h20/.ruleLMD

Next, compute the various approximations for the slope (at zero crust thickness) of Europa’s degree-two tilt factor, for the elastic-like rheology model:

approx={0,0,0,0};
approx[[1]]=First@Limit[(ΛT[((μECF[[1]])/(ρ*g*R)),νECF[[1]],0,0,0,ϵ,n]*(h20/ϵ))/.ruleEuropa,{ϵ->0}];
approx[[2]]=First@Limit[(ΛT[((μECF[[1]])/(ρ*g*R)),νECF[[1]],χECF[[1]],0,0,ϵ,n]*(h20/ϵ))/.ruleEuropa,{ϵ->0}];
approx[[3]]=First@Limit[(ΛT[((μECF[[1]])/(ρ*g*R)),νECF[[1]],χECF[[1]],δρhat,0,ϵ,n]*(h20/ϵ))/.ruleEuropa,{ϵ->0}];
approx[[4]]=First@Limit[(tilt[((μECF[[1]])/(ρ*g*R)),νECF[[1]],χECF[[1]],δρhat,0,ϵ,h20,n]*(1/ϵ))/.ruleEuropa,{ϵ->0}];
caseE=approx/.ruleLMD;

Do the same for the critical rheology model:

approx={0,0,0,0};
approx[[1]]=First@Limit[(ΛT[((μECF[[2]])/(ρ*g*R)),νECF[[2]],0,0,0,ϵ,n]*(h20/ϵ))/.ruleEuropa,{ϵ->0}];
approx[[2]]=First@Limit[(ΛT[((μECF[[2]])/(ρ*g*R)),νECF[[2]],χECF[[2]],0,0,ϵ,n]*(h20/ϵ))/.ruleEuropa,{ϵ->0}];
approx[[3]]=First@Limit[(ΛT[((μECF[[2]])/(ρ*g*R)),νECF[[2]],χECF[[2]],δρhat,0,ϵ,n]*(h20/ϵ))/.ruleEuropa,{ϵ->0}];
approx[[4]]=First@Limit[(tilt[((μECF[[2]])/(ρ*g*R)),νECF[[2]],χECF[[2]],δρhat,0,ϵ,h20,n]*(1/ϵ))/.ruleEuropa,{ϵ->0}];
caseC=approx/.ruleLMD;

Do the same for the fluid-like rheology model:

approx={0,0,0,0};
approx[[1]]=First@Limit[(ΛT[((μECF[[3]])/(ρ*g*R)),νECF[[3]],0,0,0,ϵ,n]*(h20/ϵ))/.ruleEuropa,{ϵ->0}];
approx[[2]]=First@Limit[(ΛT[((μECF[[3]])/(ρ*g*R)),νECF[[3]],χECF[[3]],0,0,ϵ,n]*(h20/ϵ))/.ruleEuropa,{ϵ->0}];
approx[[3]]=First@Limit[(ΛT[((μECF[[3]])/(ρ*g*R)),νECF[[3]],χECF[[3]],δρhat,0,ϵ,n]*(h20/ϵ))/.ruleEuropa,{ϵ->0}];
approx[[4]]=First@Limit[(tilt[((μECF[[3]])/(ρ*g*R)),νECF[[3]],χECF[[3]],δρhat,0,ϵ,h20,n]*(1/ϵ))/.ruleEuropa,{ϵ->0}];
caseF=approx/.ruleLMD;

As the instructions for Figure 7 of MW are rather long, I give instead the numerical values:

Print["ELASTIC-LIKE: LEFT PANELS OF MW-FIG 7"]
TableForm[caseE,TableHeadings->{{"Light","Mixed","Dense"},None}]
Print["CRITICAL: MIDDLE PANELS OF MW-FIG 7"]
TableForm[caseC,TableHeadings->{{"Light","Mixed","Dense"},None}]
Print["FLUID-LIKE: RIGHT PANELS OF MW-FIG 7"]
TableForm[caseF,TableHeadings->{{"Light","Mixed","Dense"},None}]

4. Static Love numbers

Example 4.1: Radial Love number of Europa as a function of the inverse frequency/inverse viscosity of the crust, for diurnal tides (MP-Fig 7A)

There is no density contrast between crust and ocean, so that the fluid-crust Love numbers are computed with the bulk density and surface gravity of the physical model.

For safety, clear variables to be replaced:

Clear[ξ,μEhat,νE,fϵ,δ1,δ2,δρhat,Λω,ϵ,n,hn0];
Clear[Rm,R,ρ,ρb,g,hn0];

First, define the replacement rules for the Europa model considered in the MP paper (MP-Tables 2, 3, and 4):

mult=10^6;(*the viscosities of the two crust layers are in a constant ratio*)
d0=20.*10^3;(*the crust is 20 km thick*)
μEicehat=((μEice)/(ρWater*RadiusEuropa*SurfaceGravityEuropa));(*nondimensional shear modulus of elastic ice*)
ruleMPfig7={ξ->((ρWater)/BulkDensityEuropa),μEhat->μEicehat,νE->νEice,fϵ->fE,δ2->mult*δ1,δρhat->0,Λω->0,ϵ->(d0/RadiusEuropa),n->2};

Next, compute the fluid-crust Love numbers for diurnal tides: elastic mantle (h20E) or rigid mantle (h20R):

h20E=hVL0[(Rm/R),(ρ/(ρb)),((μEm)/(ρb*g*R)),n]/.ruleEuropa/.{ρ->ρWater}
h20R=hRL0[(ρ/(ρb)),n]/.ruleEuropa/.{ρ->ρWater}

Then, express the radial Love number as a function of (Subscript[δ, 1]) and another variable called χfac:
I leave the possibility to include (χfac=1) or exclude (χfac=0) the compressibility factor:

Clear[h2E,h2R];
h2E[δ1_,χfac_]=hn[μbar[μEhat,fϵ,δ1,δ2],νbar[νE,fϵ,δ1,δ2],χfac*χbar[νE,fϵ,δ1,δ2],δρhat,ϵ,ξ,hn0,n]/.ruleMPfig7/.{hn0->h20E};
h2R[δ1_,χfac_]=hn[μbar[μEhat,fϵ,δ1,δ2],νbar[νE,fϵ,δ1,δ2],χfac*χbar[νE,fϵ,δ1,δ2],δρhat,ϵ,ξ,hn0,n]/.ruleMPfig7/.{hn0->h20R};

Finally, plot radial Love number for diurnal tides:

{δ1Min,δ1Max}={10^-8,10^2};
Fig1=Show[
LogLinearPlot[Re@h2E[δ1,0],{δ1,δ1Min,δ1Max},PlotStyle->{Blue,Thick}],
LogLinearPlot[Re@h2R[δ1,0],{δ1,δ1Min,δ1Max},PlotStyle->{Blue,Dotted,Thick}],
LogLinearPlot[Re@h2E[δ1,1],{δ1,δ1Min,δ1Max},PlotStyle->{Darker@Green,Thick}],
LogLinearPlot[Re@h2R[δ1,1],{δ1,δ1Min,δ1Max},PlotStyle->{Darker@Green,Dotted,Thick}],
Graphics[{Text[Style["rigid mantle",Black],{Log[0.001],1.217}]}],
Graphics[{Text[Style["elastic mantle",Black],{Log[0.001],1.243}]}],
FrameLabel->{"(Subscript[δ, 1])","Re ((Subscript[h, 2]))","Diurnal tides: χ=0 (blue), χ≠0 (green)",None},PlotRange->{1.18,1.27},AspectRatio->0.75,ImageSize->400
];
Fig2=Show[
LogLinearPlot[Im@h2E[δ1,1],{δ1,δ1Min,δ1Max},PlotStyle->{Red}],
LogLinearPlot[Im@h2R[δ1,1],{δ1,δ1Min,δ1Max},PlotStyle->{Red,Dotted,Thick}],
LogLinearPlot[Im@h2E[δ1,0],{δ1,δ1Min,δ1Max},PlotStyle->{Orange}],
LogLinearPlot[Im@h2R[δ1,0],{δ1,δ1Min,δ1Max},PlotStyle->{Orange,Dotted,Thick}],
FrameLabel->{"(Subscript[δ, 1])","Im ((Subscript[h, 2]))","Diurnal tides: χ=0 (red), χ≠0 (orange)",None},PlotRange->{-0.021,0.002},AspectRatio->0.75,ImageSize->400
];
Show[GraphicsRow[{Fig1,Fig2}]]

Example 4.2: Radial Love number of Europa as a function of the inverse frequency/inverse viscosity of the crust, for NSR tides (MP-Fig 7B)

The replacement rule is the same as in Example 4.1.

First, compute the fluid-crust Love numbers for NSR tides (due to nonsynchronous rotation);
In comparison with diurnal tides, the shear modulus of the mantle is reduced: consider either a soft mantle with μ=0.2 GPa (h20S) or a completely fluid mantle (h20F):

Clear[μEmNSR];
h20S=hVL0[(Rm/R),(ρ/(ρb)),((μEmNSR)/(ρb*g*R)),n]/.ruleEuropa/.{ρ->ρWater}/.{μEmNSR->μEmantleSoft};
h20F=hVL0[(Rm/R),(ρ/(ρb)),0,n]/.ruleEuropa/.{ρ->ρWater};

Next, express the radial Love number as a function of (Subscript[δ, 1]) and another variable called χfac:
I leave the possibility to include (χfac=1) or exclude (χfac=0) the compressibility factor:

Clear[h2S,h2F];
h2S[δ1_,χfac_]=hn[μbar[μEhat,fϵ,δ1,δ2],νbar[νE,fϵ,δ1,δ2],χfac*χbar[νE,fϵ,δ1,δ2],δρhat,ϵ,ξ,hn0,n]/.ruleMPfig7/.{hn0->h20S};
h2F[δ1_,χfac_]=hn[μbar[μEhat,fϵ,δ1,δ2],νbar[νE,fϵ,δ1,δ2],χfac*χbar[νE,fϵ,δ1,δ2],δρhat,ϵ,ξ,hn0,n]/.ruleMPfig7/.{hn0->h20F};

Finally, plot radial Love number for NSR tides:

{δ1Min,δ1Max}={10^-8,10^2};
Fig1=Show[
LogLinearPlot[Re@h2S[δ1,0],{δ1,δ1Min,δ1Max},PlotStyle->{Blue,Thick}],
LogLinearPlot[Re@h2F[δ1,0],{δ1,δ1Min,δ1Max},PlotStyle->{Blue,Dotted,Thick}],
LogLinearPlot[Re@h2S[δ1,1],{δ1,δ1Min,δ1Max},PlotStyle->{Darker@Green,Thick}],
LogLinearPlot[Re@h2F[δ1,1],{δ1,δ1Min,δ1Max},PlotStyle->{Darker@Green,Dotted,Thick}],
Graphics[{Text[Style["soft mantle",Black],{Log[0.001],1.81}]}],
Graphics[{Text[Style["fluid mantle",Black],{Log[0.001],1.965}]}],
FrameLabel->{"(Subscript[δ, 1])","Re ((Subscript[h, 2]))","NSR tides: χ=0 (blue), χ≠0 (green)",None},PlotRange->{1.74,2.06},AspectRatio->0.75,ImageSize->400
];
Fig2=Show[
LogLinearPlot[Im@h2S[δ1,1],{δ1,δ1Min,δ1Max},PlotStyle->Red],
LogLinearPlot[Im@h2F[δ1,1],{δ1,δ1Min,δ1Max},PlotStyle->{Red,Dotted,Thick}],
LogLinearPlot[Im@h2S[δ1,0],{δ1,δ1Min,δ1Max},PlotStyle->Orange],
LogLinearPlot[Im@h2F[δ1,0],{δ1,δ1Min,δ1Max},PlotStyle->{Orange,Dotted,Thick}],
FrameLabel->{"(Subscript[δ, 1])","Im ((Subscript[h, 2]))","NSR tides: χ=0 (red), χ≠0 (orange)",None},PlotRange->{-0.036,0.002},AspectRatio->0.75,ImageSize->400
];
Show[GraphicsRow[{Fig1,Fig2}]]

Example 4.3: Tidal Love numbers of Titan as functions of crust thickness for different crust/ocean density models (MW-Fig 5)

For safety, clear variables to be replaced:

Clear[Rm,R,g,ϵ,hn0,n]
Clear[ρ,ρb,δρ,δρhat]
Clear[μ,ν,χ,μEm]

First, evaluate partially the fluid-crust Love number h20 with the VL model, using the replacement rule ‘ruleTitan’;
the core-mantle layer is assumed to be elastic;
don’t forget to use the fluid-crust bulk density and gravity;
h20 now depends on the three variables (ρ, δρ, ϵ):

h20=hVL0[(Rm/R),(ρ/(ρb0[ρb,δρ,ϵ])),((μEm)/(ρb0[ρb,δρ,ϵ]*g0[g,ρb,δρ,ϵ]*R)),n]/.ruleTitan;

Next, define the gravitational and radial Love numbers (degree 2) of Titan as functions of crust thickness;
As the crust is elastic, use the replacement rule ‘ruleElastic’.
Define (Subscript[k, n])and (Subscript[h, n]) functions for the three density models (Light, Mixed, Dense), using the replacement rule ‘ruleLMD’:

Clear[knL,knM,knD,hnL,hnM,hnD]
{knL[ϵ_],knM[ϵ_],knD[ϵ_]}=kn[(μ/(ρ*g*R)),ν,χ,δρhat,ϵ,(ρ/(ρb0[ρb,δρ,ϵ])),hn0,n]/.ruleTitan/.ruleElastic/.{hn0->h20}/.ruleLMD;
{hnL[ϵ_],hnM[ϵ_],hnD[ϵ_]}=hn[(μ/(ρ*g*R)),ν,χ,δρhat,ϵ,(ρ/(ρb0[ρb,δρ,ϵ])),hn0,n]/.ruleTitan/.ruleElastic/.{hn0->h20}/.ruleLMD;

Plot the Love numbers as functions of crust thickness (MW-Fig 5):

{ϵMin,ϵMax}=100*{0,0.101};
Fig1=Show[
Plot[knL[ϵ/100],{ϵ,ϵMin,ϵMax},PlotStyle->{Black,Thick}],
Plot[knM[ϵ/100],{ϵ,ϵMin,ϵMax},PlotStyle->{Red,Thick}],
Plot[knD[ϵ/100],{ϵ,ϵMin,ϵMax},PlotStyle->{Blue,Thick}],
FrameLabel->{"Relative crust thickness [%]","(Subscript[k, 2])","Light/Mixed/Dense = black/red/blue", None},
PlotRange->{0.32,0.721},ImageSize->320
];
Fig2=Show[
Plot[hnL[ϵ/100],{ϵ,ϵMin,ϵMax},PlotStyle->{Black,Thick}],
Plot[hnM[ϵ/100],{ϵ,ϵMin,ϵMax},PlotStyle->{Red,Thick}],
Plot[hnD[ϵ/100],{ϵ,ϵMin,ϵMax},PlotStyle->{Blue,Thick}],
FrameLabel->{"Relative crust thickness [%]","(Subscript[h, 2])","Light/Mixed/Dense = black/red/blue", None},
PlotRange->{1.04,1.76},ImageSize->320
];
Show[GraphicsRow[{Fig1,Fig2}]]

Example 4.4: Tidal Love numbers of Titan in function of crust thickness: choose another fluid-crust density (not in MW)

Similar to Example 4.3, except that the mantle is infinitely rigid.
Compare the results with two different fluid-crust densities: either ρ0=ρ or ρ0=ρbar

For safety, clear variables to be replaced:

Clear[R,Rm,g,ϵ,hn0,n]
Clear[ρ,ρb,δρ,δρhat]
Clear[μ,ν,χ]

As in Example 4.3, define (Subscript[k, n])and (Subscript[h, n]) functions for the three density models, assuming that ρ0=ρ (the mantle is infinitely rigid):

h20R=hRL0[(ρ/(ρb0[ρb,δρ,ϵ])),n]/.ruleTitan;
Clear[knLR,knMR,knDR,hnLR,hnMR,hnDR]
{knLR[ϵ_],knMR[ϵ_],knDR[ϵ_]}=kn[(μ/(ρ*g*R)),ν,χ,δρhat,ϵ,(ρ/(ρb0[ρb,δρ,ϵ])),hn0,n]/.ruleTitan/.ruleElastic/.{hn0->h20R}/.ruleLMD;
{hnLR[ϵ_],hnMR[ϵ_],hnDR[ϵ_]}=hn[(μ/(ρ*g*R)),ν,χ,δρhat,ϵ,(ρ/(ρb0[ρb,δρ,ϵ])),hn0,n]/.ruleTitan/.ruleElastic/.{hn0->h20R}/.ruleLMD;

Define new formulas for (Subscript[k, n])and (Subscript[h, n]), assuming that ρ0=ρbar (fluid-crust density=density of the physical crust):

Clear[knBis,hnBis];
knBis[μhat_,ν_,χ_,ϵ_,ξ_,hn0_,n_]:=hn0*(((1+((3ξ)/(2n+1))*(hn0*((Λ[μhat,ν,ϵ,n])/(1+Λ[μhat,ν,ϵ,n])))))^(-1))-1
hnBis[μhat_,ν_,χ_,δρhat_,ϵ_,ξ_,hn0_,n_]:=hn0*(((1+(1+((3ξ)/(2n+1))*hn0)*Λ[μhat,ν,ϵ,n]+Λχ[ν,χ,ϵ,n]+Hρ[δρhat,0,ν,ϵ,ξ,hn0,n]))^(-1))

Now, define (Subscript[k, n])and (Subscript[h, n]) functions for the three density models, assuming that ρ0=ρbar (the mantle is infinitely rigid):

h20RBis=hRLL0[1-ϵ,(ρ/(ρb)),((ρbar)/(ρb)),n]/.ruleTitan;
{knLRBis[ϵ_],knMRBis[ϵ_],knDRBis[ϵ_]}=knBis[(μ/(ρ*g*R)),ν,χ,ϵ,(ρ/(ρb)),hn0,n]/.ruleTitan/.ruleElastic/.{hn0->h20RBis}/.ruleLMD;
{hnLRBis[ϵ_],hnMRBis[ϵ_],hnDRBis[ϵ_]}=hnBis[(μ/(ρ*g*R)),ν,χ,δρhat,ϵ,(ρ/(ρb)),hn0,n]/.ruleTitan/.ruleElastic/.{hn0->h20RBis}/.ruleLMD;

Plot the Love numbers as functions of crust thickness, for the two different fluid-crust densities:

{ϵMin,ϵMax}=100*{0,0.101};
Fig1=Show[
Plot[knLR[ϵ/100],{ϵ,ϵMin,ϵMax},PlotStyle->Black],
Plot[knMR[ϵ/100],{ϵ,ϵMin,ϵMax},PlotStyle->Red],
Plot[knDR[ϵ/100],{ϵ,ϵMin,ϵMax},PlotStyle->Blue],
Plot[knLRBis[ϵ/100],{ϵ,ϵMin,ϵMax},PlotStyle->{Black,Dashed}],
Plot[knMRBis[ϵ/100],{ϵ,ϵMin,ϵMax},PlotStyle->{Red,Dashed}],
Plot[knDRBis[ϵ/100],{ϵ,ϵMin,ϵMax},PlotStyle->{Blue,Dashed}],
FrameLabel->{"Relative crust thickness [%]","(Subscript[k, 2])","(ρ^0)=ρ (solid); (ρ^0)=(OverBar[ρ]) (dashed)", None},
PlotRange->{0.32,0.721}
];
Fig2=Show[
Plot[hnLR[ϵ/100],{ϵ,ϵMin,ϵMax},PlotStyle->Black],
Plot[hnMR[ϵ/100],{ϵ,ϵMin,ϵMax},PlotStyle->Red],
Plot[hnDR[ϵ/100],{ϵ,ϵMin,ϵMax},PlotStyle->Blue],
Plot[hnLRBis[ϵ/100],{ϵ,ϵMin,ϵMax},PlotStyle->{Black,Dashed}],
Plot[hnMRBis[ϵ/100],{ϵ,ϵMin,ϵMax},PlotStyle->{Red,Dashed}],
Plot[hnDRBis[ϵ/100],{ϵ,ϵMin,ϵMax},PlotStyle->{Blue,Dashed}],
FrameLabel->{"Relative crust thickness [%]","(Subscript[h, 2])","(ρ^0)=ρ (solid); (ρ^0)=(OverBar[ρ]) (dashed)", None},
PlotRange->{1.04,1.76}
];
Show[GraphicsRow[{Fig1,Fig2}]]

Conclusion:
The results obtained with the different fluid-crust densities are equivalent to order ϵ.
Comparing this figures with Figure 5 of Paper MW shows that the value of (Subscript[k, 2]) computed with ρ0=ρbar (dashed curves) is closer to the benchmark, especially for the density model ‘Mixed’ (red curve). This was expected since the (Subscript[k, 2]) formula, if ρ0=ρbar, does not include any density corrections: the density structure of the whole body, crust included, is taken into account into the fluid-crust Love number.

Example 4.5: Load Love numbers of Titan in function of crust thickness for different crust/ocean density models (not in MW)

For safety, clear variables to be replaced:

Clear[Rm,R,g,ϵ,hn0,n]
Clear[ρ,ρb,δρ,δρhat]
Clear[μ,ν,χ,χμ,μEm]

Replacement rules (the crust is supposed to be elastic) and fluid-crust Love number (it is a function of (ρ,δρ,ϵ)):

h20=hVL0[(Rm/R),(ρ/(ρb0[ρb,δρ,ϵ])),((μEm)/(ρb0[ρb,δρ,ϵ]*g0[g,ρb,δρ,ϵ]*R)),n]/.ruleTitan;

Gravity and radial load Love numbers (degree 2) of Titan in function of crust thickness;
three density models (Light, Mixed, Dense):

Clear[knLoadL,knLoadM,knLoadD,hnLoadL,hnLoadM,hnLoadD]
{knLoadL[ϵ_],knLoadM[ϵ_],knLoadD[ϵ_]}=knLoad[(μ/(ρ*g*R)),ν,χ,δρhat,ϵ,(ρ/(ρb0[ρb,δρ,ϵ])),hn0,n]/.ruleTitan/.ruleElastic/.{hn0->h20}/.ruleLMD;
{hnLoadL[ϵ_],hnLoadM[ϵ_],hnLoadD[ϵ_]}=hnLoad[(μ/(ρ*g*R)),ν,χ,(ρ*g*R)*χμ,δρhat,ϵ,(ρ/(ρb)),(ρ/(ρb0[ρb,δρ,ϵ])),hn0,n]/.ruleTitan/.ruleElastic/.{χμ->χEice/μEice}/.{hn0->h20}/.ruleLMD;

Plot as in Example 4.3:

{ϵMin,ϵMax}=100*{0,0.101};
Fig1=Show[
Plot[knLoadL[ϵ/100],{ϵ,ϵMin,ϵMax},PlotStyle->{Black,Thick}],
Plot[knLoadM[ϵ/100],{ϵ,ϵMin,ϵMax},PlotStyle->{Red,Thick}],
Plot[knLoadD[ϵ/100],{ϵ,ϵMin,ϵMax},PlotStyle->{Blue,Thick}],
FrameLabel->{"Relative crust thickness [%]","(Subscript[k, 2])'","Light/Mixed/Dense = black/red/blue", None},
PlotRange->All,ImageSize->320
];
Fig2=Show[
Plot[hnLoadL[ϵ/100],{ϵ,ϵMin,ϵMax},PlotStyle->{Black,Thick}],
Plot[hnLoadM[ϵ/100],{ϵ,ϵMin,ϵMax},PlotStyle->{Red,Thick}],
Plot[hnLoadD[ϵ/100],{ϵ,ϵMin,ϵMax},PlotStyle->{Blue,Thick}],
FrameLabel->{"Relative crust thickness [%]","(Subscript[h, 2])'","Light/Mixed/Dense = black/red/blue", None},
PlotRange->All,ImageSize->320
];
Show[GraphicsRow[{Fig1,Fig2}]]

Example 4.6: Tilt factor of Europa as a function of crust thickness (MW-Fig 8)

For safety, clear variables to be replaced:

Clear[Rm,R,g,ϵ,hn0,n]
Clear[ρ,ρb,δρ,δρhat]
Clear[μ,ν,χ,μEm]

First, evaluate partially the fluid-crust Love number h20 with the replacement rule ‘ruleEuropa’;
don’t forget to use the fluid-crust bulk density and gravity;
h20 now depends on the three variables (ρ, δρ, ϵ):

h20=hVL0[(Rm/R),(ρ/(ρb0[ρb,δρ,ϵ])),((μEm)/(ρb0[ρb,δρ,ϵ]*g0[g,ρb,δρ,ϵ]*R)),n]/.ruleEuropa;

Next, evaluate partially the radial Love number of three rheology models (not strictly necessary for the tilt factor, but it is a better approximation as explained in Section 6.3 of MW).
The results (h2E,h2C,h2F) depend on the three variables (ρ, δρ, ϵ):

h2E=hn[(μ/(ρ*g*R)),ν,χ,((δρ)/ρ),ϵ,(ρ/(ρb0[ρb,δρ,ϵ])),hn0,n]/.{hn0->h20}/.ruleEuropa/.ruleEbar;
h2C=hn[(μ/(ρ*g*R)),ν,χ,((δρ)/ρ),ϵ,(ρ/(ρb0[ρb,δρ,ϵ])),hn0,n]/.{hn0->h20}/.ruleEuropa/.ruleCbar;
h2F=hn[(μ/(ρ*g*R)),ν,χ,((δρ)/ρ),ϵ,(ρ/(ρb0[ρb,δρ,ϵ])),hn0,n]/.{hn0->h20}/.ruleEuropa/.ruleFbar;

Then, express the tilt factor (degree 2) of Europa as a function of crust thickness;
Do it for the three rheology models (Elastic-like, Critical, Fluid-like) and for the three density models (Light, Mixed, Dense):

Clear[tiltEL,tiltEM,tiltED,tiltCL,tiltCM,tiltCD,tiltFL,tiltFM,tiltFD]
{tiltEL[ϵ_],tiltEM[ϵ_],tiltED[ϵ_]}=Abs[tilt[(μ/(ρ*g*R)),ν,χ,δρhat,0,ϵ,h2E,n]/.ruleEuropa/.ruleLMD/.ruleEbar];
{tiltCL[ϵ_],tiltCM[ϵ_],tiltCD[ϵ_]}=Abs[tilt[(μ/(ρ*g*R)),ν,χ,δρhat,0,ϵ,h2C,n]/.ruleEuropa/.ruleLMD/.ruleCbar];
{tiltFL[ϵ_],tiltFM[ϵ_],tiltFD[ϵ_]}=Abs[tilt[(μ/(ρ*g*R)),ν,χ,δρhat,0,ϵ,h2F,n]/.ruleEuropa/.ruleLMD/.ruleFbar];

Finally, plot the tilt factor as in MW-Fig 8 (Light/Mixed/Dense = black/red/blue):

{ϵMin,ϵMax}=100*{0,0.1025};
{yMin,yMax}={-0.015,0.365};
Fig1=Show[
Plot[tiltEL[ϵ/100],{ϵ,ϵMin,ϵMax},PlotStyle->{Black,Thick}],
Plot[tiltEM[ϵ/100],{ϵ,ϵMin,ϵMax},PlotStyle->{Red,Thick}],
Plot[tiltED[ϵ/100],{ϵ,ϵMin,ϵMax},PlotStyle->{Blue,Thick}],
FrameLabel->{"Relative crust thickness [%]","Abs (Subscript[γ, 2])","Elastic-like", None},
PlotRange->{yMin,yMax}
];
Fig2=Show[
Plot[tiltCL[ϵ/100],{ϵ,ϵMin,ϵMax},PlotStyle->{Black,Thick}],
Plot[tiltCM[ϵ/100],{ϵ,ϵMin,ϵMax},PlotStyle->{Red,Thick}],
Plot[tiltCD[ϵ/100],{ϵ,ϵMin,ϵMax},PlotStyle->{Blue,Thick}],
FrameLabel->{"Relative crust thickness [%]","Abs (Subscript[γ, 2])","Critical",None},
PlotRange->{yMin,yMax}
];
Fig3=Show[
Plot[tiltFL[ϵ/100],{ϵ,ϵMin,ϵMax},PlotStyle->{Black,Thick}],
Plot[tiltFM[ϵ/100],{ϵ,ϵMin,ϵMax},PlotStyle->{Red,Thick}],
Plot[tiltFD[ϵ/100],{ϵ,ϵMin,ϵMax},PlotStyle->{Blue,Thick}],
FrameLabel->{"Relative crust thickness [%]","Abs (Subscript[γ, 2])","Fluid-like",None},
PlotRange->{yMin,yMax}
];
Show[GraphicsRow[{Fig1,Fig2,Fig3}]]

5. Dynamical Love numbers

Example 5.1: Dynamical tilt factor of Europa as a function of ocean thickness (MW-Fig 14)

For safety, clear variables to be replaced:

Clear[Rm,R,g,n,qω,z,ϵ,ϵO]
Clear[ρ,ρb,δρhat]
Clear[μ,ν,χ,μEm]

First, evaluate partially the fluid-crust Love number h20 with the replacement rule ‘ruleEuropa’;
It is sufficient here to approximate it by its value at zero crust thickness; the fluid-crust bulk density and gravity are thus the same as in the physical model.
h20 now depends on the variable ρ:

h20=hVL0[(Rm/R),(ρ/(ρb)),((μEm)/(ρb*g*R)),n]/.ruleEuropa;

The dynamical correction Λω depends on z, which is defined by z=(Subscript[R, m])/(Subscript[R, ϵ]) and can be written in terms of ϵ and (Subscript[ϵ, O])=(Subscript[d, O])/R, where (Subscript[d, O]) is the ocean thickness:

zϵ=FullSimplify[(Rm/(Rϵ))/.{Rm->R-d-dO,Rϵ->R-d}/.{d->ϵ*R,dO->ϵO*R}]

Express the dynamical tilt factor as a function of the relative crust thickness and of the relative ocean thickness.
Since the crust is elastic, use the replacement rule ‘ruleElastic’.
For the Light density model, use the replacement rule ‘ruleL’.

Clear[tiltDyn];
tiltDyn[ϵ_,ϵO_]=(tilt[(μ/(ρ*g*R)),ν,χ,δρhat,Λωn[qω,z,n],ϵ,h20,n]*(1/ϵ))/.ruleElastic/.ruleEuropa/.{qω->qωEuropa,z->zϵ}/.ruleL;

Assume that the crust thickness is either d=10 km or d=50 km:

{ϵ10,ϵ50}={10.,50.}*10^3/RadiusEuropa;

Plot the dynamical tilt factor of Europa as a function of ocean thickness in units of d/R (MW-Fig 14):

R0=RadiusEuropa/1000;
{Dmin,Dmax}={0,120};
{xMin,xMax}={-3,123};
{yMin,yMax}={-0.4,4.8};
Show[
Plot[tiltDyn[ϵ10,dO/R0],{dO,Dmin,Dmax},PlotStyle->{Blue,Thick}],
Plot[tiltDyn[ϵ50,dO/R0],{dO,Dmin,Dmax},PlotStyle->{Red,Thick}],
FrameLabel->{"Ocean thickness [km]","Tilt factor [in units of d/R]","d=10 km (blue), d=50 km (red)",None},
PlotRange->{{xMin,xMax},{yMin,yMax}}, AspectRatio->0.75,ImageSize->350
]

Example 5.2: Dynamical resonance (MW-Fig 15)

For safety, clear variables to be replaced:

Clear[ξ0,δρhat,ϵ,qω,z,n,ρ,g,R]
Clear[μE,νE,χE,f,δ1,δ2]

The variable z can be replaced by zϵ as in Example 5.1:

zϵ=1+((ϵO)/(-1+ϵ));

Crust and ocean have the same density (1000 kg/m3) so that ρ=ρWater and δρhat=0.
The crust is incompressible so that (Subscript[ν, E])=0.5 and (Subscript[χ, E])=0.
For the critical model, use the replacement rule ‘ruleC’.

Clear[k2Dyn,h2Dyn];
k2Dyn[ϵ_,ϵO_]=knrDyn[((μbar[μE,fϵ,δ1,δ2])/(ρ*g*R)),νE,χE,δρhat,ϵ,qω,z,ξ0,n]/.{ξ0->ρ/ρb,νE->0.5,χE->0,δρhat->0}/.{ρ->ρWater}/.ruleC/.ruleEuropa/.{qω->qωEuropa,z->zϵ};
h2Dyn[ϵ_,ϵO_]=hnrDyn[((μbar[μE,fϵ,δ1,δ2])/(ρ*g*R)),νE,χE,δρhat,ϵ,qω,z,ξ0,n]/.{ξ0->ρ/ρb,νE->0.5,χE->0,δρhat->0}/.{ρ->ρWater}/.ruleC/.ruleEuropa/.{qω->qωEuropa,z->zϵ};

Assume that the crust thickness is 10 km:

ϵ10=10.*10^3/RadiusEuropa;

Plot the dynamical resonance enhancing the Love numbers of Europa as a function of ocean thickness (MW-Fig 15);
assume that the crust thickness is 10 km:

R0=RadiusEuropa/1000;
Show[
LogLinearPlot[Abs@h2Dyn[ϵ10,dO/R0],{dO,0.01,100.},PlotStyle->{Blue,Thick}],
LogLinearPlot[Abs@k2Dyn[ϵ10,dO/R0],{dO,0.01,100.},PlotStyle->{Red,Thick}],
LogLinearPlot[-Im@k2Dyn[ϵ10,dO/R0],{dO,0.01,100.},PlotStyle->{Red,Dotted,Thick}],
FrameLabel->{"Ocean thickness [km]","Love numbers","|(Subscript[h, 2])| in blue; |(Subscript[k, 2])| in red; -Im(\!\(\*SubscriptBox[\(k\), \(\(2\)\(\\ \)\)]\)) in red & dotted",None},
PlotRange->{-0.15,1.6},AspectRatio->0.75,ImageSize->350]

Evaluate the resonant thickness (in meters) for d=10 km and d=50 km:

Clear[μE,fϵ,δ1,δ2,ϵ,n,ρ,g,R]
{ϵ10,ϵ50}={10.*^3,50.*^3}/RadiusEuropa;
Dresonant[R,qω,ξ,((μbar[μE,fϵ,δ1,δ2])/(ρ*g*R)),νbar[νE,fϵ,δ1,δ2],ϵ,n]/.{ξ->(ρ/(ρb))}/.{ρ->ρWater}/.ruleEuropa/.{qω->qωEuropa}/.ruleE/.{{ϵ->ϵ10},{ϵ->ϵ50}}

6. A simple thick shell model: the homogeneous crust model (static and dynamical)

Example 6.1: Dependence of geometrical factors on crust thickness (MP-Fig 13; MW-Table 11)

This model is useful to estimate deviations from the membrane formulas due to the finite thickness of the crust.
In this model, the (Subscript[l, 2])-(Subscript[h, 2]) and (Subscript[k, 2])-(Subscript[h, 2]) relations depend on the factors (Subscript[z, h]) and (Subscript[z, l]):
(Subscript[k, 2])+1=(1+(Subscript[z, h])(OverHat[μ])) (Subscript[h, 2])
(Subscript[l, 2])=(Subscript[z, l])(Subscript[h, 2])
The factors (Subscript[z, h]) and (Subscript[z, l]) depend on (OverHat[μ]) and (Subscript[Λ, ω]), and on five geometrical factors ((Subscript[X, a]),(Subscript[X, b]),(Subscript[X, c]),(Subscript[X, d]),(Subscript[X, e])).

The geometrical factors are given by the following functions of x (x=(R-d)/R):

Xlist={Xa[x],Xb[x],Xc[x],Xd[x],Xe[x]}

The numerators and denominators can be factorized but the resulting expressions are more complicated:

Factor[Xlist]

In the thin shell limit, the geometrical factors become (MW-Table 11):

Series[Xlist/.{x->1-ϵ},{ϵ,0,2}]

In the static limit, the membrane approximations of (Subscript[z, h]) and (Subscript[z, l]) are, at leading order,

zhLin[x_]:=(24/11)(1-x);
zlMemb=(3/11);

In the homogeneous body limit, the geometrical factors become (MW-Table 11):

Xlist/.{x->0}

Plot the factors (Subscript[z, h]) and (Subscript[z, l]) if Λω=0;
plot also the membrane approximations given by MP-115 (MP-Fig 13A):

Fig1=Show[
Plot[Xa[x],{x,0,1}],
Plot[zhLin[x],{x,0,1},PlotStyle->{Black,Dotted,Thick}],
FrameLabel->{"x","(Subscript[z, h])  (if (Subscript[Λ, ω])=0)","(Subscript[z, h]) (solid), membrane approx. (dotted)",None},PlotRange->{-0.2,4.2},AspectRatio->0.75,ImageSize->370];
Fig2=Show[
Plot[Xc[x],{x,0,1}],
Graphics[{PointSize[0.025],Black,Point[{1.,zlMemb}]}],
FrameLabel->{"x","(Subscript[z, l])  (if (Subscript[Λ, ω])=0)","(Subscript[z, l]) (solid), membrane approx. (dot)",None},PlotRange->{0.195,0.305},AspectRatio->0.75,ImageSize->370];
Show[GraphicsRow[{Fig1,Fig2}]]

Zoom on the above plots, in the range 0.8 < x < 1:

Fig1=Show[
Plot[Xa[x],{x,0.8,1},FrameLabel->{"x","(Subscript[z, h])  (if (Subscript[Λ, ω])=0)","(Subscript[z, h]) (solid), membrane approx. (dotted)",None}],
Plot[zhLin[x],{x,0.8,1},PlotStyle->{Black,Dotted,Thick}],
PlotRange->{-0.02,0.62},AspectRatio->0.75,ImageSize->370];
Fig2=Show[
Plot[Xc[x],{x,0.8,1},FrameLabel->{"x","(Subscript[z, l])  (if (Subscript[Λ, ω])=0)","(Subscript[z, l]) (solid), membrane approx. (dot)",None}],
Graphics[{PointSize[0.025],Black,Point[{1.,zlMemb}]}],
PlotRange->{0.218,0.282},AspectRatio->0.75,ImageSize->370];
Show[GraphicsRow[{Fig1,Fig2}]]