ESS4 524: Heat and Mass Flow Modeling in Earth Sciences
University of Washington
Spring 2020

[an error occurred while processing this directive]
Earth and Space Sciences ESS 524
Introduction to Heat and Mass Flow Modelling in Earth Sciences

Sample Code
This page links to sample matlab code groups on the right sidebar that illustrate ideas in class on heat and mass flow.

Matlab Basics
Demonstration of some Matlab operations and matrix manipulation.

Code Group 1: SS 2D diffusion Practice B
Steady Diffusion in 2D on a Rectangle using Patankar's Practice B (page 70) for node and volume edge positions.

Code Group 2: Transient diffusion - Stability and Accuracy
This 1D code allows you to set time-step size and time-step mixing parameter "alpha" to explore linear computational instability.

More Sample Codes

Transient Heat Flow - I

2-D transient diffusion with implicit time stepping.

  • fvm_soln_transient_bdry_nodes.m sets up parameters a_P etc and includes time-stepping loop.
  • includes a (kludged) variable mixing factor "0<=theta<=1" to allow exploration of implicit, Crank-Nicolson, and explicit schemes.
  • uses same old "solver.m" to solve matrix equation at each time step.
Ice Cap Growth -

Because ice deformation rate depends on surface slope, the surface evolution can be cast as a transient nonlinear diffusion problem for the surface topography. The source term is the net snow accumulation or melt at each location.

  • 2-D transient diffusion with implicit time stepping.

    SS Centered-difference Advection

    This code solves steady advective-diffusion in 1-D using a central-difference representation of advection. This method can have negative coefficients when F=F/D>2. Don't use it for real problems!

    • 1 row of finite volumes
    • zero flux out transverse sides
    • specified values at top and bottom

    m-file diffusion_SS_1D_adv.m runs several different values of nz, number of finite volumes in z direction. Different Peclet numbers are associated with each dz. Some are >2, some are <2.

    Execution "Pauses" after each nz is selected. Press any key to continue.

    SS 1-D Advective Diffusion, Power-Law scheme

    Solves steady 1-D advective diffusion for a series of models with various Peclet numbers and grid sizes. and for either central-difference scheme (bad) or power-law scheme (good).

    Hits "pause" after each plot - press any key to resume.

    SS 2-D Advective Diffusion

    This code puts finite-volume edges on boundaries, but also includes extra nodes on those boundaries to facilitate incorporation of either phi-value or phi-gradient Boundary Conditions at every boundary node. (Patankar Practice B, p. 69).

    Patankar's Power-law advective scheme is used.

    SS Illustration of False Diffusion in 2-D

    SS 2-D Adv-Diff code above is used to run a sequence of models illustrating false diffusion when strong flow is not aligned with coordinate axes.

    Transient 2-D Advective Diffusion

    By incorporating minor changes to the SS 2-D Advective Diffusion code above, this code solves transient problems.

    phi-Dependent Coefficients

    This variation on 2-D advective diffusive code solves a steady 1-D problem where the coefficient Gamma depends on solution phi

    Internal BC

    Sometimes you may want to fix the value of phi at volumes remote from the actual domain boundary. Following Patankar (1980) p. 145, this can be done with the source term. The example uses the 2-D transient advective diffusive code, with a block of volumes held at a constant value of phi.

    SH Wave Propagation

    By converting the first tme derivative into a second time derivative, the diffusion equation can be transformed into a wave equation, applicable to SH waves traveling through the Earth.

    phi becomes displacement u, and Gamma becomes shear modulus.

    1D Stability Analysis

    This code finds wavenumber transfer functions for 1D transient diffusion, for specified kappa, dx, and dt. It then carries out a corresponding 1D time-domain finite difference simulation. Users can see how the transfer functions are useful.

    One-point Transient Response

    Code evaluates the response of a single node when adjacent points are held fixed. Solutions are shown for 3 levels of discretization:  

    • no discretization
    • only space discretized
    • space and time discretized


One-point Transient Response


Linear Computational Instability
same old solver
Notes - a'comin'

SS 1-D Advective Diffusion, central difference
SS 1-D Advective Diffusion, power-law scheme

SS 2-D Advective Diffusion

False Diffusion

Transient 2-D advective diffusion

Ice-cap model

phi-Dependent coefficients

Internal Boundary Conditions

SH Wave propagation