# Exercise 2b: Modeling steady-state heat flow in the lithosphere

Elco Luijendijk

November 2019

<elco.luijendijk@geo.uni-goettingen.de>

## Objectives

* Learn to model heat flow in the crust and lithosphere using Python
* Learn what typical geothermal gradients in the lithosphere look like and how they are affected by thermal parameters and the thickness of the lithosphere


**Deadline**: **17 Jan 2020**. Hand in a version of your jupyter notebook and a short word document with answers to the assignments and the result figures of your numerical model. 

**Grading**: Each assignment is 1 point, for a total of 5 points.

Do not hesitate to ask questions if you get stuck anywhere. You can reach me by email <eluijen@gwdg.de> or pass by at my office, room 122 in the Structural Geology dept.

*Good luck !*

## Introduction
In this exercise we will update the notebook that you completed in exercise 2a to model steady-state heat flow instead of groundwater flow. Heat conduction and groundwater flow are governed by very similar diffusion laws. We will adjust the parameters and boundary conditions of your groundwater model to instead calculate steady-state geothermal gradients in the lithosphere. 

## Rewrite your notebook to model heat flow

The heat flow equation that we will use in this exercise is a combination of Fouriers's law:

\begin{equation}
    q = - K \dfrac{\partial T}{\partial x}
\end{equation}

and the heat balance equation:

\begin{equation}
    \dfrac{\partial q}{\partial x} = W
\end{equation}

Note the similarity of Fourier's law and Darcy's law. Instead of hydraulic head (*h*) we use temperature (*T*), and instead of hydraulic conductivity we use thermal conductivity (both denoted by *K*). The source term in this case is not recharge or groundwater pumping, but the generation of heat by the decay of radioactive elements. 

Copy the fast version of the diffusion function (``solve_steady_state_diffusion_eq_faster``) from your notebook of exercise 2a into this notebook below.

Change the function to model heat flow model by changing all the variables called ``h``, ``h_new``, ``h_old`` to ``T``, ``T_new``, ``T_old``.


```python
# importing external python modules
# numpy for working with arrays:
import numpy as np
# and matplotlib to make nice looking figures
import matplotlib.pyplot as pl
```


```python
# copy your finite difference diffusion function here...



```

**Assignment 1** Write down the finite difference approximation of the steady-state heat flow equation. This is the heat flow version of equation 10 in your handout of exercise 1. What units do heat flow (*q*) and heat production (*W*) have?

## Assigning new parameters

Next we have to assign new parameter values to our notebook. We will model the geothermal gradient in the lithosphere. We will increase the size of our model (variable ``L``) to a 100 km, which is a value representative of the thickness of the continental lithosphere. To change this, copy the code block where you assign parameters in exercise 2a below and change the line assigning ``L`` in your notebook like this:

~~~~python
L = 100000
~~~~

For thermal conductivity we will use an initial value of 2.5 W m^-1 K^-1. Our model will be one dimensional, so the thickness parameter ``b`` should be equal to 1. And for the grid cell size of our numerical model we will use a value of 200 m (``dx = 200``). For now we will leave the source term zero, so ``W_array[:] = 0``.


```python
# copy the code blocks where you assign parameters here....
```

## Adjust the boundary conditions

In contrast to the groundwater model, the right hand side / bottom of our model domain is not insulated. Instead there is a heat flux to the lithosphere from the astenosphere. The base of the lithosphere is usually defined as a thermal boundary, and this boundary is commonly assumed to be a temperature of 1300 °C. This means that we can define the right hand boundary condition as a specified temperature of 1300 °C. Change the following line in the ``solve_steady_state_diffusion_eq_faster`` function:

~~~~python
T_new[-1] = ...
~~~~

to 

~~~~python
T_new[-1] = 1300.0
~~~~

This makes sure that the last node in our model always has a temperature of 1300 °C. For the top boundary we can use an average global surface temperature of 10 °C.

## Calculate a steady-state geothermal gradient

With the newly assigned variables and boundary condition we are ready to run our numerical model and calculate an average steady-state geothermal gradient in the lithosphere. 

