diff --git a/include/aspect/parameters.h b/include/aspect/parameters.h
index 201fb425efa..72568bc8b92 100644
--- a/include/aspect/parameters.h
+++ b/include/aspect/parameters.h
@@ -208,6 +208,7 @@ namespace aspect
           hydrostatic_compression,
           reference_density_profile,
           implicit_reference_density_profile,
+          fully_compressible,
           incompressible,
           projected_density_field,
           ask_material_model
@@ -228,6 +229,8 @@ namespace aspect
             return Formulation::MassConservation::reference_density_profile;
           else if (input == "implicit reference density profile")
             return Formulation::MassConservation::implicit_reference_density_profile;
+          else if (input == "fully compressible")
+            return Formulation::MassConservation::fully_compressible;
           else if (input == "incompressible")
             return Formulation::MassConservation::incompressible;
           else if (input == "projected density field")
diff --git a/include/aspect/simulator/assemblers/interface.h b/include/aspect/simulator/assemblers/interface.h
index fa918170c0e..962b56e1c6a 100644
--- a/include/aspect/simulator/assemblers/interface.h
+++ b/include/aspect/simulator/assemblers/interface.h
@@ -174,7 +174,8 @@ namespace aspect
 
           std::vector<Tensor<1,dim>>          phi_u;
           std::vector<Tensor<1,dim>>          velocity_values;
