[aspect-devel] Continental extension model

[John Naliboff]

Hi Mohamed,

An example continental extension input file is attached, which will be 
submitted as a cookbook example in the next few days.

In this example, extension is driven by prescribed outflow on the sides 
and inflow at the base.  The rheology is dislocation creep with 
different flow laws for the upper/lower crust and mantle. Internal 
friction angle (20 degrees) (20 MPa) are within the range of commonly 
used values.  No strain-weakening, but I have a pull request open that 
implements this.

Most extension problems include the asthenosphere, while this model only 
goes down to 100 km.  As such, I would only use this type of setup for 
studying early stage extension.  Any problems examining extension from 
start to breakup should extend to at least 150 km. In this case, one 
might alter the simple boundary conditions I've prescribed.

Hope this helps and let me know you if you have any questions about the 
input file.

Cheers,
John




*************************************************
John Naliboff
Assistant Project Scientist, CIG
Earth & Planetary Sciences Dept., UC Davis

On 10/10/2016 12:59 AM, Mohamed Gouiza wrote:
>
> Hi John,
>
> I looked for the continental extension model that you showed in the 
> online workshop last month, but couldn’t find it in tests/ folder.
>
> I’ve been running several extension models with the visco-plastic 
> material model and I am interested in knowing the visco-plastic law 
> parameters that you used and how do you prescribe the boundary 
> conditions: is the extension rate defined the same way as in the 
> crustal deformation example in the cookbook by Cedric?
>
> Thank you
>
> -------------------------------------------------
> Mohamed Gouiza, Research Fellow
> Basin Structure Group, Institute of Applied Geosciences
> University of Leeds, School of Earth and Environment
>
> Leeds,  LS2 9JT, UK
>
> M.Gouiza at leeds.ac.uk <mailto:M.Gouiza at leeds.ac.uk>
> +44 7985 782073
> -------------------------------------------------
>
>
>
> _______________________________________________
> Aspect-devel mailing list
> Aspect-devel at geodynamics.org
> http://lists.geodynamics.org/cgi-bin/mailman/listinfo/aspect-devel

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.geodynamics.org/pipermail/aspect-devel/attachments/20161010/862263e2/attachment-0001.html>
-------------- next part --------------
#####  Global parameters
set Dimension                              = 2
set Start time                             = 0
set End time                               = 5e6
set Use years in output instead of seconds = true
set Linear solver tolerance                = 1e-7
set Nonlinear solver scheme                = iterated Stokes
set Nonlinear solver tolerance             = 1e-4
set Max nonlinear iterations               = 10
set Number of cheap Stokes solver steps    = 0
set CFL number                             = 0.5
set Output directory                       = output
set Timing output frequency                = 1
set Pressure normalization                 = no

#### Parameters describing the model

# Model geometry (400x100 km, 2 km spacing)
subsection Geometry model
  set Model name = box

  subsection Box
    set X repetitions = 200
    set Y repetitions = 50
    set X extent      = 400e3
    set Y extent      = 100e3
  end
end

# Mesh refinement specifications
subsection Mesh refinement
  set Initial adaptive refinement        = 0
  set Initial global refinement          = 0
  set Time steps between mesh refinement = 0
end

# Boundary classifications
subsection Model settings
  set Include adiabatic heating               = false
  set Include shear heating                   = false
  set Fixed temperature boundary indicators   = bottom, top
  set Prescribed velocity boundary indicators = left x: function, right x:function, bottom y:function
  set Free surface boundary indicators        = top
end

# Advecting the free surface vertically rather than
# in the surface normal direction can result in a
# more stable mesh when the deformation is large
subsection Free surface
  set Surface velocity projection = vertical
end

# Velocity on boundaries characterized by functions
subsection Boundary velocity model
  subsection Function
    set Variable names      = x,y
    set Function constants  = m=0.0025, year=1
    set Function expression =  if (x<200e3 , -1*m/year, 1*m/year) ; 0.00125
  end
end


# Number and name of compositional fields
subsection Compositional fields
  set Number of fields = 4
  set Names of fields = upper, lower, mantle, seed
end


# Spatial domain of different compositional fields
subsection Compositional initial conditions
  set Model name = function
  subsection Function
    set Variable names      = x,y
    set Function expression = if(y>=80.e3, 1, 0); \
                              if(y<80.e3 && y>=70.e3, 1, 0); \
                              if(y<70.e3 && y>-100.e3,1, 0); \ 
                              if(y<68.e3 && y>60.e3 && x>=198.e3 && x<=202.e3 , 1, 0);
  end
