In the fluid limit, the shell above the ocean becomes fluid, so that the
relevant variables within the shell are now $y_{5}^{J}$ and $y_{7}^{J}$. These
variables depend on one unknown constant and satisfy one boundary condition at
the surface:

$y_{7}^{J}(R_{s})=0\,.$ (B.4)

This boundary condition is obtained by substituting either Eq. (A.4) or Eq.
(A.5) into Eq. (B.3). The unknown constant is thus equal to zero, so that the
fluid variables vanish within the ocean and fluid shell:

$y_{5}^{J}(r)=y_{7}^{J}(r)=0\,.$ (B.5)

Physically, this means that the fluid cannot support surface and internal
loads: from the gravitational point of view, the loads are completely
compensated. Although $y_{6}^{I}$ is discontinuous at the shell-ocean boundary
(Eq. (A.7)), $y_{7}^{I}$ is defined so as to be continuous and vanishes
through the fluid shell and the ocean.

The condition that $y_{5}^{J}=0$ implies with Eq. (A.9) that the gravitational
fluid Love numbers cancel the direct effect of the load:

$\left(k_{s}^{L\circ},k_{o}^{L\circ},k_{s}^{I\circ},k_{o}^{I\circ}\right)=\left(-1,-x^{n},-x^{n+1},-1\right).$
(B.6)

At the surface, the fluid constraint combined with the condition
$y_{5}^{J}(r)=0$ yields

$\rho_{s}g_{s}\,y_{1}^{J}(R_{s})=y_{2}^{J}(R_{s})\,,$ (B.7)

from which the radial fluid Love numbers can be computed by application of the
surface boundary conditions (Eqs. (A.4)-(A.5)):

$\displaystyle h_{s}^{L\circ}$ $\displaystyle=$
$\displaystyle-\frac{2n+1}{3\,\xi_{s1}}\,,$ (B.8) $\displaystyle
h_{s}^{I\circ}$ $\displaystyle=$ $\displaystyle 0\,.$ (B.9)

Eq. (B.8) can be found in Wu and Peltier [1982] for a homogeneous model, and
with a restriction to high harmonic degrees for an inhomogeneous model; its
general validity was stated without proof by Greff-Lefftz et al. [2010].

At the shell-ocean boundary, the fluid constraint combined with the condition
$y_{5}^{J}(r)=0$ yields

$\displaystyle\rho_{s}g_{s}\,y_{1}^{J}(R_{o}^{+})=y_{2}^{J}(R_{o}^{+})\,,$
(B.10)
$\displaystyle\rho_{o}g_{s}\,y_{1}^{J}(R_{o}^{-})=y_{2}^{J}(R_{o}^{-})\,.$
(B.11)

For a surface load, these equations are only consistent with the continuity of
$y_{1}^{J}$ and $y_{2}^{J}$ if the two variables vanish. In other words, the
shell and ocean are unperturbed below the compensated surface load. Thus, the
radial fluid Love number vanishes at the shell-ocean boundary:

$h_{o}^{L\circ}=0\,.$ (B.12)

For an internal load, Eqs. (B.10)-(B.11) combined with the continuity of
$y_{1}^{I}$ and the discontinuity of $y_{2}^{I}$ (Eq. (A.6)) yield

$h_{o}^{I\circ}=-\frac{2n+1}{3\left(\xi_{o1}-\xi_{s1}\right)}\,\frac{R_{s}}{R_{o}}\,.$
(B.13)

##  Appendix C Gravity-deformation relations in a 3-layer body

Consider an incompressible body stratified in layers of homogeneous density
(the rheology may be depth-dependent). In that case, the gravitational
perturbation at an arbitrary radius can be computed as a function of the
radial deformations of interfaces with a density contrast. For a 3-layer body
(shell/ocean/core), Eqs. (7)-(8) yield

$\displaystyle\Gamma_{sn}$ $\displaystyle=$ $\displaystyle
g_{s}\left(\xi_{sn}\,H_{sn}+\Delta\xi_{n}\,x^{n+2}\,H_{on}+\Delta\xi^{\prime}_{n}\,y^{n+2}\,H_{cn}\right),$
(C.1) $\displaystyle\Gamma_{on}$ $\displaystyle=$ $\displaystyle
g_{s}\left(\xi_{sn}\,x^{n}\,H_{sn}+\Delta\xi_{n}\,x\,H_{on}+\Delta\xi^{\prime}_{n}\,x\,(y/x)^{n+2}\,H_{cn}\right),$
(C.2) $\displaystyle\Gamma_{cn}$ $\displaystyle=$ $\displaystyle
g_{s}\left(\xi_{sn}\,y^{n}\,H_{sn}+\Delta\xi_{n}\,x\,(y/x)^{n}\,H_{on}+\Delta\xi^{\prime}_{n}\,y\,H_{cn}\right).$
(C.3)

Using Eqs. (14)-(16) and Eqs. (17)-(19), I can write these relations in terms
of Love numbers (even better, in terms of deviatoric Love numbers):

$\displaystyle\hat{k}_{s}^{J}$ $\displaystyle=$
$\displaystyle\xi_{sn}\,\hat{h}_{s}^{J}+\Delta\xi_{n}\,x^{n+2}\,\hat{h}_{o}^{J}+\Delta\xi^{\prime}_{n}\,y^{n+2}\,\hat{h}_{c}^{J}\,,$
(C.4) $\displaystyle\hat{k}_{o}^{J}$ $\displaystyle=$
$\displaystyle\xi_{sn}\,x^{n}\,\hat{h}_{s}^{J}+\Delta\xi_{n}\,x\,\hat{h}_{o}^{J}+\Delta\xi^{\prime}_{n}\,x\,(y/x)^{n+2}\,\hat{h}_{c}^{J}\,,$
(C.5) $\displaystyle\hat{k}_{c}^{J}$ $\displaystyle=$
$\displaystyle\xi_{sn}\,y^{n}\,\hat{h}_{s}^{J}+\Delta\xi_{n}\,x\,(y/x)^{n}\,\hat{h}_{o}^{J}+\Delta\xi^{\prime}_{n}\,y\,\hat{h}_{c}^{J}\,.$
(C.6)

If the core is infinitely rigid and $n\geq 2$, $\hat{h}_{c}^{J}=0$ so that
$(\hat{k}_{s}^{J},\hat{k}_{o}^{J})$ do not depend on the properties of the
core. If the core is fluid, the condition that the core surface is an
equipotential implies that the shape of the core follows the local geoid:

$H_{cn}=\frac{\Gamma_{cn}}{g_{c}}\,\,\Leftrightarrow\,\,\hat{k}_{c}^{J}=\gamma_{c}\,\hat{h}_{c}^{J}\,,$
(C.7)

where $\gamma_{c}=g_{c}/g_{s}$ (see Eq. (5)).

If the core is elastic and incompressible, the relation between $k_{c}^{J}$
and $h_{c}^{J}$ can be derived with the propagator matrix approach. The 6 non-
singular radial functions $y_{1...6}(r)$ depending on 3 unknowns are
propagated from the center to the core-ocean boundary, where two constraints
are applied: the free-slip condition, $y_{4}(R_{c})=0$, and the fluid
constraint, $y_{2}(R_{c})=\rho_{o}(g_{o}y_{1}(R_{c})-y_{5}(R_{c}))$ [Saito,
1974]. These equations can be solved to yield a relation between the
gravitational perturbation $\hat{k}_{c}^{J}$ and the radial deformation
$\hat{h}_{c}^{J}$ of the core-ocean boundary. If the core is homogeneous with
shear modulus $\mu_{\rm c}$, the resulting relation can be written as

$\hat{k}_{c}^{J}=\left(\gamma_{c}+f_{n}\,\frac{\bar{\mu}_{\rm
c}}{y\,\Delta\xi^{\prime}_{1}}\right)\hat{h}_{c}^{J}\,,$ (C.8)

where $\bar{\mu}_{\rm c}=\mu_{\rm c}/(\rho_{b}g_{s}R_{s})$ is the reduced
shear modulus of the core and $\Delta\xi^{\prime}_{1}$ is the core-ocean
density ratio contrast given by Eq. (6). The function $f_{n}$ is defined by

$f_{n}=\frac{2\left(n-1\right)\left(2n^{2}+4n+3\right)}{n\left(2n+1\right)}\,.$
(C.9)

At degree one, $f_{1}=0$ so that Eq. (C.8) reduces to Eq. (C.7): the core is
rigidly translated and its surface is equipotential.

Solving Eq. (C.6) and Eq. (C.8) for $\hat{h}_{c}^{J}$, I get

$\hat{h}_{c}^{J}=\frac{\xi_{sn}\,y^{n}\,\hat{h}_{s}^{J}+\Delta\xi_{n}\,x\,(y/x)^{n}\,\hat{h}_{o}^{J}}{\gamma_{c}+f_{n}\,\bar{\mu}_{\rm
c}/(y\,\Delta\xi^{\prime}_{1})-y\,\Delta\xi^{\prime}_{n}}\,.$ (C.10)

Substituting the result into Eqs. (C.4)-(C.5), I write the gravitational
perturbations in terms of the radial deformations at the shell boundaries:

$\displaystyle\hat{k}_{s}^{J}$ $\displaystyle=$
$\displaystyle\xi_{sn}\left(1+K_{n}\right)\hat{h}_{s}^{J}+\Delta\xi_{n}\,x^{n+2}\left(1+K_{n}\,x^{-2n-1}\right)\hat{h}_{o}^{J}\,,$
(C.11) $\displaystyle\hat{k}_{o}^{J}$ $\displaystyle=$
$\displaystyle\left(1+K_{n}\,x^{-2n-1}\right)\left(\xi_{sn}\,x^{n}\,\hat{h}_{s}^{J}+\Delta\xi_{n}\,x\,\hat{h}_{o}^{J}\right),$
(C.12)

where

$K_{n}=\frac{\Delta\xi^{\prime}_{n}\,y^{2n+2}}{\gamma_{c}+f_{n}\,\bar{\mu}_{\rm
c}/(y\,\Delta\xi^{\prime}_{1})-y\,\Delta\xi^{\prime}_{n}}\,.$ (C.13)

In the main text, these equations are written more compactly as Eq. (46) with
the coefficients of Table 3. Instead of repeating this equation here, I give
the corresponding relations for non-deviatoric Love numbers:

$\left(\begin{array}[]{ll}k_{s}^{L}&k_{s}^{I}\\\
k_{o}^{L}&k_{o}^{I}\end{array}\right)=\left(\begin{array}[]{ll}a&b\\\
c&b\,x^{-n-1}\end{array}\right)\left(\begin{array}[]{ll}h_{s}^{L}&h_{s}^{I}\\\
h_{o}^{L}&h_{o}^{I}\end{array}\right)+K_{n}\left(\begin{array}[]{ll}1&x^{-n}\\\
x^{-n-1}&x^{-2n-1}\end{array}\right).$ (C.14)

In the limit $\bar{\mu}_{\rm s}\rightarrow\infty$, the shell becomes
infinitely rigid so that its deformation tends to zero: $h_{j}^{J}\rightarrow
0$. Eq. (C.14) then implies that $k_{s}^{J}\rightarrow K_{n}$. In words,
$K_{n}$ measures the gravitational perturbation (due to deformed layers below
the shell) in the limit of a non-deformable shell.

In the limit $\bar{\mu}_{\rm c}\rightarrow\infty$, the core becomes infinitely
rigid and $K_{n}\rightarrow 0$. The rigid core constraint, $K_{n}=0$, only
make sense if $n\geq 2$. At degree one, $f_{1}=0$ so that Eq. (C.13) yields,
whatever the rheology of the core,