-          std::vector<double>                  velocity_divergence;
+          std::vector<double>                 velocity_divergence;
+          std::vector<double>                 old_pressure_values;
           std::vector<Tensor<1,dim>>          temperature_gradients;
 
           /**
diff --git a/include/aspect/simulator/assemblers/stokes.h b/include/aspect/simulator/assemblers/stokes.h
index 6ffe27fd8f3..1c3dec2af58 100644
--- a/include/aspect/simulator/assemblers/stokes.h
+++ b/include/aspect/simulator/assemblers/stokes.h
@@ -89,6 +89,22 @@ namespace aspect
                 internal::Assembly::CopyData::CopyDataBase<dim> &data_base) const override;
     };
 
+    /**
+     * This class assembles the compressibility term of the Stokes equation
+     * that is caused by the compressibility in the mass conservation equation
+     * without approximation. This term is expressed as
+     * $-\nabla\cdot\mathbf{u} = \frac{\beta(p - p^{old})}{\Delta t}$.
+     */
+    template <int dim>
+    class StokesCompressibleMassConservationTerm : public Assemblers::Interface<dim>,
+      public SimulatorAccess<dim>
+    {
+      public:
+        void
+        execute(internal::Assembly::Scratch::ScratchBase<dim>   &scratch_base,
+                internal::Assembly::CopyData::CopyDataBase<dim> &data_base) const override;
+    };
+
     /**
      * This class assembles the compressibility term of the Stokes equation
      * that is caused by the compressibility in the mass conservation equation.
diff --git a/source/simulator/assemblers/interface.cc b/source/simulator/assemblers/interface.cc
index b824393e2c3..2fc1ec2f4df 100644
--- a/source/simulator/assemblers/interface.cc
+++ b/source/simulator/assemblers/interface.cc
@@ -139,6 +139,7 @@ namespace aspect
           phi_u (stokes_dofs_per_cell, numbers::signaling_nan<Tensor<1,dim>>()),
           velocity_values (quadrature.size(), numbers::signaling_nan<Tensor<1,dim>>()),
           velocity_divergence(quadrature.size(), numbers::signaling_nan<double>()),
+          old_pressure_values(quadrature.size(), numbers::signaling_nan<double>()),
           temperature_gradients (quadrature.size(), numbers::signaling_nan<Tensor<1,dim>>()),
           face_material_model_inputs(face_quadrature.size(), n_compositional_fields),
           face_material_model_outputs(face_quadrature.size(), n_compositional_fields),
@@ -163,6 +164,7 @@ namespace aspect
           phi_u (scratch.phi_u),
           velocity_values (scratch.velocity_values),
           velocity_divergence (scratch.velocity_divergence),
+          old_pressure_values (scratch.old_pressure_values),
           temperature_gradients (scratch.temperature_gradients),
           face_material_model_inputs(scratch.face_material_model_inputs),
           face_material_model_outputs(scratch.face_material_model_outputs),
diff --git a/source/simulator/assemblers/stokes.cc b/source/simulator/assemblers/stokes.cc
index 164a7a7913e..92912a4a215 100644
--- a/source/simulator/assemblers/stokes.cc
+++ b/source/simulator/assemblers/stokes.cc
@@ -211,6 +211,26 @@ namespace aspect
                                               * JxW;
                   }
             }
+
+          // If the fully compressible formulation is chosen, then the (P,P)
+          // block of the system matrix is not zero, which should also be
+          // added to the preconditioner matrix (H. C. Elman et al., 2014,
+          // Eq.[8.11]).
+          if (this->get_parameters().formulation_mass_conservation ==
+              Parameters<dim>::Formulation::MassConservation::fully_compressible)
+            {
+              const double beta_over_dt
+                = this->get_timestep_number() > 0 ?
+                  scratch.material_model_outputs.compressibilities[q] / this->get_timestep() :
+                  0.0;
+
+              for (unsigned int i=0; i<stokes_dofs_per_cell; ++i)
+                for (unsigned int j=0; j<stokes_dofs_per_cell; ++j)
+                  data.local_matrix(i, j) += ( beta_over_dt *
+                                               pressure_scaling * pressure_scaling *
+                                               (scratch.phi_p[i] * scratch.phi_p[j])
+                                             ) * JxW;
+            }
         }
     }
 
@@ -568,6 +588,75 @@ namespace aspect
 
 
 
+    template <int dim>
+    void
+    StokesCompressibleMassConservationTerm<dim>::
+    execute (internal::Assembly::Scratch::ScratchBase<dim>   &scratch_base,
+             internal::Assembly::CopyData::CopyDataBase<dim> &data_base) const
+    {
+      Assert(this->get_parameters().formulation_mass_conservation ==
+             Parameters<dim>::Formulation::MassConservation::fully_compressible,
+             ExcInternalError());
+
+      internal::Assembly::Scratch::StokesSystem<dim> &scratch = dynamic_cast<internal::Assembly::Scratch::StokesSystem<dim>&> (scratch_base);
+      internal::Assembly::CopyData::StokesSystem<dim> &data = dynamic_cast<internal::Assembly::CopyData::StokesSystem<dim>&> (data_base);
+
+      if (!scratch.rebuild_stokes_matrix)
+        return;
+
+      const Introspection<dim> &introspection = this->introspection();
+      const FiniteElement<dim> &fe = this->get_fe();
+      const unsigned int stokes_dofs_per_cell = data.local_dof_indices.size();
+      const unsigned int n_q_points = scratch.finite_element_values.n_quadrature_points;
+      const double pressure_scaling = this->get_pressure_scaling();
+      const double time_step = this->get_timestep();
+
+      for (unsigned int q=0; q<n_q_points; ++q)
+        {
+          for (unsigned int i=0, i_stokes=0; i_stokes<stokes_dofs_per_cell; /*increment at end of loop*/)
+            {
+              if (introspection.is_stokes_component(fe.system_to_component_index(i).first))
+                {
+                  scratch.phi_p[i_stokes] = scratch.finite_element_values[introspection.extractors.pressure].value(i, q);
+                  ++i_stokes;
+                }
+              ++i;
+            }
+
+          const double beta_over_dt
+            = this->get_timestep_number() > 0 ?
+              scratch.material_model_outputs.compressibilities[q] / time_step :
+              0.0;
+
+          const double thermal_alpha
+            = scratch.material_model_outputs.thermal_expansion_coefficients[q];
+
+          const double JxW = scratch.finite_element_values.JxW(q);
+
+          for (unsigned int i=0; i<stokes_dofs_per_cell; ++i)
+            {
+              data.local_rhs(i) -= (
+                                     pressure_scaling *
+                                     (
+                                       // pressure part
+                                       scratch.old_pressure_values[q] * beta_over_dt
+                                       // thermal part
+                                       + thermal_alpha *
+                                       (scratch.velocity_values[q] * scratch.temperature_gradients[q])
+                                     ) * scratch.phi_p[i]
+                                   ) * JxW;
+
+              for (unsigned int j=0; j<stokes_dofs_per_cell; ++j)
+                data.local_matrix(i, j) -= ( beta_over_dt *
+                                             pressure_scaling * pressure_scaling *
+                                             (scratch.phi_p[i] * scratch.phi_p[j])
+                                           ) * JxW;
+            }
+        }
+    }
+
+
+
     template <int dim>
     void
     StokesReferenceDensityCompressibilityTerm<dim>::
