# 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. ). 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. 178 illustrates the model setup.

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 modelling 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].