$K_{1}=\frac{\Delta\xi^{\prime}_{1}}{\xi_{o1}}\,y^{3}=\frac{1-\xi_{s1}-x^{3}\Delta\xi_{1}}{\xi_{o1}}\,.$
(C.15)

This value does not depend on the shear modulus of the core because the core
is rigidly translated.

##  Appendix D Love numbers if homogeneous shell

Consider an incompressible 3-layer body with a homogeneous shell (the ocean
and the core are not necessarily homogeneous). This problem can be solved
analytically with the propagator matrix method (see Appendix A.4). Since the
shell is homogeneous, the shell propagator matrix and its inverse appear each
once in the elastic-gravitational equations. The Love numbers can thus be
written as the ratio of two degree-2 polynomials in the shear modulus of the
shell.

###  D.1 Radial Love numbers

Radial Love numbers can thus be written in the following generic form:

$h_{j}^{J}=\frac{A_{j}^{J}+B_{j}^{J}\mu_{\rm s}+C_{j}^{J}\mu_{\rm
s}^{2}}{D+E\mu_{\rm s}+F\mu_{\rm s}^{2}}\,,$ (D.1)

where $j=(s,o)$ and $J=(L,I)$. The coefficients $(D,E,F)$ are common to all
Love numbers (radial, tangential, or gravitational), whatever the type of load
(applied at the surface or at the shell-ocean boundary) or the location where
they are evaluated (surface or shell-ocean boundary). The degree of the
polynomial in the denominator is related to the number of relaxation modes of
a viscoelastic shell (see companion paper). All coefficients except
$A_{j}^{J}$ are polynomials in the interior parameters
$(x,y,\xi_{s1},\xi_{o1},\bar{\mu}_{\rm c})$.

The fluid limit ($h_{j}^{J}\rightarrow h_{j}^{J\circ}$, see Table 2) and the
limit of an infinitely rigid shell ($h_{j}^{J}\rightarrow 0$) imply that

$\displaystyle A^{J}_{j}$ $\displaystyle=$ $\displaystyle
D\,h^{J\circ}_{j}\,,$ (D.2) $\displaystyle C_{j}^{J}$ $\displaystyle=$
$\displaystyle 0\,.$ (D.3)

Two of the identities required for the $\mu$-invariance of zero-deflection
isostasy, Eqs. (94)-(95), imply that

$\displaystyle A_{s}^{L}\,B_{o}^{I}+A_{o}^{I}\,B_{s}^{L}$ $\displaystyle=$
$\displaystyle(E/D)\,A_{s}^{L}\,A_{o}^{I}\,,$ (D.4) $\displaystyle
B_{s}^{L}\,B_{o}^{I}-B_{o}^{L}\,B_{s}^{I}$ $\displaystyle=$
$\displaystyle(F/D)\,A_{s}^{L}\,A_{o}^{I}\,,$ (D.5)

while Eq. (93) is automatically satisfied because
$h_{o}^{L}/h_{s}^{I}=B_{o}^{L}/B_{s}^{I}$ does not depend on $\mu_{\rm s}$.

Deviatoric radial Love numbers read

$\hat{h}_{j}^{J}=\frac{\hat{A}_{j}^{J}+\hat{B}_{j}^{J}\mu_{\rm
s}+\hat{C}_{j}^{J}\mu_{\rm s}^{2}}{D+E\mu_{\rm s}+F\mu_{\rm s}^{2}}\,,$ (D.6)

where

$\displaystyle\hat{A}_{j}^{J}$ $\displaystyle=$ $\displaystyle 0\,,$ (D.7)
$\displaystyle\hat{B}_{j}^{J}$ $\displaystyle=$ $\displaystyle
B_{j}^{J}-E\,(A_{j}^{J}/D)\,,$ (D.8) $\displaystyle\hat{C}_{j}^{J}$
$\displaystyle=$ $\displaystyle C_{j}^{J}-F\,(A_{j}^{J}/D)\,.$ (D.9)

Using Eqs. (D.2) and (D.4), I can write the coefficients $\hat{B}_{j}^{J}$ as

$\left(\hat{B}^{L}_{s},\hat{B}^{L}_{o},\hat{B}^{I}_{s},\hat{B}^{I}_{o}\right)=\left(-\frac{h_{s}^{L\circ}}{h_{o}^{I\circ}}\,B_{o}^{I},B^{L}_{o},B^{I}_{s},-\frac{h_{o}^{I\circ}}{h_{s}^{L\circ}}\,B_{s}^{L}\right).$
(D.10)

Partial derivatives of radial Love numbers (deviatoric or not, see Eq. (30))
with respect to $\mu_{\rm s}$ and evaluated at $\mu_{\rm s}=0$ are equal to

$\dot{h}^{J}_{j}=\hat{B}_{j}^{J}/D\,.$ (D.11)

###  D.2 Gravitational Love numbers

Gravitational Love numbers are written in the same generic form as the radial
Love numbers, except that the coefficients are denoted by a ‘prime’
superscript:

$k_{j}^{J}=\frac{A^{\prime J}_{j}+B^{\prime J}_{j}\mu_{\rm s}+C^{\prime
J}_{j}\mu_{\rm s}^{2}}{D+E\mu_{\rm s}+F\mu_{\rm s}^{2}}\,.$ (D.12)

The fluid limit (see Table 2) and the limit of an infinitely rigid shell (see
after Eq. (C.14)) imply that

$\displaystyle A^{\prime J}_{j}$ $\displaystyle=$ $\displaystyle
D\,k^{J\circ}_{j}\,,$ (D.13) $\displaystyle\left(C^{\prime L}_{s},C^{\prime
L}_{o},C^{\prime I}_{s},C^{\prime I}_{o}\right)$ $\displaystyle=$
$\displaystyle\left(1,x^{-n-1},x^{-n},x^{-2n-1}\right)K_{n}\,F\,.$ (D.14)

Deviatoric gravitational Love numbers read

$\hat{k}_{j}^{J}=\frac{\hat{A}^{\prime J}_{j}+\hat{B}^{\prime J}_{j}\mu_{\rm
s}+\hat{C}^{\prime J}_{j}\mu_{\rm s}^{2}}{D+E\mu_{\rm s}+F\mu_{\rm s}^{2}}\,,$
(D.15)

where $\hat{A}^{\prime J}_{j}=0$. The coefficients $(\hat{B}^{\prime
J}_{j},\hat{C}^{\prime J}_{j})$ are related to $(B^{\prime J}_{j},C^{\prime
J}_{j})$ by equations similar to Eqs. (D.8)-(D.9) and to
$(\hat{B}^{J}_{j},\hat{C}^{J}_{j})$ by Eq. (46).

The Saito-Molodensky relation (Eq. (20)) gives the following (equivalent)
constraints:

$\displaystyle B^{I}_{s}-B^{\prime I}_{s}$ $\displaystyle=$ $\displaystyle
x\left(\gamma_{o}\,B^{L}_{o}-B^{\prime L}_{o}\right)\,,$ (D.16)
$\displaystyle\hat{B}^{I}_{s}-\hat{B}^{\prime I}_{s}$ $\displaystyle=$
$\displaystyle x\left(\gamma_{o}\,\hat{B}^{L}_{o}-\hat{B}^{\prime
L}_{o}\right)\,.$ (D.17)

###  D.3 Rigid core, fluid core, and thin shell limits

If the core is elastic, the analytical expressions for the Love number
coefficients are quite long. More manageable formulas are obtained by relating
the elastic core model to simpler models in which the core is infinitely
rigid, fluid-like, or point-like. Since the core is homogeneous and extends to
the centre of the body, the propagator matrix for the core is used only once
and its inverse is not needed. As a result, the numerator and denominator of
the Love number coefficients depend linearly on the shear modulus of the core.
Denoting generically the coefficients ($B_{j}^{J}$ and so forth) of the
elastic core model by the symbol $Z_{e}$, I can thus write

$Z_{e}=Z_{f}+p\,\bar{\mu}_{\rm c}\,Z_{r}\,,$ (D.18)

where $Z_{f}$ and $Z_{r}$ are the values of the coefficient $Z_{e}$ in the
fluid core limit ($\bar{\mu}_{\rm c}\rightarrow 0$) and in the rigid core
limit ($\bar{\mu}_{\rm c}\rightarrow\infty$), respectively. The factor of
proportionality $p$ is common to all coefficients, so that it simplifies in
the Love numbers when taking the rigid core limit. Moreover, the fluid core
model should reduce to the rigid core model in the point-core limit
($y\rightarrow 0$) because the gravity field outside a rigid core is
indistinguishable from the one produced by a point-like core. Thus, I can
write

$Z_{f}=q\,Z_{r}+Z_{f0}\,,$ (D.19)

where $q$ is a factor of proportionality common to all coefficients. $Z_{f0}$
is the remainder function which vanishes in the point-core limit. Solving
explicitly the elastic core model, I get

$\displaystyle p$ $\displaystyle=$ $\displaystyle
2\left(n-1\right)\left(2n^{2}+4n+3\right)x^{2n+1}\,y^{4}\,,$ (D.20)
$\displaystyle q$ $\displaystyle=$ $\displaystyle
2n\left(n-1\right)x^{2n+1}(1-\xi_{s1}-x^{3}\Delta\xi_{1})^{2}\,.$ (D.21)

In this way, the lengthy coefficients $Z_{e}$ can be expressed in terms of the
shorter coefficients $Z_{r}$ and $\tilde{Z}_{f}$.

In the thin shell limit, the limiting values of Love numbers (Eqs. (64)-(65))
imply that the generic coefficients tend to

$\lim_{\varepsilon\rightarrow 0}\frac{(D,F,B^{J}_{j},B^{\prime
J}_{j})}{E}=\left(0,0,-\frac{1}{\xi_{on}},-1\right).$ (D.22)

The identity (D.5) then entails that the following ratios of $B_{j}^{J}$
coefficients are equal up to order $\varepsilon$ included:

$B_{s}^{L}/B_{o}^{L}=B_{s}^{I}/B_{o}^{I}+{\cal O}(\varepsilon^{2})\,.$ (D.23)

Combined with Eq. (101), this result shows that elastic isostasy in a thin
incompressible homogeneous shell does not depend on the choice of boundary
conditions (represented by the isostatic family parameter), so that ZDI and
MSI prescriptions have the same thin shell limit.

##  Appendix E Integral identities and variation of Love numbers

###  E.1 Integral of motion

Let $y^{J}_{i}$ and $y^{K}_{i}$ ($i=1...6$) be two solutions of the equations
of motion, with $(J,K$) denoting the type of load (and the the boundary
conditions). Consider the following bilinear form:

$Y^{JK}=r^{2}\left(y^{J}_{1}\,y^{K}_{2}-\delta_{n}\,y^{J}_{3}\,y^{K}_{4}+\frac{y^{J}_{5}\,y^{K}_{6}}{4\pi{}G}\right),$
(E.1)

where $\delta_{n}=-n(n+1)$. In the solid layers (shell and core), the
following quantities are integrals of motion [Takeuchi and Saito, 1972; Okubo
and Saito, 1983; Tobie et al., 2005]:

$\displaystyle\left[Y^{JK}\right]_{0}^{R_{c}}$ $\displaystyle=$
$\displaystyle\int_{0}^{R_{c}}{\cal L}^{JK}_{s}\,dr\,,$ (E.2)
$\displaystyle\left[Y^{JK}\right]_{R_{o}}^{R_{s}}$ $\displaystyle=$
$\displaystyle\int_{R_{o}}^{R_{s}}{\cal L}^{JK}_{s}\,dr\,,$ (E.3)

where

