diff --git a/Documentation/bibliography/other.bib b/Documentation/bibliography/other.bib
index f62df77bd36ee13f0866ce974acd5f6dde2abf8a..2b2ab112fc09a4bb2d56da3e3eae80967bb86ac6 100644
--- a/Documentation/bibliography/other.bib
+++ b/Documentation/bibliography/other.bib
@@ -165,3 +165,13 @@
   publisher={Springer}
 }
 
+@article{rutqvist2011implementation,
+  title={Implementation of the Barcelona Basic Model into TOUGH--FLAC for simulations of the geomechanical behavior of unsaturated soils},
+  author={Rutqvist, Jonny and Ijiri, Yuji and Yamamoto, Hajime},
+  journal={Computers \& Geosciences},
+  volume={37},
+  number={6},
+  pages={751--762},
+  year={2011},
+  publisher={Elsevier}
+}
diff --git a/MaterialLib/MPL/Properties/SwellingStress/LinearSaturationSwellingStress.h b/MaterialLib/MPL/Properties/SwellingStress/LinearSaturationSwellingStress.h
index 6e0f25350aa22f170582775c0266f10ee2f163f4..2837f5e807ed0a4a2e6119c7940f966c2d92e108 100644
--- a/MaterialLib/MPL/Properties/SwellingStress/LinearSaturationSwellingStress.h
+++ b/MaterialLib/MPL/Properties/SwellingStress/LinearSaturationSwellingStress.h
@@ -23,14 +23,17 @@ class Phase;
  *   Clay materials like bentonite have a high swelling capacity in dry state,
  * and their swelling property can be described by this model.
  *
- *   The model takes the form
+ *  The original model was proposed in \cite rutqvist2011implementation
+ *  (equations (39) and (40) on pages 758--759). With a simplification of the
+ *  parameters of the original formula and introducing a constraint to avoid
+ *  shrinkage stress when saturation is below the initial saturation, the model
+ *  takes the form
  *    \f[  {\mathbf{\sigma}}^{\text{sw}} =
  * {\alpha}_{\text{sw}}  (S-S_0) \mathbf{I}, \, \forall S \in
  * [S_0, S_\text{max}] \f]
  *  where
- *  \f${\alpha}_{\text{sw}}\f$ is a coefficient,  and \f$S_0\f$ is the
- *  initial saturation, and \f$S_{\text{max}}\f$ is the maximum
- * saturation.
+ *  \f${\alpha}_{\text{sw}}\f$ is a coefficient, and \f$S_0\f$ is the
+ *  initial saturation, and \f$S_{\text{max}}\f$ is the maximum saturation.
  *  The coefficient gives the swelling stress at full saturation, which can be
  *  computed as
  *    \f[ {\alpha}_{\text{sw}} =