In [1]:
import underworld as uw
import glucifer
import numpy
import matplotlib.pyplot as plt
import mpi4py

uw.matplotlib_inline()
plt.ion()

comm = mpi4py.MPI.COMM_WORLD
rank = comm.Get_rank()
size = comm.Get_size()

Temperature Advection Due to Groundwater Flow

Underworld can efficiently solve the diffusion-advection equation and can set up as many simultaneous instances of it in a model as we would like. Heat-flow and groundwater flow are both described by this equation, so we can set these up and couple them, to model temperature advection by groundwater.

Temperature Diffusion-Advection Equation

$\frac{\partial T} {\partial t} = \kappa \nabla^2 T + G$

where: $\kappa = \frac{k}{C \rho}$

($\kappa$ - thermal diffusivity, $T$ - temperature, $t$ - time, $G$ - groundwater advection term, $C$ - rock heat capacity, $\rho$ - rock density)

Groundwater Flow Equation

$ \nabla\cdot q = 0 $

$ q = \kappa_H \left(\nabla H + \rho_w g\right) $ where: $\kappa_H = \frac{k_H}{\mu S}$

In this particular example, we'll leave out the gravitational term ($\rho_w g$) and drive flow with just a pressure gradient.

($\kappa_H$ - hydraulic diffusivity, $H$ - groundwater pressure, $k_H$ - hydraulic permeability, $\mu$ - water viscosity, $S$ - specific storage, $g$ - gravitational accelleration, $\rho_w$ - density of water)

Coupling

$G = -\nabla T \cdot \frac{\rho_w C_w}{\rho C} q$

($C_w$ - heat capacity of water)

Coupling is implemented by handing an 'effective velocity field', $\frac{\rho_w C_w}{\rho C} q$ to the temperature diffusion-advection solver.

Set up mesh and fields

In [2]:
elementType = "Q1"
resX = 128
resY = 64
mesh = uw.mesh.FeMesh_Cartesian( elementType = (elementType), 
                                 elementRes  = (resX, resY), 
                                 minCoord    = (0, -1.), 
                                 maxCoord    = (2., 0.)) 

temperatureField    = mesh.add_variable( nodeDofCount=1 )
temperatureDotField    = mesh.add_variable( nodeDofCount=1 )

gwPressureField    = mesh.add_variable( nodeDofCount=1 )
hydraulicDiffusivityField    = mesh.add_variable( nodeDofCount=1 )
gwPressureDotField    = mesh.add_variable( nodeDofCount=1 )

velocityField    = mesh.add_variable( nodeDofCount=2 )

domainVolume = (mesh.maxCoord[1] - mesh.minCoord[1]) * (mesh.maxCoord[0] - mesh.minCoord[0])
averageTemperature = uw.utils.Integral(fn = temperatureField / domainVolume,mesh= mesh,integrationType="volume")

Set up the types of boundary conditions

In [3]:
iWalls = mesh.specialSets["MinI_VertexSet"] + mesh.specialSets["MaxI_VertexSet"]
jWalls = mesh.specialSets["MinJ_VertexSet"] + mesh.specialSets["MaxJ_VertexSet"]

temperatureBC = uw.conditions.DirichletCondition( variable        = temperatureField, 
                                               indexSetsPerDof = ( jWalls) )

gwPressureBC = uw.conditions.DirichletCondition( variable        = gwPressureField, 
                                               indexSetsPerDof = ( jWalls) )

Set initial conditions and the values of the boundary conditions

In [4]:
#Lower groundwater pressure boundary condition
maxgwpressure = 1.

yCoordFn = uw.function.input()[1]

initialFn = -1. * yCoordFn
temperatureField.data[:] = initialFn.evaluate(mesh)

initialFn = -1. * maxgwpressure * yCoordFn
gwPressureField.data[:] = initialFn.evaluate(mesh)
In [5]:
figMaterial = glucifer.Figure( figsize=(800,400), title="Initial Temperature Field" )
figMaterial.append( glucifer.objects.Surface(mesh,temperatureField ))
if size == 1:
    figMaterial.show()