${\cal L}^{JK}_{s}=\mu\,{\cal H}^{JK}_{\mu}+\kappa\,{\cal
H}^{JK}_{\kappa}+\rho\,{\cal H}^{JK}_{\rho}+{\cal H}^{JK}_{0}\,,$ (E.4)

in which $(\mu,\kappa,\rho)$ are the depth-dependent values of the shear
modulus, bulk modulus, and density. The functions ${\cal H}^{JK}_{\mu}$ and
${\cal H}^{JK}_{\kappa}$, defined by Eqs. (116)-(117), were independently
obtained by integrating the strain energy density over the volume (see Eqs.
(114)-(115)): $(\mu\,{\cal H}^{JK}_{\mu}+\kappa\,{\cal H}^{JK}_{\kappa})/2$ is
actually the angular average of the strain energy density. For the sake of
completeness, the other functions are given by

$\displaystyle{\cal H}^{JK}_{\rho}$ $\displaystyle=$
$\displaystyle(n+1)r\left(y^{J}_{1}\,y^{K}_{5}+y^{J}_{5}\,y^{K}_{1}\right)+r\delta_{n}\left(y^{J}_{3}\,y^{K}_{5}+y^{J}_{5}\,y^{K}_{3}\right)$
(E.6)
$\displaystyle-\,4g_{r}r\,y^{J}_{1}\,y^{K}_{1}-g_{r}r\,\delta_{n}\left(y^{J}_{1}\,y^{K}_{3}+y^{J}_{1}\,y^{K}_{3}\right),$
$\displaystyle{\cal H}^{JK}_{0}$ $\displaystyle=$
$\displaystyle\frac{1}{4\pi{}G}\,r^{2}\,y^{J}_{6}\,y^{K}_{6}\,.$ (E.7)

In the fluid layer (i.e. the ocean), the function $y_{3}$ (tangential
displacement potential) is undetermined in the static limit. Eq. (E.4) must
then be replaced by

$\left[Y^{JK}\right]_{R_{c}}^{R_{o}}=\int_{R_{c}}^{R_{o}}{\cal
L}^{JK}_{f}\,dr+\left[\frac{r^{2}}{\rho_{o}g_{r}}\,y^{J}_{2}\,y^{K}_{2}\right]_{R_{c}}^{R_{o}},$
(E.8)

where $y^{J}_{4}=y^{K}_{4}=0$ in the LHS (no shear stress in a fluid layer),
while the integrand in the RHS is given by

${\cal L}^{JK}_{f}=2(n-1)\,\frac{\rho
r}{g}\,y^{J}_{5}\,y^{K}_{5}+\frac{r^{2}}{4\pi{}G}\,y^{J}_{7}\,y^{K}_{7}\,,$
(E.9)

in which $y^{J}_{7}=y^{J}_{6}+(4\pi{}G/g_{r})y^{J}_{2}$. Adding Eqs.
(E.2)-(E.3) and Eq. (E.8) yields Eq. (12) of Okubo and Saito [1983] modified
to allow for discontinuities at the shell-ocean boundary::

$\left[Y^{JK}\right]_{0}^{R_{c}}+\left[Y^{JK}\right]_{R_{c}}^{R^{-}_{o}}+\left[Y^{JK}\right]_{R^{+}_{o}}^{R_{s}}=\int_{0}^{R_{s}}{\cal
L}^{JK}\,dr+\left[\frac{r^{2}}{\rho_{o}g_{r}}\,y^{J}_{2}\,y^{K}_{2}\right]_{R_{c}}^{R_{o}}\,,$
(E.10)

where ${\cal L}^{JK}={\cal L}^{JK}_{s}$ in the solid layers and ${\cal
L}^{JK}={\cal L}^{JK}_{f}$ in the fluid layer.

###  E.2 Saito-Molodensky relation

The integral of motion implies a relation between two sets of solutions of the
equations of motion. Since the RHS of Eq. (E.10) is symmetric under the
exchange of $J$ and $K$,

$\left[Y^{JK}-Y^{KJ}\right]_{0}^{R_{c}}+\left[Y^{JK}-Y^{KJ}\right]_{R_{c}}^{R_{o}}+\left[Y^{JK}-Y^{KJ}\right]_{R_{o}}^{R_{s}}=0\,.$
(E.11)

Applying the boundary and discontinuity conditions for the surface load
($J=L$) and the internal load ($K=I$), I obtain a relation between the surface
load and internal load Love numbers:

$k_{s}^{I}-h_{s}^{I}=x\,\Big{(}k_{o}^{L}-\gamma_{o}\,h_{o}^{L}\Big{)}\,,$
(E.12)

where $\gamma_{o}=g_{o}/g_{s}$. This is a new example of the Saito-Molodensky
relations already known between tidal, surface load, and shear Love numbers
[Molodensky, 1977; Saito, 1978].

###  E.3 Variation of Love numbers

Applying variational calculus on the integral of motion (Eq. (E.10)), Okubo
and Saito [1983] determine how variations of $(\mu,\kappa,\rho)$ affect the
tidal, surface load and shear Love numbers. I will only consider here
variations with respect to the elastic parameters (the density does not
change). The equations of Okubo and Saito [1983] must be slightly modified to
allow for discontinuities at the shell-ocean boundary. Define

$Z^{JK}=r^{2}\left(\Big{(}y^{J}_{1}\,\delta y^{K}_{2}-y^{J}_{2}\,\delta
y^{K}_{1}\Big{)}-\delta_{n}\Big{(}y^{J}_{3}\,\delta
y^{K}_{4}-y^{J}_{4}\,\delta y^{K}_{3}\Big{)}+\frac{y^{J}_{5}\,\delta
y^{K}_{6}-y^{J}_{6}\,\delta y^{K}_{5}}{4\pi{}G}\right),$ (E.13)

where $\delta_{n}=-n(n+1)$. Then Eq. (16) of Okubo and Saito [1983] becomes

$\left[Z^{JK}\right]_{0}^{R_{c}}+\left[Z^{JK}\right]_{R_{c}}^{R^{-}_{o}}+\left[Z^{JK}\right]_{R^{+}_{o}}^{R_{s}}=\int_{R_{o}}^{R_{s}}\left({\cal
H}^{JK}_{\mu}\,\delta\mu+{\cal H}^{JK}_{\kappa}\,\delta\kappa\right)dr\,.$
(E.14)

Applying the boundary and discontinuity conditions for the surface load
($J=L,I$) and the internal load ($K=L,I$), I obtain

$\displaystyle\delta h_{s}^{L}-\delta k_{s}^{L}$ $\displaystyle=$
$\displaystyle\chi^{-1}\int_{d}\left({\cal H}^{LL}_{\mu}\,\delta\mu+{\cal
H}^{LL}_{\kappa}\,\delta\kappa\right)dr\,,$ (E.15) $\displaystyle\delta
h_{s}^{I}-\delta k_{s}^{I}$ $\displaystyle=$
$\displaystyle\chi^{-1}\int_{d}\left({\cal H}^{LI}_{\mu}\,\delta\mu+{\cal
H}^{LI}_{\kappa}\,\delta\kappa\right)dr\,,$ (E.16) $\displaystyle
x\,\Big{(}\gamma_{o}\,\delta h_{o}^{L}-\delta k_{o}^{L}\Big{)}$
$\displaystyle=$ $\displaystyle\chi^{-1}\int_{d}\left({\cal
H}^{IL}_{\mu}\,\delta\mu+{\cal H}^{IL}_{\kappa}\,\delta\kappa\right)dr\,,$
(E.17) $\displaystyle x\,\Big{(}\gamma_{o}\,\delta h_{o}^{I}-\delta
k_{o}^{I}\Big{)}$ $\displaystyle=$ $\displaystyle\chi^{-1}\int_{d}\left({\cal
H}^{II}_{\mu}\,\delta\mu+{\cal H}^{II}_{\kappa}\,\delta\kappa\right)dr\,.$
(E.18)

where $\chi=(2n+1)R_{s}/(4\pi G)$. As ${\cal H}^{JK}_{\mu}$ and ${\cal
H}^{JK}_{\kappa}$ are symmetric under the exchange of $J$ and $K$, the LHS of
the second and third equations must be identical: $\delta h_{s}^{I}-\delta
k_{s}^{I}=x\,\Big{(}\gamma_{o}\,\delta h_{o}^{L}-\delta k_{o}^{L}\Big{)}$,
which also follows from the Saito-Molodensky relation (Eq. (E.12)). In the
main text, these identities are expressed in terms of partial derivatives of
Love numbers (Eqs. (122)-(123)).

##  Appendix F MSI $\mu$-invariance and fluid limit

###  F.1 MSI $\mu$-invariance

For an incompressible body with a homogeneous shell, ZDI $\mu$-invariance
implies MSI $\mu$-invariance. I will prove this statement by showing that MSI
shape ratio is $\mu$-invariant if the conditions on Love numbers required for
ZDI $\mu$-invariance (Eqs. (93)-(95)) are satisfied. The MSI shape ratio is
given by

$S^{\rm MSI}_{n}=\frac{E_{os,p}+\beta\,E_{oo,p}}{E_{ss,p}+\beta\,E_{so,p}}\,.$
(F.1)

This formula is equivalent to Eq. (131) because $E_{os,p}=E_{so,p}$, but we
will see that there is a benefit to keep these coefficients formally distinct
when proving $\mu$-invariance. The essence of the proof consists in expressing
stress-energy coefficients in terms of partial derivatives of Love numbers,
eliminating gravitational Love numbers with gravity-deformation relations, and
finally factoring out derivatives with the conditions for ZDI
$\mu$-invariance. The resulting formula depends on the same products of Love
numbers as in zero deflection isostasy.

For the purpose of the proof, the stress/energy coefficients are written as

$E_{jk,p}=\frac{1}{N^{2}}\left(\hat{h}_{j}^{I}\left(\hat{h}_{k}^{I}\,E^{LL}_{p}-\hat{h}_{k}^{L}\,E^{LI}_{p}\right)+\hat{h}_{j}^{L}\left(\hat{h}_{k}^{L}\,E^{II}_{p}-\hat{h}_{k}^{I}\,E^{IL}_{p}\right)\right).$
(F.2)

This formula is equivalent to Eq. (126) because $E^{IL}_{p}=E^{LI}_{p}$. Using
Eqs. (122)-(123), I write the quantities within the interior brackets as

$\displaystyle\hat{h}_{k}^{I}\,E^{LL}_{p}-\hat{h}_{k}^{L}\,E^{LI}_{p}$
$\displaystyle=$
$\displaystyle\bar{\chi}\left(\hat{h}_{k}^{I}\,\partial_{\mu}\Big{(}h_{s}^{L}-k_{s}^{L}\Big{)}-\hat{h}_{k}^{L}\,\partial_{\mu}\Big{(}h_{s}^{I}-k_{s}^{I}\Big{)}\right),$
(F.3) $\displaystyle\hat{h}_{k}^{L}\,E^{II}_{p}-\hat{h}_{k}^{I}\,E^{IL}_{p}$
$\displaystyle=$
$\displaystyle\bar{\chi}\,x\left(\hat{h}_{k}^{L}\,\partial_{\mu}\Big{(}\gamma_{o}\,h_{o}^{I}-k_{o}^{I}\Big{)}-\hat{h}_{k}^{I}\,\partial_{\mu}\Big{(}\gamma_{o}\,h_{o}^{L}-k_{o}^{L}\Big{)}\right)\,,$
(F.4)

where common constants are gathered into

$\bar{\chi}=\frac{p}{2}\,(\mu_{\rm
s})^{p}\,(g_{s}R_{s})^{2}\,\frac{(2n+1)R_{s}}{4\pi G}\,.$ (F.5)

Eliminating the gravitational Love numbers with Eq. (46), I get

