diff --git a/ProcessLib/ThermoMechanics/Tests.cmake b/ProcessLib/ThermoMechanics/Tests.cmake
index 42980f4ea455fa50c5114b8302525a105a65a09c..46b4e57ca20f72aa4bfaa65eb82dd0f6b3d25744 100644
--- a/ProcessLib/ThermoMechanics/Tests.cmake
+++ b/ProcessLib/ThermoMechanics/Tests.cmake
@@ -102,3 +102,17 @@ AddTest(
DIFF_DATA
SimpleAxisymmetricCreepWithAnalyticSolution.vtu SimpleAxisymmetricCreepWithAnalyticalSolution_pcs_0_ts_1000_t_100.000000.vtu analytic_strain_yy epsilon_yy 1e-10 1e-4
)
+
+AddTest(
+ NAME ThermoMechanics_CreepAfterExcavation
+ PATH ThermoMechanics/CreepBGRa/CreepAfterExcavation
+ EXECUTABLE ogs
+ EXECUTABLE_ARGS CreepAfterExcavation.prj
+ WRAPPER time
+ TESTER vtkdiff
+ REQUIREMENTS NOT (OGS_USE_LIS OR OGS_USE_MPI)
+ DIFF_DATA
+ ExpectedCreepAfterExcavation_pcs_0_ts_61_t_4320000.000000.vtu CreepAfterExcavation_pcs_0_ts_61_t_4320000.000000.vtu sigma_xx sigma_xx 1e-16 1e-11
+ ExpectedCreepAfterExcavation_pcs_0_ts_61_t_4320000.000000.vtu CreepAfterExcavation_pcs_0_ts_61_t_4320000.000000.vtu sigma_yy sigma_yy 1e-16 1e-11
+ ExpectedCreepAfterExcavation_pcs_0_ts_61_t_4320000.000000.vtu CreepAfterExcavation_pcs_0_ts_61_t_4320000.000000.vtu displacement displacement 1e-16 1e-9
+)
diff --git a/Tests/Data/ThermoMechanics/CreepBGRa/CreepAfterExcavation/CreepAfterExcavation.prj b/Tests/Data/ThermoMechanics/CreepBGRa/CreepAfterExcavation/CreepAfterExcavation.prj
new file mode 100644
index 0000000000000000000000000000000000000000..de8f4d0112986216915464022f995436f04c3d3e
--- /dev/null
+++ b/Tests/Data/ThermoMechanics/CreepBGRa/CreepAfterExcavation/CreepAfterExcavation.prj
@@ -0,0 +1,262 @@
+<?xml version="1.0" encoding="ISO-8859-1"?>
+<OpenGeoSysProject>
+ <mesh>mesh.vtu</mesh>
+ <geometry>geometry.gml</geometry>
+ <search_length_algorithm>
+ <type>fixed</type>
+ <value>1e-1</value>
+ </search_length_algorithm>
+ <processes>
+ <process>
+ <name>ThermoMechanics</name>
+ <type>THERMO_MECHANICS</type>
+ <integration_order>2</integration_order>
+ <constitutive_relation>
+ <type>CreepBGRa</type>
+ <youngs_modulus>E</youngs_modulus>
+ <poissons_ratio>nu</poissons_ratio>
+ <a>2.0833333333333334e-06</a>
+ <n>4.9</n>
+ <sigma0>1e6</sigma0>
+ <q>54000</q>
+ <nonlinear_solver>
+ <maximum_iterations>1000</maximum_iterations>
+ <error_tolerance>2e-8</error_tolerance>
+ </nonlinear_solver>
+ </constitutive_relation>
+ <reference_solid_density>rho_sr</reference_solid_density>
+ <linear_thermal_expansion_coefficient>alpha</linear_thermal_expansion_coefficient>
+ <specific_heat_capacity>cs</specific_heat_capacity>
+ <thermal_conductivity>lambda</thermal_conductivity>
+ <process_variables>
+ <displacement>displacement</displacement>
+ <temperature>temperature</temperature>
+ </process_variables>
+ <secondary_variables>
+ <secondary_variable type="static" internal_name="sigma_xx" output_name="sigma_xx"/>
+ <secondary_variable type="static" internal_name="sigma_yy" output_name="sigma_yy"/>
+ <secondary_variable type="static" internal_name="sigma_zz" output_name="sigma_zz"/>
+ <secondary_variable type="static" internal_name="sigma_xy" output_name="sigma_xy"/>
+ <secondary_variable type="static" internal_name="epsilon_xx" output_name="epsilon_xx"/>
+ <secondary_variable type="static" internal_name="epsilon_yy" output_name="epsilon_yy"/>
+ <secondary_variable type="static" internal_name="epsilon_zz" output_name="epsilon_zz"/>
+ <secondary_variable type="static" internal_name="epsilon_xy" output_name="epsilon_xy"/>
+ </secondary_variables>
+ <specific_body_force>0 -9.81</specific_body_force>
+ </process>
+ </processes>
+ <time_loop>
+ <processes>
+ <process ref="ThermoMechanics">
+ <nonlinear_solver>basic_newton</nonlinear_solver>
+ <convergence_criterion>
+ <type>DeltaX</type>
+ <norm_type>NORM2</norm_type>
+ <abstol>1e-11</abstol>
+ </convergence_criterion>
+ <time_discretization>
+ <type>BackwardEuler</type>
+ </time_discretization>
+ <output>
+ <variables>
+ <variable>displacement</variable>
+ <variable>temperature</variable>
+ <variable>sigma_xx</variable>
+ <variable>sigma_yy</variable>
+ <variable>sigma_zz</variable>
+ <variable>sigma_xy</variable>
+ <variable>epsilon_xx</variable>
+ <variable>epsilon_yy</variable>
+ <variable>epsilon_zz</variable>
+ <variable>epsilon_xy</variable>
+ </variables>
+ </output>
+ <time_stepping>
+ <type>FixedTimeStepping</type>
+ <t_initial>0</t_initial>
+ <t_end>43200e+2</t_end>
+ <timesteps>
+ <pair>
+ <repeat>1</repeat>
+ <delta_t>86400.0e-3</delta_t>
+ </pair>
+ <pair>
+ <repeat>10</repeat>
+ <delta_t>86400.0e-2</delta_t>
+ </pair>
+ <pair>
+ <repeat>500</repeat>
+ <delta_t>86400.0</delta_t>
+ </pair>
+ </timesteps>
+ </time_stepping>
+ </process>
+ </processes>
+ <output>
+ <type>VTK</type>
+ <prefix>CreepAfterExcavation</prefix>
+ <timesteps>
+ <pair>
+ <repeat>10</repeat>
+ <each_steps>2</each_steps>
+ </pair>
+ <pair>
+ <repeat>1</repeat>
+ <each_steps>10</each_steps>
+ </pair>
+ </timesteps>
+ </output>
+ </time_loop>
+
+ <parameters>
+ <parameter>
+ <name>E</name>
+ <type>MeshElement</type>
+ <field_name>E</field_name>
+ </parameter>
+ <parameter>
+ <name>nu</name>
+ <type>MeshElement</type>
+ <field_name>nu</field_name>
+ </parameter>
+ <parameter>
+ <name>rho_sr</name>
+ <type>MeshElement</type>
+ <field_name>rho</field_name>
+ </parameter>
+ <parameter>
+ <name>alpha</name>
+ <type>Constant</type>
+ <value>4e-5</value>
+ </parameter>
+ <parameter>
+ <name>cs</name>
+ <type>Constant</type>
+ <value>0.</value>
+ </parameter>
+ <parameter>
+ <name>lambda</name>
+ <type>Constant</type>
+ <value>5</value>
+ </parameter>
+ <parameter>
+ <name>temperature_ic</name>
+ <type>MeshNode</type>
+ <field_name>T0</field_name>
+ </parameter>
+ <parameter>
+ <name>temperature_t</name>
+ <type>Constant</type>
+ <value>310</value>
+ </parameter>
+ <parameter>
+ <name>temperature_circ</name>
+ <type>Constant</type>
+ <value>300</value>
+ </parameter>
+ <parameter>
+ <name>temperature_b</name>
+ <type>Constant</type>
+ <value>320</value>
+ </parameter>
+ <parameter>
+ <name>displacement0</name>
+ <type>Constant</type>
+ <values>0 0</values>
+ </parameter>
+ <parameter>
+ <name>u_bc</name>
+ <type>Constant</type>
+ <value>0.</value>
+ </parameter>
+ <parameter>
+ <name>s_top</name>
+ <type>Constant</type>
+ <value>-1.0e+7</value>
+ </parameter>
+ </parameters>
+ <process_variables>
+ <process_variable>
+ <name>displacement</name>
+ <components>2</components>
+ <order>1</order>
+ <initial_condition>displacement0</initial_condition>
+ <boundary_conditions>
+ <boundary_condition>
+ <geometrical_set>geo</geometrical_set>
+ <geometry>left</geometry>
+ <type>Dirichlet</type>
+ <component>0</component>
+ <parameter>u_bc</parameter>
+ </boundary_condition>
+ <boundary_condition>
+ <geometrical_set>geo</geometrical_set>
+ <geometry>right</geometry>
+ <type>Dirichlet</type>
+ <component>0</component>
+ <parameter>u_bc</parameter>
+ </boundary_condition>
+ <boundary_condition>
+ <geometrical_set>geo</geometrical_set>
+ <geometry>bottom</geometry>
+ <type>Dirichlet</type>
+ <component>1</component>
+ <parameter>u_bc</parameter>
+ </boundary_condition>
+ <boundary_condition>
+ <geometrical_set>geo</geometrical_set>
+ <geometry>top</geometry>
+ <type>Neumann</type>
+ <component>1</component>
+ <parameter>s_top</parameter>
+ </boundary_condition>
+ </boundary_conditions>
+ </process_variable>
+ <process_variable>
+ <name>temperature</name>
+ <components>1</components>
+ <order>1</order>
+ <initial_condition>temperature_ic</initial_condition>
+ <boundary_conditions>
+ <boundary_condition>
+ <geometrical_set>geo</geometrical_set>
+ <geometry>top</geometry>
+ <type>Dirichlet</type>
+ <parameter>temperature_t</parameter>
+ </boundary_condition>
+ <boundary_condition>
+ <geometrical_set>geo</geometrical_set>
+ <geometry>bottom</geometry>
+ <type>Dirichlet</type>
+ <parameter>temperature_b</parameter>
+ </boundary_condition>
+ <boundary_condition>
+ <geometrical_set>geo</geometrical_set>
+ <geometry>circle</geometry>
+ <type>Dirichlet</type>
+ <parameter>temperature_circ</parameter>
+ </boundary_condition>
+ </boundary_conditions>
+ </process_variable>
+ </process_variables>
+ <nonlinear_solvers>
+ <nonlinear_solver>
+ <name>basic_newton</name>
+ <type>Newton</type>
+ <max_iter>100</max_iter>
+ <linear_solver>general_linear_solver</linear_solver>
+ </nonlinear_solver>
+ </nonlinear_solvers>
+ <linear_solvers>
+ <linear_solver>
+ <name>general_linear_solver</name>
+ <lis>-i bicgstab -p jacobi -tol 1e-11 -maxiter 10000</lis>
+ <eigen>
+ <solver_type>BiCGSTAB</solver_type>
+ <precon_type>DIAGONAL</precon_type>
+ <max_iteration_step>10000</max_iteration_step>
+ <error_tolerance>1e-16</error_tolerance>
+ </eigen>
+ </linear_solver>
+ </linear_solvers>
+</OpenGeoSysProject>
diff --git a/Tests/Data/ThermoMechanics/CreepBGRa/CreepAfterExcavation/ExpectedCreepAfterExcavation_pcs_0_ts_61_t_4320000.000000.vtu b/Tests/Data/ThermoMechanics/CreepBGRa/CreepAfterExcavation/ExpectedCreepAfterExcavation_pcs_0_ts_61_t_4320000.000000.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..0a9e194a5b667c2710e83779b8d4547f44bc4e51
--- /dev/null
+++ b/Tests/Data/ThermoMechanics/CreepBGRa/CreepAfterExcavation/ExpectedCreepAfterExcavation_pcs_0_ts_61_t_4320000.000000.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:8e5c6736ac05a9b66ca8bc4f9d771453bbfabe1530b2465353019e64d4fadae3
+size 185031
diff --git a/Tests/Data/ThermoMechanics/CreepBGRa/CreepAfterExcavation/geometry.gml b/Tests/Data/ThermoMechanics/CreepBGRa/CreepAfterExcavation/geometry.gml
new file mode 100644
index 0000000000000000000000000000000000000000..f3afc61452c2a2d5fff4bc838560358eb8ed6f4c
--- /dev/null
+++ b/Tests/Data/ThermoMechanics/CreepBGRa/CreepAfterExcavation/geometry.gml
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:81381efca2c71882c910c6c489467b8f7c70dfb267cd454c7d1b2a665221bcd0
+size 7944
diff --git a/Tests/Data/ThermoMechanics/CreepBGRa/CreepAfterExcavation/mesh.vtu b/Tests/Data/ThermoMechanics/CreepBGRa/CreepAfterExcavation/mesh.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..87f1d669d17a5abfb739cfe324dfe14dcf54b3d9
--- /dev/null
+++ b/Tests/Data/ThermoMechanics/CreepBGRa/CreepAfterExcavation/mesh.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:9e6693047520f77b89ff0f43e63229fe54efaed674290d93c1ad38bf05022196
+size 176019
diff --git a/web/content/docs/benchmarks/creep-after-excavation-bgra/CreepAfterExcavation.md b/web/content/docs/benchmarks/creep-after-excavation-bgra/CreepAfterExcavation.md
new file mode 100644
index 0000000000000000000000000000000000000000..5b72eacec498bccc4c3811ffccd308a415118ca6
--- /dev/null
+++ b/web/content/docs/benchmarks/creep-after-excavation-bgra/CreepAfterExcavation.md
@@ -0,0 +1,127 @@
++++
+date = "2018-08-10T11:46:45+01:00"
+title = "Creep in a rectangle domain with hole under pressure on its top surface"
+weight = 171
+project = "ThermoMechanics/CreepBGRa/CreepAfterExcavation/CreepAfterExcavation.prj"
+author = "Wenqing Wang"
+
+[menu]
+ [menu.benchmarks]
+ parent = "thermo-mechanics"
+
++++
+
+{{< data-link >}}
+
+-------------------------------------------------------------------------
+Following up the benchmark about the BGRa creep model described on the
+[this page](https://www.opengeosys.org/docs/benchmarks/creepbgra/creepbrga/),
+this example represents the creep in the near field of
+drift in the deep rock salt after excavation. The domain and the
+geometry are shown in the following figure:
+<figure>
+ <img src="../mesh.png" alt="Mesh and Geometry" id="fig_6"
+ style="height:400px;width:490px;">
+ <figcaption>Mesh and Geometry.</figcaption>
+</figure>
+
+The domain has two material groups, which are highlighted by different
+colors. The material group that is in the
+top part of the domain represents a cap rock type, while the other
+material group is for rock salt. The material properties of the rocks are given in the following table:
+
+---------------------- ---------- ----------- ------------
+| |Cap rock | Rock salt| Unit|
+| --- |:---------:| -----:|----:|
+|Density | 2000 | 2170 | kg/m$^{3}$|
+|Young’s Modulus | 7.0 | 7.65 | GPa |
+|Poisson ratio | 0.3 | 2.7 | - |
+|Thermal conductivity | 5 | 5 | W/(mK) |
+---------------------- ---------- ----------- ------------
+
+The parameters of the BGRa creep model are $A=0.18\, \mbox{d}^{-1}$,
+$m=5$, $Q=54 \mbox{ kJ/mol}$.
+
+
+The width
+and the height of of the domain are 300 m and 340 m, respectively. The
+height of the cap rock portion is 40 m. The drift to be excavated has a
+radius of 50 m.
+
+Here we consider the creep of the rock after excavation. Therefore, we
+assume a steady state of stress and temperature of excavation at the
+beginning of the current simulation. For this assumption, the boundary
+conditions are given as:
+
+- top boundary: $\tau_x = \sigma_x=0$, $\tau_y=\sigma_y=70$ MPa,
+ $T=310$ K.
+
+- two lateral boundaries: normal displacement is fixed, and no heat
+ flux.
+
+- bottom boundary: normal displacement is fixed, and $T=320$ K.
+
+- circle of drift surface: traction free
+ (${ \mathbf\sigma}\cdot \mathbf n = 0$), and $T=300$ K as excavation
+ conditions.
+
+The boundary conditions of the mechanical process lead to a distinct
+stress field after excavation. Besides, heat capacity is neglected in
+order to pose a steady state temperature field after excavation. The
+initial stresses are zero and the initial temperature is linearly
+distributed from top to bottom with the boundary values of temperature.
+
+The time step sizes of the simulation are: One step of 0.001 day, 10 steps of 0.1 day, and the
+remaining steps of 1 day.
+
+The following three figures are plotted by using the results of the simulation of 1000 days creep.
+The three figures display the distribution of horizontal and vertical stresses at times of
+108 days, 409 days and 1000 days, respectively.
+
+In theses three figures,
+the left sub-figure show the time variations of horizontal and vertical stresses
+ at a position just close to the top of the drift, and the
+green vertical line in it marks the time of the displayed stress field.
+
+<figure>
+ <img src="../stress_xx_yy_20.png" alt="Stress distribution at the time of 109 days." id="fig_2">
+ <figcaption>Stress distribution at the time of 109 days.</figcaption>
+</figure>
+
+<figure>
+ <img src="../stress_xx_yy_50.png" alt="Stress distribution at the time of 409 days." id="fig_3">
+ <figcaption>Stress distribution at the time of 409 days.</figcaption>
+</figure>
+
+<figure>
+ <img src="../stress_xx_yy_110.png" alt="Stress distribution at the time of 1000 days." id="fig_4">
+ <figcaption>Stress distribution at the time of 1000 days.</figcaption>
+</figure>
+
+
+The above three figures show that the absolute value of the horizontal stress increase
+signification in the areas above and beneath the drift due to the creep,
+and the change of the vertical stress is slow compared to that of the
+horizontal stress.
+
+The following figure shows the strain distribution at the
+end of the simulation time at 1000 days.
+<figure>
+ <img src="../strain.png" alt="Strain distribution at the time of 1000 days." id="fig_5">
+ <figcaption>Strain distribution at the time of 1000 days.</figcaption>
+</figure>
+
+
+The steady-state temperature distribution is displayed in the following figure
+<figure>
+ <img src="../T.png" alt="Temperature distribution at the time of 1000 days." id="fig_6">
+ <figcaption>Temperature distribution at the time of 1000 days.</figcaption>
+</figure>
+
+## Note:
+For the automatic benchmarking, the time duration of creep is reduced to 50 days in order to reduce the run time.
+ If one wants to test this benchmark for 1000 days' creep, please change the end time in the tag of `<t_end> `
+in the project file as
+ ```
+<t_end>1000</t_end>
+```
\ No newline at end of file
diff --git a/web/content/docs/benchmarks/creep-after-excavation-bgra/T.png b/web/content/docs/benchmarks/creep-after-excavation-bgra/T.png
new file mode 100644
index 0000000000000000000000000000000000000000..3ad5f8b981098fb74e5dfda508f009763501d546
--- /dev/null
+++ b/web/content/docs/benchmarks/creep-after-excavation-bgra/T.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:a360404be4d5764b85c0bead5ba92d92faaa9542a12957f72d9d706b4b4e39c6
+size 384898
diff --git a/web/content/docs/benchmarks/creep-after-excavation-bgra/mesh.png b/web/content/docs/benchmarks/creep-after-excavation-bgra/mesh.png
new file mode 100644
index 0000000000000000000000000000000000000000..9b4029da7e61d2fac8e848b5ece92b690ec3a93e
--- /dev/null
+++ b/web/content/docs/benchmarks/creep-after-excavation-bgra/mesh.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:c2095fc571a43c62a3eb1816d8aaf0ead9f4c66f1974a544fda5ced0d2beed49
+size 326201
diff --git a/web/content/docs/benchmarks/creep-after-excavation-bgra/strain.png b/web/content/docs/benchmarks/creep-after-excavation-bgra/strain.png
new file mode 100644
index 0000000000000000000000000000000000000000..f0e960169e3b6baa242ca08e5f069fd020fd12b5
--- /dev/null
+++ b/web/content/docs/benchmarks/creep-after-excavation-bgra/strain.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:ff2fb210d6f7d543082f490c55ab265e27095efef30fe60f6dd596df0d87ce67
+size 153593
diff --git a/web/content/docs/benchmarks/creep-after-excavation-bgra/stress_xx_yy_110.png b/web/content/docs/benchmarks/creep-after-excavation-bgra/stress_xx_yy_110.png
new file mode 100644
index 0000000000000000000000000000000000000000..d192fbc4e577658358db571ff8bd1aba005aaff3
--- /dev/null
+++ b/web/content/docs/benchmarks/creep-after-excavation-bgra/stress_xx_yy_110.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:abde5d089d6edc8847dc425c8feec31515832449415bc6bb5f8a414d05dd632a
+size 870467
diff --git a/web/content/docs/benchmarks/creep-after-excavation-bgra/stress_xx_yy_20.png b/web/content/docs/benchmarks/creep-after-excavation-bgra/stress_xx_yy_20.png
new file mode 100644
index 0000000000000000000000000000000000000000..0323c62ccc22eae7c819e42ac8e5a474d32b2900
--- /dev/null
+++ b/web/content/docs/benchmarks/creep-after-excavation-bgra/stress_xx_yy_20.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:3ea312ed37bd6c593ca4b4a7f4a7b1d7d064c83af53bc70c675eadd2517adef2
+size 847877
diff --git a/web/content/docs/benchmarks/creep-after-excavation-bgra/stress_xx_yy_50.png b/web/content/docs/benchmarks/creep-after-excavation-bgra/stress_xx_yy_50.png
new file mode 100644
index 0000000000000000000000000000000000000000..e3cc98ff1c410a7007d97db00dda36ff8b16a6eb
--- /dev/null
+++ b/web/content/docs/benchmarks/creep-after-excavation-bgra/stress_xx_yy_50.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:f50ecaaeaaf49a8e6ac371571f55ab0278df5b27a15b43c015ebdf7827b93c32
+size 858276