Skip to main content
Geosciences LibreTexts

3.5: Isostasy Jupyter Notebook

  • Page ID
    11678
  • Instructions

    • To use this page, read the text between each of the python code windows, then press RUN to execute the code in the box.
    • It is necessary to run each of the boxes in order.
    • Don't click restart in the cells unless you come back to the top of the page and start over.
    • If you modify the code, this modification will not remain after you leave or refresh this page.
    • You must run simulations/cells that require user input to completion. Make sure you see the Message informing you a cell is complete before continuing to the next cell. Here's why you must do this:
      • These cells may not be re-run until the simulation is completed, otherwise no output will be generated.
      • All following cells will not run until the user input cell is completed.
    • For Dropdown Menus: You do not need to rerun the cells with the menu to change your selection. However, you must rerun all of the following cells to implement the change.

    Be patient, some times it can take 1-2 minutes for the juypter kernel to start.

    Interactive examples of the principle of isostasy starting with an iceberg, and then moving to more complicated examples of the crust and the mantle, and a crust with multiple thicknesses.

    Key Questions: Consider these as you work your way through this page.

    Iceberg in Water Example

    • For a given thickness, does the iceberg extend father down into the water, or extend higher above the water?
    • How does the height of the iceberg relative to its depth change between thicknesses of 10 m, 60 m, and 100 m? (Does the height/depth ratio increase, decrease, or remain the same?) What combination of parameters in the equation determines this ratio? 

    Crust in Mantle "Ocean" Example

    • Run the simulation for a few different crustal thicknesses (1, 20, and 40 km): does the height/depth ratio increase, decrease, or remain the same? 
    • Compare how the height/depth ratio depends on layer thickness for this example with the iceberg example. Why are they the same or different?

    Continent Isostasy Example

    • As the oceanic thickness increases, does the elevation of the continent increase, decrease, or remain the same? What combination of parameters controls this relationship?
    • Run the simulation at a few different oceanic crustal thicknesses (1, 5, and 20 km), does the height/depth ratio of the continental crust increase, decrease, or remain the same? In other words, how does the thickness of the oceanic crust affect the height/depth ratio of the continental crust?
    • What happens when you use an oceanic crustal thickness above 22 km? Why does this happen?

     

    First, let's import the required packages and define constants for this problem. rho_w is the density of seawater at the surface and rho_i is the density of ice. 

    import numpy as np
    import matplotlib.pyplot as plt
    
    rho_w = 1028 #kg/m^3
    rho_i = 917 #kg/m^3
    
    print('Done, Move on to the next cell.')

    Iceberg in Water

    We can use Archimedes' Principle to calculate the height and depth of the iceberg by setting the weight of the iceberg equal to the weight of the displaced water. Let's define that function here first, then run a simulation to see what happens as the thickness of the iceberg changes. The T in the function stands for the thickness of the iceberg. x is the horizontal distance range that will be considered, for simplicity the width of the iceberg will remain constant at 20 m. After the function, there is a plotting routine that will update the plot every time you enter a new thickness. After you run the cell, enter a thickness between 1 - 100 m and consider the key questions.

    def iceberg_isos(rho_w,rho_i,T):
        depth = (rho_i / rho_w) * T
        height = T - depth
        return -depth, height
    
    handle = display(None, display_id=True)
    x = np.arange(-25,25) # Horizontal distance in the Graph
    
    bbox = dict(boxstyle ="round", fc = '1') 
    fig, ax = plt.subplots(figsize = (9,7))
    ax.hlines(0,-25,-10,color='blue')
    ax.hlines(0,10,25,color='blue')
    ax.fill_betweenx((0,-100),-25,-10,color = 'blue')
    ax.fill_betweenx((0,-100),10,25,color = 'blue')
    ax.set_xlim(-25,25)
    ax.set_title('Iceberg Isostasy')
    ax.set_ylabel('Elevation (m)')
    ax.set_xlabel('Distance from the Center of the Iceberg (m)')
    ax.set_ylim(-100,15)
    ax.annotate('Air',(18,3))
    print('Enter a Iceberg Thickness between 1 and 100 m: ')
    T = float(input())
    while T< 1 or T > 100:
        print('Invalid Thickness, please enter a number between 1 and 100 km')
        T = float(input())
    while T >= 1 and T <= 100:
        print('Now Plotting:', T, 'm')
        depth, height = iceberg_isos(rho_w,rho_i,T)
        iceberg_bottom = ax.hlines(depth,-10,10,color = 'black')
        iceberg_top = ax.hlines(height,-10,10,color = 'black')
        iceberg_left = ax.vlines(-10,depth,height,color='black')
        iceberg_right = ax.vlines(10,depth,height,color='black')
        iceberg_fill = ax.fill_betweenx((depth,height),-10,10,color='lightsteelblue',label = 'Iceberg')
        below_iceberg = ax.fill_betweenx((-100,depth),-10,10,color='blue',label = 'Water/Ocean')
        depth_label = ax.annotate('Depth = '+'{:.2f}'.format(-depth)+' m',(-5,depth/2),bbox=bbox)
        height_label = ax.annotate('Height = '+'{:.2f}'.format(height)+' m',(-23,8),bbox=bbox)
        ax.legend(bbox_to_anchor=(1.3, 0.5))
        handle.update(fig)
        iceberg_bottom.remove()
        iceberg_top.remove()
        iceberg_left.remove()
        iceberg_right.remove()
        iceberg_fill.remove()
        below_iceberg.remove()
        depth_label.remove()
        height_label.remove()
        print('Enter a number outside of the provided range to exit the simulation when ready.')
        print('Enter a Iceberg Thickness between 1 and 100 m to continue: ')
        T = float(input())
    plt.close()
    print('The Simulation is Complete, Consider the Key Questions')

    Crust in a Mantle "Ocean"

    Now, let's consider an example more closely related to geophysics. Here, we replace the water with mantle rock, treating the mantle like an ocean. We then replace the ice with crustal rock. (Think of the first crust to crystallize in the magma ocean that covered the early Earth). To be more realistic with the appropriate length scales, we also move from meters to kilometers in this example. In this scenario, the math is exactly the same as above, so we will reuse our iceberg_isos function from above, replacing rho_i with rho_crust, and rho_w with rho_mantle. 

    rho_crust = 2850 # kg/m^3
    rho_mantle = 4000 # kg/m^3
    
    handle = display(None, display_id=True)
    x = np.arange(-25,25) # Horizontal distance in the Graph
    
    
    bbox = dict(boxstyle ="round", fc = '1') 
    fig, ax = plt.subplots(figsize = (9,7))
    ax.hlines(0,-25,-10,color='red')
    ax.hlines(0,10,25,color='red')
    ax.fill_betweenx((0,-50),-25,-10,color = 'red')
    ax.fill_betweenx((0,-50),10,25,color = 'red')
    ax.set_xlim(-25,25)
    ax.set_title('Crust-Mantle Ocean Isostasy')
    ax.set_ylabel('Elevation (km)')
    ax.set_xlabel('Distance from the Center of the Crustal Block (km)')
    ax.set_ylim(-40,20)
    ax.annotate('Air',(18,3))
    print('Enter a Crustal Thickness between 1 and 40 km: ')
    T = float(input())
    while T< 1 or T > 40:
        print('Invalid Thickness, please enter a number between 1 and 40 km')
        T = float(input())
    while T >= 1 and T <= 40:
        print('Now Plotting:',T,'km')
        T *= 1000
        depth, height = iceberg_isos(rho_mantle,rho_crust,T)
        height /= 1000
        depth /= 1000
        crust_bottom = ax.hlines(depth,-10,10,color = 'black')
        crust_top = ax.hlines(height,-10,10,color = 'black')
        crust_left = ax.vlines(-10,depth,height,color='black')
        crust_right = ax.vlines(10,depth,height,color='black')
        crust_fill = ax.fill_betweenx((depth,height),-10,10,color='tan',label = 'Crust')
        below_crust = ax.fill_betweenx((-50,depth),-10,10,color='red',label = 'Mantle')
        depth_label = ax.annotate('Depth = '+'{:.2f}'.format(-depth)+' km',(-5,depth/2),bbox=bbox)
        height_label = ax.annotate('Height = '+'{:.2f}'.format(height)+' km',(-23,15),bbox=bbox)
        ax.legend(bbox_to_anchor=(1.3, 0.5))
        handle.update(fig)
        crust_bottom.remove()
        crust_top.remove()
        crust_left.remove()
        crust_right.remove()
        crust_fill.remove()
        below_crust.remove()
        depth_label.remove()
        height_label.remove()
        print('Enter a number outside of the provided range to exit the simulation when ready.')
        print('Enter a Crustal Thickness between 1 and 40 km: ')
        T = float(input())
    plt.close()
    print('The Simulation is Complete, Consider the Key Questions')

    Continent Isostasy

    Now let's build on your understanding from the two above examples to move to a complex example that resembles a more realistic situation. This time we will consider a setup that has air, ocean, lithospheric mantle, and oceanic crust in the outer columns. We will fix the continental crust thickness and the thickness of the ocean while varying the thickness of the oceanic crust. The first portion of this cell defines the known densities and thicknesses. The next part defines a new isostasy function that will take in all of the known thicknesses and densities, and calculates the two unknown thicknesses: air and mantle. The equations within this function can be found in section 3.6 Isostasy. The rest of the cell is simply a plotting routine. Go ahead and run the cell, and then answer the key questions relating to this example.

    Note: in the present-day earth the thickness of oceanic crust formed at a mid-ocean ridge is remarkably constant, and is about 7.5 km, with a range of 1-10 km. However, this thickness is determined by melting processes and so it could be different in the past or future when the mantle was/is hotter or colder. It could also be different on different planet with a different bulk composition.

    The average elevation of the continents above sea level is about 800m, or 0.8km. You should find a number close to that if you use the average oceanic thickness of 7.5km.

    rho_air = 1 # kg/m^3
    rho_mantle = 3300 # kg/m^3
    rho_crust = 2900 # kg/m^3
    rho_oc = 3000 # kg/m^3
    h_c = 35 # km
    h_w = 4 # km
    rho_water = 1030 # kg/m^3
    
    def isostasy(rho_a,rho_m,rho_c,rho_oc,h_c,h_oc,rho_w,h_w):
        h_a = ((h_c*(rho_m - rho_c)) - (h_w*(rho_m - rho_w)) - (h_oc*(rho_m - rho_oc))) / (rho_m - rho_a)
        h_m = h_c - h_oc - h_a - h_w
        return h_a, h_m
    
    print('Please enter the thickness of the oceanic crust between 1 and 40 km')
    h_oc = float(input())
    while h_oc < 1 or h_oc > 40:
        print('Invalid thickness, enter a number between 1 and 40 km')
        h_oc = float(input())
    
    
    handle = display(None, display_id=True)
    fig, ax = plt.subplots(figsize = (9,7))
    ax.set_xlim(-25,25)
    ax.set_ylim(-50,20)
    ax.annotate('Air',(18,14))
    ax.set_title('Continental Isostasy')
    ax.set_xlabel('Distance from the Center of the Continent (km)')
    ax.set_ylabel('Elevation (km)')
    while h_oc >= 1 and h_oc <= 40:
        print('Now plotting a oceanic crust thickness of',h_oc,'km')
        h_a,h_m, = isostasy(rho_air,rho_mantle,rho_crust,rho_oc,h_c,h_oc,rho_water,h_w)
        center_crust = ax.fill_betweenx((h_a-h_c,h_a),-10,10,color='tan',label = 'Crust')
        water_right = ax.fill_betweenx((0,-h_w),10,25,color='blue',label = 'Water/Ocean')
        water_left = ax.fill_betweenx((0,-h_w),-25,-10,color='blue')
        o_crust_left = ax.fill_betweenx((-h_w,-h_w-h_oc),-25,-10,color='tan')
        o_crust_right = ax.fill_betweenx((-h_w,-h_w-h_oc),10,25,color='tan')
        mantle_left = ax.fill_betweenx((-h_w-h_oc,-1000),-25,-10,color='red',label = 'Lithospheric Mantle')
        mantle_right = ax.fill_betweenx((-h_w-h_oc,-1000),10,25,color='red')
        mantle_center = ax.fill_betweenx((h_a-h_c,-1000),-10,10,color='red')
        outline_top = ax.hlines(h_a,-10,10,color='black')
        outline_bottom = ax.hlines(h_a-h_c,-10,10)
        outline_right = ax.vlines(-10,h_a-h_c,h_a)
        outline_left = ax.vlines(10,h_a-h_c,h_a)
        outline_left_oc = ax.hlines(-h_w,-25,-10)
        outline_right_oc = ax.hlines(-h_w,10,25)
        outline_bot_left_oc = ax.hlines(-h_w-h_oc,-25,-10)
        outline_bot_right_oc = ax.hlines(-h_w-h_oc,10,25)
        outline_oc_outer_right = ax.vlines(25,-h_w,-h_w-h_oc)
        outline_oc_outer_left = ax.vlines(-25,-h_w,-h_w-h_oc)
        if h_a < 0:
          top_ocean = ax.fill_betweenx((0,h_a),-10,10,color='blue')
          right_outline = ax.vlines(10,h_a-h_c,-h_w-h_oc)
          left_outline = ax.vlines(-10,h_a-h_c,-h_w-h_oc)
        depth_label = ax.annotate('Depth = '+'{:.2f}'.format(-(h_a-h_c))+' km',(-5,(h_a-h_c)/2),bbox=bbox)
        height_label = ax.annotate('Height = '+'{:.2f}'.format(h_a)+' km',(-23,15),bbox=bbox)
        thickness_label = ax.annotate('h_oc = '+'{:.2f}'.format(h_oc)+' km',(-23,(-h_w-h_oc/2)),bbox=bbox)
        ax.legend(bbox_to_anchor=(1.3, 0.5))
        handle.update(fig)
        if h_a < 0:
          top_ocean.remove()
          right_outline.remove()
          left_outline.remove()
        center_crust.remove()
        water_right.remove()
        water_left.remove()
        o_crust_left.remove()
        o_crust_right.remove()
        mantle_left.remove()
        mantle_right.remove()
        mantle_center.remove()
        outline_top.remove()
        outline_bottom.remove()
        outline_right.remove()
        outline_left.remove()
        outline_left_oc.remove()
        outline_right_oc.remove()
        outline_bot_left_oc.remove()
        outline_bot_right_oc.remove()
        outline_oc_outer_right.remove()
        outline_oc_outer_left.remove()
        depth_label.remove()
        height_label.remove()
        thickness_label.remove()
        print('You may exit the simulation when ready by entering a number outside the range')
        print('Please enter the thickness of the oceanic crust between 1 and 40 km')
        h_oc = float(input())
    plt.close()
    
    print('The simulation is complete. Consider the key questions.')

    Key Questions: Consider these as you work your way through this page.

    Iceberg in Water Example

    • For a given thickness, does the iceberg extend father down into the water, or extend higher above the water?
    • How does the height of the iceberg relative to its depth change between thicknesses of 10 m, 60 m, and 100 m? (Does the height/depth ratio increase, decrease, or remain the same?) What combination of parameters in the equation determines this ratio? 

    Crust in Mantle "Ocean" Example

    • Run the simulation for a few different crustal thicknesses (1, 20, and 40 km): does the height/depth ratio increase, decrease, or remain the same? 
    • Compare how the height/depth ratio depends on layer thickness for this example with the iceberg example. Why are they the same or different? 

    Continent Isostasy Example

    • As the oceanic thickness increases, does the elevation of the continent increase, decrease, or remain the same? What combination of parameters controls this relationship?
    • Run the simulation at a few different oceanic crustal thicknesses (1, 5, and 20 km), does the height/depth ratio of the continental crust increase, decrease, or remain the same? In other words, how does the thickness of the oceanic crust affect the height/depth ratio of the continental crust?
    • What happens when you use an oceanic crustal thickness above 22 km? Why does this happen?
    • Was this article helpful?