The solitary wave benchmark

The solitary wave benchmark#

This section was contributed by Juliane Dannberg and is based on a section in Dannberg and Heister [2016] by Juliane Dannberg and Timo Heister.

One of the most widely used benchmarks for codes that model the migration of melt through a compacting and dilating matrix is the propagation of solitary waves (e.g. [Keller et al., 2013, Schmeling, 2000, Simpson and Spiegelman, 2011]). The benchmark is intended to test the accuracy of the solution of the two-phase flow equations as described in Calculations with melt transport and makes use of the fact that there is an analytical solution for the shape of solitary waves that travel through a partially molten rock with a constant background porosity without changing their shape and with a constant wave speed. Here, we follow the setup of the benchmark as it is described in Barcilon and Richter [1986], which considers one-dimensional solitary waves.

The model features a perturbation of higher porosity with the amplitude \(A \phi_0\) in a uniform low-porosity (\(\phi=\phi_0\)) background. Due to its lower density, melt migrates upwards, dilating the solid matrix at its front and compacting it at its end.

Assuming constant shear and compaction viscosities and using a permeability law of the form

\[k_\phi = k_0 \phi^3,\]

implying a Darcy coefficient

\[K_D(\phi) = \frac{k_0}{\eta_f} \phi^3 ,\]

and the non-dimensionalization

\[\begin{split}\begin{aligned} x &= \delta x' && \text{ with the compaction length } \delta = \sqrt{K_D(\phi_0)(\xi + \frac{4}{3}\eta)} , \\ \phi &= \phi_0 \phi ' && \text{ with the background porosity } \phi_0 , \\ (\mathbf u_s, \mathbf u_f) &= u_0 (\mathbf u_s, \mathbf u_f)' && \text{ with the separation flux } \phi_0 u_0 = K_D(\phi_0) \Delta\rho g , \end{aligned}\end{split}\]

the analytical solution for the shape of the solitary wave can be written in implicit form as:

\[\begin{aligned} x(\phi) &= \pm (A + 0.5) \left[ -2 \sqrt{A-\phi} + \frac{1}{\sqrt{A-1}} \ln \frac{\sqrt{A-1} - \sqrt{A-\phi}}{\sqrt{A-1} + \sqrt{A-\phi}} \right] \end{aligned}\]

and the phase speed \(c\), scaled back to physical units, is \(c = u_0 (2A+1)\). This is only valid in the limit of small porosity \(\phi_0 \ll 1\). Fig. 208 illustrates the model setup.

../../../../../_images/setup.svg

Fig. 208 Setup of the solitary wave benchmark. The domain is 400 m high and includes a low porosity (ϕ = 0.001) background with an initial perturbation (ϕ = 0.1). The solid density is 3300 kg/m3 and the melt density is 2500 kg/m3. We apply the negative phase speed of the solitary wave us =  − cez as velocity boundary condition, so that the wave will stay at its original position while the background is moving.#

The parameter file and material model for this setup can be found in benchmarks/solitary_wave/solitary_wave.prm and benchmarks/solitary_wave/solitary_wave.cc. The most relevant sections are shown in the following paragraph.

# Listing of Parameters
# ---------------------
# Set up the solitary wave benchmark
# (Barcilon & Richter, 1986; Simpson & Spiegelman, 2011;
# Keller et al., 2013; Schmeling, 200

set Additional shared libraries            = ./libsolitary_wave.so

# A non-linear solver has to be used for models with melt migration
set Nonlinear solver scheme                = iterated Advection and Stokes
set Max nonlinear iterations               = 10
set Nonlinear solver tolerance             = 1e-5

# The end time is chosen in such a way that the solitary wave travels
# approximately 5 times its wavelength during the model time.
set End time                               = 6e6

# To model melt migration, there has to be a compositional field with
# the name 'porosity'.
subsection Compositional fields
  set Number of fields = 1
  set Names of fields = porosity
end

# Enable modeling of melt migration in addition to the advection of
# solid material.
subsection Melt settings
  set Include melt transport                  = true
end

######### Parameters for the porosity field ########################

# We use the initial conditions and material model from the
# solitary wave plugin and choose a wave with an amplitude of
# 0.01 and a background porosity of 0.001.
subsection Initial composition model
  set Model name = Solitary wave initial condition

  subsection Solitary wave initial condition
    set Offset = 200
    set Read solution from file = true
    set Amplitude = 0.01
    set Background porosity = 0.001
  end
end

subsection Material model
  set Model name = Solitary Wave
end

# As material is flowing in, we prescribe the porosity at the
# upper and lower boundary.
subsection Boundary composition model
  set List of model names = initial composition
  set Fixed composition boundary indicators   = 2,3
end

# As we know that our solution does not have any steep gradients
# we can use a low stabilization to avoid too much diffusion.
subsection Discretization
  subsection Stabilization parameters
    set beta  = 0.001
  end
end

######### Model geometry ##############################################

# Our domain is a pseudo-1D-profile 400 m in height, but only a few elements wide
subsection Geometry model
  set Model name = box

  subsection Box
    set X extent  = 10
    set Y extent  = 400
    set Y repetitions = 40
  end
end

######### Velocity boundary conditions ################################

# We apply the phase speed of the wave here, so that it always stays in the
# same place in our model. The phase speed is c = 5.25e-11 m/s, but we have
# to convert it to m/years using the same conversion that is used internally
# in ASPECT: year_in_seconds = 60*60*24*365.2425.
subsection Boundary velocity model
  set Tangential velocity boundary indicators = 0,1
  set Prescribed velocity boundary indicators = 2:function, 3:function

  subsection Function
    set Function expression = 0;-1.65673998e-4
  end
end

# Postprocessor for the error calculation
subsection Postprocess
  set List of postprocessors = solitary wave statistics
end

subsection Solver parameters
  subsection Stokes solver parameters
    set Linear solver tolerance = 1e-10
  end
end

The benchmark uses a custom model to generate the initial condition for the porosity field as specified by the analytical solution, and its own material model, which includes the additional material properties needed by models with melt migration, such as the permeability, melt density and compaction viscosity. The solitary wave postprocessor compares the porosity and pressure in the model to the analytical solution, and computes the errors for the shape of the porosity, shape of the compaction pressure and the phase speed. We apply the negative phase speed of the solitary wave as a boundary condition for the solid velocity. This changes the reference frame, so that the solitary wave stays in the center of the domain, while the solid moves downwards. The temperature evolution does not play a role in this benchmark, so all temperature and heating-related parameters are disabled or set to zero.

And extensive discussion of the results and convergence behavior can be found in Dannberg and Heister [2016].