1. Thermal convection
Download
Here you can download the complete parameter file and the matplotlib scriptfor the animation.
Parameter file line by line
The parameter file contains all parameters that can be set. Not all of them are needed for every setting, therefore here only the active ones are highlighted.
Output files settings
We start with the directory name
# --- Name of the directory with results ---
name = "convection"
The directory protection will be turned off, meaning that when lauched repeatedly, the results will be overwritten each time
# Protection from overwriting the directory above
protect_directory = False
We choose the time units millions of years and output every 250 kyr.
# --- In what units the time will be? ---
# 1.0 - seconds or nondimensional, or yr, kyr or Myr
time_units = Myr
# --- String representation of time_units ---
time_units_string = "Myr"
# --- Output method for Paraview, HDF5 and tracers ---
output_frequency = ["time", 250*kyr]
As there will be no tracers, their saving will be turned off and the properties and their headers may be empty
# --- Save tracers into files? ---
save_tracers = False
# --- What properties to save? Must be one of the following (order does not matter):
Tracers_Output = []
# --- Headers for the columns in the text file for tracers
Tracers_header = []
To visualize the results in Paraview, we will save temperature and velocity. For initial condition, we will save only temperature.
# --- What functions to write into Paraview and HDF5 file ---
Paraview_Output = ["temperature", "velocity"]
Paraview_Output_Ini = ["temperature"]
To the statistics text file, time (in the units specified above, i.e., Myr) and the duration of the time step will be saved. We also choose corresponding headers.
# --- What values to print in a text file every time step? ---
stat_output = [ "time", "timestep"]
# --- Headers for the columns in the text file
stat_header = ["Time (h)", "dt (s)"]
Reloading options
As we define our initial condition and not reload, the following options need to be set.
The number at reload_HDF5_step and reload_tracers_step does not matter, but needs
to be defined.
# --- Name of the directory from which the data will be reloaded ---
reload_name = ""
# --- Whether or not to load HDF5 data from data_reload_name/HDF5/data.h5 ---
reload_HDF5 = False
# --- Time stamp of the HDF5 file which will be loaded ---
reload_HDF5_step = 0
# --- What functions to read from HDF5 file when reloading ---
reload_HDF5_functions = []
reload_tracers = False
reload_tracers_step = 0 # Second column in data_timestamp.dat
# --- Whether to reset time when reloading ---
restart_time = False
Simulation duration
The simulation will be stopped at 50 Myr.
# --- Criterion for ending the simulation, e.g. ["time", 1*Myr] or ["step", 1000] ---
termination_condition = ["time", 50*Myr]
Tracers options
For this case, anything related to tracers is irrelevant, therefore we can leave any values here
sm_thickness = 500.0 # m
om_thickness = 500.0 # m
integration_method = "RK4" # Euler, RK2, RK4
tracers_per_cell = 20
weight_tracers_by_ratio = False
Mesh geometry
The computational domain will be 100 km thick and 200 km long. The mesh will be created during the inicialization, using crossed geometry of the elements. The basic elements will have the longest side of 5 km and the mesh will be refined in the bottom half.
# --- Mesh height ---
height = 100e3 # m
# --- Mesh length ---
length = 200e3 # m
# --- Whether to read an external mesh ---
loading_mesh = False
# --- Mesh to be loaded ---
mesh_name = ""
# --- Basic mesh resolution if not loading mesh ---
z_div = 20
# --- Number of nodes in horizontal direction ---
x_div = int(z_div*(length/height)) # (keeps aspect ratio 1)
# --- Method of dividing basic squares into mesh triangle elements. ---
triangle_types = "crossed" # crossed, left, right, left/right, right/left
# --- Rectangular mesh refinement ---
refinement = [0, length, 0, height/2]
Material composition
The ice will consist of pure ice, therefore the list materials will remain empty.
Other options can be left as they are.
materials = []
empty_cells_allowed = False
empty_cells_composition = 0
empty_cells_region = []
Stokes problem settings
For the discretization of the velocity, we choose Mini elements The time step will be computed from the highest velocity and smallest element side in the domain, and assigned after solving the Stokes problem. The CFL coefficient will be 0.5. Since the rheology will be Newtonian (only temperature-dependent viscosity), the variables related to Picard iterations can be disregarded.
stokes_elements = "Mini"
time_step_position = "stokes" #stokes (right after Stokes problem) or "end" (at the end of the time loop)
time_step_strategy = "domain"
# --- The CFL parameter ---
cfl = 0.5
# --- Maximum number of Stokes solver Picard iterations ---
Picard_iter_max = 25
# --- Minimum relative error in velocity field ---
Picard_iter_error = 1e-3
error_type = "maximum" # "maximum or integrated"
For the boundary conditions we prescribe free slip everywhere. The mesh will not be moving, therefore the mesh displacement settings do not matter.
# Boundary conditions for velocity (free_slip, no_slip, free surface, velocity, velocity_x, velocity_y)
BC_Stokes_problem = [["free_slip"], # top boundary (1)
["free_slip"], # bottom boundary (2)
["free_slip"], # left boundary (3)
["free_slip"]] # right boundary (4)
mesh_displacement_laplace = "full" #"full_laplace" or "z_only"
# Method for correcting velocity field in case of multiple free surface conditions
stokes_null = "boundary"
Thermal solver settings
Since we solve thermal convection, the heat transfer equation needs to be solved. The nonlinearity is rather weak (\(k\) and \(c_p\)) are temperature dependent, therefore linear solver can be used. At the surface and ice-water boudnary, we prescribe temperature, while the lateral boundaries are adiabatic. The reference temperature is set to 273 K. Parameters related to radiative boundary condition will not be used.
# --- Solve heat transfer? ---
solve_energy_problem = True
# --- Whether the heat transfer should be solved with nonlinear solver ---
nonlinear_heat_equation = False
# --- Boundary condition for heat transfer equation ---
BC_heat_transfer = [["temp", 90.0], # top boundary (1)
["temp", 265.0], # bottom boundary (2)
["heat_flux", 0.0], # left boundary (3)
["heat_flux", 0.0]] # right boundary (4)
# --- Reference temperature ---
T_ref = 273.0 # K
# --- Insolation parametes ---
emis = 0.97 # Ice emissivity
albedo = 0.67
dist_AU = 5.2 # Object's distance from Sun in AU
insolation = insol_1AU/dist_AU**2
In the absence of tidal heating, the melting is not expected to occur. As an initial condition, a conductive profile will be solved and perturbed by a cosine wave.
# --- Melting inside the shell ---
# --- Whether generate partial melt if the temperature of the solid reaches T_melt ---
internal_melting = False
# --- Melting temperature of the solid ---
T_melt = 270.0 # K
# --- Tidal heating ---
tidal_dissipation = False
initial_tidal_dissipation = False
heating_model = "Maxwell" # Maxwell, Andrade or none
H_max = 4e-6 # W m^{-3}
# Andrade parameters
alpha_and = 0.2
# --- Find conductive initial condition ---
init_cond_profile = True
# --- Cosine perturbation of initial temperature ---
cos_perturbation = True
# --- Thermal perturbation amplitude ---
perturb_ampl = 5 # K
# --Number of half-cosine waves in lateral direction ---
perturb_freq = 1
Rheology
To mimic Titan’s ice shell, we prescribe corresponding gravity and ice physical properties.
Since the density is directly given by a temperature-dependent formula
in m_material_properties.py, the thermal expansivity here is inactive.
# --- Gravity acceleratuin ---
g = 1.3 # m s^-2
# --- Density of the domain solid ---
rho_s = 920.0 # kg/m3
# --- Density of the liquid below the domain ---
rho_l = 1000.0 # kg/m3
# --- Density of the melt produced by the solid ---
rho_m = rho_l # kg/m3
alpha_exp = 1.6e-4 # K^-1
Solving for a viscous flow only, the plasticity and elasticity will not be taken account. Since there is a free slip boundary condition at the bottom boundary, the phase transition and its parameters are inactive. The top surface is immobile as well, therefore adaptive smoothing is inactive.
# --- Rheology ---
elasticity = False
plasticity = False
# --- Phase transition at the bottom boundary ----
phase_transition = False
# --- Strength of the phase transition at the ice-water boundary ---
DAL_factor = 1e-1 # W/m3
# --- Latent heat of the material ---
Lt = 334.0e3 # J/kg
# --- Adaptive topography diffusion and topography diffusion factor ---
adaptive_smoothing = False
topo_diff = 1e-8*time_scaling
The viscosity will be onpy temperature-dependent, with viscosity at melting temperature \(10^{15}\) Pa s, activation energy 60 kJ/mol, and cut-off value \(10^{23}\) Pa s. The grain size and parameters for plasticity are inactive.
# --- VISCOSITY ---
viscosity_type = "temp-dep" # constant, temp-dep, GK_2001 (Goldsby and Kohlstedt, 2001) or composition
# --- Parameters for constant / temperature-dependent viscosity ---
eta_0 = 1e15 # Pa.s
# --- Activation energy ---
Q_activ = 60e3 # J/mol
# --- Grain size ---
d_grain = 1.0e-3
# --- Upper cut-off viscosity ---
eta_max = 1e23
# --- Lower cut-off value for plastic viscosity ---
eta_min_plast = eta_max/1e6
# # --- Elasticity ---
# If elasticity is off, we can compute the new viscosity directly from the strain rate (it will be faster)
# However, the stress formula + VEP iterations would work too
stress_iter_error = 1e-4
# --- Plasticity ---
# --- Angle of internal friction in degrees ---
int_friction_angle = 16.0
int_friction_angle2 = 0.0
# --- Cohesion of an undamaged material ---
cohesion_strong = 1e6
# --- Cohesion of a fully damaged material ---
cohesion_weak = 0.0
yield_stress_max, yield_stress_min = 1e8, 0.1
# --- Value of the plastic strain at which the material starts to accumulate damage ---
eps_strong = 0
# --- Critical value of the plastic strain beyond which the material is considered as fully damaged ---
eps_weak = 0.2
# --- Whether the accumulated plastic strain decreases in time ---
healing = False
# --- Characteristic time scale for the microscopic healing processes ---
healing_timescale = 300*kyr