Set up swarm

In [6]:
swarm         = uw.swarm.Swarm( mesh=mesh )
swarmLayout   = uw.swarm.layouts.PerCellGaussLayout( swarm=swarm, gaussPointCount=1 )
swarm.populate_using_layout( layout=swarmLayout )

Assign materials to particles

Here we are deciding on the distribution of different materials. The benchmark is for a large region which extends from the bottom to the top of the domain, in this case on the left hand half of the model. Once you are happy with the benchmark, feel free to change the model geometries!

In [7]:
materialIndex  = swarm.add_variable( dataType="int",    count=1 )

materialPorous        = 0
materialImpermeable   = 1

xCoordFn = uw.function.input()[0]
yCoordFn = uw.function.input()[1]

conditions = [ (yCoordFn > 0., materialPorous),
                ( xCoordFn < 1. , materialPorous),
               ( True , materialImpermeable )]


materialIndex.data[:]  = uw.function.branching.conditional( conditions ).evaluate(swarm)

Visualise the materials we have just assigned

In [8]:
figMaterial = glucifer.Figure( figsize=(800,400), title="Initial Material Distribution" )
figMaterial.append( glucifer.objects.Points(swarm, materialIndex, pointSize=4.0) )
figMaterial.append( glucifer.objects.Mesh(mesh))
if size == 1:
    figMaterial.show()

Assign material properties

Here we have material properties which are chosen relative to the thermal diffusion time-scale. Leave the thermal diffusivity as 1, but you can change the maximum hydraulic diffusivity to control their relative time-scale (along with the chosen pressure gradient above).

The storage capacity effectively controls the ratio of groundwater pressure diffusion to groundwater flow. The ratio of $\rho_w c_w / \left( \rho c \right)$ controls the rate of energy transfer from the water to the rock matrix. If you turn on the gravity term, $\rho_w$ also affects the gravity term in the ground-water flow equation ... denser water will sink in the absence of an upward pressure decrease.

In [9]:
maxHydDiff = 1.
hydraulicDiffusivityMap = { materialPorous : maxHydDiff, 
                 materialImpermeable : 1e-2}
hydraulicDiffusivityMapFn = uw.function.branching.map( fn_key = materialIndex, mapping = hydraulicDiffusivityMap )

#Hydraulic storage capacity
Storage = 1.

# coeff is equivalent to rho_water / rho_rock * c_water / c_rock
coeff = 1.

thermalDiff = 1.
thermalDiffusivityMap = { materialPorous : thermalDiff, 
                 materialImpermeable : thermalDiff}
thermalDiffusivityMapFn = uw.function.branching.map( fn_key = materialIndex, mapping = thermalDiffusivityMap )

Setup groundwater equations

In [10]:
gwadvDiff = uw.systems.SteadyStateHeat( temperatureField = gwPressureField, fn_diffusivity = hydraulicDiffusivityMapFn,conditions=[gwPressureBC])
gwadvDiff = uw.systems.SteadyStateDarcyFlow( velocityField=velocityField, pressureField = gwPressureField, 
                                            fn_diffusivity=hydraulicDiffusivityMapFn, fn_bodyforce=(0.,0.),
                                           voronoi_swarm=swarm, conditions=gwPressureBC)

gwsolver = uw.systems.Solver(gwadvDiff)

Gradients - needed for calculating ground-water flow and advection

In [11]:
tempGrad = temperatureField.fn_gradient

Ground-water Velocity

We solve for the groundwater pressure diffusion first, as it temperature-independent. There is no source term, so we can use the steady-state solver.

In [12]:
gwsolver.solve()
In [13]:
figMaterial = glucifer.Figure( figsize=(800,400), title="Initial ground-water pressure and velocity vectors" )
figMaterial.append( glucifer.objects.Surface(mesh,gwPressureField ))
figMaterial.append(glucifer.objects.VectorArrows(mesh,velocityField,scaling=5e-2,arrowHead=0.3))
if size == 1:                   
    figMaterial.show()

