[web/doc] upload Matlab script for semianalytical solution

web/content/docs/benchmarks/thermal-two-phase-flow/TCE-diffusion/AH2010.m

+% 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)');

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}$.