diff --git a/include/aspect/material_model/reactive_fluid_transport.h b/include/aspect/material_model/reactive_fluid_transport.h
index 599b3f1bf21..4b72839e4f5 100644
--- a/include/aspect/material_model/reactive_fluid_transport.h
+++ b/include/aspect/material_model/reactive_fluid_transport.h
@@ -179,6 +179,11 @@ namespace aspect
         std::vector<double> csat_sediment_poly_coeffs {-0.150662, 0.301807, 1.01867};
         std::vector<double> Td_sediment_poly_coeffs {2.83277, -24.7593, 85.9090, 524.898};
 
+        // The polynomials breakdown above certain pressures, 10 GPa for peridotite, 26 GPa for gabbro, 16 GPa for MORB,
+        // and 50 GPa for sediment. These cutoff pressures were determined by extending the pressure range in Tian et al. (2019)
+        // and observing where the maximum allowed water contents jump towards infinite values.
+        const std::vector<double> pressure_cutoffs {10, 26, 16, 50};
+
         std::vector<std::vector<double>> devolatilization_enthalpy_changes {LR_peridotite_poly_coeffs, LR_gabbro_poly_coeffs, \
                                                                              LR_MORB_poly_coeffs, LR_sediment_poly_coeffs
                                                                             };
diff --git a/source/material_model/reactive_fluid_transport.cc b/source/material_model/reactive_fluid_transport.cc
index da5a7622253..b384cf909fa 100644
--- a/source/material_model/reactive_fluid_transport.cc
+++ b/source/material_model/reactive_fluid_transport.cc
@@ -56,10 +56,6 @@ namespace aspect
     tian_equilibrium_bound_water_content (const MaterialModel::MaterialModelInputs<dim> &in,
                                           unsigned int q) const
     {
-      // Pressure, which must be in GPa for the parametrization, or GPa^-1
-      const double pressure = in.pressure[q]<=0 ? 1e-12 : in.pressure[q]/1.e9;
-      const double inverse_pressure = std::pow(pressure, -1);
-
       // Create arrays that will store the values of the polynomials at the current pressure
       std::vector<double> LR_values(4);
       std::vector<double> csat_values(4);
@@ -71,6 +67,12 @@ namespace aspect
       // Td polynomials are defined in equations 15, B3, B11, and B19.
       for (unsigned int i = 0; i<devolatilization_enthalpy_changes.size(); ++i)
         {
+          // Pressure, which must be in GPa for the parametrization, or GPa^-1. The polynomials for each lithology
+          // breaks down above certain pressures, make sure that we cap the pressure just before this break down.
+          // Introduce minimum pressure to avoid a division by 0.
+          const double minimum_pressure = 1e-12;
+          const double pressure = std::min(std::max(minimum_pressure, in.pressure[q]/1.e9), pressure_cutoffs[i]);
+          const double inverse_pressure = 1.0/pressure;
           for (unsigned int j = 0; j<devolatilization_enthalpy_changes[i].size(); ++j)
             {
               LR_values[i] += devolatilization_enthalpy_changes[i][j] * std::pow(inverse_pressure, devolatilization_enthalpy_changes[i].size() - 1 - j);
diff --git a/tests/tian_MORB/screen-output b/tests/tian_MORB/screen-output
index 62d7e789b7c..f156d5429ce 100644
--- a/tests/tian_MORB/screen-output
+++ b/tests/tian_MORB/screen-output
@@ -31,8 +31,8 @@ Number of degrees of freedom: 989 (162+73+162+25+81+81+81+81+81+81+81)
    Rebuilding Stokes preconditioner...
    Solving Stokes system... 0+0 iterations.
    Solving fluid velocity system... 1 iterations.
-      Relative nonlinear residuals (temperature, compositional fields, Stokes system): 1.54103e-16, 0, 1.39664e-16, 0, 0, 8.51473e-17, 0, 4.64648e-10
-      Relative nonlinear residual (total system) after nonlinear iteration 2: 4.64648e-10
+      Relative nonlinear residuals (temperature, compositional fields, Stokes system): 1.54103e-16, 0, 1.39664e-16, 0, 0, 8.51473e-17, 0, 4.64647e-10
+      Relative nonlinear residual (total system) after nonlinear iteration 2: 4.64647e-10
 
 
    Postprocessing:
@@ -76,7 +76,7 @@ Number of degrees of freedom: 989 (162+73+162+25+81+81+81+81+81+81+81)
    Rebuilding Stokes preconditioner...
    Solving Stokes system... 0+0 iterations.
    Solving fluid velocity system... 1 iterations.
-      Relative nonlinear residuals (temperature, compositional fields, Stokes system): 2.5218e-12, 1.52682e-10, 5.02896e-12, 0, 0, 1.07956e-16, 0, 1.87227e-09
+      Relative nonlinear residuals (temperature, compositional fields, Stokes system): 2.5218e-12, 1.52682e-10, 5.02898e-12, 0, 0, 1.07956e-16, 0, 1.87227e-09
       Relative nonlinear residual (total system) after nonlinear iteration 3: 1.87227e-09
 
 
diff --git a/tests/tian_gabbro/screen-output b/tests/tian_gabbro/screen-output
index 9fc3c376235..d9e23709225 100644
--- a/tests/tian_gabbro/screen-output
+++ b/tests/tian_gabbro/screen-output
@@ -31,8 +31,8 @@ Number of degrees of freedom: 989 (162+73+162+25+81+81+81+81+81+81+81)
    Rebuilding Stokes preconditioner...
    Solving Stokes system... 0+0 iterations.
    Solving fluid velocity system... 1 iterations.
-      Relative nonlinear residuals (temperature, compositional fields, Stokes system): 1.54103e-16, 0, 1.39664e-16, 0, 8.51473e-17, 0, 0, 4.64648e-10
-      Relative nonlinear residual (total system) after nonlinear iteration 2: 4.64648e-10
+      Relative nonlinear residuals (temperature, compositional fields, Stokes system): 1.54103e-16, 0, 1.39664e-16, 0, 8.51473e-17, 0, 0, 4.64647e-10
+      Relative nonlinear residual (total system) after nonlinear iteration 2: 4.64647e-10
 
 
    Postprocessing:
@@ -76,8 +76,8 @@ Number of degrees of freedom: 989 (162+73+162+25+81+81+81+81+81+81+81)
    Rebuilding Stokes preconditioner...
    Solving Stokes system... 0+0 iterations.
    Solving fluid velocity system... 1 iterations.
-      Relative nonlinear residuals (temperature, compositional fields, Stokes system): 2.771e-12, 1.64474e-10, 4.67296e-12, 0, 9.75727e-17, 0, 0, 2.23762e-09
-      Relative nonlinear residual (total system) after nonlinear iteration 3: 2.23762e-09
+      Relative nonlinear residuals (temperature, compositional fields, Stokes system): 2.77098e-12, 1.64474e-10, 4.67294e-12, 0, 9.75727e-17, 0, 0, 2.23761e-09
+      Relative nonlinear residual (total system) after nonlinear iteration 3: 2.23761e-09
 
 
    Postprocessing:
diff --git a/tests/tian_peridotite/screen-output b/tests/tian_peridotite/screen-output
index 9b539c1aa20..b89016b8b78 100644
--- a/tests/tian_peridotite/screen-output
+++ b/tests/tian_peridotite/screen-output
@@ -31,8 +31,8 @@ Number of degrees of freedom: 989 (162+73+162+25+81+81+81+81+81+81+81)
    Rebuilding Stokes preconditioner...
    Solving Stokes system... 0+0 iterations.
    Solving fluid velocity system... 1 iterations.
-      Relative nonlinear residuals (temperature, compositional fields, Stokes system): 1.54103e-16, 0, 1.39664e-16, 8.51473e-17, 0, 0, 0, 4.64648e-10
-      Relative nonlinear residual (total system) after nonlinear iteration 2: 4.64648e-10
+      Relative nonlinear residuals (temperature, compositional fields, Stokes system): 1.54103e-16, 0, 1.39664e-16, 8.51473e-17, 0, 0, 0, 4.64647e-10
+      Relative nonlinear residual (total system) after nonlinear iteration 2: 4.64647e-10
 
 
    Postprocessing:
@@ -63,7 +63,7 @@ Number of degrees of freedom: 989 (162+73+162+25+81+81+81+81+81+81+81)
    Rebuilding Stokes preconditioner...
    Solving Stokes system... 7+0 iterations.
    Solving fluid velocity system... 1 iterations.
-      Relative nonlinear residuals (temperature, compositional fields, Stokes system): 2.181e-05, 0.000574823, 0.000145666, 1.02675e-16, 0, 0, 0, 0.000630417
+      Relative nonlinear residuals (temperature, compositional fields, Stokes system): 2.181e-05, 0.000574823, 0.000145666, 9.16123e-17, 0, 0, 0, 0.000630417
       Relative nonlinear residual (total system) after nonlinear iteration 2: 0.000630417
 
    Solving temperature system... 4 iterations.
@@ -76,7 +76,7 @@ Number of degrees of freedom: 989 (162+73+162+25+81+81+81+81+81+81+81)
    Rebuilding Stokes preconditioner...
    Solving Stokes system... 4+0 iterations.
    Solving fluid velocity system... 1 iterations.
-      Relative nonlinear residuals (temperature, compositional fields, Stokes system): 4.66691e-08, 2.35667e-06, 3.02875e-07, 1.25585e-16, 0, 0, 0, 2.19117e-06
+      Relative nonlinear residuals (temperature, compositional fields, Stokes system): 4.66691e-08, 2.35667e-06, 3.02875e-07, 1.22071e-16, 0, 0, 0, 2.19117e-06
       Relative nonlinear residual (total system) after nonlinear iteration 3: 2.35667e-06
 
 
diff --git a/tests/tian_sediment/screen-output b/tests/tian_sediment/screen-output
index b1131d74d48..96aa664d87e 100644
--- a/tests/tian_sediment/screen-output
+++ b/tests/tian_sediment/screen-output
@@ -31,8 +31,8 @@ Number of degrees of freedom: 989 (162+73+162+25+81+81+81+81+81+81+81)
    Rebuilding Stokes preconditioner...
    Solving Stokes system... 0+0 iterations.
    Solving fluid velocity system... 1 iterations.
-      Relative nonlinear residuals (temperature, compositional fields, Stokes system): 1.54103e-16, 0, 1.39664e-16, 0, 0, 0, 8.51473e-17, 4.64648e-10
-      Relative nonlinear residual (total system) after nonlinear iteration 2: 4.64648e-10
+      Relative nonlinear residuals (temperature, compositional fields, Stokes system): 1.54103e-16, 0, 1.39664e-16, 0, 0, 0, 8.51473e-17, 4.64647e-10
+      Relative nonlinear residual (total system) after nonlinear iteration 2: 4.64647e-10
 
 
    Postprocessing:
@@ -76,7 +76,7 @@ Number of degrees of freedom: 989 (162+73+162+25+81+81+81+81+81+81+81)
    Rebuilding Stokes preconditioner...
    Solving Stokes system... 0+0 iterations.
    Solving fluid velocity system... 1 iterations.
-      Relative nonlinear residuals (temperature, compositional fields, Stokes system): 1.77101e-12, 1.37697e-10, 2.93736e-12, 0, 0, 0, 1.04501e-16, 1.53067e-09
+      Relative nonlinear residuals (temperature, compositional fields, Stokes system): 1.77098e-12, 1.37697e-10, 2.93735e-12, 0, 0, 0, 1.04501e-16, 1.53067e-09
       Relative nonlinear residual (total system) after nonlinear iteration 3: 1.53067e-09