end

# Temperature boundary conditions
subsection Boundary temperature model
  set Model name = box

  subsection Box
    set Bottom temperature = 1573
    set Left temperature   =  273
    set Right temperature  =  273
    set Top temperature    =  273

  end
end

# Initial temperature field
# Typical continental geotherm based equations 4-6 from Chapman 1986 (Geol. Soc. Lon.)
# The initial constraints are:
#   Layer Surface Temperature - upper crust (ts1) = 273 K; 
#                               mantle      (ts3) = 823 K;  
#   Model Base Temperature - (tb) = 1573 K;
#   Heat Production - upper crust (A) = 1.5e-6 W/m^3; 
#   Thermal Conductivity - upper crust (k1) = 2.5 (W/(m K)); 
#                          lower crust (k2) = 2.5 (W/(m K)); 
#                          mantle      (k3) = 3.3 (W/(m K));
# To satisfy these constraints, the following values are required:
#   Layer Surface Heat Flow - upper crust (qs1) = 0.065357 W/m^2; 
#                             lower crust (qs2) = 0.035357 W/m^2; 
#                             mantle      (qs3) = 0.035357 W/m^2;
#   Temperature - base of upper crust (ts2) = 681.5714
subsection Initial conditions
  set Model name = function
  subsection Function
    set Variable names = x,y
    set Function constants = h=100e3,ts1=273,ts2=681.5714,ts3=823., \
                                     k1=2.5,k2=2.5,k3=3.3,A=1.5e-6, \
                             qs1=0.0653571,qs2=0.035357,qs3=0.035357,qb3=0.035357
    set Function expression = if( (h-y)<=20.e3, \
                                  ts1 + (qs1/k1)*(h-y) - (A*(h-y)*(h-y))/(2.0*k1), \
                                  if((h-y)>20.e3 && (h-y)<=30.e3, ts2 + (qs2/k2)*(h-y-20.e3),\
                                  ts3 + (qs3/k3)*(h-y-30.e3)));
  end
end

# Internal heating (only internal heating for crust)
# Assuming a very large decay time, the heat production (W/m^3) is effectively:
#   radioactive_heating_rate (W/kg) * density (kg/m^3)
# For a reference crustal density of 2800 kg/m^3 and 1.5e-6 W/m^3 heating production,
# the heating rate (assuming a single element) is 1.5e-6/2800. = 5.357142857e-10.
subsection Heating model
  set Model name = radioactive decay
  subsection Radioactive decay
    set Number of elements    = 1
    set Heating rates                 = 5.357142857e-10
    set Half decay times              = 1.e20
    set Initial concentrations mantle = 0.0
    set Initial concentrations crust  = 1.0
    set Crust defined by composition  = true
    set Crust composition number      = 0
  end
end

# Material model
subsection Material model
  set Model name = visco plastic

  subsection Visco Plastic

    set Reference temperature = 293
    set Minimum strain rate = 1.e-20
    set Reference strain rate = 1.e-16
    set Minimum viscosity = 1e18
    set Maximum viscosity = 1e26
    set Reference viscosity = 1e22

    set Thermal diffusivities = 1.333333e-6,1.190476e-6,1.149425e-6,1.333333e-6,1.333333e-6
    set Heat capacities = 750.,750.,750.,750.,750.
    set Densities = 3300,2800,2900,3300,3300
    set Thermal expansivities = 2e-5,2e-5,2e-5,2e-5,2e-5

    set Viscosity averaging scheme = harmonic
    set Viscous flow law = dislocation

    set Prefactors for dislocation creep = 6.52e-16,8.57e-28,7.13e-18,6.52e-16,7.13e-18
    set Stress exponents for dislocation creep = 3.5,4.0,3.0,3.5,3.0
    set Activation energies for dislocation creep = 530.e3,223.e3,345.e3,530.e3,345.e3
    set Activation volumes for dislocation creep = 18.e-6,0.,0.,18.e-6,0.

    set Angles of internal friction = 20.,20.,20.,20.,20.
    set Cohesions = 20.e6,20.e6,20.e6,20.e6,20.e6

  end
end

# Gravity model
subsection Gravity model
  set Model name = vertical

  subsection Vertical
    set Magnitude = 9.81
  end
end


# Post processing
subsection Postprocess
  set List of postprocessors = velocity statistics, basic statistics, temperature statistics, visualization
  subsection Visualization
    set List of output variables = density, viscosity, strain rate, error indicator
    set Output format = vtu
    set Time between graphical output = 1e6
    set Interpolate output = true
  end
end