$\displaystyle\hat{h}_{k}^{I}\,E^{LL}_{p}-\hat{h}_{k}^{L}\,E^{LI}_{p}$
$\displaystyle=$
$\displaystyle\bar{\chi}\left(\left(a-1\right)Y_{ks}+b\,Y_{ko}\right),$ (F.6)
$\displaystyle\hat{h}_{k}^{L}\,E^{II}_{p}-\hat{h}_{k}^{I}\,E^{IL}_{p}$
$\displaystyle=$
$\displaystyle\bar{\chi}\,x\left(-c\,Y_{ks}+d\,Y_{ko}\right)\,,$ (F.7)

where

$Y_{jk}=\hat{h}_{j}^{L}\,\partial_{\mu}h_{k}^{I}-\hat{h}_{j}^{I}\,\partial_{\mu}h_{k}^{L}\,.$
(F.8)

The conditions for ZDI $\mu$-invariance (Eqs. (93)-(95)) are equivalent to

$\displaystyle\partial_{\mu}\Big{(}\hat{h}_{s}^{I}-\hat{h}_{o}^{I}\,h_{s}^{L}/\hat{h}_{o}^{L}\big{)}$
$\displaystyle=$ $\displaystyle 0\,,$ (F.9)
$\displaystyle\partial_{\mu}\left(\hat{h}_{o}^{L}-\hat{h}_{s}^{L}\,h_{o}^{I}/\hat{h}_{s}^{I}\right)$
$\displaystyle=$ $\displaystyle 0\,,$ (F.10)
$\displaystyle\partial_{\mu}\left(\Big{(}h_{o}^{I\circ}\,\hat{h}_{s}^{L}-h_{s}^{L\circ}\,\hat{h}_{o}^{I}\Big{)}/\hat{h}_{s}^{I}\right)$
$\displaystyle=$ $\displaystyle 0\,.$ (F.11)

Recall that $h_{o}^{L}=\hat{h}_{o}^{L}$ and $h_{s}^{I}=\hat{h}_{s}^{I}$ while
$h_{s}^{L}=h_{s}^{L\circ}+\hat{h}_{s}^{L}$ and
$h_{o}^{I}=h_{o}^{I\circ}+\hat{h}_{o}^{I}$. The first equation results from
combining Eq. (93) with Eq. (94), by taking the product of the arguments of
the derivatives; the second one is identical to Eq. (95); the third one
results from subtracting Eq. (95) from Eq. (94). Expanding Eqs. (F.9)-(F.11)
and rearranging terms, I obtain three relations between the four different
$Y_{jk}$:

$\displaystyle Y_{os}/h_{s}^{L}$ $\displaystyle=$ $\displaystyle
Y_{oo}/h_{o}^{L}\,,$ (F.12) $\displaystyle Y_{so}/h_{o}^{I}$ $\displaystyle=$
$\displaystyle Y_{ss}/h_{s}^{I}\,,$ (F.13) $\displaystyle Y_{oo}/h_{o}^{L}$
$\displaystyle=$ $\displaystyle-(cx/b)\,Y_{ss}/h_{s}^{I}\,.$ (F.14)

For the last equation, I relied on the fact that
$h^{I\circ}_{o}/h^{L\circ}_{s}=cx/b$. With these relations, I can express Eqs.
(F.3)-(F.4) in terms of only one $Y_{jk}$, for example $Y_{ss}$, and
substitute the results into Eq. (F.2). Using in addition the Saito-Molodensky
relation (Eq. (48)), I write the stress/energy coefficients as

$\displaystyle E_{js,p}$ $\displaystyle=$
$\displaystyle\Big{(}-c\,X_{js}+d\,X_{jo}\Big{)}\left(x\,\frac{\bar{\chi}}{N^{2}}\,\frac{Y_{ss}}{h^{I}_{s}}\right),$
(F.15) $\displaystyle E_{jo,p}$ $\displaystyle=$
$\displaystyle\frac{c}{b}\,\Big{(}(a-1)\,X_{js}+b\,X_{jo}\Big{)}\left(x\,\frac{\bar{\chi}}{N^{2}}\,\frac{Y_{ss}}{h^{I}_{s}}\right),$
(F.16)

where

$X_{jk}=\hat{h}_{j}^{L}\,h_{k}^{I}-\hat{h}_{j}^{I}\,h_{k}^{L}\,.$ (F.17)

Substituting Eqs. (F.15)-(F.16) into (Eq. (F.1), I simplify the derivatives of
Love numbers out of the shape ratio:

$S^{\rm
MSI}_{n}=\frac{-c\left(\beta\left(1-a\right)+b\right)X_{os}+b\left(\beta\,c+d\right)X_{oo}}{-c\left(\beta\left(1-a\right)+b\right)X_{ss}+b\left(\beta\,c+d\right)X_{so}}\,.$
(F.18)

A quick comparison with the ZDI shape ratio (Eq. (92)) shows that that both
formulas depend on the same combinations of Love numbers $X_{jk}$. Therefore
the conditions for ZDI $\mu$-invariance also imply MSI $\mu$-invariance, at
least for an incompressible body with a homogeneous shell. The added benefit
mentioned at the start of the section is that the ZDI and MSI shape ratios are
actually identical if the isostatic family parameters $\alpha$ and $\beta$ are
such that

$\alpha=-\frac{b}{c}\,\frac{\beta\,c+d}{\beta\left(1-a\right)+b}\,.$ (F.19)

In the main text (using the result of Appendix F.2), I derive this duality in
a simpler way by using the fluid limit under the assumption of
$\mu$-invariance.

Although not required for the computation of the shape ratio, it is
interesting to know the dependence on $\mu$ of the stress/energy coefficients.
It is found by using the generic form of Love numbers appropriate to a
homogeneous shell (Appendix. D). Using the identities (D.4)-(D.5) and (D.10),
one can check that

$\displaystyle N$ $\displaystyle=$
$\displaystyle\frac{h_{s}^{L\circ}\,h_{o}^{I\circ}\,F\mu_{s}^{2}}{D+E\mu_{\rm
s}+F\mu_{\rm s}^{2}}\,,$ (F.20) $\displaystyle X_{js}$ $\displaystyle=$
$\displaystyle-h_{s}^{L\circ}\,\frac{\hat{B}_{j}^{I}\,\mu_{\rm s}}{D+E\mu_{\rm
s}+F\mu_{\rm s}^{2}}\,,$ (F.21) $\displaystyle X_{jo}$ $\displaystyle=$
$\displaystyle h_{o}^{I\circ}\,\frac{\hat{B}_{j}^{L}\,\mu_{\rm s}}{D+E\mu_{\rm
s}+F\mu_{\rm s}^{2}}\,,$ (F.22) $\displaystyle\frac{Y_{ss}}{h_{s}^{I}}$
$\displaystyle=$ $\displaystyle h_{s}^{L\circ}\,\frac{F\mu_{\rm
s}}{D+E\mu_{\rm s}+F\mu_{\rm s}^{2}}\,.$ (F.23)

The stress/energy coefficients (Eqs. (F.15)-(F.16)) become

$\displaystyle E_{js,p}$ $\displaystyle=$
$\displaystyle\frac{\bar{\chi}}{\mu_{\rm
s}^{2}}\,\frac{b\,\hat{B}_{j}^{I}+d\,x\,\hat{B}_{j}^{L}}{h_{s}^{L\circ}\,h_{o}^{I\circ}\,F}\,,$
(F.24) $\displaystyle E_{jo,p}$ $\displaystyle=$
$\displaystyle\frac{\bar{\chi}}{\mu_{\rm
s}^{2}}\,\frac{(1-a)\,\hat{B}_{j}^{I}+c\,x\,\hat{B}_{j}^{L}}{h_{s}^{L\circ}\,h_{o}^{I\circ}\,F}\,.$
(F.25)

Since $\bar{\chi}\sim\mu_{\rm s}^{p}$, the stress/energy coefficients are
proportional to $(\mu_{\rm s})^{p-2}$. Observe that the substitution of
$X_{jk}$ (Eqs. (F.21)-(F.22)) into the general ZDI shape ratio (Eq. (92)),
taking into account Eq. (D.10), confirms the formula derived in the fluid
limit (Eq. (101)).

###  F.2 MSI fluid limit

Consider an incompressible body with a homogeneous shell. Under the assumption
of $\mu$-invariance, the MSI loading ratio can be computed by taking the fluid
limit of Eq. (130):

$\zeta_{n}^{\circ}=\frac{\left(\dot{h}_{s}^{I}\,H^{LL}_{0}-\dot{h}_{s}^{L}\,H^{LI}_{0}\right)+\beta\left(\dot{h}_{o}^{I}\,H^{LL}_{0}-\dot{h}_{o}^{L}\,H^{LI}_{0}\right)}{\left(\dot{h}_{s}^{L}\,H^{II}_{0}-\dot{h}_{s}^{I}\,H^{LI}_{0}\right)+\beta\left(\dot{h}_{o}^{L}\,H^{II}_{0}-\dot{h}_{o}^{I}\,H^{LI}_{0}\right)}\,,$
(F.26)

where $H^{JK}_{0}$ results from the fluid limit of Eq. (121),

$H^{JK}_{0}=\frac{1}{\chi}\left[\int_{d}\,{\cal
H}^{JK}_{\mu}\,dr\right]_{\mu=0}\,,$ (F.27)

in which $\chi=(2n+1)R_{s}/(4\pi G)$. The expressions within brackets can be
transformed into derivatives of Love numbers with Eqs. (122)-(123).
Gravitational Love numbers are then eliminated with Eq. (46). The resulting
expressions are all proportional to
$X_{so}^{0}=\dot{h}_{s}^{L}\,\dot{h}_{o}^{I}-\dot{h}_{s}^{I}\,\dot{h}_{o}^{L}$:

$\displaystyle\dot{h}_{s}^{I}\,H^{LL}_{0}-\dot{h}_{s}^{L}\,H^{LI}_{0}$
$\displaystyle=$
$\displaystyle\dot{h}_{s}^{I}\left(\dot{h}_{s}^{L}-\dot{k}_{s}^{L}\right)-\dot{h}_{s}^{L}\left(\dot{h}_{s}^{I}-\dot{k}_{s}^{I}\right)\,\,=\,\,b\,X_{so}^{0}\,,$
$\displaystyle\dot{h}_{o}^{I}\,H^{LL}_{0}-\dot{h}_{o}^{L}\,H^{LI}_{0}$
$\displaystyle=$
$\displaystyle\dot{h}_{o}^{I}\left(\dot{h}_{s}^{L}-\dot{k}_{s}^{L}\right)-\dot{h}_{o}^{L}\left(\dot{h}_{s}^{I}-\dot{k}_{s}^{I}\right)\,\,=\,\,\left(1-a\right)X_{so}^{0}\,,$
$\displaystyle\dot{h}_{s}^{L}\,H^{II}_{0}-\dot{h}_{s}^{I}\,H^{LI}_{0}$
$\displaystyle=$
$\displaystyle\dot{h}_{s}^{L}\,x\left(\gamma_{o}\,\dot{h}_{o}^{I}-\dot{k}_{o}^{I}\right)-\dot{h}_{s}^{I}\,x\left(\gamma_{o}\,\dot{h}_{o}^{L}-\dot{k}_{o}^{L}\right)\,\,=\,\,d\,x\,X_{so}^{0}\,,$
$\displaystyle\dot{h}_{o}^{L}\,H^{II}_{0}-\dot{h}_{o}^{I}\,H^{LI}_{0}$
$\displaystyle=$
$\displaystyle\dot{h}_{o}^{L}\,x\left(\gamma_{o}\,\dot{h}_{o}^{I}-\dot{k}_{o}^{I}\right)-\dot{h}_{o}^{I}\,x\left(\gamma_{o}\,\dot{h}_{o}^{L}-\dot{k}_{o}^{L}\right)\,\,=\,\,c\,x\,X_{so}^{0}\,,$
(F.28)

with the coefficients $(a,b,c,d)$ of Table 3. In the fluid limit, the MSI
loading ratio is thus given by

$\zeta_{n}^{\circ}=\frac{1}{x}\,\frac{b+\beta\left(1-a\right)}{d+\beta\,c}\,.$
(F.29)

##  Appendix G Isostatic ratios if three incompressible homogeneous layers

Consider an incompressible body with three homogeneous layers. The ZDI
topographic ratio, given by Eq. (106), is a linear fractional transform of
$\alpha$ (see Eq. (45)). The shape ratio and compensation factor are related
to this formula by the $S_{n}$-$T_{n}$ and $F_{n}$-$T_{n}$ relations (Eqs.
(49) and (52)). Since these relations, as well as Eq. (106), are linear
fractional transformations, they can be composed through the product of their
associated matrices (Eq. (45))). The MSI isostatic ratios are obtained by
considering $\alpha_{\rm MSI}$ as a linear fractional transformation of
$\beta$ (Eq. (138)). In the following, $(a,b,c,d)$ are the parameters defined
in Table 3 in terms of the interior structure and of the harmonic degree $n$.
The dependence of $P_{n}$, $Q_{n}$, and $Q^{\prime}_{n}$ on $x$ is implicit.