@@ -979,6 +1068,7 @@ namespace aspect
   template class StokesCompressiblePreconditioner<dim>; \
   template class StokesIncompressibleTerms<dim>; \
   template class StokesCompressibleStrainRateViscosityTerm<dim>; \
+  template class StokesCompressibleMassConservationTerm<dim>; \
   template class StokesReferenceDensityCompressibilityTerm<dim>; \
   template class StokesImplicitReferenceDensityCompressibilityTerm<dim>; \
   template class StokesIsentropicCompressionTerm<dim>; \
diff --git a/source/simulator/assembly.cc b/source/simulator/assembly.cc
index 5142237c083..b28e71105a7 100644
--- a/source/simulator/assembly.cc
+++ b/source/simulator/assembly.cc
@@ -121,6 +121,12 @@ namespace aspect
         assemblers->stokes_system.push_back(
           std::make_unique<aspect::Assemblers::StokesProjectedDensityFieldTerm<dim>>());
       }
+    else if (parameters.formulation_mass_conservation ==
+             Parameters<dim>::Formulation::MassConservation::fully_compressible)
+      {
+        assemblers->stokes_system.push_back(
+          std::make_unique<aspect::Assemblers::StokesCompressibleMassConservationTerm<dim>>());
+      }
     else
       AssertThrow(false,
                   ExcMessage("Unknown mass conservation equation approximation. There is no assembler"
@@ -606,7 +612,12 @@ namespace aspect
     if (parameters.formulation_mass_conservation == Parameters<dim>::Formulation::MassConservation::hydrostatic_compression)
       scratch.finite_element_values[introspection.extractors.temperature].get_function_gradients(current_linearization_point,
           scratch.temperature_gradients);
-
+    if (parameters.formulation_mass_conservation == Parameters<dim>::Formulation::MassConservation::fully_compressible)
+      scratch.finite_element_values[introspection.extractors.pressure].get_function_values(old_solution, scratch.old_pressure_values);
+    if (parameters.formulation_mass_conservation == Parameters<dim>::Formulation::MassConservation::fully_compressible ||
+        parameters.formulation_mass_conservation == Parameters<dim>::Formulation::MassConservation::hydrostatic_compression)
+      scratch.finite_element_values[introspection.extractors.temperature].get_function_gradients(current_linearization_point,
+          scratch.temperature_gradients);
 
     const bool use_reference_density_profile = (parameters.formulation_mass_conservation == Parameters<dim>::Formulation::MassConservation::reference_density_profile)
                                                || (parameters.formulation_mass_conservation == Parameters<dim>::Formulation::MassConservation::implicit_reference_density_profile);
diff --git a/source/simulator/core.cc b/source/simulator/core.cc
index aaa31fbd920..81f4ccd025d 100644
--- a/source/simulator/core.cc
+++ b/source/simulator/core.cc
@@ -999,10 +999,13 @@ namespace aspect
                 coupling[x.pressure][x.velocities[d]] = DoFTools::always;
               }
 
-            // For equal-order interpolation, we need a stabilization term
-            // in the bottom right of Stokes matrix. Make sure we have the
-            // necessary entries.
-            if (parameters.use_equal_order_interpolation_for_stokes == true)
+            // For equal-order interpolation and/or fully compressible Stokes
+            // equation, the bottom right of Stokes matrix is nonzero.
+            // Make sure we have the necessary entries.
+            if (parameters.use_equal_order_interpolation_for_stokes == true
+                ||
+                parameters.formulation_mass_conservation ==
+                Parameters<dim>::Formulation::MassConservation::fully_compressible)
               coupling[x.pressure][x.pressure] = DoFTools::always;
           }
       }
diff --git a/source/simulator/parameters.cc b/source/simulator/parameters.cc
index 630301c8fce..1ec99281034 100644
--- a/source/simulator/parameters.cc
+++ b/source/simulator/parameters.cc
@@ -663,7 +663,7 @@ namespace aspect
                          ":::");
 
       prm.declare_entry ("Mass conservation", "ask material model",
-                         Patterns::Selection ("incompressible|isentropic compression|hydrostatic compression|"
+                         Patterns::Selection ("incompressible|fully compressible|isentropic compression|hydrostatic compression|"
                                               "reference density profile|implicit reference density profile|"
                                               "projected density field|"
                                               "ask material model"),