**Assignment 2** Copy the line where you run the diffusion function (``h = solve_steady....``) from exercise 2a. Run your new heat flow model. Try to experiment with the number of iterations. What number of iterations do you need (approximately) to reach a steady-state geothermal gradient?


```python
# copy the line that calls the diffusion function here:



```

## Adding realistic thermal parameters

In exercise 1 we have already have set up an array for the source term (``W``), which we can now use to vary heat production in the lithosphere. As you have learned in the Python tutorial you can perform operations like assigning numbers to sections of arrays like this:

~~~~python
W_array[10:30] = 1e-6
~~~~

This assign a heat production of 1 x 10-6^ W m^-3 from node 10 to 30. Note: this is an example, do not insert this line into your script yet. There is an even more convenient way to assign values based on their depth. We can also assign values to parts of arrays like this:

~~~~python
W_array[x < 10000] = 1e-6
~~~~

This will change the heat flow values in the upper 10000 m of you model domain. For example, to assign numbers to a section running from 10000 to 20000 m, try the following:

~~~~python
W_array[(x >= 10000) & (x < 20000)] = 1e-6
~~~~

Note that the sign ``>=`` means larger than or equal to a number. 

Look up heat production for the upper crust, lower crust and the mantle from Cloetingh et al. (2005), Table 1. Use the reference at the bottom of this notebook and look up the paper on google scholar. Follow the examples above to assign heat production values for these three sections of the lithosphere. For the depth of the upper crust we can use a value of 12 km and for the depth of the lower crust 35 km. Try to place the new lines of code where you define heat production *after* the line where you calculate the first steady-state geothermal gradient (ie., after first call of the ``solve_steady_state_diffusion_eq function``). Next, add a line to calculate the new geothermal gradient including heat production:

~~~~python
T2 = solve_steady_state_diffusion_eq_faster(dx, K, W_array, T0)
~~~~

This makes sure that you store the new more exact geothermal gradient in a new variable ``T2``, which you can then compare to the old geothermal gradient ``T``.

Note that thermal conductivity also varies in the crust and in the mantle. However, the solution we derived for steady-state groundwater flow and heat flow equations assumes that *K* is constant, since we moved *K* out of the derivative. See the handout of exercise 1, equations 5 and 6. Therefore with our current simplified model we cannot model a variable thermal conductivity.

**Assignment 3** Run your new model with improved thermal parameters. Plot both the old and new temperature curve in the same panel. Try to explain the shape of the new geothermal gradient in a few words or sentences, why is there a curvature?

## Calculate heat flow

We can use Fourier's law to calculate heat flow, using the values of the geothermal gradient that you just calculated. Use the following lines to calculate the gradient and the heat flow in the lithosphere:


```python
T_gradient = (T2[1:] - T2[:-1]) / dx
q = T_gradient * K
```

## Add temperature and heat flow to a figure

We will try to make a figure that shows both the change in temperature and heat flow in the lithosphere. First copy the code block that generates a figure from exercise 2a below:


```python
# copy code block that makes a figure here:


```

First try to change the command where you plot temperature so that temperature is shown on the x-axis and depth (x) on the y-axis, instead of the other way around.

Next, we will try to set up a figure with two panels, one for temperature and one for heat flow. Add one more panel to the figure to show temperature and heat flow side by side. The following line creates a new figure containing two panels side by side.

~~~~python
fig, panels = pl.subplots(1, 2, figsize=(10, 6))
~~~~

*Replace* the current line that creates the new figure with this new line. The variable ``panels`` is now a list containing two panels. Each panel can be called using either ``panels[0]`` or ``panels[1]``. 
Replace all existing lines of code where you use ``panel.`` with ``panels[0].``. Then add some lines to plot heat flow in the second panel (``panels[1]``):

~~~~python
x_mid = (x[1:] + x[:-1]) / 2.0 
panels[1].plot(q, x_mid, color='black')
~~~~

we can also flip the y-axis to make sure the surface (x=0) is at the top of the figure and not the bottom:

~~~~python
panels[0].set_ylim(100000, 0)
panels[1].set_ylim(100000, 0)
~~~~

**Assignment 4** Calculate heat flow in the lithosphere and make a figure of temperature and heat flow. What percentage of the heat flow at the surface is supplied by the mantle? (hint: compare *q* at x=0 with *q* at x=mantle depth). 