The ZDI shape ratio results from the composition of
$S_{n}=LFT(\gamma_{o}T_{n})$ with $T_{n}=T_{n}^{\rm ZDI}=LFT(\alpha)$:

$\displaystyle S_{n}^{\rm ZDI}$ $\displaystyle=$ $\displaystyle
LFT\left[\left(\begin{array}[]{cc}1-a&c\\\
b&d\end{array}\right)\left(\begin{array}[]{cc}-\xi_{s1}\,P_{n}&-\xi_{s1}\,Q_{n}\\\
\Delta\xi_{1}\,Q^{\prime}_{n}&\Delta\xi_{1}\,x^{2}P_{n}\end{array}\right)\right](\alpha)$
(G.5) $\displaystyle=$
$\displaystyle\frac{\alpha\left(-\left(1-a\right)\xi_{s1}\,P_{n}+c\,\Delta\xi_{1}\,Q^{\prime}_{n}\right)-\left(1-a\right)\xi_{s1}\,Q_{n}+c\,\Delta\xi_{1}\,x^{2}P_{n}}{\alpha\left(-b\,\xi_{s1}\,P_{n}+d\,\Delta\xi_{1}\,Q^{\prime}_{n}\right)-b\,\xi_{s1}\,Q_{n}+d\,\Delta\xi_{1}\,x^{2}P_{n}}\,.$
(G.6)

The ZDI compensation factor results from the composition of
$\xi_{sn}F_{n}=LFT(\gamma_{o}T_{n})$ with $T_{n}=T_{n}^{\rm ZDI}=LFT(\alpha)$:

$\displaystyle\xi_{sn}F_{n}^{\rm ZDI}$ $\displaystyle=$ $\displaystyle
LFT\left[\left(\begin{array}[]{cc}b&a\,d+b\,c\\\
b&d\end{array}\right)\left(\begin{array}[]{cc}-\xi_{s1}\,P_{n}&-\xi_{s1}\,Q_{n}\\\
\Delta\xi_{1}\,Q^{\prime}_{n}&\Delta\xi_{1}\,x^{2}P_{n}\end{array}\right)\right](\alpha)$
(G.11) $\displaystyle=$
$\displaystyle\frac{\alpha\left(-b\,\xi_{s1}\,P_{n}+\left(ad+bc\right)\Delta\xi_{1}\,Q^{\prime}_{n}\right)-b\,\xi_{s1}\,Q_{n}+\left(ad+bc\right)\Delta\xi_{1}\,x^{2}P_{n}}{\alpha\left(-b\,\xi_{s1}\,P_{n}+d\,\Delta\xi_{1}\,Q^{\prime}_{n}\right)-b\,\xi_{s1}\,Q_{n}+d\,\Delta\xi_{1}\,x^{2}P_{n}}\,.$
(G.12)

The MSI topographic ratio results from the composition of $T_{n}^{\rm
ZDI}=LFT(\alpha)$ with $\alpha=LFT(\beta)$:

$\displaystyle\gamma_{o}\,T_{n}^{\rm MSI}$ $\displaystyle=$ $\displaystyle
LFT\left[\left(\begin{array}[]{cc}-\xi_{s1}\,P_{n}&-\xi_{s1}\,Q_{n}\\\
\Delta\xi_{1}\,Q^{\prime}_{n}&\Delta\xi_{1}\,x^{2}P_{n}\end{array}\right)\left(\begin{array}[]{cc}-b\,c&-b\,d\\\
\left(1-a\right)c&b\,c\end{array}\right)\right](\beta)$ (G.17)
$\displaystyle=$
$\displaystyle-\frac{\xi_{s1}}{\Delta\xi_{1}}\,\frac{\beta\,c\left(-b\,P_{n}+\left(1-a\right)Q_{n}\right)+b\left(-d\,P_{n}+c\,Q_{n}\right)}{\beta\,c\left(-b\,Q^{\prime}_{n}+\left(1-a\right)x^{2}\,P_{n}\right)+b\left(c\,x^{2}\,P_{n}-d\,Q^{\prime}_{n}\right)}\,.$
(G.18)

The MSI shape ratio results from the composition of
$S_{n}=LFT(\gamma_{o}T_{n})$ with $T_{n}=T_{n}^{\rm ZDI}=LFT(\alpha)$ and
$\alpha=LFT(\beta)$, which is a bit too long to expand here:

$S_{n}^{\rm MSI}=LFT\left[\left(\begin{array}[]{cc}1-a&c\\\
b&d\end{array}\right)\left(\begin{array}[]{cc}-\xi_{s1}\,P_{n}&-\xi_{s1}\,Q_{n}\\\
\Delta\xi_{1}\,Q^{\prime}_{n}&\Delta\xi_{1}\,x^{2}P_{n}\end{array}\right)\left(\begin{array}[]{cc}-b\,c&-b\,d\\\
\left(1-a\right)c&b\,c\end{array}\right)\right](\beta)\,.$ (G.19)

The MSI compensation factor results from the composition of
$\xi_{sn}F_{n}=LFT(\gamma_{o}T_{n})$ with $T_{n}=T_{n}^{\rm ZDI}=LFT(\alpha)$
and $\alpha=LFT(\beta)$:

$\xi_{sn}F_{n}^{\rm MSI}=LFT\left[\left(\begin{array}[]{cc}b&a\,d+b\,c\\\
b&d\end{array}\right)\left(\begin{array}[]{cc}-\xi_{s1}\,P_{n}&-\xi_{s1}\,Q_{n}\\\
\Delta\xi_{1}\,Q^{\prime}_{n}&\Delta\xi_{1}\,x^{2}P_{n}\end{array}\right)\left(\begin{array}[]{cc}-b\,c&-b\,d\\\
\left(1-a\right)c&b\,c\end{array}\right)\right](\beta)\,.$ (G.20)

