# Widget to show equilibrium profile of an ice sheet


```python
import matplotlib.pyplot as plt
import numpy as np
from ipywidgets import interactive
import ipywidgets as widgets
```

Let's derive the equation first. Let's assume a horizontal ice sheet bed at $z=0$. In this case the basal shear stress, $\tau_b$, balances the driving stress (Cuffey & Paterson eq. 8.132):

\begin{equation}
\tau_b = -\rho_i g h \frac{dh}{dx},
\end{equation}

where $\rho_i$ is the ice density, $g$ is the acceleration due to gravity, $h$ is the ice thickness, and $x$ is the distance from the center of the ice sheet. If we assume perfect plasticity, the ice-sheet profile will adjust so that $\tau_b$ equals the yield stress, $\tau_0$ at every point:

\begin{equation}
\tau_0 = -\rho_i g h \frac{dh}{dx}.
\end{equation}

We ultimately want an equation for $h(x)$, so let's rearrange and integrate:

\begin{eqnarray}
\tau_0 &=& -\rho_i g h \frac{dh}{dx} \\
hdh &=&  -\frac{\tau_0}{\rho_i g}dx \\
\frac{h^2}{2} &=& -\frac{\tau_0}{\rho_i g}x + C
\end{eqnarray}

We need a boundary condition to solve for the integration constant $C$. In this case, let's say the ice sheet has length $L$, meaning at $x=L$, $h=0$.

\begin{eqnarray}
\frac{h^2}{2} &=& -\frac{\tau_0}{\rho_i g}x + C \\
\frac{0^2}{2} &=& -\frac{\tau_0}{\rho_i g}L + C \\
0 &=& -\frac{\tau_0}{\rho_i g}L + C  \\
C &=& \frac{\tau_0}{\rho_i g}L
\end{eqnarray}

Plug that constant back into the equation:
\begin{eqnarray}
\frac{h^2}{2} &=& -\frac{\tau_0}{\rho_i g}x + \frac{\tau_0}{\rho_i g}L \\
h^2 &=& \frac{2\tau_0}{\rho_i g}(L-x)\\
h(x) &=& \sqrt{\frac{2\tau_0}{\rho_i g}(L-x)}
\end{eqnarray}


```python
rho = 918 # ice density in kg/m3
g = 9.8 # acceleration due to gravity in m/s2

# define a function that takes an ice sheet length (in km) and yield stress (in kPa)
# and it calculates and plots the equilibrium profile 
def h(length, yieldstress):
    xvec = np.arange(0,length) # vector of x locations
    hvec = np.sqrt(2*yieldstress*1000/(rho*g) * (length*1000 - xvec*1000)) # calc height profile
    plt.fill(np.concatenate((xvec,np.flip(xvec))),
             np.concatenate((hvec,xvec*0)))
    plt.xlim(0,5000)
    plt.ylim(0,14000)
    plt.xlabel('distance [km]')
    plt.ylabel('height [m]')
    plt.show()

# make ice-sheet length slider bar
l = widgets.IntSlider(description="length (km)", 
                      min = 0, max = 5000, value = 3000, step=100)
# make yield stress slider bar
tau = widgets.IntSlider(description="\N{greek small letter tau}$_o$ (kPa)", 
                        min = 0, max = 150, value = 100, step=5)

#create and show the widget
img = interactive(h, length=l, yieldstress=tau)
display(img)
```


    interactive(children=(IntSlider(value=3000, description='length (km)', max=5000, step=100), IntSlider(value=10…



```python

```

f0099c67d872ac96d79435a436ff94395b53fc06
ipynb
Jupyter Notebook
glaciers/eq_profile.ipynb
mrsiegfried/class_widgets
8d2b54cc87406cf759c4179fb5a1a227d50cc24c
BSD-2-Clause
glaciers/eq_profile.ipynb
mrsiegfried/class_widgets
8d2b54cc87406cf759c4179fb5a1a227d50cc24c
BSD-2-Clause
glaciers/eq_profile.ipynb
mrsiegfried/class_widgets
8d2b54cc87406cf759c4179fb5a1a227d50cc24c
BSD-2-Clause
Qwen/Qwen-72B
1. YES
2. YES