Setup temperature advection-diffusion solver

In [14]:
tempDiff = uw.systems.SteadyStateHeat( temperatureField, fn_diffusivity=thermalDiffusivityMapFn,
            conditions=temperatureBC,
           fn_heating = uw.function.math.dot(-1.*coeff * velocityField, tempGrad))

tempSolver = uw.systems.Solver(tempDiff)

For a ground-water velocity $q$, this the following is an analytic solution to the 1D temperature diffusion-advection equation: $\frac{d^2 T}{dy^2} + q'\frac{dT}{dy} = 0$

where: $q' = \frac{\rho_w c_w}{\rho c} \frac{q}{\kappa}$

The solution, assuming $H=0$ and $T=0$ at the surface and $H = H_{max}$ and $T=1$ at the base of the domain, is:

$T(y) = \alpha \left(e^{q' y} - 1 \right)\ \ \ $ where: $\alpha = \left( e^{-q'}-1 \right)^{-1}$ (you can check it's actually a solution)

We can use this to benchmark an effectively 1D part of the domain.

In [15]:
def analyticsol(y):
    qp = coeff * maxHydDiff / thermalDiff * Storage * (maxgwpressure)
    
    if qp == 0.:
        return -y
    else:
        return (numpy.exp(qp*y) - 1.) / (numpy.exp(-qp) -1. )

Now solve for the temperature.

Because the heat-source term (G) depends on the temperature gradient, the solution is non-linear. For our simple geometry, it should converge after one iteration though.

In [16]:
tempSolver.solve(nonLinearIterate=True)

Here is our temperature field. Does it look perturbed?

In [17]:
figMaterial = glucifer.Figure( figsize=(800,400),title="Temperature Field" )
figMaterial.append( glucifer.objects.Surface(mesh,temperatureField ))
if size == 1:
    figMaterial.show()

figMaterial.save_image("TemperatureField.png")
Out[17]:
'TemperatureField.png'

Plot profiles of temperature with depth at the left wall, centre and right wall. The walls should agree with the analytic solutions with and without advection ... unless you're playing around with the geometry. Then they're a useful reference.

In [18]:
n = 50
arrY = numpy.linspace(-1.,0.,n)
arrT = numpy.zeros(n)

if rank == 0:
    plt.clf()

for xpos in [0.,1.,2.]:
    arrT = numpy.zeros(n)
    for i in range(n):
        arrT[i] = temperatureField.evaluate_global((xpos,arrY[i]))

    if rank == 0:
        plt.plot(arrT,arrY,label="x=%i" %xpos)


# Analytic Solution
if rank == 0:
    arrY = numpy.linspace(-1.,0.,10)
    arrAnalytic = numpy.zeros(10)
    for i in range(10):
        arrAnalytic[i] = analyticsol(arrY[i])

    plt.scatter(arrAnalytic,arrY,label="Analytic Advection",lw=0,c="Blue")
    plt.scatter(-1. * arrY,arrY, label="Analytic Non-Advect",lw=0,c="Black")

    plt.legend(loc=0,frameon=False)
    plt.xlabel('Temperature')
    plt.ylabel('Depth')

    plt.xlim(0,1)
    plt.ylim(-1,0)
    plt.title('Temperature Profiles')
    if size == 1:
        plt.show()

    plt.savefig("Geotherms.pdf")

<matplotlib.figure.Figure at 0x115fb40d0>

Now go feel free to go back and make some changes:

  1. Lower the maximum groundwater pressure until the perturbed geotherm is colder than the linear geotherm
  2. Remove the low permeability material completely (i.e. give everything a high permeability) and check that the analaytical solution almost exact. Then explore how wide the high permeablity zone has to be for it to be effectively 1D
  3. Go to where the materials are set up and make some more complicated geometries.
  4. Set up a horizontal groundwater pressure gradient, make sure the thermal diffusivity is homogeneous over the domain and turn off gravity. You will need to alter the boundary conditions. If done correctly, the temperature field should not be perturbed.
  5. What does the ratio of hydraulic diffusivity to storage control?