## References

  * Airy [1855] Airy, G. B. (1855), On the computation of the effect of the attraction of mountain-masses, as disturbing the apparent astronomical latitude of stations in geodetic surveys, Phil. Trans. R. Soc., 145, 101–104, doi:10.1098/rstl.1855.0003.
  * Banerdt et al. [1982] Banerdt, W. B., R. S. Saunders, R. J. Phillips, and N. H. Sleep (1982), Thick shell tectonics on one-plate planets - Applications to Mars, J. Geophys. Res., 87, 9723–9733, doi:10.1029/JB087iB12p09723.
  * Beuthe [2013] Beuthe, M. (2013), Spatial patterns of tidal heating, Icarus, 223, 308–329, doi:10.1016/j.icarus.2012.11.020.
  * Beuthe [2015] Beuthe, M. (2015), Tidal Love numbers of membrane worlds: Europa, Titan, and Co., Icarus, 258, 239–266, doi:10.1016/j.icarus.2015.06.008.
  * Beuthe [2020a] Beuthe, M. (2020a), Equivalence between elastic, viscoelastic, and viscous isostasy, in Europlanet Science Congress 2020, vol. 14, pp. EPSC2020–863, doi:10.5194/epsc2020-863.
  * Beuthe [2020b] Beuthe, M. (2020b), IsostasyWithLove: Mathematica notebooks (version 1.0), Zenodo, doi:10.5281/zenodo.4297495.
  * Beuthe et al. [2016] Beuthe, M., A. Rivoldini, and A. Trinh (2016), Enceladus’s and Dione’s floating ice shells supported by minimum stress isostasy, Geophys. Res. Lett., 43, 10,088–10,096, doi:10.1002/2016GL070650.
  * Choblet et al. [2017] Choblet, G., G. Tobie, C. Sotin, M. Běhounková, O. Čadek, F. Postberg, and O. Souček (2017), Powering prolonged hydrothermal activity inside Enceladus, Nature Astronomy, 1, 841–847, doi:10.1038/s41550-017-0289-8.
  * Dahlen [1982] Dahlen, F. A. (1982), Isostatic geoid anomalies on a sphere, J. Geophys. Res., 87, 3943–3947, doi:10.1029/JB087iB05p03943.
  * Farrell [1972] Farrell, W. E. (1972), Deformation of the Earth by surface loads, Rev. Geophys., 10, 761–797, doi:10.1029/RG010i003p00761.
  * Garland [1965] Garland, G. D. (1965), The Earth’s shape and gravity, 183 pp., Pergamon Press, New York, doi:10.1016/C2013-0-05382-0.
  * Greff-Lefftz and Legros [1997] Greff-Lefftz, M., and L. Legros (1997), Some remarks about the degree-one deformation of the Earth, Geophys. J. Int., 131, 699–723, doi:10.1111/j.1365-246X.1997.tb06607.x.
  * Greff-Lefftz et al. [2010] Greff-Lefftz, M., L. Métivier, and J. Besse (2010), Dynamic mantle density heterogeneities and global geodetic observables, Geophys. J. Int., 180, 1080–1094, doi:10.1111/j.1365-246X.2009.04490.x.
  * Heiskanen and Vening Meinesz [1958] Heiskanen, W. A., and F. A. Vening Meinesz (1958), The Earth and its gravity field, McGraw-Hill, New York.
  * Hemingway and Matsuyama [2017] Hemingway, D. J., and I. Matsuyama (2017), Isostatic equilibrium in spherical coordinates and implications for crustal thickness on the Moon, Mars, Enceladus, and elsewhere, Geophys. Res. Lett., 44(15), 7695–7705, doi:10.1002/2017GL073334.
  * Hinderer and Legros [1989] Hinderer, J., and H. Legros (1989), Elasto-gravitational deformation, relative gravity changes and earth dynamics, Geophys. J. Int., 97, 481–495, doi:10.1111/j.1365-246X.1989.tb00518.x.
  * Jaeger et al. [2007] Jaeger, J. C., N. G. W. Cook, and R. W. Zimmerman (2007), Fundamentals of Rock Mechanics, 4th edition, Blackwell.
  * Jeffreys [1932] Jeffreys, H. (1932), On the stresses in the Earth’s crust required to support surface inequalities. (Second paper.), Geophys. Suppl. MNRAS, 3, 60–69, doi:10.1111/j.1365-246X.1932.tb00400.x.
  * Jeffreys [1943] Jeffreys, H. (1943), The stress-differences in the Earth’s shell, Geophys. Suppl. MNRAS, 5, 71–89, doi:10.1111/j.1365-246X.1943.tb00344.x.
  * Jeffreys [1959] Jeffreys, H. (1959), The Earth, 4th Ed., Cambridge University Press, Cambridge.
  * Kaula [1963] Kaula, W. M. (1963), Elastic models of the mantle corresponding to variations in the external gravity field, J. Geophys. Res., 68, 4967–4978, doi:10.1029/JZ068i017p04967.
  * Kaula [1968] Kaula, W. M. (1968), An Introduction to Planetary Physics: The Terrestrial Planets, John Wiley, New York.
  * Lambeck [1980] Lambeck, K. (1980), Estimates of stress differences in the crust from isostatic considerations, J. Geophys. Res., 85, 6397–6402, doi:10.1029/JB085iB11p06397.
  * Lambert [1930] Lambert, W. D. (1930), The form of the geoid on the hypothesis of complete isostatic compensation, Bull. Géod., 26, 98–106, doi:0.1007/BF03030028.
  * Liao et al. [2020] Liao, Y., F. Nimmo, and J. A. Neufeld (2020), Heat production and tidally driven fluid flow in the permeable core of Enceladus, J. Geophys. Res., 125(9), e06209, doi:10.1029/2019JE006209.
  * Melosh [2011] Melosh, H. J. (2011), Planetary Surface Processes, Cambridge University Press, Cambridge, doi:10.1017/CBO9780511977848.
  * Molodensky [1977] Molodensky, S. M. (1977), Relation between Love numbers and load factors, Izv., Phys. Solid Earth, 13, 147–149.
  * Okubo and Saito [1983] Okubo, S., and M. Saito (1983), Partial derivative of Love numbers, Bull. Géod., 57(1-4), 167–179, doi:10.1007/BF02520923.
  * Roberts [2015] Roberts, J. H. (2015), The fluffy core of Enceladus, Icarus, 258, 54–66, doi:10.1016/j.icarus.2015.05.033.
  * Sabadini and Vermeersen [2004] Sabadini, R., and B. Vermeersen (2004), Global Dynamics of the Earth, Kluwer Academic Publishers, Dordrecht.
  * Saito [1974] Saito, M. (1974), Some problems of static deformation of the earth, J. Phys. Earth, 22, 123–140, doi:10.4294/jpe1952.22.123.
  * Saito [1978] Saito, M. (1978), Relationship between tidal and load Love numbers, J. Phys. Earth, 26, 13–16, doi:10.4294/jpe1952.26.13.
  * Takeuchi and Saito [1972] Takeuchi, H., and M. Saito (1972), Seismic surface waves, in Methods in Computational Physics, vol. 1, edited by Bolt, B.A., pp. 217–295, Academic Press, New York, doi:10.1016/B978-0-12-460811-5.50010-6.
  * Tobie et al. [2005] Tobie, G., A. Mocquet, and C. Sotin (2005), Tidal dissipation within large icy satellites: Applications to Europa and Titan, Icarus, 177, 534–549, doi:10.1016/j.icarus.2005.04.006.
  * Travis and Schubert [2015] Travis, B. J., and G. Schubert (2015), Keeping Enceladus warm, Icarus, 250, 32–42, doi:10.1016/j.icarus.2014.11.017.
  * Trinh [2019] Trinh, A. (2019), Modelling the interior of Enceladus: A combined view from gravity, topography, and libration, Ph.D. thesis, Université catholique de Louvain, Louvain-la-Neuve, Belgium, http://hdl.handle.net/2078.1/216772.
  * Turcotte and Schubert [2014] Turcotte, D. L., and G. Schubert (2014), Geodynamics, 3rd ed., 626 pp., Cambridge University Press, Cambridge, doi:10.1017/CBO9780511843877.
  * Čadek et al. [2019] Čadek, O., O. Souček, and M. Běhounková (2019), Is Airy isostasy applicable to icy moons?, Geophys. Res. Lett., 46(24), 14,299–14,306, doi:10.1029/2019GL085903.
  * Vening Meinesz [1946] Vening Meinesz, F. A. (1946), The indirect isostatic or Bowie reduction and the equilibrium figure of the Earth, Bull. Géod., 1, 33–107, doi:10.1007/BF02519025.
  * Wieczorek [2015] Wieczorek, M. (2015), Gravity and topography of the terrestrial planets, in Treatise on Geophysics, vol. 10, edited by G. Schubert, second ed., chap. 5, pp. 153–193, Elsevier, Oxford, doi:10.1016/B978-0-444-53802-4.00169-X.
  * Wu and Peltier [1982] Wu, P., and W. R. Peltier (1982), Viscous gravitational relaxation, Geophys. J. R. astr. Soc., 70(2), 435–485, doi:10.1111/j.1365-246X.1982.tb04976.x.


# QoS-Aware Load Balancing in Wireless Networks using Clipped Double
Q-Learning

Pedro Enrique Iturria-Rivera, , Melike Erol-Kantarci  School of Electrical
Engineering and Computer ScienceUniversity of Ottawa
Ottawa, Canada
{pitur008<EMAIL_ADDRESS>
###### Abstract

In recent years, long-term evolution (LTE) and 5G NR (5th Generation New
Radio) technologies have showed great potential to utilize Machine Learning
(ML) algorithms in optimizing their operations, both thanks to the
availability of fine-grained data from the field, as well as the need arising
from growing complexity of networks. The aforementioned complexity sparked
mobile operators’ attention as a way to reduce the capital expenditures
(CAPEX) and the operational (OPEX) expenditures of their networks through
network management automation (NMA). NMA falls under the umbrella of Self-
Organizing Networks (SON) in which 3GPP has identified some challenges and
opportunities in load balancing mechanisms for the Radio Access Networks
(RANs). In the context of machine learning and load balancing, several studies
have focused on maximizing the overall network throughput or the resource
block utilization (RBU). In this paper, we propose a novel Clipped Double
Q-Learning (CDQL)-based load balancing approach considering resource block
utilization, latency and the Channel Quality Indicator (CQI). We compare our
proposal with a traditional handover algorithm and a resource block
utilization based handover mechanism. Simulation results reveal that our
scheme is able to improve throughput, latency, jitter and packet loss ratio in
comparison to the baseline algorithms.

###### Index Terms:

load balancing, double clipped q-learning, wireless networks.

##  I Introduction

The exponential increase in mobile network usage alongside the growing data
consumption demand by several use cases such as AR/VR (Augmented
Reality/Virtual Reality) and video streaming, have led wireless networks to
evolve in order to satisfy the pressing requirements. More specifically,
video-dominated applications related to video streaming, video conferencing
and high quality buffered media have faced a great spike due the recent impact
of the COVID-19 pandemic. Some studies showed that after the lockdown, an
increase of 215-285% in VoIP and videoconferencing traffic and a 20-40%
increment in streaming and web video consumption have been observed [1].
Besides, the bandwidth requirement of the aforementioned use cases, they also
have tight delay requirements. The trend for multimedia usage is expected to
continue after the pandemic and mobile network operators will shoulder a
majority of the load. This calls for optimized use of resources at the RAN,
the transport network and the core.

LTE (Long Term Evolution) and its successor 5G (5th Generation) NR (New Radio)
are able to support self-optimization functionalities which 3GPP has
identified in the context of SON. These also include load balancing. Load
balancing for the RAN, involves handing over UEs (User Equipments) to less
occupied base stations. Naturally, a base station’s load would be in relation
to the number of UEs that are associated to it and the traffic demand of these
UEs.

As part of the handover mechanism, each UE in the network will send periodical
measurement reports to its respective serving cell. In practical terms, the
serving cell will send via a Radio Resource Control (RRC) message the
indication of what type of measurements each connected UE must gather and
consequently report of the vicinity cells and of itself. Then, BSs will look
for handover opportunities by verifying some possible events to trigger or not
the handover procedure. Some of the events described by 3GPP are: [2]:

  * •

A2: Serving cell Reference Signal Received Quality (RSRQ) becomes worse than
threshold

  * •

A3: Neighbour Reference Signal Received Power (RSRP) becomes better than
serving cell

  * •

A4: Neighbour cell RSRQ becomes better than threshold

In the context of handover, there are two well-known mechanisms. These are
A2-A4 handover where two conditions must be satisfied corresponding to the
events A2 and A4: $i)$ to trigger handover the serving cell of the UE must
fall below certain RSRQ serving cell threshold and $ii)$ the handover is only
performed if the difference between the best neighbor and the serving cell
RSRQ is greater than certain predefined neighbor cell offset. The A3 handover
algorithm or “strongest cell handover algorithm” is a simpler approach where
handover is triggered for the UE to the best cell in the measurement report.
When the best cell in terms of RSRP is selected, this value must be greater
than the current serving cell by an hysteresis value and must be maintained
for a time (TTT, time to trigger) in order to avoid ping-pong effect. Handover
using the previous mechanism is illustrated in Fig. 1.

Figure 1: Illustration of A3 handover algorithm.

In this paper, we address the load balancing problem in wireless networks
where we seek to balance resource block utilization and QoS metrics such as
latency and throughput. To do so, we use clipped double Q-learning which
tackles the overestimation issues presented by previous deterministic action
space RL algorithms. We consider a centralized approach where our agent can
choose the best cell individual offset (CIO) per BSs that maximizes the key
performance indicators of the network. Such a centralized learning can be
utilized in Cloud RAN (C-RAN). We compare our results with the classical A3
Handover algorithm and a simplistic resource block utilization-based baseline
handover algorithm. Our results show an improvement in the overall performance
of the network in terms of throughput, latency, jitter, and packet loss ratio
when the proposed algorithm is used.

This paper is organized as follows. Section II presents literature related to
load balancing in wireless networks. In section III, the description of the
system model and the motivation of this paper is introduced. Section IV
provides a description of the proposed scheme. Section V depicts the
performance evaluation and comparison with the baseline algorithms. Finally,
section VI concludes the paper.

##  II Related work

There are several works in the literature related to load balancing in RANs
due its significant impact on saving resources and maximizing the network
performance. In most of the cases, the latter problem is tackled by modifying
the existing handover (HO) strategies either by modifying some HO parameters
or by constantly tracking the network Key Performance Indicators (KPIs). [3,
4].

Recently, machine learning has been used for load balancing as well. In [5],
the authors proposed a 5G handover algorithm based on Q-learning for
optimizing the handover mechanism by the selection of data link beams and
access beams in 5G cellular networks. In [6] the authors presented a
supervised learning solution based on deep learning by considering the
variation of the SINR (Signal to Noise Ratio) to calculate the probability of
Radio Link Failure (RLF) and based on this metric, they performed handover to
the cell that is less prone to experience RLF. In [7] the authors presented an
RL-based mobility load balancing (MLB) algorithm addressed to deal with ultra-
dense networks. The proposal consisted of a two-layer architecture with the
first layer in charge of building small clusters and the second layer where in
each intra-cluster, an MLB algorithm is executed to obtain the optimal HO
parameters by minimizing the resource block utilization. In [8] the authors
presented an RL-based load balancing algorithm for LTE where maximizing the
instantaneous throughput of the overall network is the main objective.

In summary, previous works do not consider end-to-end delay and the channel
quality in the optimization process of HO parameters. Different than other
papers, we consider additional QoS metrics alongside the resource block
utilization seeking to balance both while maximizing the network performance.

##  III System model

We consider a network consisting of a set $\Gamma$ of size $M_{T}$ of base
stations (BS). The network serves a set of $\Psi$ of size $N_{T}$ stationary
mobile users deployed randomly around the BSs.

