<jupyter_start><jupyter_text>> **Note**: This is a version of the practical without the code. I will use this during class and fill in the code live. I encourage students to type along with me since it can help with understanding and retention of the content. You will always have the full version as a reference. Cooling of the oceanic lithosphere**Goal:** Implement the half-space cooling model and compare its predictions to observations of heat flow, lithosphere age, and bathymetry.The **half-space cooling model** simulates the cooling of the oceanic lithosphere as it moves away from the mid-ocean ridges where it's formed (due to  seafloor spreading).Figure caption: Sketch of the half-space cooling model setup. The lithosphere-asthenosphere boundary (LAB) is marked by an isotherm. The ridge is at the origin and the plates spread symmetrically with constant velocity $u$.The cooling happens by vertical condution of heat to the water layer and can be modelled by the 1D diffusion equation:$$\dfrac{\partial^2 T}{\partial z^2}= \dfrac{1}{\alpha} \dfrac{\partial T}{\partial t}$$The solution to this equation, given the starting condition of the temperature being equal to the asthenosphere temperature $T(z, t=0) = T_a$ and boundary condition of the temperature at the surface being constant $T(z=0, t) = T_0$, is:$$T(z, t) = T_0 + (T_a - T_0)\ \text{erf}\left(\dfrac{z}{2\sqrt{\alpha t}}\right)$$in which $\alpha$ is the thermal diffusivity of the lithosphere, $t$ is the age of the lithosphere, $z$ is the depth, and $\text{erf}$ is the error function.From this equation, we can apply Fourier's law to get **heat flow at the surface**:$$\ q_0(t) = k \dfrac{T_a - T_0}{\sqrt{\pi \alpha t}}$$where $k$ is the thermal conductivity of the lithosphere.We can also apply the conditions for isostatic equilibrium (since we know that the oceanic lithosphere is largely in equilibrium from gravity disturbance data) to predict **bathymetry** from the half-space cooling model:$$w(t) = w_r + \dfrac{2 \rho_a \alpha_V (T_a - T_0)}{\rho_a - \rho_w}\sqrt{\dfrac{\alpha t}{\pi}}$$in which $w_r$ is the bathymetric depth at the ridge, $\rho_a$ is the density of the asthenosphere, $\rho_w$ is the density of water, and $\alpha_V$ is the coefficient of thermal expansion of the lithosphere.**Our job now is to code the equations for bathymetry and surface heat flow in Python**.  The data**Download:** [global-geophysical-data.nc](https://github.com/leouieda/lithosphere/raw/main/data/global-geophysical-data.nc) (place it in the **same folder** as this notebook)The data grids are stored in a [netCDF](https://en.wikipedia.org/wiki/NetCDF) file, which is a very common standard for storing and distributing gridded data. It's used throughout the Earth, ocean, and climate sciences and can be read an manipulated with most software platforms. It contains global grids of gravity, topography, and heat flow. All grids were interpolated onto the same resolution of 0.2 degrees (a compromise between resolution and size of the file) so that we can easily compare and use all of the data together.* **Gravity**: generated from the [EIGEN-6C4 spherical harmonic model](https://doi.org/10.5880/icgem.2015.1) through the [ICGEM webservice](http://icgem.gfz-potsdam.de/home). It includes two grids: the gravity disturbance calculated with the WGS84 ellipsoid and the Bouguer gravity disturbance calculated using a topography density of 2670 kg/m³. The data for both grids was generated on top of the Earth's surface (the water surface in the oceans and topography on land).* **Topography and bathymetry**: a downsampled version of [ETOPO1](https://doi.org/10.7289/V5C8276M) and includes a version smoothed using a Gaussian filter with 1 degree width (~100 km).* **Heat flow**: an upsampled version of the compilation by [Lucazeau (2019)](https://doi.org/10.1029/2019GC008389) (originally 0.5 degree resolution).* **Lithosphere age (oceanic)**: a downsampled version of the grid by [M<NAME> al. (2008)](https://doi.org/10.1029/2007GC001743) (originally 6 arc-minutes).  Import the required librariesLoad the required software to load the data, make maps, and perform calculations.<jupyter_code># For arrays and linear algebra
import numpy as np
# To make figures and maps
import matplotlib.pyplot as plt
# To load and manipulate grids
import xarray as xr<jupyter_output><empty_output><jupyter_text>Load the data gridsOnce again, we'll use xarray to load the data from the netCDF file. For this practical, we'll use a subsection of the grid that overlaps with a mid-ocean ridge in the South Pacific (a large section around the island of [Rapa Nui](https://en.wikipedia.org/wiki/Easter_Island)). We can now plot the observed heat flow, bathymetry, and lithospheric age to see what they look like. The main mid-ocean ridge is cutting across the area from Southwest to Northeast, marked by several transform faults. You may also notice that there is a separate smaller ridge coming from the East and (almost) connecting with the main ridge to form a triple junction. This is clearer in the heat flow map below. Another useful visualization is to make a point plot of lithosphere age versus bathymetry and heat flow. This is the relationship that we'll be predicting from our half-space cooling model so let's plot it and see what the data look like. Heat flow decays very rapidily with age while the bathymetry decays more slowly. There is also significant spread in the bathymetry, which is expected because of all the seamounts and small islands that dot the oceans. These features aren't formed from the processes at mid-ocean ridges and aren't supported by isostatic equilibrium so we don't expect them to comform to our model. To remove these effects, we can downsample the grids (reduce the number of points) by taking the median value of the data in blocks of a certain size.> **Important**: This is the recommended way of downsampling gridded data. You **cannot just take every N points** from the grid. This would lead to aliasing effects in the final product. Even this is a bit of an approximation since it doesn't take into account variations is distance between points due to the latitude. But since this is a small area the effects will be negligible. Now we can plot the downsampled version of the data and get a much cleaner trend in both heat flow and bathymetry.  Heat flow from half-space coolingNow we can create a function that calculates heat flow at the surface of lithosphere (the bottom the ocean) predicted by the half-space cooling model. As a reminder, this is the equation that we want to calculate:$$q_0(t) = k \dfrac{T_a - T_0}{\sqrt{\pi \alpha t}}$$An important thing to take into account are the **units of the data and parameters**. In this case, we'll use the following as inputs and outputs:| Name | Variable | Unit ||--:|:--:|:--|| Surface heat flow | $$q_0$$ | $$\frac{mW}{m^2}$$ | | Thermal conductivity | $$k$$ | $$\frac{W}{m \cdot K}$$ | | Thermal diffusivity | $$\alpha$$ | $$\frac{mm^2}{s}$$ | | Age of the lithosphere | $$t$$ | $$myr = 10^6\ \text{years}$$ | | Asthenosphere temperature | $$T_a$$ | $$K$$ | | Surface temperature | $$T_0$$ | $$K$$ | This means that our function is going to have to do some unit conversions in order to match units of all inputs and return heat flow in $\frac{mW}{m^2}$.<jupyter_code>def halfspace_heat_flow(
    age, 
    surface_temperature, 
    asthenosphere_temperature, 
    diffusivity, 
    conductivity,
):
    """
    Heat flow from the half space cooling model.
    """
   



    return<jupyter_output><empty_output><jupyter_text>Now we can make a prediction of heat flow and compare it against our data. In this case, let's assume:|Input parameters||--||$$T_0 = 273\ K$$||$$T_a = 1600\ K$$||$$\alpha = 1\ \frac{mm^2}{s}$$||$$k = 3.3\ \frac{W}{m \cdot K}$$|  Bathymetry from half-space coolingNow let's predict bathymetry and see how well it fits our data. Here is the equation for bathymetry again:$$w(t) = w_r + \dfrac{2 \rho_a \alpha_V (T_a - T_0)}{\rho_a - \rho_w} \sqrt{\dfrac{\alpha t}{\pi}}$$In this case, the units we'll consider are:| Name | Variable | Unit ||--:|:--:|:--|| Bathymetry | $$w$$ | $$m$$ | | Ridge depth | $$w_r$$ | $$m$$ | | Density of water | $$\rho_w$$ | $$\frac{kg}{m^3}$$ | | Density of asthenosphere | $$\rho_a$$ | $$\frac{kg}{m^3}$$ | | Thermal diffusivity | $$\alpha$$ | $$\frac{mm^2}{s}$$ | | Age of the lithosphere | $$t$$ | $$myr = 10^6\ \text{years}$$ | | Asthenosphere temperature | $$T_a$$ | $$K$$ | | Surface temperature | $$T_0$$ | $$K$$ | | Thermal expansion coefficient | $$\alpha_V$$ | $$K^{-1}$$ |  ---- **Your turn!**In groups:1. Finish the code for the `halfspace_bathymetry` function. 1. Use your function to predict observations using the input parameters in the table below.1. Plot your model predictions against the observed data.|Input parameters||--||$$w_r = 2500\ m$$||$$\rho_w = 1000\ \frac{kg}{m^3}$$||$$\rho_a = 3300\ \frac{kg}{m^3}$$||$$T_0 = 273\ K$$||$$T_a = 1600\ K$$||$$\alpha = 1\ \frac{mm^2}{s}$$||$$\alpha_V = 3\cdot 10^{-5}\ K^{-1}$$|Discuss the following points based on the results you've obtained so far:* How well does the model fit both observations?* Is the fit equal across all lithosphere ages?* What does this say about the model and the way that the lithosphere changes with age?Take notes and prepare to share a summary of your discussion with the class.As a bonus (if you have time), select data from a different mid-ocean ridge system and check if the results are compatible. In other words, are these results universal or specific to this spreading center in the Pacific? For example, these coordinates will get a slice of the Atlantic: longitude between -30 and 5 and latitude between -50 and 0.<jupyter_code>def halfspace_bathymetry(
    age, 
    ridge_depth, 
    asthenosphere_density, 
    water_density, 
    surface_temperature, 
    asthenosphere_temperature, 
    thermal_expansion, 
    diffusivity,
):
    """
    Bathymetry from the half-space cooling model.
    """
    # Get everything into compatible units first
    
 

    
    bathymetry = (
        
        
        
    )
    # -1 because the equation gives us thickness of the water layer
    return -1 * bathymetry<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10353559949352294, 'fraction_numerical': 0.0161683062238239, 'mean_word_length': 4.345132743362832, 'pattern_counts': {'":': 0, '<': 14, '<?xml version=': 0, '>': 16, 'https://': 7, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '3', 'dir': 'jupyter-structured-clean-dedup', 'id': '107929', 'license': 'CC-BY-4.0', 'n_tokens_mistral': 3112, 'n_tokens_neox': 2809, 'n_words': 1388, 'path': 'practicals/practical2-empty.ipynb', 'repo_name': 'leouieda/lithosphere', 'size': '10269'}