diff --git a/web/content/docs/benchmarks/thermal-two-phase-flow/TCE-diffusion/AH2010.m b/web/content/docs/benchmarks/thermal-two-phase-flow/TCE-diffusion/AH2010.m
new file mode 100644
index 0000000000000000000000000000000000000000..0b20a31b19ef7a61b85d56097aa4d4132bd77067
--- /dev/null
+++ b/web/content/docs/benchmarks/thermal-two-phase-flow/TCE-diffusion/AH2010.m
@@ -0,0 +1,40 @@
+% This MATLAB script is used to calculate the steady-state solution of
+% diffusion across the capillary fringe. It is based on the method
+% suggested by Atteia and Ho``hener (2010).
+
+% Atteia, O., & Ho``hener, P. (2010). Semianalytical model predicting
+% transfer of volatile pollutants from groundwater to the soil surface.
+% Environmental science & technology, 44(16), 6228-6232.
+
+% Author: Boyan Meng
+% Email: boyan107@gmail.com
+
+clear;
+% Saturation profile at steady-state.
+z = [1	0.980568	0.961883	0.943917	0.926642	0.910031	0.89406	0.878702	0.863935	0.849736	0.836084	0.822956	0.810333	0.798196	0.786526	0.775304	0.764514	0.754139	0.744163	0.734571	0.725348	0.716479	0.707952	0.699752	0.691868	0.684287	0.676998	0.669989	0.66325	0.656769	0.650538	0.644547	0.638786	0.633247	0.627921	0.622799	0.617875	0.61314	0.608587	0.604209	0.6	0.596172	0.59253	0.589056	0.585734	0.58255	0.579491	0.576546	0.573706	0.57096	0.568301	0.565722	0.563214	0.560773	0.558392	0.556066	0.55379	0.551559	0.549369	0.547216	0.545096	0.543006	0.540942	0.5389	0.536879	0.534874	0.532883	0.530903	0.528931	0.526964	0.525	0.523036	0.521069	0.519097	0.517117	0.515126	0.513121	0.5111	0.509058	0.506994	0.504904	0.502784	0.500631	0.498441	0.49621	0.493934	0.491608	0.489227	0.486786	0.484278	0.481699	0.47904	0.476294	0.473454	0.470509	0.46745	0.464266	0.460944	0.45747	0.453828	0.45	0.445888	0.441627	0.437211	0.432635	0.427893	0.422979	0.417886	0.41261	0.407141	0.401475	0.395603	0.389517	0.383212	0.376677	0.369906	0.362888	0.355617	0.348081	0.340273	0.332181	0.323795	0.315106	0.306101	0.29677	0.2871	0.27708	0.266696	0.255935	0.244785	0.233229	0.221255	0.208846	0.195988	0.182663	0.168854	0.154545	0.139717	0.124351	0.108428	0.0919271	0.0748279	0.0571084	0.0387463	0.0197182	0];
+S_w = [0.00634606	0.00725321	0.00828275	0.00945034	0.010773	0.0122694	0.0139604	0.0158688	0.0180196	0.0204401	0.02316	0.0262117	0.02963	0.0334521	0.0377181	0.0424703	0.0477534	0.0536138	0.0601001	0.0672615	0.075148	0.0838096	0.093295	0.103651	0.11492	0.127141	0.140347	0.154563	0.169803	0.186072	0.203364	0.221657	0.240916	0.261092	0.282118	0.303916	0.326391	0.349439	0.372943	0.396778	0.420821	0.443602	0.466042	0.488095	0.509719	0.530879	0.551545	0.571695	0.591309	0.610375	0.628884	0.646829	0.664209	0.681023	0.697275	0.71297	0.728113	0.742713	0.756779	0.770321	0.78335	0.795875	0.80791	0.819465	0.830553	0.841184	0.851372	0.861126	0.870459	0.879381	0.887903	0.896036	0.903789	0.911172	0.918195	0.924865	0.931193	0.937186	0.942851	0.948197	0.953231	0.95796	0.962391	0.96653	0.970384	0.97396	0.977263	0.980301	0.983081	0.985608	0.987891	0.989937	0.991754	0.993351	0.994737	0.995925	0.996924	0.997747	0.998409	0.998925	0.99931	0.999585	0.999755	0.99984	0.999867	0.999873	0.999878	0.999884	0.999889	0.999894	0.9999	0.999905	0.99991	0.999916	0.999921	0.999926	0.999932	0.999937	0.999942	0.999947	0.999952	0.999957	0.999961	0.999965	0.99997	0.999974	0.999977	0.999981	0.999984	0.999987	0.999989	0.999991	0.999993	0.999995	0.999996	0.999997	0.999998	0.999999	0.999999	1	1	1	1	1	1	1];
+depth = 1 - z;
+phi = 0.38; % porosity
+D_a = 8.3e-6; % TCE diffusion coefficient in air 
+D_w = 9.1e-10; % TCE diffusion coefficient in water 
+theta_a = phi * (1-S_w); % volumetric air content
+theta_w = phi * S_w; % volumetric water content
+alpha_z = 1e-3; % vertical transverse dispersivity
+v = 0; % seepage velocity
+H = 0.38; % Henry constant of TCE (non-dimensional)
+N_G = 101325/8.31446/298.15; % molar density of ideal gas at 25 dC
+m = 0.8; % van Genuchten paramter
+k_rel_wet = sqrt(S_w).*(1-(1-S_w.^(1/m)).^m).^2; % relative permeability of wetting phase
+A = D_a*theta_a.^(10/3)/phi/phi + (D_w*theta_w.^(10/3)/phi/phi + alpha_z*v*k_rel_wet)/H;
+F = 1./A;
+I = zeros (1,110);
+for i=2:110
+    ddepth = depth(i)-depth(i-1);
+    I(i) = I(i-1) + (F(i)+F(i-1))/2 * ddepth;
+end
+c_rel = I./I(end); % relative TCE concentration in the gas phase 
+d2WT = depth(110) - depth; % distance to water table (+ means above)
+semilogx(c_rel, d2WT(1:110));
+ylim([0 0.6]);
+xlabel('Rel. TCE (g) conc. (-)');
+ylabel('Distance above water table (m)');
diff --git a/web/content/docs/benchmarks/thermal-two-phase-flow/TCE-diffusion/index.md b/web/content/docs/benchmarks/thermal-two-phase-flow/TCE-diffusion/index.md
index d88d648569c19cdc5f4bb02b2c7a3408f9f49859..f74a3cb3cddc7f23dfe3284debb03971d5e1912a 100644
--- a/web/content/docs/benchmarks/thermal-two-phase-flow/TCE-diffusion/index.md
+++ b/web/content/docs/benchmarks/thermal-two-phase-flow/TCE-diffusion/index.md
@@ -42,7 +42,7 @@ By performing numerical integration, we have
 \begin{equation}
     x^c_a=c+J\cdot I(z)\ \mathrm{with}\ I(z)=\int\frac{1}{A(z)}dz.
 \end{equation}
-where $I(z)$ is integrated w.r.t. $z$ in the unsaturated zone. Since at steady state, the diffusive flux should be uniform at any depth, the integral constants $c$ and $J$ can be solved by the upper and lower values of $x^c_a$ which are given as Dirichlet boundary conditions.
+where $I(z)$ is integrated w.r.t. $z$ in the unsaturated zone. Since at steady state, the diffusive flux should be uniform at any depth, the integral constants $c$ and $J$ can be solved by the upper and lower values of $x^c_a$ which are given as Dirichlet boundary conditions. A [MATLAB script](AH2010.m) is provided for the above solution.
 
 Figure 2 compares the normalized gas phase TCE concentration profiles given by the numerical model and the above semianalytical solution at steady state. Figure 3 shows the absolute and relative errors between the two curves in Fig. 2. The maximum relative error is less than 1\% except one point. Thus, a very good agreement was obtained. It is interesting to note the sharp decline of the TCE concentration immediately above the water table, which is due to the several orders of magnitude difference between $D_{0a}$ and $D_{0w}$.