According the Downlink (DL) bandwidth configuration chosen, $B$ MHz, let us
define $N_{RB}^{DL}$ as the number of resource blocks available to be assigned
by the Media Access Control (MAC) scheduler. The Channel and QoS Aware (CQA)
scheduler is chosen as MAC scheduler. This scheduler assigns resource blocks
by prioritizing users with greater Head of Line (HOL) delay and maximizing MAC
layer throughput[9]. In this paper, we consider one Resource Block Group (RBG)
($RBG=N_{RB}^{DL}/K$) as the smallest resource unit where $K\in[1,2,3,4]$
according the DL bandwidth used. In addition, a centralized approach is
assumed where a central agent is able to monitor the BSs and UEs Key
Performance Indicators (KPIs). This can be conveniently applicable to C-RANs,
as mentioned before. The agent is capable of altering the CIO of each $i$ cell
based on the proposed machine learning algorithm. The CIO of each cell denoted
as $\o_{i}\in[\o_{min},\o_{max}]$ dB and is modified in order to trigger the
handover algorithm among cells. For A3 handover algorithm, if a user is served
by some cell $i$, it will start a handover to cell $j$ request through the
$X2$ interface if the following condition holds:

$RSPR_{j}+\o_{j\rightarrow i}>Hys+RSPR_{i}+\o_{i\rightarrow j}$ (1)

where $RSPR_{j}$ and $RSPR_{i}$ are the measured RSRP values in dB of the
serving cell and the neighbor cell. $Hys$ is the hysteresis value used to
avoid ping-pong scenarios. $\o_{j\rightarrow i}$ and $\o_{i\rightarrow j}$ are
the cell individual offsets. The values are independent for each cell.

##  IV Clipped Double Q-Learning Based Load Balancing

In the proposed approach, an agent observes the environment parameters, such
as CQI, packet delay and resource block utilization per BS. The agent’s
actions consist of modifying the BS’s individual CIO values in order to
maximize the agent’s reward. In the following subsection, we present an
overview of the Clipped Double Q-Learning (CDQL) and then formally define our
solution.

###  IV-A Clipped Double Q-Learning

In this work, we use a state-of-the-art deterministic action space RL
algorithm named Clipped Double Q-Learning. This algorithm is presented in [10]
as part of TD3 (Twin Delayed Deep Deterministic Policy Gradient Algorithm),
which builds on the Deep Deterministic Policy Gradient algorithm (DDPG) [11].
It may be seen as a continuation of the work of [12] with its Double Deep
Q-Learning algorithm (DDQN). The aforementioned algorithms suffered from
Q-value overestimation due the usage of the $argmax$ operator for selecting
the maximum Q-value. CDQL tackles overestimation issues by following the
strategy of having two neural networks that learn at the same time meanwhile
the reward is calculated based on the minimum Q-value of such networks hence,
reducing overestimation. The pseudo-code for CDQL is presented in Algorithm 1.

Figure 2: Learning performance for varying number of UEs and for the CDQL
algorithm.

###  IV-B Action space selection

The actions of our agent are defined as a vector with the CIO values assigned
to each BS. For CDQL the set of actions will be deterministicly predefined by
the possible permutations of the set of predefined values.

For CDQL algorithm, the size of the action space will be
${}^{l}\\!P_{k}=\frac{l!}{(l-k)!}$ where $l$ is the size of the set of the
possible CIO values that can take each BS and $k$ will be equivalent to
$M_{T}$ or the amount of BSs in the network. Thus,

$A(t)=\begin{bmatrix}\o_{1}(t)&\o_{2}(t)&...&\o_{M_{T}}(t)\end{bmatrix}$ (2)

The individual CIO value defined as $\o_{i}$ will be lower and upper bounded
by two predefined values $\o_{min}$ and $\o_{max}$ as described as follows:
$\o_{i}(t)\in[\o_{min},\o_{max}]$.

Initialize $Q_{\theta_{1}}$,$Q_{\theta_{2}}$ and target networks with random
weights, replay buffer $J$, $\epsilon$, $\epsilon_{decay}$ and
$\epsilon_{min}$ ;

foreach  _environment step_ do

Observe state $s_{t}$;

Select $a_{t}\sim\pi(a_{t},s_{t})$ if $\epsilon\geq x\sim\emph{U}(0,1)$
otherwise select random action $a_{t}$ ;

Execute $a_{t}$ and observe next state $s_{t+1}$ and reward
$r_{t}=R(s_{t},a_{t})$;

Store $(s_{t},a_{t},r_{t},s_{t+1})$ in replay buffer $J$;

foreach  _update _ do

sample $e_{t}=(s_{t},a_{t},r_{t},s_{t+1})\sim J$;

Compute $Y_{t}^{CDQL}$:

$Y_{t}^{CDQL}=r_{t}+\gamma
min_{i=1,2}Q_{\theta_{i}}(s_{t+1},argmax_{a^{\prime}}$
$Q_{\theta_{i}}(s_{s+1},a^{\prime}))$;

Perform gradient descent step on
$(Y_{t}^{CDQL}-Q_{\theta_{i}}(s_{t},a_{t}))^{2}$;

Update target networks parameters;

$\theta_{i}^{\prime}\leftarrow\tau*\theta_{i}+(1-\tau)*\theta_{i}^{\prime}$ ;

Update $epsilon$ if $\epsilon>\epsilon_{min}$;

$\epsilon*=\epsilon_{decay}$

end foreach

end foreach

Algorithm 1 Clipped Double Q-Learning

###  IV-C State space selection

The state space is composed by two terms. The first one corresponds to the
number of attached UEs’ ratio in each BS and the second term is the Resource
Block Utilization (RBU) vector. Thus, the state $S(t)$ will be represented by
the concatenation of both metrics as:

$S(t)=\begin{bmatrix}\mathbf{U(t)}&\mathbf{P(t)}\end{bmatrix}$ (3)

where $\mathbf{U}(t)$ corresponds to a $M_{T}$-length vector comprised by the
ratio of UEs attached to each cell $i$.

$\mathbf{U}(t)=\begin{bmatrix}u_{1}(t)&...&u_{M_{T}}(t)\end{bmatrix}$ (4)

where $\Psi_{\Gamma_{i}}$ is the total of UEs attached to cell $i$ and
$u_{i}=\dfrac{\Psi_{\Gamma_{i}}}{N_{T}}$

Additionally, $\mathbf{P}(t)$ corresponds to a $M_{T}$-length vector conformed
by the RBU for each BS at time $t$ as:

$\mathbf{P}(t)=\begin{bmatrix}p_{1}(t)&...&p_{M}(t)\end{bmatrix}$ (5)

Given the fact that resource allocation in each cell is reported every TTI,
which is a smaller time frame compared with the observable time interval of
our agent (1s), we consider the resource block utilization in each TTI as a
discrete random variable. Thus, we can model the resource block utilization
during the observation time $p_{i}$ as the expected value of the resource
block utilization for each TTI $p_{tti}$.

$p_{i}=\mathbb{E}[p_{tti}]$ (6)

###  IV-D Reward

The total reward function is calculated based on Quality of Service (QoS)
parameters and the load of each BS in terms resource block utilization as
follows:

$R_{T}(t)=w_{1}*R_{D}(t)+w_{2}*R_{RB}(t)+w_{3}*R_{CQI}(t)$ (7)

Here $R_{D}\in[-1,1]$ is the reward related to delay constrains,
$R_{RB}\in[-1,1]$ is the reward related to resource allocation usage and
$R_{CQI}\in[-1,1]$ corresponds to a reward measuring the quality of the
modulation used by the users in the network. Finally, $w_{1},w_{2},w_{3}$ are
the weights for each individual reward.

Figure 3: Performance metrics of CDQL, A3 and ReBuHA algorithms. (a)
Throughput, (b) Delay, (c) Jitter and (d) PLR

The first component of the proposed reward can be defined as:

$R_{D}(t)\mathrel{\ensurestackMath{\stackon[1pt]{=}{\scriptstyle\Delta}}}\frac{1}{N_{T}}\sum_{i=1}^{N_{T}}\mathbb{1}\cdot\delta_{i}(t)$
(8)

and,

$\delta_{i}(t)=\begin{cases}-1&\text{if $\exists\Psi_{i}\notin\Psi_{C}$ ,}\\\
\daleth(D_{avg})&\textit{otherwise}\end{cases}$ (9)

where $\Psi_{i}$ represents the $i$th UE and $\Psi_{C}$, $\Psi_{C}\subset\Psi$
is the set of connected UEs to any BS $\in\Gamma$. On the other hand,
$\daleth(D_{avg})$ is a sigmoid function defined as follows:

$\daleth(D_{avg})=1+\dfrac{c}{1+e^{-o(D_{avg}-\mathcal{F})}}$ (10)

Here $c$ establishes the upper bound of the slope, $o$ adjusts the slope of
the sigmoid and $\mathcal{F}=2/3*PDB$ controls the target packet delay. $PDB$
is the Packet Delay Budget which according the type of traffic used in the
network.

The objective of the first term is to reward each UE’s average latency based
on the PDB of a defined packet type. As can be seen, we penalize the cases if
a UE gets disconnected from the network as a result of the load balancing
decision.

The second component of the proposed reward can be defined as:

$R_{RB}(t)\mathrel{\ensurestackMath{\stackon[1pt]{=}{\scriptstyle\Delta}}}1+\dfrac{c}{1+e^{-a(max(P(t))-\mathcal{D})}}$
(11)

where $P(t)$ corresponds to a $M_{T}$-length vector composed of the resource
block utilization of each BS at time t. The term $max(P(t))$ allows selecting
the most loaded BS to either penalize or to reward based on a predefined
threshold defined as $\gamma_{RB}$. As the maximum value of $P(t)$ decreases,
higher reward is obtained. $c$ establishes the upper bound of the slope, $a$
adjusts the slope of the sigmoid as in equation 11 and $\mathcal{D}$ controls
the target data utilization.

Finally, the third component is defined by:

$R_{CQI}(t)\mathrel{\ensurestackMath{\stackon[1pt]{=}{\scriptstyle\Delta}}}\frac{1}{N_{T}}\sum_{i=1}^{N_{T}}\mathbb{1}\cdot\Omega_{i}(t)$
(12)

where,

$\Omega_{i}(t)=\begin{cases}-1&\text{if $\Psi_{i_{CQI}}(t)<6$,}\\\ 0&\text{if
$7\leq\Psi_{i_{CQI}}(t)\leq 9$ ,}\\\ 1&\textit{otherwise}\end{cases}$ (13)

$\Psi_{i_{CQI}}$ corresponds to the CQI measurement of the $i$th UE. The
objective of this term is to reward or penalize each UE’s CQI based on the
fact that higher CQI will be translated in a higher modulation scheme and thus
a more efficient usage of the resource block allocation.

Figure 4: Average BS resource block utilization (RBU) and UE per BS ratio
after convergence of each strategy. (a) 30 UEs, (b) 35 UEs, (c) 40 UEs, (d) 45
UEs, (e) 50 UEs

.

###  IV-E Baseline: Resource Block Utilization based Handover Algorithm
(ReBUHA)

In this work, we consider a baseline algorithm named ReBUHA in addition to the
classical A3 handover algorithm to perform a comparison of our results. ReBUHA
algorithm replaces the A3 handover which mainly relies on received power (see
eq. 1) with a resource allocation awareness method to trigger the handover
procedure.

Data: RBU information of all BSs

Result: Resource allocation based handover

if  _all $p_{i}(t)\in\mathbf{P}(t):p_{i}(t)>\gamma_{RB}$_ then

break;

else

foreach  _$p_{j\neq i}(t)\in\mathbf{P}(t)$ _ do

Check if UE $\Psi_{i}$ is attached to an overloaded BS;

