# Using a free surface in a model with a crust#

This section was contributed by William Durkin.

This cookbook is a modification of the previous example that explores changes in the way topography develops when a highly viscous crust is added. In this cookbook, we use a material model in which the material changes from low viscosity mantle to high viscosity crust at $$z = z_j = \texttt{jump height}$$, i.e., the piecewise viscosity function is defined as

\begin{split}\begin{aligned} \eta(z) = \left\{ \begin{matrix} \eta_U & \text{for}\ z > z_j, \\ \eta_L & \text{for}\ z \le z_j. \end{matrix} \right. \end{aligned}\end{split}

where $$\eta_U$$ and $$\eta_L$$ are the viscosities of the upper and lower layers, respectively. This viscosity model can be implemented by creating a plugin that is a small modification of the simpler material model (from which it is otherwise simply copied). We call this material model “SimplerWithCrust.” In particular, what is necessary is an evaluation function that looks like this:

template <int dim>
void
SimplerWithCrust<dim>::
evaluate(const typename Interface<dim>::MaterialModelInputs &in,
typename Interface<dim>::MaterialModelOutputs &out ) const
{
for (unsigned int i=0; i<in.n_evaluation_points(); ++i)
{
const double z = in.position[i][1];

if (z>jump_height)
out.viscosities[i] = eta_U;
else
out.viscosities[i] = eta_L;

out.densities[i] = reference_density * (1.0 - thermal_expansion_coefficient * (in.temperature[i] - reference_temperature));
out.thermal_expansion_coefficients[i] = thermal_expansion_coefficient;
out.specific_heat[i] = reference_specific_heat;
out.thermal_conductivities[i] = thermal_conductivity;
out.compressibilities[i] = 0.0;
out.entropy_derivative_pressure[i] = 0.0;
out.entropy_derivative_temperature[i] = 0.0;
for (unsigned int c=0; c<in.composition[i].size(); ++c)
out.reaction_terms[i][c] = 0.0;
}
}


Additional changes make the new parameters Jump height, Lower viscosity, and Upper viscosity available to the input parameter file, and corresponding variables available in the class and used in the code snippet above. The entire code can be found in cookbooks/free_surface_with_crust/plugin/simpler_with_crust.cc. Refer to The idea of plugins and the SimulatorAccess and Introspection classes for more information about writing and running plugins.

The following changes are necessary compared to the input file from the cookbook shown in Using a free surface to include a crust:

• Load the plugin implementing the new material model:

# This model loads an additional shared library, which provides another
# material model that is not part of the main ASPECT code, but instead
# lies in in a subdirectory of this cookbook folder.
set Additional shared libraries            = ./plugin/libsimpler_with_crust.so

• Declare values for the new parameters:

# Because we load the additional shared library at the top, we here have the
# option to choose the new material model that is described inside that library
# called 'simpler with crust'. This model has additional input options for the
# two layers with different viscosities and the depth of the transition between
# the two layers.
subsection Material model
set Model name = simpler with crust

subsection Simpler with crust model
set Reference density             = 3300
set Reference specific heat       = 1250
set Reference temperature         = 0.0
set Thermal conductivity          = 1.0
set Thermal expansion coefficient = 4e-5
set Lower viscosity               = 1.e20
set Upper viscosity               = 1.e23
set Jump height                   = 170.e3
end
end


Note that the height of the interface at 170km is interpreted in the coordinate system in which the box geometry of this cookbook lives. The box has dimensions $$500\text{ km}\times 200\text{ km}$$, so an interface height of 170km implies a depth of 30km.

The entire parameter file is located in cookbooks/free_surface_with_crust/free_surface_with_crust.prm.

Running this input file generates a crust that is 30 km thick and 1000 times as viscous as the lower layer. Fig. 57 shows that adding a crust to the model causes the maximum topography to both decrease and occur at a later time. Heat flows through the system primarily by advection until the temperature anomaly reaches the base of the crustal layer (approximately at the time for which Fig. 57 shows the temperature profile). The crust’s high viscosity reduces the temperature anomaly’s velocity substantially, causing it to affect the surface topography at a later time. Just as the cookbook shown in Using a free surface, the topography returns to zero after some time.