Note: make sure that your model has reached steady state, you may need to run significantly more iterations once you add heat production.

## Heat flow in different geological settings

The thickness of the lithosphere exerts a strong influence on geothermal gradients, because it is a thermal boundary. We will calculate representative geothermal gradients for three geological settings with a different lithosperic thickness: 1) a craton, 2) average continental crust and 3) oceanic crust and pick a fourth location that your are interested in (the Goettingen campus, mount everest, antarctica, etc...). 
Pick each of these settings and look up values for the lithosphere thickness for a location of choice. For oceanic lithosphere look up thickness in Conrad et al. (2006). Refer to Artemieva et al. (2006) for estimates of the thickness of the continental lithosphere. Estimate the thickness of the crust using Reguzzoni et al. (2013), Fig. 12. You can find these references below.

The thermal parameters of oceanic and continental lithosphere also vary. For the continental lithosphere you can use the parameters that you have already assigned. The heat generation in the oceanic lithosphere is relatively low, on average 0.5 x 10^-6 W m^-3 (Allen & Allen 2005, note: this reference is not available online). The average thermal conductivity (*K*) of the oceanic lithosphere is 3.14 W m^-1 K^-1 (Stein 1992). 

**Assignment 5** Model heat flow for the four different geological settings, make a figure of temperature and heat flow and compare the modeled heat flow with a database of global heat flow (Davies & Davies 2010): http://www.heatflow.und.edu/. Describe in words why the heat flow and geothermal gradients vary between the different settings, and why the models match the global heat flow data or why there is a difference.

## References

Allen, P.A., Allen, J.R., 2005. Basin analysis: principles and applications. Blackwell publishing, Oxford.

Artemieva, I.M., 2006. Global thermal model TC1 for the continental lithosphere: Implications for lithosphere secular evolution. Tectonophysics 416, 245–277.

Cloetingh, S.A.P.L., Ziegler, P.A., Beekman, F., Andriessen, P.A.M., Hardebol, N., Dèzes, P., 2005. Intraplate deformation and 3D rheological structure of the Rhine Rift System and adjacent areas of the northern Alpine foreland. Int. J. Earth Sci. 94, 758–778.

Conrad, C.P., Lithgow-Bertelloni, C., 2006. Influence of continental roots and asthenosphere on plate-mantle coupling. Geophys. Res. Lett. 33, 2–5. doi:10.1029/2005GL025621

Davies, J.H., Davies, D.R., 2010. Earth’s surface heat flux. Solid Earth 1, 5–24. doi:10.5194/se-1-5-2010

Reguzzoni, M., Sampietro, D., Sansò, F., 2013. Global moho from the combination of the CRUST2.0 model and GOCE data. Geophys. J. Int. 195, 222–237. doi:10.1093/gji/ggt247

Stein, C.A., Stein, S., 1992. A model for the global variation in oceanic depth and heat flow with lithospheric age. Nature 356, 133–135. doi:10.1038/359123a0

(note use google scholar to find any of these papers, https://scholar.google.com/)


```python

```

77a36bec0128acbd0d27070bf35e650bb87b4e00
ipynb
Jupyter Notebook
exercises/exercise_2_steady_state_flow_models/.ipynb_checkpoints/exercise_2b-checkpoint.ipynb
ElcoLuijendijk/fluids_in_the_crust
c2cadb0a91e9f9ed62094ac5e796168fef0d1a3e
CC-BY-4.0
2021-01-12T19:08:16.000Z
2021-01-13T14:27:42.000Z
exercises/exercise_2_steady_state_flow_models/.ipynb_checkpoints/exercise_2b-checkpoint.ipynb
ElcoLuijendijk/fluids_in_the_crust
c2cadb0a91e9f9ed62094ac5e796168fef0d1a3e
CC-BY-4.0
exercises/exercise_2_steady_state_flow_models/.ipynb_checkpoints/exercise_2b-checkpoint.ipynb
ElcoLuijendijk/fluids_in_the_crust
c2cadb0a91e9f9ed62094ac5e796168fef0d1a3e
CC-BY-4.0
Qwen/Qwen-72B
1. YES
2. YES