if  _$p_{j}(t) <\gamma_{RB}$ _ then

trigger handover to cell $j\in\Gamma$;

else

break;

end if

end foreach

end if

Algorithm 2 ReBUHA algorithm

As seen in Algorithm 2 $\mathbf{P}(t)$ is $M_{T}$-size vector where each
element corresponds to the ratio of RBs used in the $m^{th}$ BS at the time t.
Every $t=1s$, the algorithm will use $\mathbf{P}(t)$ information and look for
handover opportunities in a centralized way. $\gamma_{RB}$ corresponds to a
resource block utilization threshold defined as the ratio of usage in which
the handover algorithm will look for opportunities for handover.

##  V Performance Evaluation

###  V-A Simulation Setting

TABLE I: Network settings

Parameter | Value
---|---
Inter-site Distance | 720 m
$M_{T}$ | 3
$N_{T}$ | 30,35,40,45,50
Center Frequency | 2 MHz
System Bandwidth | 5 MHz (25 resource blocks)
Pathloss Model | Log Distance Propagation Loss Model
| 95 + 27 $log_{10}(distance[km])$
BS antenna height | 30 m
UE antenna height | 1.5 m
Max Tx power | 20 dbm
MAC scheduler | CQA scheduler
User distribution | Stationary and uniformly distributed
Traffic Model | Conversational video (live streaming) and Poisson
| Packet payload size = 250 Bytes
| Interval = 10 ms
| Packet delay budget = 150 ms
Handover algorithm | A3-event based
| Time to trigger = 8 ms
| Hysteresis = 2 dBm

The simulations are performed by using the discrete network simulator ns-3
[13]. In Table I and Table II we provide the settings utilized in our
simulations and the RL parameters, respectively. Each BS is positioned with an
inter-site distance of 720 meters. Five different scenarios are tested under
the proposed algorithms with 30, 35, 40, 45 and 50 UEs. Each scenario is
designed by distributing a percentage of the total users on the edge of each
cell in a random disc and the rest uniformly allocated throughout the coverage
of the middle BS. We initialize the simulations by attaching all UEs to the BS
that sits in between 2 BSs and then according to the policy of each strategy,
handover is triggered or not. The traffic is a mixture of 20 UEs using CBR,
and the rest following Poisson arrivals. For the case of Poisson traffic, we
use a small payload of 32 bytes with a traffic load of 0.1 Mbps, meanwhile for
the CBR which emulates video traffic we use a larger payload of 250 bytes with
an interval of 10 ms. Simulation results are collected by averaging 15
simulations per scenario. Each simulation consists of 150 episodes with 50
iterations per episode. OpenAI Gym is used as the interface between ns-3 and
our agent [14]. Lastly, we consider each parameter of our reward of equal
importance, thus the weights of equation 7 are equal to 1.

TABLE II: RL environment settings

Parameter | Value
---|---
Number of iterations/episode | 50
Number of episodes | 150
Gym environment step time | 1s
Batch size | 32
CDQL | $\o$ set $\\{-9,-6,-3,0,3,6,9\\}$ dBm
| $w_{1},w_{2},w_{3}=1$
|  $\mathcal{F}=2/3*PDB$ where $PDB=150$ ms, $c=-2$, $o=75$
|  $\gamma_{RB}=0.6$, $a=20$
| Optimizer : Adam $(1e-3)$
| Number of hidden layers $(N_{h})=$2
| Loss function : Huber Loss $(clipdelta=1.0)$
| Update target model type : Polyak averaging
|  $\gamma=0.95$, $\epsilon=1.0$, $\epsilon_{min}=0.001$,
$\epsilon_{decay}=0.995$

###  V-B Simulation Results

To assess the performance of our proposed scheme, we present throughput,
delay, jitter, packet loss ratio (PLR) with 90% confidence interval as well as
the convergence of the machine learning algorithm.

Figure 2 shows the learning performance achieved by CDQL for the simulated
scenarios. The trend shows how the reward value per episode converges in all
cases and the converged reward value becomes lower as the number of UEs
increase. This is due that some of the KPIs that are tracked by our agent are
affected by the increment of the number of users. Figure 3. presents
throughput, end-to-end delay, jitter, and packet loss ratio (PLR). In Figure
3(a) it can be seen that our algorithm achieves an average improvement of 6.1%
and 9.5% in throughput in comparison with the A3 and ReBuHa algorithms,
respectively. Similarly, the other figures (b, c, d and e) show a noticeable
improvement with respect to the baselines with a gain of 49.8% and 52.9% in
terms of delay, 55% and 51% in terms of jitter and 34% and 55.2% in terms of
PLR for the cases of the A3 and ReBuHa algorithms, respectively. Note that,
the delay results of ReBuHa algorithm decreases at 45 UEs in comparison with
A3. This is because, after 45 UEs the algorithm triggers its handover
procedure as the middle BS surpasses the resource block utilization threshold
predefined by the value of $\gamma_{RB}$ . However, metrics such as PLR
continue to worsen reassuring our thesis that in high traffic scenarios a
closed track of QoS metrics is needed. Fig. 4 presents the components of the
vectors $\mathbf{P}$ and $\mathbf{U}$ at the end of the simulated scenarios.
It can be seen how the A3 handover algorithm does not trigger in any of the
scenarios by keeping all the UEs attached in the middle BS. The latter occurs
because none of the UEs comply with the event trigger condition of the A3
algorithm. Furthermore, it is noticeable that in Figure 4 (d) the ReBuHa
algorithm starts distributing the load in agreement with the behavior
described in the previous figure. For our proposed algorithm it is observable
that it distributes the UEs over the BSs based on choosing the ”best” CIO
value per BS that will maximize our agent’s reward objective function. Note
that the proposed algorithm for (d) and (e) does not meet the target resource
block utilization goal which is established by $\gamma_{RB}$. This behavior
can be explained based on the needed minimization of not only the resource
block utilization in the network but also improvement of QoS metrics as well.

Figure 5: Learning performance of CDQL for 30 UEs and 0%, 10%, 20% and 30% of
mobile users with UE’s speed of 20 m/s.  Figure 6: Different performance
metrics of CDQL and A3 algorithms for different mobility percentages and UE’s
speed of 20 m/s. (a) Throughput, (b) Delay, (c) Jitter and (d) PLR

Additionally, we present the performance of our scheme under mobility. We
consider 30 UEs where 10%, 20% and 30% of the total UEs in the network are
mobile and they use random walk with a speed of 20 m/s. Figure 5 shows that,
similar to the non-mobility scenario (represented as 0%), our scheme converges
for mobile scenarios in a similar fashion. Furthermore, we show in Figure 6
throughput, end-to-end delay, jitter, and PLR for non-mobile and mobile
scenarios. Note that we only presented A3 Algorithm and CDQL since for such
number of users the ReBuHa algorithm behave identical as the A3. In terms of
throughput and PLR both CDQL and A3 algorithm perform similar. For delay, the
proposed CDQL offers slightly lower latency than A3 algorithm. More
specifically, our scheme is able to achieve an improvement of 64% respecting
delay with no considerable difference in the other KPIs. Finally, it is worth
to mention that our scheme shows its real potential when the number of users
increases, in other words, when the scarcity of resources increases. We chose
the A3 best performance scenario (30 UEs) to show the steady behavior of our
scheme.

##  VI Conclusions

In this paper, we presented a Clipped Double Q-Learning strategy that performs
load balancing with awareness of QoS metrics. As main difference from previous
works, our RL method uses a state-of-the-art algorithm and a QoS-aware load
balancing approach enhancing the overall KPI metrics such as throughput,
delay, jitter, and packet delivery ratio. We compared our proposed scheme with
two baselines: the traditional A3 handover algorithm and a resource block
utilization based handover scheme, named ReBuHa. The results showed an average
improvement up to 6.1% and 9.5% in terms of throughput, 49.8% and 52.9% in
terms of delay, 55% and 51% in terms of jitter and 34% and 55.2% in terms of
PLR, in comparison to the A3 and ReBuHa algorithms, respectively.
Additionally, we evaluated the performance of our scheme under mobility. The
results revealed the importance of performing load balancing while maintaining
latency and CQI metrics.

##  VII Acknowledgment

This research is supported by the 5G ENCQOR program and Ciena.

## References

  * [1] A. Lutu, D. Perino, M. Bagnulo, E. Frias-Martinez, and J. Khangosstar, “A Characterization of the COVID-19 Pandemic Impact on a Mobile Network Operator Traffic,” in _Proceedings of the ACM SIGCOMM Internet Measurement Conference, IMC_ , 2020.
  * [2] 3GPP, “TS 36.331: Radio Resource Control (RRC); Protocol specification (Release 15),” _3Gpp_ , 2018.
  * [3] R. Ahmad, E. A. Sundararajan, N. E. Othman, and M. Ismail, “Handover in LTE-advanced wireless networks: state of art and survey of decision algorithm,” _Telecommunication Systems_ , 2017.
  * [4] M. Tayyab, X. Gelabert, and R. Jantti, “A Survey on Handover Management: From LTE to NR,” 2019.
  * [5] V. Yajnanarayana, H. Ryden, and L. Hevizi, “5G Handover using Reinforcement Learning,” in _2020 IEEE 3rd 5G World Forum, 5GWF 2020 - Conference Proceedings_ , 2020.
  * [6] Z.-H. Huang, Y.-L. Hsu, P.-K. Chang, and M.-J. Tsai, “Efficient Handover Algorithm in 5G Networks using Deep Learning,” in _GLOBECOM 2020 - 2020 IEEE Global Communications Conference_. IEEE, dec 2020, pp. 1–6. [Online]. Available: https://ieeexplore.ieee.org/document/9322618/
  * [7] Y. Xu, W. Xu, Z. Wang, J. Lin, and S. Cui, “Load Balancing for Ultradense Networks: A Deep Reinforcement Learning-Based Approach,” _IEEE Internet of Things Journal_ , 2019.
  * [8] K. Attiah, K. Banawan, A. Gaber, A. Elezabi, K. Seddik, Y. Gadallah, and K. Abdullah, “Load Balancing in Cellular Networks: A Reinforcement Learning Approach,” in _2020 IEEE 17th Annual Consumer Communications and Networking Conference, CCNC 2020_ , 2020.
  * [9] B. Bojovic and N. Baldo, “A new channel and QoS aware scheduler to enhance the capacity of voice over LTE systems,” in _2014 IEEE 11th International Multi-Conference on Systems, Signals and Devices, SSD 2014_ , 2014\.
  * [10] S. Fujimoto, H. Van Hoof, and D. Meger, “Addressing Function Approximation Error in Actor-Critic Methods,” in _35th International Conference on Machine Learning, ICML 2018_ , 2018.
  * [11] T. P. Lillicrap, J. J. Hunt, A. Pritzel, N. Heess, T. Erez, Y. Tassa, D. Silver, and D. Wierstra, “Continuous control with deep reinforcement learning,” in _4th International Conference on Learning Representations, ICLR 2016 - Conference Track Proceedings_ , 2016.
  * [12] H. Van Hasselt, A. Guez, and D. Silver, “Deep reinforcement learning with double Q-Learning,” in _30th AAAI Conference on Artificial Intelligence, AAAI 2016_ , 2016.
  * [13] N. Baldo, M. Miozzo, M. Requena-Esteso, and J. Nin-Guerrero, “An open source product-oriented LTE network simulator based on ns-3,” in _MSWiM’11 - Proceedings of the 14th ACM International Conference on Modeling, Analysis, and Simulation of Wireless and Mobile Systems_ , 2011.
  * [14] P. Gawłowicz and A. Zubow, “ns3-gym: Extending OpenAI gym for networking research,” 2018.