{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Tutorial 8: Subduction Model\n",
" \n",
"Rebecca Farrington rebecca.farrington (at) unimelb.edu.au"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This two dimensional subduction model has a dense, high viscosity 3 layered plate overlying a lower viscosity mantle. The upper and lower plate layers have a visco-plastic rheology, yielding under large stresses. The middle, core layer has a viscous only rheology, maintaining strength during bending. The top 1000 km of the mantle is included, the upper & lower mantle is partitioned with a viscosity contrast of 100x at 600 km depth. The velocity boundary conditions on the domain are period side, free-slip top and no-slip bottom wall. \n",
"\n",
"**References**\n",
"\n",
"1. OzBench, M.; Regenauer-Lieb, K.; Stegman, D. R.; Morra, G.; Farrington, R.; Hale, A.; May, D. A.; Freeman, J.; Bourgouin, L.; Mühlhaus, H. & Moresi, L. A model comparison study of large-scale mantle-lithosphere dynamics driven by subduction. Physics of the Earth and Planetary Interiors, 2008, 171, 224-234. [OzBench, 2008](http://www.sciencedirect.com/science/article/pii/S0031920108002318)"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"loaded rc file /opt/venv/lib/python3.7/site-packages/UWGeodynamics/uwgeo-data/uwgeodynamicsrc\n"
]
}
],
"source": [
"import UWGeodynamics as GEO\n",
"from UWGeodynamics import visualisation as vis"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"u = GEO.UnitRegistry"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Scaling"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"model_length = 1000. * u.kilometer\n",
"refDensity = 3300. * u.kilogram / u.meter**3\n",
"refViscosity = 1.4e19 * u.pascal * u.second\n",
"\n",
"KL = model_length\n",
"KM = refDensity * KL**3\n",
"Kt = 1.0 / (refViscosity / KM * KL)\n",
"\n",
"GEO.scaling_coefficients[\"[length]\"] = KL\n",
"GEO.scaling_coefficients[\"[time]\"] = Kt\n",
"GEO.scaling_coefficients[\"[mass]\"]= KM"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Model Geometry"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"Model = GEO.Model(elementRes=(288, 128), \n",
" minCoord=(0. * u.kilometer, 400. * u.kilometer), \n",
" maxCoord=(4000. * u.kilometer, 1000. * u.kilometer),\n",
" periodic = [True, False],\n",
" gravity=(0.0, -10. * u.meter / u.second**2))"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"Model.outputDir=\"outputs_tutorial8\""
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"upperMantleShape = GEO.shapes.Layer2D(top=Model.top, bottom=400. * u.kilometer)\n",
"upperSlabShape = GEO.shapes.Polygon([(1200. * u.kilometer, 1000. * u.kilometer),\n",
" (3400. * u.kilometer, 1000. * u.kilometer),\n",
" (3350. * u.kilometer, 975. * u.kilometer),\n",
" (1200. * u.kilometer, 975. * u.kilometer),\n",
" (1020. * u.kilometer, 900. * u.kilometer),\n",
" (1020. * u.kilometer, 925. * u.kilometer)])\n",
"lowerSlabShape = GEO.shapes.Polygon([(1200. * u.kilometer, 925. * u.kilometer),\n",
" (3250. * u.kilometer, 925. * u.kilometer),\n",
" (3200. * u.kilometer, 900. * u.kilometer),\n",
" (1200. * u.kilometer, 900. * u.kilometer),\n",
" (1020. * u.kilometer, 825. * u.kilometer),\n",
" (1020. * u.kilometer, 850. * u.kilometer)])\n",
"coreSlabShape = GEO.shapes.Polygon([(1200. * u.kilometer, 975. * u.kilometer),\n",
" (3350. * u.kilometer, 975. * u.kilometer),\n",
" (3250. * u.kilometer, 925. * u.kilometer),\n",
" (1200. * u.kilometer, 925. * u.kilometer),\n",
" (1020. * u.kilometer, 850. * u.kilometer),\n",
" (1020. * u.kilometer, 900. * u.kilometer)])"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"upperMantle = Model.add_material(name=\"Upper Mantle\", shape=upperMantleShape)\n",
"upperSlab = Model.add_material(name=\"Upper Slab\", shape=upperSlabShape)\n",
"lowerSlab = Model.add_material(name=\"Lower Slab\", shape=lowerSlabShape)\n",
"coreSlab = Model.add_material(name=\"Core Slab\", shape=coreSlabShape)"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"Fig = vis.Figure(figsize=(1200,400))\n",
"Fig.Points(Model.swarm, Model.materialField, colours='white green red purple blue', fn_size=2.0, discrete=True)\n",
"Fig.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Densities"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [],
"source": [
"upperMantle.density = 3300 * u.kilogram / u.metre**3\n",
"upperSlab.density = 3380 * u.kilogram / u.metre**3\n",
"lowerSlab.density = 3380 * u.kilogram / u.metre**3\n",
"coreSlab.density = 3380 * u.kilogram / u.metre**3"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Viscosities"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [],
"source": [
"upperMantle.viscosity = 1.4e19 * u.pascal * u.second\n",
"upperSlab.viscosity = 200 * 1.4e19 * u.pascal * u.second\n",
"lowerSlab.viscosity = 200 * 1.4e19 * u.pascal * u.second\n",
"coreSlab.viscosity = 20000 * 1.4e19 * u.pascal * u.second"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Plasticity"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [],
"source": [
"upperSlab.plasticity = GEO.VonMises(cohesion=48. * u.megapascal)\n",
"lowerSlab.plasticity = GEO.VonMises(cohesion=48. * u.megapascal)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Elasticity"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [],
"source": [
"coreSlab.elasticity = GEO.Elasticity(shear_modulus=4.*u.gigapascal, observation_time=20000. * u.year)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Boundary conditions"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
""
]
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"Model.set_velocityBCs(bottom=[0.,0.], top=[None, 0.])"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"Fig = vis.Figure(figsize=(1200,400))\n",
"Fig.Points(Model.swarm, GEO.dimensionalise(Model.viscosityField, u.pascal * u.second), fn_size=2.0, logScale=True)\n",
"Fig.show()"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [],
"source": [
"Model.solver.set_inner_method(\"lu\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"Model.run_for(50000.* u.years)"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.9.6"
}
},
"nbformat": 4,
"nbformat_minor": 2
}