diff --git a/CHANGELOG.md b/CHANGELOG.md index 674c8815..1f25150b 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,5 +1,14 @@ # Unreleased +**Changed:** + +- All `cond*_` functions are refcatored and moved to `+conductivity` folder in + [#189](https://github.com/EcoExtreML/STEMMUS_SCOPE/pull/189) + +**Fixed:** + +- issue [#181](https://github.com/EcoExtreML/STEMMUS_SCOPE/issues/181) + # [1.3.0](https://github.com/EcoExtreML/STEMMUS_SCOPE/releases/tag/1.3.0) - 22 Jun 2023 diff --git a/run_model_on_snellius/exe/STEMMUS_SCOPE b/run_model_on_snellius/exe/STEMMUS_SCOPE index f7a8fbff..af863234 100755 Binary files a/run_model_on_snellius/exe/STEMMUS_SCOPE and b/run_model_on_snellius/exe/STEMMUS_SCOPE differ diff --git a/src/+conductivity/+hydraulicConductivity/calculateDTheta.m b/src/+conductivity/+hydraulicConductivity/calculateDTheta.m new file mode 100644 index 00000000..d214616b --- /dev/null +++ b/src/+conductivity/+hydraulicConductivity/calculateDTheta.m @@ -0,0 +1,43 @@ +function dtheta = calculateDTheta(subRoutine, heat_term, theta_s, theta_r, theta_m, gamma_hh, theta_ll, theta_l, theta_uu, theta_u, hh, h, hh_frez, h_frez, phi_s, lamda, alpha, n, m) + + % get soil constants + SoilConstants = io.getSoilConstants(); + + switch subRoutine + case 0 + % for DTheta_LLh, heat_term = hh + % for DTheta_UUh, heat_term = hh + hh_frez + A = (-theta_r) / (abs(heat_term) * log(abs(SoilConstants.hd / SoilConstants.hm))) * (1 - (1 + abs(alpha * heat_term)^n)^(-m)); + B = alpha * n * m * (theta_m - gamma_hh * theta_r); + C = ((1 + abs(alpha * heat_term)^n)^(-m - 1)); + D = (abs(alpha * heat_term)^(n - 1)); + dtheta = A - B * C * D; + if (hh + hh_frez) <= SoilConstants.hd + dtheta = 0; + end + case 1 + % heat_term = hh or hh + hh_frez + A = (theta_m - theta_r) * alpha * n; + B = abs(alpha * heat_term)^(n - 1) * (-m); + C = (1 + abs(alpha * heat_term)^n)^(-m - 1); + dtheta = A * B * C; + case 2 + A = theta_ll - theta_l; + B = hh + hh_frez - h - h_frez; + dtheta = A / B; + case 3 + A = theta_uu - theta_u; + B = hh - h; + dtheta = A / B; + case 4 + dtheta = theta_s / phi_s * (hh / phi_s)^(-1 * lamda - 1); + case 5 + % this case differs with case 1 only regarding `theta_m` and `theta_s` + % see issue 181, item 7 + % heat_term = hh or hh + hh_frez + A = (theta_s - theta_r) * alpha * n; + B = abs(alpha * heat_term)^(n - 1) * (-m); + C = (1 + abs(alpha * heat_term)^n)^(-m - 1); + dtheta = A * B * C; + end +end diff --git a/src/+conductivity/+hydraulicConductivity/calculateDTheta_UUh.m b/src/+conductivity/+hydraulicConductivity/calculateDTheta_UUh.m new file mode 100644 index 00000000..ef5514dd --- /dev/null +++ b/src/+conductivity/+hydraulicConductivity/calculateDTheta_UUh.m @@ -0,0 +1,65 @@ +function dtheta_uuh = calculateDTheta_UUh(theta_uu, theta_m, theta_ll, gamma_hh, SoilVariables, VanGenuchten) + theta_s = VanGenuchten.Theta_s; + theta_r = VanGenuchten.Theta_r; + alpha = VanGenuchten.Alpha; + n = VanGenuchten.n; + m = VanGenuchten.m; + + theta_u = SoilVariables.Theta_U; + theta_l = SoilVariables.Theta_L; + hh = SoilVariables.hh; + h = SoilVariables.h; + h_frez = SoilVariables.h_frez; + hh_frez = SoilVariables.hh_frez; + phi_s = SoilVariables.Phi_s; + lamda = SoilVariables.Lamda; + + % get model settings + ModelSettings = io.getModelSettings(); + + heatTerm = hh + hh_frez; + if ModelSettings.SWCC == 1 + if ModelSettings.SFCC == 1 + if hh >= -1 + if hh_frez >= 0 + dtheta_uuh = 0; + else + if ModelSettings.Thmrlefc + subRoutine = 0; + dtheta_uuh = conductivity.hydraulicConductivity.calculateDTheta(subRoutine, heatTerm, theta_s, theta_r, theta_m, gamma_hh, theta_ll, theta_l, theta_uu, theta_u, hh, h, hh_frez, h_frez, phi_s, lamda, alpha, n, m); + elseif abs(hh - h) < 1e-3 + subRoutine = 1; + dtheta_uuh = conductivity.hydraulicConductivity.calculateDTheta(subRoutine, heatTerm, theta_s, theta_r, theta_m, gamma_hh, theta_ll, theta_l, theta_uu, theta_u, hh, h, hh_frez, h_frez, phi_s, lamda, alpha, n, m); + else + subRoutine = 2; + dtheta_uuh = conductivity.hydraulicConductivity.calculateDTheta(subRoutine, heatTerm, theta_s, theta_r, theta_m, gamma_hh, theta_ll, theta_l, theta_uu, theta_u, hh, h, hh_frez, h_frez, phi_s, lamda, alpha, n, m); + end + end + elseif ModelSettings.Thmrlefc + subRoutine = 0; + dtheta_uuh = conductivity.hydraulicConductivity.calculateDTheta(subRoutine, heatTerm, theta_s, theta_r, theta_m, gamma_hh, theta_ll, theta_l, theta_uu, theta_u, hh, h, hh_frez, h_frez, phi_s, lamda, alpha, n, m); + end + elseif hh >= -1e-6 || hh <= -1e7 + dtheta_uuh = 0; + elseif ModelSettings.Thmrlefc || abs(hh - h) < 1e-3 + subRoutine = 1; + dtheta_uuh = conductivity.hydraulicConductivity.calculateDTheta(subRoutine, heatTerm, theta_s, theta_r, theta_m, gamma_hh, theta_ll, theta_l, theta_uu, theta_u, hh, h, hh_frez, h_frez, phi_s, lamda, alpha, n, m); + else + subRoutine = 3; + dtheta_uuh = conductivity.hydraulicConductivity.calculateDTheta(subRoutine, heatTerm, theta_s, theta_r, theta_m, gamma_hh, theta_ll, theta_l, theta_uu, theta_u, hh, h, hh_frez, h_frez, phi_s, lamda, alpha, n, m); + end + else + if hh >= phi_s || hh <= -1e7 + dtheta_uuh = 0; + elseif ModelSettings.Thmrlefc + subRoutine = 4; + dtheta_uuh = conductivity.hydraulicConductivity.calculateDTheta(subRoutine, heatTerm, theta_s, theta_r, theta_m, gamma_hh, theta_ll, theta_l, theta_uu, theta_u, hh, h, hh_frez, h_frez, phi_s, lamda, alpha, n, m); + elseif abs(hh - h) < 1e-3 + subRoutine = 5; + dtheta_uuh = conductivity.hydraulicConductivity.calculateDTheta(subRoutine, heatTerm, theta_s, theta_r, theta_m, gamma_hh, theta_ll, theta_l, theta_uu, theta_u, hh, h, hh_frez, h_frez, phi_s, lamda, alpha, n, m); + else + subRoutine = 3; + dtheta_uuh = conductivity.hydraulicConductivity.calculateDTheta(subRoutine, heatTerm, theta_s, theta_r, theta_m, gamma_hh, theta_ll, theta_l, theta_uu, theta_u, hh, h, hh_frez, h_frez, phi_s, lamda, alpha, n, m); + end + end +end diff --git a/src/+conductivity/+hydraulicConductivity/calculateGamma_hh.m b/src/+conductivity/+hydraulicConductivity/calculateGamma_hh.m new file mode 100644 index 00000000..92c18643 --- /dev/null +++ b/src/+conductivity/+hydraulicConductivity/calculateGamma_hh.m @@ -0,0 +1,12 @@ +function gamma_hh = calculateGamma_hh(hh) + % get soil constants + SoilConstants = io.getSoilConstants(); + + if abs(hh) >= abs(SoilConstants.hd) + gamma_hh = 0; + elseif abs(hh) >= abs(SoilConstants.hm) + gamma_hh = log(abs(SoilConstants.hd) / abs(hh)) / log(abs(SoilConstants.hd) / abs(SoilConstants.hm)); + else + gamma_hh = 1; + end +end diff --git a/src/+conductivity/+hydraulicConductivity/calculateKL_h.m b/src/+conductivity/+hydraulicConductivity/calculateKL_h.m new file mode 100644 index 00000000..64b420a6 --- /dev/null +++ b/src/+conductivity/+hydraulicConductivity/calculateKL_h.m @@ -0,0 +1,23 @@ +function kl_h = calculateKL_h(mu_w, se, Ks, m) + + % load Constants + Constants = io.define_constants(); + + MU_WN = Constants.MU_W0 * exp(Constants.MU1 / (8.31441 * (20 + 133.3))); + CKT = MU_WN / mu_w; + if se == 0 + kl_h = 0; + else + kl_h = CKT * Ks * (se^(0.5)) * (1 - (1 - se^(1 / m))^m)^2; + end + if kl_h <= 1E-20 + kl_h = 1E-20; + end + if isnan(kl_h) == 1 + kl_h = 0; + warning('\n case "isnan(kl_h) == 1", set "kl_h = 0" \r'); + end + if ~isreal(kl_h) + warning('\n case "~isreal(kl_h)", dont know what to do! \r'); + end +end diff --git a/src/+conductivity/+hydraulicConductivity/calculateSe.m b/src/+conductivity/+hydraulicConductivity/calculateSe.m new file mode 100644 index 00000000..08a98e9a --- /dev/null +++ b/src/+conductivity/+hydraulicConductivity/calculateSe.m @@ -0,0 +1,57 @@ +function se = calculateSe(theta_ll, gamma_hh, SoilVariables) + % get model settings + ModelSettings = io.getModelSettings(); + + % get soil constants + SoilConstants = io.getSoilConstants(); + + hh_frez = SoilVariables.hh_frez; + POR = SoilVariables.POR; + hh = SoilVariables.hh; + phi_s = SoilVariables.Phi_s; + se = SoilVariables.Se; + seCalculated = theta_ll / POR; + + if ModelSettings.SWCC == 1 + if ModelSettings.SFCC == 1 + if hh >= -1 + if hh_frez >= 0 + se = 1; + elseif Thmrlefc + if (hh + hh_frez) <= SoilConstants.hd + se = 0; + else + se = seCalculated; + end + end + elseif ModelSettings.Thmrlefc + if (hh + hh_frez) <= SoilConstants.hd + se = 0; + else + se = seCalculated; + end + end + elseif theta_ll <= 0.06 || hh <= -1e7 + se = 0; + else + se = seCalculated; + end + else + if hh + hh_frez <= -1e7 || hh <= -1e7 + se = 0; + elseif hh + hh_frez >= phi_s + se = 1; + else + se = seCalculated; + end + end + + if se > 1 + se = 1; + elseif se < 0 + se = 0; + end + if isnan(se) == 1 + warning('\n case "isnan(se) == 1" happens. Dont know what to do! \r'); + end +end diff --git a/src/+conductivity/+hydraulicConductivity/calculateTheta.m b/src/+conductivity/+hydraulicConductivity/calculateTheta.m new file mode 100644 index 00000000..01fc9f13 --- /dev/null +++ b/src/+conductivity/+hydraulicConductivity/calculateTheta.m @@ -0,0 +1,25 @@ +function theta = calculateTheta(subRoutine, theta_m, heat_term, gamma_hh, theta_s, theta_r, lamda, phi_s, alpha, n, m) + % get soil constants + SoilConstants = io.getSoilConstants(); + hd = SoilConstants.hd; + + switch subRoutine + case 0 + % heat_term = hh or hh + hh_frez + theta = gamma_hh * theta_r + (theta_m - gamma_hh * theta_r) / (1 + abs(alpha * heat_term)^n)^m; + if heat_term <= hd + theta = 0; + end + case 1 + % heat_term = hh or hh + hh_frez + theta = theta_s * (heat_term / phi_s)^(-1 * lamda); + if heat_term <= -1e7 + theta = Theta_r; + elseif heat_term >= phi_s + theta = theta_s; + end + case 2 + % heat_term = hh or hh + hh_frez + theta = Theta_r + (theta_s - Theta_r) / (1 + abs(alpha * heat_term)^n)^m; + end +end diff --git a/src/+conductivity/+hydraulicConductivity/calculateTheta_II.m b/src/+conductivity/+hydraulicConductivity/calculateTheta_II.m new file mode 100644 index 00000000..125f191c --- /dev/null +++ b/src/+conductivity/+hydraulicConductivity/calculateTheta_II.m @@ -0,0 +1,23 @@ +function theta_ii = calculateTheta_II(tt, xcap, hh, Theta_II) + + % get model settings + ModelSettings = io.getModelSettings(); + + Tf1 = 273.15 + 1; + Tf2 = 273.15 - 3; + + theta_ii = Theta_II; + + if hh <= -1e7 + theta_ii = 0; + end + if ModelSettings.SWCC == 1 && ModelSettings.SFCC ~= 1 + if tt + 273.15 > Tf1 + theta_ii = 0; + elseif tt + 273.15 >= Tf2 + theta_ii = 0.5 * (1 - sin(pi() * (tt + 273.15 - 0.5 * Tf1 - 0.5 * Tf2) / (Tf1 - Tf2))) * xcap; + else + theta_ii = xcap; + end + end +end diff --git a/src/+conductivity/+hydraulicConductivity/calculateTheta_LL.m b/src/+conductivity/+hydraulicConductivity/calculateTheta_LL.m new file mode 100644 index 00000000..a580dfbd --- /dev/null +++ b/src/+conductivity/+hydraulicConductivity/calculateTheta_LL.m @@ -0,0 +1,54 @@ +function theta_ll = calculateTheta_LL(theta_uu, theta_ii, theta_m, gamma_hh, SoilVariables, VanGenuchten) + % get model settings + ModelSettings = io.getModelSettings(); + + % load Constants + Constants = io.define_constants(); + + hh = SoilVariables.hh; + hh_frez = SoilVariables.hh_frez; + phi_s = SoilVariables.Phi_s; + lamda = SoilVariables.Lamda; + + theta_s = VanGenuchten.Theta_s; + theta_r = VanGenuchten.Theta_r; + alpha = VanGenuchten.Alpha; + n = VanGenuchten.n; + m = VanGenuchten.m; + + heatTerm = hh + hh_frez; + + % calculate theta_ll + theta_ll = SoilVariables.Theta_LL; + if ModelSettings.SWCC == 1 + if ModelSettings.SFCC == 1 + if hh >= -1 && hh_frez >= 0 + theta_ll = theta_s; + elseif ModelSettings.Thmrlefc + subRoutine = 0; + theta_ll = conductivity.hydraulicConductivity.calculateTheta(subRoutine, theta_m, heatTerm, gamma_hh, theta_s, theta_r, lamda, phi_s, alpha, n, m); + end + else + if hh >= -1e-6 && TT + 273.15 > (273.15 + 1) + theta_ll = theta_s; + elseif hh <= -1e7 + theta_ll = theta_r; + elseif TT + 273.15 > (273.15 + 1) + subRoutine = 2; + theta_ll = conductivity.hydraulicConductivity.calculateTheta(subRoutine, theta_m, heatTerm, gamma_hh, theta_s, theta_r, lamda, phi_s, alpha, n, m); + else + theta_ll = theta_uu - theta_ii * Constants.RHOI / Constants.RHOL; + end + if theta_ll <= 0.06 + theta_ll = 0.06; + end + end + else + if hh >= phi_s + subRoutine = 1; + theta_ll = conductivity.hydraulicConductivity.calculateTheta(subRoutine, theta_m, heatTerm, gamma_hh, theta_s, theta_r, lamda, phi_s, alpha, n, m); + elseif hh <= -1e7 + theta_ll = theta_r; + end + end +end diff --git a/src/+conductivity/+hydraulicConductivity/calculateTheta_UU.m b/src/+conductivity/+hydraulicConductivity/calculateTheta_UU.m new file mode 100644 index 00000000..8d4348a3 --- /dev/null +++ b/src/+conductivity/+hydraulicConductivity/calculateTheta_UU.m @@ -0,0 +1,43 @@ +function theta_uu = calculateTheta_UU(theta_m, gamma_hh, SoilVariables, VanGenuchten) + + hh = SoilVariables.hh; + phi_s = SoilVariables.Phi_s; + lamda = SoilVariables.Lamda; + + theta_s = VanGenuchten.Theta_s; + theta_r = VanGenuchten.Theta_r; + alpha = VanGenuchten.Alpha; + n = VanGenuchten.n; + m = VanGenuchten.m; + + % get model settings + ModelSettings = io.getModelSettings(); + + % calculate theta_uu + if ModelSettings.SWCC == 1 + if ModelSettings.SFCC == 1 + if hh >= -1 + theta_uu = theta_s; + elseif ModelSettings.Thmrlefc + if gamma_hh == 0 + theta_uu = 0; + else + subRoutine = 0; + theta_uu = conductivity.hydraulicConductivity.calculateTheta(subRoutine, theta_m, hh, gamma_hh, theta_s, theta_r, lamda, phi_s, alpha, n, m); + end + end + else + if hh >= -1e-6 + theta_uu = theta_s; + elseif hh <= -1e7 + theta_uu = theta_r; + else + subRoutine = 2; + theta_uu = conductivity.hydraulicConductivity.calculateTheta(subRoutine, theta_m, hh, gamma_hh, theta_s, theta_r, lamda, phi_s, alpha, n, m); + end + end + else + subRoutine = 1; + theta_uu = conductivity.hydraulicConductivity.calculateTheta(subRoutine, theta_m, hh, gamma_hh, theta_s, theta_r, lamda, phi_s, alpha, n, m); + end +end diff --git a/src/+conductivity/+hydraulicConductivity/calculateTheta_m.m b/src/+conductivity/+hydraulicConductivity/calculateTheta_m.m new file mode 100644 index 00000000..98e45a11 --- /dev/null +++ b/src/+conductivity/+hydraulicConductivity/calculateTheta_m.m @@ -0,0 +1,15 @@ +function theta_m = calculateTheta_m(gamma_hh, VanGenuchten, POR) + + theta_s = VanGenuchten.Theta_s; + theta_r = VanGenuchten.Theta_r; + alpha = VanGenuchten.Alpha; + n = VanGenuchten.n; + m = VanGenuchten.m; + + theta_m = gamma_hh * theta_r + (theta_s - gamma_hh * theta_r) * (1 + abs(alpha * (-1))^n)^m; + if theta_m >= POR + theta_m = POR; + elseif theta_m <= theta_s + theta_m = theta_s; + end +end diff --git a/src/+conductivity/+hydraulicConductivity/calcuulateDTheta_LLh.m b/src/+conductivity/+hydraulicConductivity/calcuulateDTheta_LLh.m new file mode 100644 index 00000000..eaf543c6 --- /dev/null +++ b/src/+conductivity/+hydraulicConductivity/calcuulateDTheta_LLh.m @@ -0,0 +1,57 @@ +function dtheta_llh = calcuulateDTheta_LLh(dtheta_uuh, theta_m, theta_uu, theta_ll, gamma_hh, SoilVariables, VanGenuchten) + theta_s = VanGenuchten.Theta_s; + theta_r = VanGenuchten.Theta_r; + alpha = VanGenuchten.Alpha; + n = VanGenuchten.n; + m = VanGenuchten.m; + + theta_u = SoilVariables.Theta_U; + hh = SoilVariables.hh; + h = SoilVariables.h; + hh_frez = SoilVariables.hh_frez; + phi_s = SoilVariables.Phi_s; + lamda = SoilVariables.Lamda; + theta_l = SoilVariables.Theta_L; + h_frez = SoilVariables.h_frez; + + % get model settings + ModelSettings = io.getModelSettings(); + heatTerm = hh + hh_frez; + if ModelSettings.SWCC == 1 + if ModelSettings.SFCC == 1 + if hh >= -1 + dtheta_llh = 0; + elseif ModelSettings.Thmrlefc + if gamma_hh == 0 + dtheta_llh = 0; + else + subRoutine = 0; + dtheta_llh = conductivity.hydraulicConductivity.calculateDTheta(subRoutine, hh, theta_s, theta_r, theta_m, gamma_hh, theta_ll, theta_l, theta_uu, theta_u, hh, h, hh_frez, h_frez, phi_s, lamda, alpha, n, m); + end + elseif abs(hh - h) < 1e-3 + subRoutine = 1; + dtheta_llh = conductivity.hydraulicConductivity.calculateDTheta(subRoutine, hh, theta_s, theta_r, theta_m, gamma_hh, theta_ll, theta_l, theta_uu, theta_u, hh, h, hh_frez, h_frez, phi_s, lamda, alpha, n, m); + else + subRoutine = 3; + dtheta_llh = conductivity.hydraulicConductivity.calculateDTheta(subRoutine, hh, theta_s, theta_r, theta_m, gamma_hh, theta_ll, theta_l, theta_uu, theta_u, hh, h, hh_frez, h_frez, phi_s, lamda, alpha, n, m); + end + elseif hh <= -1e7 || theta_ll <= 0.06 + dtheta_llh = 0; + else + dtheta_llh = dtheta_uuh; + end + else + if hh <= -1e7 || hh + hh_frez <= -1e7 || hh + hh_frez >= phi_s + dtheta_llh = 0; + elseif ModelSettings.Thmrlefc + subRoutine = 4; + dtheta_llh = conductivity.hydraulicConductivity.calculateDTheta(subRoutine, heatTerm, theta_s, theta_r, theta_m, gamma_hh, theta_ll, theta_l, theta_uu, theta_u, hh, h, hh_frez, h_frez, phi_s, lamda, alpha, n, m); + elseif abs(hh - h) < 1e-3 + subRoutine = 5; + dtheta_llh = conductivity.hydraulicConductivity.calculateDTheta(subRoutine, heatTerm, theta_s, theta_r, theta_m, gamma_hh, theta_ll, theta_l, theta_uu, theta_u, hh, h, hh_frez, h_frez, phi_s, lamda, alpha, n, m); + else + subRoutine = 2; + dtheta_llh = conductivity.hydraulicConductivity.calculateDTheta(subRoutine, heatTerm, theta_s, theta_r, theta_m, gamma_hh, theta_ll, theta_l, theta_uu, theta_u, hh, h, hh_frez, h_frez, phi_s, lamda, alpha, n, m); + end + end +end diff --git a/src/+conductivity/+hydraulicConductivity/fixHeat.m b/src/+conductivity/+hydraulicConductivity/fixHeat.m new file mode 100644 index 00000000..4dcb8d97 --- /dev/null +++ b/src/+conductivity/+hydraulicConductivity/fixHeat.m @@ -0,0 +1,22 @@ +function [hh, hh_frez] = fixHeat(hh, hh_frez, Phi_s) + + % get model settings + ModelSettings = io.getModelSettings(); + + if ModelSettings.SWCC == 1 + if ModelSettings.SFCC ~= 1 + if hh > -1e-6 + hh = -1e-6; + elseif hh < -1e7 + hh = -1e7; + end + end + else + if hh >= Phi_s + hh = Phi_s; + elseif hh <= -1e7 + hh = -1e7; + hh_frez = -1e-6; + end + end +end diff --git a/src/+conductivity/calculateGasConductivity.m b/src/+conductivity/calculateGasConductivity.m new file mode 100644 index 00000000..ccca06b3 --- /dev/null +++ b/src/+conductivity/calculateGasConductivity.m @@ -0,0 +1,38 @@ +function k_g = calculateGasConductivity(InitialValues, TransportCoefficient, VanGenuchten, SoilVariables) + %{ + This is to calculate the intrinsic permeability of soil for gas flow. + Scanlon, B. R. (2000), Soil gas movement in unsaturated systems, in + Handbook of Soil Science, edited by M. E. Sumner, pp. A277-A319, CRC + Press, Boca Raton, Fla. + Zeng, Y., Su, Z., Wan, L. and Wen, i.: Numerical analysis of + air-water-heat flow in unsaturated soil: Is it necessary to consider + airflow in land surface models?, i. Geophys. Res. Atmos., 116(D20), + 20107, doi:10.1029/2011JD015835, 2011. + %} + + % get model settings + ModelSettings = io.getModelSettings(); + + % load Constants + Constants = io.define_constants(); + + k_g = InitialValues.k_g; + m = VanGenuchten.m; + + for i = 1:ModelSettings.NL + for j = 1:ModelSettings.nD + MN = i + j - 1; + Sa = SoilVariables.Theta_g(i, j) / SoilVariables.POR(i); + if ModelSettings.SWCC == 1 + k_g(i, j) = SoilVariables.Ks(i) * TransportCoefficient.MU_W(i, j) * (1 - Sa^0.5) * (1 - (1 - (1 - Sa^(1 / m(i))))^m(i))^2 / (Constants.g * Constants.RHOL); + else + k_g(i, j) = 0; + end + if ModelSettings.Soilairefc == 1 + k_g(i, j) = (SoilVariables.Ks(i) * TransportCoefficient.MU_W(i, j) * (1 - Sa^0.5) * (1 - (1 - (1 - Sa^(1 / m(i))))^m(i))^2 / (Constants.g * Constants.RHOL)) * 10^(-1 * SoilVariables.Imped(MN) * SoilVariables.Ratio_ice(i, j)); + else + k_g(i, j) = 0; + end + end + end +end diff --git a/src/+conductivity/calculateGasDispersivity.m b/src/+conductivity/calculateGasDispersivity.m new file mode 100644 index 00000000..4bbea955 --- /dev/null +++ b/src/+conductivity/calculateGasDispersivity.m @@ -0,0 +1,45 @@ +function GasDispersivity = calculateGasDispersivity(InitialValues, SoilVariables, P_gg, k_g) + %{ + This is to calculate the gas phase longitudinal dispersivity. + Zeng, Y., Su, Z., Wan, L. and Wen, i.: Numerical analysis of + air-water-heat flow in unsaturated soil: Is it necessary to consider + airflow in land surface models?, i. Geophys. Res. Atmos., 116(D20), + 20107, doi:10.1029/2011JD015835, 2011. + %} + + % get model settings + ModelSettings = io.getModelSettings(); + + % get model constants + Constants = io.define_constants(); + + V_A = InitialValues.V_A; + D_Vg = InitialValues.D_Vg; + Beta_g = InitialValues.Beta_g; + DPgDZ = InitialValues.DPgDZ; + Alpha_Lg = InitialValues.Alpha_Lg; + + for i = 1:ModelSettings.NL + for j = 1:ModelSettings.nD + Sa = 1 - SoilVariables.Se(i, j); + Beta_g(i, j) = (k_g(i, j) / Constants.MU_a); + Alpha_Lg(i, j) = 0.078 * (13.6 - 16 * Sa + 3.4 * Sa^5) * 100; + end + + DPgDZ(i) = (P_gg(i + 1) - P_gg(i)) / ModelSettings.DeltZ(i); + Alpha_LgBAR(i) = (Alpha_Lg(i, 1) + Alpha_Lg(i, 2)) / 2; + + Beta_gBAR = (Beta_g(i, 1) + Beta_g(i, 2)) / 2; + V_A(i) = -Beta_gBAR * DPgDZ(i); + if SoilVariables.KaT_Switch == 1 + D_Vg(i) = Alpha_LgBAR(i) * abs(V_A(i)); + else + D_Vg(i) = 0; + end + end + + GasDispersivity.D_Vg = D_Vg; + GasDispersivity.V_A = V_A; + GasDispersivity.Beta_g = Beta_g; + GasDispersivity.DPgDZ = DPgDZ; +end diff --git a/src/+conductivity/calculateHydraulicConductivity.m b/src/+conductivity/calculateHydraulicConductivity.m new file mode 100644 index 00000000..b6c8d2c4 --- /dev/null +++ b/src/+conductivity/calculateHydraulicConductivity.m @@ -0,0 +1,139 @@ +function SoilVariables = calculateHydraulicConductivity(SoilVariables, VanGenuchten, KIT, L_f) + %{ + This is to calculate the hydraulic conductivity of soil, based on + hydraulic conductivity models (like VG and others). + van Genuchten, M. T. (1980), A closed-form equation for predicting the + hydraulic conductivity of unsaturated soils, Soil Sci. Soc. Am. J., 44, + 892-898, + Zeng, Y., Su, Z., Wan, L. and Wen, J.: Numerical analysis of + air-water-heat flow in unsaturated soil: Is it necessary to consider + airflow in land surface models?, J. Geophys. Res. Atmos., 116(D20), + 20107, doi:10.1029/2011JD015835, 2011. + %} + + % get model settings + ModelSettings = io.getModelSettings(); + + % load Constants + Constants = io.define_constants(); + + function [sliced] = sliceVector(structure, lengthX, ix) + %{ + slice fields of a structure if they are vectors. + %} + + for fieldName = fieldnames(structure)' + fieldData = structure.(fieldName{1}); + if isvector(fieldData) && length(fieldData) >= (lengthX) + sliced.(fieldName{1}) = fieldData(ix); + else + sliced.(fieldName{1}) = fieldData; + end + end + end + + function [sliced] = sliceMatrix(structure, lengthX, lengthY, ix, iy) + %{ + slice fields of a structure if they are matrices. + %} + + for fieldName = fieldnames(structure)' + fieldData = structure.(fieldName{1}); + fieldDataSize = size(fieldData); + if ismatrix(fieldData) + if fieldDataSize(1) >= lengthX && fieldDataSize(2) >= lengthY + sliced.(fieldName{1}) = fieldData(ix, iy); + elseif fieldDataSize(1) >= lengthY && fieldDataSize(2) >= lengthX + sliced.(fieldName{1}) = fieldData(iy, ix); + else + sliced.(fieldName{1}) = fieldData; + end + else + sliced.(fieldName{1}) = fieldData; + end + end + end + + for i = 1:ModelSettings.NL + % slice vectors + VG = sliceVector(VanGenuchten, ModelSettings.NL, i); + + for j = 1:ModelSettings.nD + MN = i + j - 1; + + SV = SoilVariables; + + % slice + % only for these variables, i is used, for the rest MN. + SV.POR = SoilVariables.POR(i); + SV.XCAP = SoilVariables.XCAP(i); + SV.Ks = SoilVariables.Ks(i); + + lengthX = ModelSettings.NL + j - 1; + SV = sliceVector(SV, lengthX, MN); + SV = sliceMatrix(SV, lengthX, ModelSettings.nD, MN, j); + + [hh, hh_frez] = conductivity.hydraulicConductivity.fixHeat(SV.hh, SV.hh_frez, SV.Phi_s); + SV.hh = hh; + SV.hh_frez = hh_frez; + + Gamma_hh = conductivity.hydraulicConductivity.calculateGamma_hh(SV.hh); + Theta_m = conductivity.hydraulicConductivity.calculateTheta_m(Gamma_hh, VG, SV.POR); + Theta_UU = conductivity.hydraulicConductivity.calculateTheta_UU(Theta_m, Gamma_hh, SV, VG); + + % circular calculation of Theta_II! See issue 181, item 3 + % Theta_II is soil ice content, + % Theta_LL is liquid water content, + % Theta_UU is the total water content before soil freezing. The + % 'Theta_UU' is set as saturation. + Theta_II = conductivity.hydraulicConductivity.calculateTheta_II(SV.TT, SV.XCAP, SV.hh, SV.Theta_II); + Theta_LL = conductivity.hydraulicConductivity.calculateTheta_LL(Theta_UU, Theta_II, Theta_m, Gamma_hh, SV, VG); + Theta_II = (Theta_UU - Theta_LL) * Constants.RHOL / Constants.RHOI; % ice water contentTheta_II + + DTheta_UUh = conductivity.hydraulicConductivity.calculateDTheta_UUh(Theta_UU, Theta_m, Theta_LL, Gamma_hh, SV, VG); + DTheta_LLh = conductivity.hydraulicConductivity.calcuulateDTheta_LLh(DTheta_UUh, Theta_m, Theta_UU, Theta_LL, Gamma_hh, SV, VG); + Se = conductivity.hydraulicConductivity.calculateSe(Theta_LL, Gamma_hh, SV); + + % Ratio_ice used in Condg_k_g.m + if Theta_UU ~= 0 + Ratio_ice = Constants.RHOI * Theta_II / (Constants.RHOL * Theta_UU); % ice ratio + else + Ratio_ice = 0; + end + if KIT + % The calculation of MU_W repeated in CondL_Tdisp and used in Condg_k_g + if SV.TT < -20 + MU_W = 3.71e-2; + elseif SV.TT > 150 + MU_W = 1.81e-3; + else + MU_W = Constants.MU_W0 * exp(Constants.MU1 / (8.31441 * (SV.TT + 133.3))); + end + + KL_h = conductivity.hydraulicConductivity.calculateKL_h(MU_W, Se, SV.Ks, VG.m); + + if Gamma_hh ~= 1 + KfL_h = KL_h * 10^(-1 * SV.Imped * Ratio_ice); % hydraulic conductivity for freezing soil + else + KfL_h = KL_h * 10^(-1 * SV.Imped * Ratio_ice); % hydraulic conductivity for freezing soil + end + KfL_T = helpers.heaviside1(SV.TT_CRIT - (SV.TT + ModelSettings.T0)) * L_f * 1e4 / (Constants.g * (ModelSettings.T0)); % thermal consider for freezing soil + else + KL_h = 0; + KfL_h = 0; + KfL_T = 0; + end + SoilVariables.KL_h(i, j) = KL_h; + SoilVariables.KfL_h(i, j) = KfL_h; + SoilVariables.KfL_T(i, j) = KfL_T; + SoilVariables.Theta_LL(i, j) = Theta_LL; + SoilVariables.DTheta_LLh(i, j) = DTheta_LLh; + SoilVariables.Theta_II(i, j) = Theta_II; + SoilVariables.Theta_UU(i, j) = Theta_UU; + SoilVariables.DTheta_UUh(i, j) = DTheta_UUh; + SoilVariables.Se(i, j) = Se; + SoilVariables.Ratio_ice(i, j) = Ratio_ice; + SoilVariables.Gamma_hh(MN) = Gamma_hh; + end + end +end diff --git a/src/+conductivity/calculateSoilThermalProperites.m b/src/+conductivity/calculateSoilThermalProperites.m new file mode 100644 index 00000000..596271a7 --- /dev/null +++ b/src/+conductivity/calculateSoilThermalProperites.m @@ -0,0 +1,147 @@ +function [ETCON, EHCAP, TETCON, EfTCON, ZETA] = calculateSoilThermalProperites(InitialValues, ThermalConductivity, SoilVariables, VanGenuchten, DRHOVT, L, RHOV) + %{ + This is used to calculate Heat Capacity and Thermal Conductivity. + %} + + HCAP = ThermalConductivity.HCAP; + SF = ThermalConductivity.SF; + TCA = ThermalConductivity.TCA; + GA1 = ThermalConductivity.GA1; + GA2 = ThermalConductivity.GA2; + GB1 = ThermalConductivity.GB1; + GB2 = ThermalConductivity.GB2; + HCD = ThermalConductivity.HCD; + ZETA0 = ThermalConductivity.ZETA0; + CON0 = ThermalConductivity.CON0; + PS1 = ThermalConductivity.PS1; + PS2 = ThermalConductivity.PS2; + TCON_dry = ThermalConductivity.TCON_dry; + TPS1 = ThermalConductivity.TPS1; + TPS2 = ThermalConductivity.TPS2; + FEHCAP = ThermalConductivity.FEHCAP; + TCON0 = ThermalConductivity.TCON0; + TCON_s = ThermalConductivity.TCON_s; + + XWILT = SoilVariables.XWILT; + XK = SoilVariables.XK; + TT = SoilVariables.TT; + POR = SoilVariables.POR; + Theta_LL = SoilVariables.Theta_LL; + Theta_V = SoilVariables.Theta_V; + Theta_II = SoilVariables.Theta_II; + XSOC = SoilVariables.XSOC; + + D_A = InitialValues.D_A; + + Theta_s = VanGenuchten.Theta_s; + + % load Constants + Constants = io.define_constants(); + + % get model settings + ModelSettings = io.getModelSettings(); + + MN = 0; + TCON(1) = 1.37e-3 * 4.182; + for i = 1:ModelSettings.NL + for j = 1:ModelSettings.nD + MN = i + j - 1; + XXX = Theta_LL(i, j); + XII = Theta_II(i, j) * Constants.RHOI / Constants.RHOL; % MODIFIED 201905 + if Theta_LL(i, j) < XK(i) + XXX = XK(i); + end + + if XXX < XWILT(i) + SF(2) = GA2 + GB2(i) * XXX; + else + SF(2) = GA1 + GB1(i) * XXX; + end + D_A(MN) = 0.229 * (1 + TT(MN) / 273)^1.75; + TCON(2) = TCA + D_A(MN) * L(MN) * DRHOVT(MN); % TCA+D_A(MN)*L(MN)*DRHOVT(MN); % Revised from ""(D_V(i,j)*Eta(i,j)+D_Vg(i))*L(MN)*DRHOVT(MN) + TCON(6) = 5.2e-3 * 4.182; % ice thermal conductivity + SF(6) = 0.125; + TARG(2) = TCON(2) / TCON(1) - 1; + TARG(6) = TCON(6) / TCON(1) - 1; + GRAT(2) = 0.667 / (1 + TARG(2) * SF(2)) + 0.333 / (1 + TARG(2) * (1 - 2 * SF(2))); + GRAT(6) = 0.667 / (1 + TARG(6) * SF(6)) + 0.333 / (1 + TARG(6) * (1 - 2 * SF(6))); + ETCON(i, j) = (PS1(i) + XXX * TCON(1) + (POR(i) - XXX - XII) * GRAT(2) * TCON(2) + XII * GRAT(6) * TCON(6)) / (PS2(i) + XXX + (POR(i) - XXX - XII) * GRAT(2) + XII * GRAT(6)); + ZETA(i, j) = GRAT(2) / (GRAT(2) * (POR(i) - XXX) + XXX + PS2(i)); % ita_T enhancement factor + ZETA(i, j) = GRAT(2) / (GRAT(2) * (POR(i) - XXX - XII) + XXX + PS2(i) + XII * GRAT(6)); + if Theta_LL(i, j) == XXX + EHCAP(i, j) = HCD(i) + HCAP(1) * Theta_LL(i, j); + EHCAP(i, j) = EHCAP(i, j) + (0.448 * RHOV(MN) * 4.182 + HCAP(2)) * Theta_V(i, j) + HCAP(6) * XII; % The Calorie should be converted as i + else + ZETA(i, j) = ZETA0(i) + (ZETA(i, j) - ZETA0(i)) * Theta_LL(i, j) / XXX; + ETCON(i, j) = CON0(i) + (ETCON(i, j) - CON0(i)) * Theta_LL(i, j) / XXX; + EHCAP(i, j) = HCD(i) + HCAP(1) * Theta_LL(i, j); + EHCAP(i, j) = EHCAP(i, j) + (0.448 * RHOV(MN) * 4.182 + HCAP(2)) * Theta_V(i, j) + HCAP(6) * XII; % The Calorie should be converted as i + end + %%%%%% simpliefied Vries method by Tian 2016 + TSF(2) = 0.333 * (1 - (POR(i) - XXX - XII) / POR(i)); + TSF(6) = 0.333 * (1 - XII / POR(i)); + TSF(7) = 0.5; % shape factor for organic matter + D_A(MN) = 0.229 * (1 + TT(MN) / 273)^1.75; + TCON(2) = TCA + D_A(MN) * L(MN) * DRHOVT(MN); % TCA+D_A(MN)*L(MN)*DRHOVT(MN); % Revised from ""(D_V(i,j)*Eta(i,j)+D_Vg(i))*L(MN)*DRHOVT(MN) + TCON(6) = 5.2e-3 * 4.182; % ice conductivity + TCON(7) = 0.25e-2; % soil organic matter conductivity + TARG(2) = TCON(2) / TCON(1) - 1; + TARG(6) = TCON(6) / TCON(1) - 1; + TARG(7) = TCON(7) / TCON(1) - 1; + TTGRAT(2) = 0.667 / (1 + TARG(2) * TSF(2)) + 0.333 / (1 + TARG(2) * (1 - 2 * TSF(2))); + TTGRAT(6) = 0.667 / (1 + TARG(6) * TSF(6)) + 0.333 / (1 + TARG(6) * (1 - 2 * TSF(6))); + TTGRAT(7) = 0.667 / (1 + TARG(7) * TSF(7)) + 0.333 / (1 + TARG(7) * (1 - 2 * TSF(7))); + if XSOC(i) == 0 + TETCON(i, j) = (TPS1(i) + XXX * TCON(1) + (POR(i) - XXX - XII) * TTGRAT(2) * TCON(2) + XII * TTGRAT(6) * TCON(6)) / (TPS2(i) + XXX + (POR(i) - XXX - XII) * TTGRAT(2) + XII * TTGRAT(6)); + else + TETCON(i, j) = (TPS1(i) + XXX * TCON(1) + (POR(i) - XXX - XII) * TTGRAT(2) * TCON(2) + XII * TTGRAT(6) * TCON(6) + XSOC(i) * TTGRAT(7) * TCON(7)) / (TPS2(i) + XXX + (POR(i) - XXX - XII) * TTGRAT(2) + XII * TTGRAT(6) + XSOC(i) * TTGRAT(7)); + end + TZETA(i, j) = TTGRAT(2) / (TTGRAT(2) * (POR(i) - XXX) + XXX + TPS2(i)); % ita_T enhancement factor + if Theta_LL(i, j) == XXX + EHCAP(i, j) = HCD(i) + HCAP(1) * Theta_LL(i, j); + EHCAP(i, j) = EHCAP(i, j) + (0.448 * RHOV(MN) * 4.182 + HCAP(2)) * Theta_V(i, j) + HCAP(6) * XII; % The Calorie should be converted as i + else + ZETA(i, j) = ZETA0(i) + (ZETA(i, j) - ZETA0(i)) * Theta_LL(i, j) / XXX; + TETCON(i, j) = TCON0(i) + (TETCON(i, j) - TCON0(i)) * Theta_LL(i, j) / XXX; + EHCAP(i, j) = HCD(i) + HCAP(1) * Theta_LL(i, j); + EHCAP(i, j) = EHCAP(i, j) + (0.448 * RHOV(MN) * 4.182 + HCAP(2)) * Theta_V(i, j) + HCAP(6) * XII; % The Calorie should be converted as i + end + %%%%%%%%%%%%%%%%%% Farouki thermal parameter method %%%%%%%%%%% + %%%%%%%%%%%%%%%%% change 21 July 2017 switch of thermal conductivity%%%%%% + if ModelSettings.ThermCond == 4 %|| ThermCond==2 + % FEHCAP(i,j)=(2.128*Theta_sa(i)+2.385*Theta_cl(i))/(Theta_sa(i)+Theta_cl(i))*1e6; %i m-3 K-1 + % TCON_s(i)=(8.8*Theta_sa(i)+2.92*Theta_cl(i))/(Theta_sa(i)+Theta_cl(i)); % W m-1 K-1 + FFEHCAP(i, j) = (FEHCAP(i) * (1 - POR(i)) + 1000 * 4217.7 * Theta_LL(i, j) + 917 * 2117.27 * XII) * 1e-6; % The Calorie should be converted as i + EHCAP(i, j) = FFEHCAP(i, j); + end + %%%%%%%%%%%%%%%%%% Johansen thermal conductivity method %%%%%%%%%%% + TCON_qtz = 7.7; + TCON_o = 2.0; + TCON_L = 0.57; % Theta_qtz(i)=0.47; % thermal conductivities of soil quartz, other soil particles and water; unit W m-1 K-1 + if ModelSettings.ThermCond == 2 + TCON_s(i) = TCON_qtz^(Theta_qtz(i)) * TCON_o^(1 - Theta_qtz(i)); % solid soil thermal conductivity Unit W m-1 K-1 + end + TCON_i(i) = 2.2; + %%%%% Peters-Lidard et al. 1998 frozen soil thermal conductivity method %%%%%%% + Satr(i, j) = (Theta_LL(i, j) + XII) / Theta_s(i); + if Theta_II(i, j) <= 0 + if Theta_LL(i, j) / POR(i) > 0.1 + K_e(i, j) = log10(Theta_LL(i, j) / POR(i)) + 1; % Kersten coefficient, weighting dry and wet thermal conductivity + elseif Theta_LL(i, j) / POR(i) > 0.05 + K_e(i, j) = 0.7 * log10(Theta_LL(i, j) / POR(i)) + 1; % Kersten coefficient, weighting dry and wet thermal conductivity + else + K_e(i, j) = 0; + end + Coef_k = 1.9; % [4.6 3.55 1.9; 1.7 0.95 0.85]; + K_e1(i, j) = Coef_k * (Theta_LL(i, j) / POR(i)) / (1 + (Coef_k - 1) * (Theta_LL(i, j) / POR(i))); % Kersten coefficient, weighting dry and wet thermal conductivity + ALPHA = 0.9; % [1.05,0.9,0.58]; + K_e2(i, j) = exp(ALPHA * (1 - (Theta_LL(i, j) / POR(i))^(ALPHA - 1.33))); % Lu and Ren 2007; ALPHA=1.05,0.9,0.58 + TCON_sat(i, j) = TCON_s(i)^(1 - Theta_s(i)) * TCON_L^(Theta_s(i)); % saturated soil thermal conductivity Unit W m-1 K-1 + EfTCON(i, j) = K_e(i, j) * (TCON_sat(i, j) - TCON_dry(i)) + TCON_dry(i); + else + K_e(i, j) = (Theta_LL(i, j) + Theta_II(i, j)) / POR(i); % +Theta_II(i,j) + TCON_sat(i, j) = TCON_s(i)^(1 - Theta_s(i)) * TCON_i(i)^(Theta_s(i) - Theta_LL(i, j)) * TCON_L^(Theta_LL(i, j)); % saturated soil thermal conductivity Unit W m-1 K-1 + EfTCON(i, j) = K_e(i, j) * (TCON_sat(i, j) - TCON_dry(i)) + TCON_dry(i); + end + end + end diff --git a/src/+conductivity/calculateThermalConductivityCapacity.m b/src/+conductivity/calculateThermalConductivityCapacity.m new file mode 100644 index 00000000..cd9bf097 --- /dev/null +++ b/src/+conductivity/calculateThermalConductivityCapacity.m @@ -0,0 +1,55 @@ +function ThermalConductivityCapacity = calculateThermalConductivityCapacity(InitialValues, ThermalConductivity, SoilVariables, VanGenuchten, DRHOVT, L, RHOV) + %{ + This is to calculate thermal conductivity and thermal capacity. + Chung, S. O., and R. Horton (1987), Soil heat and water flow with a + partial surface mulch, Water Resour. Res., 23(12), 2175-2186, + Zeng, Y., Su, Z., Wan, L. and Wen, J.: Numerical analysis of + air-water-heat flow in unsaturated soil: Is it necessary to consider + airflow in land surface models?, J. Geophys. Res. Atmos., 116(D20), + 20107, doi:10.1029/2011JD015835, 2011. + %} + + % get model settings + ModelSettings = io.getModelSettings(); + + % load Constants + Constants = io.define_constants(); + + Theta_LL = SoilVariables.Theta_LL; + c_unsat = InitialValues.c_unsat; + Lambda_eff = InitialValues.Lambda_eff; + RHO_bulk = ThermalConductivity.RHO_bulk; + + if ModelSettings.ThmrlCondCap == 1 + [ETCON, EHCAP, TETCON, EfTCON, ZETA] = conductivity.calculateSoilThermalProperites(InitialValues, ThermalConductivity, SoilVariables, VanGenuchten, DRHOVT, L, RHOV); + end + + for i = 1:ModelSettings.NL + for j = 1:ModelSettings.nD + if ModelSettings.ThmrlCondCap == 1 + if ModelSettings.ThermCond == 1 + Lambda_eff(i, j) = ETCON(i, j); + elseif ModelSettings.ThermCond == 2 + Lambda_eff(i, j) = EfTCON(i, j) / 100; + elseif ModelSettings.ThermCond == 3 + Lambda_eff(i, j) = TETCON(i, j); + elseif ModelSettings.ThermCond == 4 + Lambda_eff(i, j) = EfTCON(i, j) / 100; + end + if Lambda_eff(i, j) <= 0 + Lambda_eff(i, j) = 0.0008; + elseif Lambda_eff(i, j) > 0.02 + Lambda_eff(i, j) = 0.02; + end + c_unsat(i, j) = EHCAP(i, j); + else + MN = i + j - 1; + Lambda_eff(i, j) = Constants.Lambda1 + Constants.Lambda2 * Theta_LL(i, j) + Constants.Lambda3 * Theta_LL(i, j)^0.5; + c_unsat(i, j) = 837 * RHO_bulk / 1000 + Theta_LL(i, j) * Constants.c_L + Theta_g(i, j) * (RHODA(MN) * Constants.c_a + RHOV(MN) * Constants.c_V); + end + end + end + ThermalConductivityCapacity.c_unsat = c_unsat; + ThermalConductivityCapacity.Lambda_eff = Lambda_eff; + ThermalConductivityCapacity.ZETA = ZETA; +end diff --git a/src/+conductivity/calculateTransportCoefficient.m b/src/+conductivity/calculateTransportCoefficient.m new file mode 100644 index 00000000..17b399d6 --- /dev/null +++ b/src/+conductivity/calculateTransportCoefficient.m @@ -0,0 +1,74 @@ +function TransportCoefficient = calculateTransportCoefficient(InitialValues, SoilVariables, VanGenuchten, Delt_t) + %{ + This is to calculate the transport coefficient for absorbed liquid flow + due to temperature gradient. + Groenevelt, P. H., and B. D. Kay (1974), On the interaction of water and + heat transport in frozen and unfrozen soils: II. The liquid phase, Soil + Sci. Soc. Am. Proc., 38, 400-404, + Zeng, Y., Su, Z., Wan, L. and Wen, i.: Numerical analysis of + air-water-heat flow in unsaturated soil: Is it necessary to consider + airflow in land surface models?, i. Geophys. Res. Atmos., 116(D20), + 20107, doi:10.1029/2011JD015835, 2011. + %} + + % Tortuosity Factor is a reverse of the tortuosity. In "L_WT", tortuosity + % should be used. That is why "f0" is in the numerator. Note: f0 in L_WT + % has been changed as 1.5. The Rv in MU_W should be 8.31441 i/mol.K. + + TransportCoefficient.W = InitialValues.W; + TransportCoefficient.WW = InitialValues.WW; + TransportCoefficient.MU_W = InitialValues.MU_W; + TransportCoefficient.D_Ta = InitialValues.D_Ta; + + % get model settings + ModelSettings = io.getModelSettings(); + + % load Constants + Constants = io.define_constants(); + + for i = 1:ModelSettings.NL + for j = 1:ModelSettings.nD + MN = i + j - 1; + if ModelSettings.W_Chg == 0 + W = 0; + WW = 0; + WARG = SoilVariables.Theta_LL(i, j) * 10^7 / ModelSettings.SSUR; + if WARG < 80 + W = Constants.W0 * exp(-WARG); + WW = Constants.W0 * exp(-WARG); + end + else + W = -0.2932 * (SoilVariables.h(MN)) / 1000; + WW = -0.2932 * (SoilVariables.hh(MN)) / 1000; + if W > 2932 + W = 2932; + end + if WW > 2932 + WW = 2932; + end + end + f0 = SoilVariables.Theta_g(i, j)^(7 / 3) / VanGenuchten.Theta_s(i)^2; + + % SSUR and RHO_bulk could also be set as an array to consider more + % than one soil type; + H_W = Constants.RHOL * WW * (SoilVariables.Theta_LL(i, j) - SoilVariables.Theta_L(i, j)) / ((ModelSettings.SSUR / Constants.RHO_bulk) * Delt_t); + if SoilVariables.TT(MN) < -20 + MU_W = 3.71e-2; + elseif SoilVariables.TT(MN) > 150 + MU_W = 1.81e-3; + else + MU_W = Constants.MU_W0 * exp(Constants.MU1 / (8.31441 * (SoilVariables.TT(MN) + 133.3))); + end + L_WT = f0 * 1e7 * 1.5550e-13 * SoilVariables.POR(i) * H_W / (Constants.b * MU_W); + if SoilVariables.KLT_Switch == 1 + D_Ta = L_WT / (Constants.RHOL * (SoilVariables.TT(MN) + 273.15)); + else + D_Ta = 0; + end + TransportCoefficient.W(i, j) = W; + TransportCoefficient.WW(i, j) = WW; + TransportCoefficient.D_Ta(i, j) = D_Ta; + TransportCoefficient.MU_W(i, j) = MU_W; + end + end +end diff --git a/src/+conductivity/calculateVaporVariables.m b/src/+conductivity/calculateVaporVariables.m new file mode 100644 index 00000000..6240e709 --- /dev/null +++ b/src/+conductivity/calculateVaporVariables.m @@ -0,0 +1,58 @@ +function VaporVariables = calculateVaporVariables(InitialValues, SoilVariables, VanGenuchten, ThermalConductivityCapacity, TT) + %{ + This is to calculate vapor diffusivity, vapor dispersivity, and vapor + enhancement factor. + Zeng, Y., Su, Z., Wan, L. and Wen, i.: Numerical analysis of + air-water-heat flow in unsaturated soil: Is it necessary to consider + airflow in land surface models?, i. Geophys. Res. Atmos., 116(D20), + 20107, doi:10.1029/2011JD015835, 2011. + Yu, L., Zeng, Y. and Su, Z.: STEMMUS-UEB v1.0.0: integrated modeling of + snowpack and soil water and energy transfer with three complexity levels + of soil physical processes, Geosci. Model Dev, 14, 7345-7376, + doi:10.5194/gmd-14-7345-2021, 2021. + %} + + Theta_LL = SoilVariables.Theta_LL; + Theta_g = SoilVariables.Theta_g; + POR = SoilVariables.POR; + XK = SoilVariables.XK; + Theta_UU = SoilVariables.Theta_UU; + Theta_s = VanGenuchten.Theta_s; + + D_V = InitialValues.D_V; + Eta = InitialValues.Eta; + D_A = InitialValues.D_A; + + % get model settings + ModelSettings = io.getModelSettings(); + + for i = 1:ModelSettings.NL + for j = 1:ModelSettings.nD + MN = i + j - 1; + if ModelSettings.ThmrlCondCap + if Theta_UU(i, j) < XK(i) + EnhnLiqIsland = POR(i); + else + EnhnLiqIsland = Theta_g(i, j) * (1 + Theta_UU(i, j) / (POR(i) - XK(i))); + end + + f0 = Theta_g(i, j)^(7 / 3) / Theta_s(i)^2; + D_A(MN) = 0.229 * (1 + TT(MN) / 273)^1.75; + + if SoilVariables.DVT_Switch == 1 + Eta(i, j) = ThermalConductivityCapacity.ZETA(i, j) * EnhnLiqIsland / (f0 * Theta_g(i, j)); + else + Eta(i, j) = 0; + end + else + f0 = Theta_g(i, j)^0.67; + D_A(MN) = 1e4 * 2.12 * 10^(-5) * (1 + TT(MN) / 273.15)^2; + Eta(i, j) = 8 + 3 * Theta_LL(i, j) / Theta_s(i) - 7 * exp(-((1 + 2.6 / sqrt(ModelSettings.fc)) * Theta_LL(i, j) / Theta_s(i))^3); + end + D_V(i, j) = f0 * Theta_g(i, j) * D_A(MN); + end + end + + VaporVariables.D_V = D_V; + VaporVariables.Eta = Eta; +end diff --git a/src/+init/applySoilHeteroWithInitialFreezing.m b/src/+init/applySoilHeteroWithInitialFreezing.m index c77965ab..07e97ae5 100644 --- a/src/+init/applySoilHeteroWithInitialFreezing.m +++ b/src/+init/applySoilHeteroWithInitialFreezing.m @@ -13,7 +13,7 @@ SoilVariables.SAVEh(i) = SoilVariables.h(i); SoilVariables.SAVEhh(i) = SoilVariables.hh(i); - [SoilVariables.Gama_h, SoilVariables.Gama_hh] = init.updateGamaH(i, SoilVariables); + [SoilVariables.Gamma_h, SoilVariables.Gamma_hh] = init.updateGammaH(i, SoilVariables); if ModelSettings.Thmrlefc == 1 SoilVariables.TT(i) = SoilVariables.T(i); diff --git a/src/+init/defineSoilVariables.m b/src/+init/defineSoilVariables.m index 64828842..5cc7a4f5 100644 --- a/src/+init/defineSoilVariables.m +++ b/src/+init/defineSoilVariables.m @@ -13,7 +13,8 @@ 'P_g', 'P_gg', 'h', 'T', 'TT', 'h_frez', 'Theta_L', ... 'Theta_LL', 'Theta_V', 'Theta_g', 'Se', 'KL_h', ... 'DTheta_LLh', 'KfL_T', 'Theta_II', 'Theta_I', ... - 'Theta_UU', 'Theta_U', 'TT_CRIT', 'KfL_h', 'DTheta_UUh' + 'Theta_UU', 'Theta_U', 'TT_CRIT', 'KfL_h', 'DTheta_UUh', ... + 'Ratio_ice' }; for field = soil_fields SoilVariables.(field{1}) = InitialValues.(field{1}); diff --git a/src/+init/updateGamaH.m b/src/+init/updateGamaH.m deleted file mode 100644 index 26002264..00000000 --- a/src/+init/updateGamaH.m +++ /dev/null @@ -1,22 +0,0 @@ -function [Gama_h, Gama_hh] = updateGamaH(i, SoilVariables) - - % get soil constants - SoilConstants = io.getSoilConstants(); - - hd = SoilConstants.hd; - hm = SoilConstants.hm; - hh = SoilVariables.hh; - - if abs(hh(i)) >= abs(hd) - Gama_h(i) = 0; - Gama_hh(i) = Gama_h(i); - elseif abs(hh(i)) >= abs(hm) - % Gama_h(i)=1-log(abs(hh(i)))/log(abs(hm)); - Gama_h(i) = log(abs(hd) / abs(hh(i))) / log(abs(hd) / abs(hm)); - Gama_hh(i) = Gama_h(i); - else - Gama_h(i) = 1; - Gama_hh(i) = Gama_h(i); - end - -end diff --git a/src/+init/updateGammaH.m b/src/+init/updateGammaH.m new file mode 100644 index 00000000..fe012cbe --- /dev/null +++ b/src/+init/updateGammaH.m @@ -0,0 +1,22 @@ +function [Gamma_h, Gamma_hh] = updateGammaH(i, SoilVariables) + + % get soil constants + SoilConstants = io.getSoilConstants(); + + hd = SoilConstants.hd; + hm = SoilConstants.hm; + hh = SoilVariables.hh; + + if abs(hh(i)) >= abs(hd) + Gamma_h(i) = 0; + Gamma_hh(i) = Gamma_h(i); + elseif abs(hh(i)) >= abs(hm) + % Gamma_h(i)=1-log(abs(hh(i)))/log(abs(hm)); + Gamma_h(i) = log(abs(hd) / abs(hh(i))) / log(abs(hd) / abs(hm)); + Gamma_hh(i) = Gamma_h(i); + else + Gamma_h(i) = 1; + Gamma_hh(i) = Gamma_h(i); + end + +end diff --git a/src/+io/getModelSettings.m b/src/+io/getModelSettings.m index 5560f754..951458cc 100644 --- a/src/+io/getModelSettings.m +++ b/src/+io/getModelSettings.m @@ -49,6 +49,7 @@ % Other settings ModelSettings.rwuef = 1; ModelSettings.rroot = 1.5 * 1e-3; + ModelSettings.SFCC = 1; % Arraies for calculating boundary flux; ModelSettings.SAVE = zeros(3, 3, 3); diff --git a/src/CondL_T.m b/src/CondL_T.m deleted file mode 100644 index 6a8ccde6..00000000 --- a/src/CondL_T.m +++ /dev/null @@ -1,16 +0,0 @@ -function [KL_T] = CondL_T(NL) - - MN = 0; - for ML = 1:NL - for ND = 1:2 - MN = ML + ND - 1; - % if KLT_Switch==1 - KL_T(ML, ND) = 0; % KL_h(ML,ND)*((hh(MN)*GWT)/Gamma0)*(-0.1425-4.76*10^(-4)*TT(MN)); %(50+2.75*TT(MN))/((50+2.75*20));% - % else - % KL_T(ML,ND)=0; - % end - end - end - - %%%%%%%% Unit of KL_T is determined by KL_h, which is subsequently %%%%%%%% - %%%%%%%% determined by Ks set at the beginning. %%%%%%%%%%%%%%%%%%%%%%%%%%% diff --git a/src/CondL_Tdisp.m b/src/CondL_Tdisp.m deleted file mode 100644 index 746dcc31..00000000 --- a/src/CondL_Tdisp.m +++ /dev/null @@ -1,58 +0,0 @@ -function [W, WW, MU_W, D_Ta] = CondL_Tdisp(InitialValues, POR, Theta_LL, Theta_L, SSUR, RHOL, TT, Theta_s, h, hh, W_Chg, NL, nD, Delt_t, Theta_g, KLT_Switch) - %% - - W = InitialValues.W; - WW = InitialValues.WW; - MU_W = InitialValues.MU_W; - D_Ta = InitialValues.D_Ta; - - % load Constants - Constants = io.define_constants(); - - MN = 0; - for ML = 1:NL - J = ML; - for ND = 1:nD - MN = ML + ND - 1; - if W_Chg == 0 - W(ML, ND) = 0; - WW(ML, ND) = 0; - WARG = Theta_LL(ML, ND) * 10^7 / SSUR; - if WARG < 80 - W(ML, ND) = Constants.W0 * exp(-WARG); - WW(ML, ND) = Constants.W0 * exp(-WARG); - end - else - W(ML, ND) = -0.2932 * (h(MN)) / 1000; % 0;% +h_frez(MN) %%% J.g^-1---Original J.Kg^-1, now is divided by 1000. - WW(ML, ND) = -0.2932 * (hh(MN)) / 1000; % 0;%+hh_frez(MN) - if W(ML, ND) >= 2932 - W(ML, ND) = 2932; - end - if WW(ML, ND) >= 2932 - WW(ML, ND) = 2932; - end - end - f0(ML, ND) = Theta_g(ML, ND)^(7 / 3) / Theta_s(J)^2; % Theta_g(ML,ND)^0.67; - H_W(ML, ND) = RHOL * WW(ML, ND) * (Theta_LL(ML, ND) - Theta_L(ML, ND)) / ((SSUR / Constants.RHO_bulk) * Delt_t); % 1e3; % 1e-4J cm-2---> g s-2 ; SSUR and RHO_bulk could also be set as an array to consider more than one soil type; - if TT(MN) < -20 - MU_W(ML, ND) = 3.71e-2; % L_WT(ML,ND)=0; - elseif TT(MN) > 150 - MU_W(ML, ND) = 1.81e-3; - else - MU_W(ML, ND) = Constants.MU_W0 * exp(Constants.MU1 / (8.31441 * (TT(MN) + 133.3))); - end - L_WT(ML, ND) = f0(ML, ND) * 1e7 * 1.5550e-13 * POR(J) * H_W(ML, ND) / (Constants.b * MU_W(ML, ND)); % kgm^-1s^-1 --> 10 g.cm^-1.s^-1; J.cm^-2---> kg.m^2.s^-2.cm^-2--> 1e7g.cm^2.s^-2.cm^-2 - - if KLT_Switch == 1 - D_Ta(ML, ND) = L_WT(ML, ND) / (RHOL * (TT(MN) + 273.15)); % 0; %0;%0; % - else - D_Ta(ML, ND) = 0; - end - end - end - %% Tortuosity Factor is a reverse of the tortuosity. In "L_WT", tortuosity should be used. That is why "f0" is in the numerator.%%% - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% NOTICE: f0 in L_WT has been changed as 1.5 %%%%%%%%%%% - %%%%%% kg.m^2.s^-2.cm^-2.kg.m^-3 --> 1e7g.cm^2.s^-2.cm^-2.g.cm^-3 - %%%%%% Unit of L_WT IS (kgm^-1s^-1)=-------------------------- cm^2=;(1.5548e-013 cm^2); Converting meter to centimeter here by multipling UnitC - %%%%%% m. kg.m^-1.s^-1 --> cm. g.cm^-1.s^-1 - %%%%%% Please note that the Rv in MU_W should be 8.31441 J/mol.K. %%%%%%%% diff --git a/src/CondL_h.m b/src/CondL_h.m deleted file mode 100644 index 0ade467c..00000000 --- a/src/CondL_h.m +++ /dev/null @@ -1,368 +0,0 @@ -function [Theta_LL, Se, KfL_h, KfL_T, DTheta_LLh, hh, hh_frez, Theta_UU, DTheta_UUh, Theta_II, KL_h] = CondL_h(SoilVariables, Theta_r, Theta_s, Alpha, hh, hh_frez, h_frez, n, m, Ks, NL, Theta_L, h, KIT, TT, Thmrlefc, POR, SWCC, Theta_U, XCAP, Phi_s, RHOI, RHOL, Lamda, Imped, L_f, g, T0, TT_CRIT, Theta_II, KfL_h, KfL_T, KL_h, Theta_UU, Theta_LL, DTheta_LLh, DTheta_UUh, Se) - - % load Constants - Constants = io.define_constants(); - - % get soil constants for StartInit - SoilConstants = io.getSoilConstants(); - hd = SoilConstants.hd; - hm = SoilConstants.hm; - - Gama_hh = SoilVariables.Gama_hh; - - SFCC = 1; - MN = 0; - for ML = 1:NL - % J=IS(ML); - J = ML; % hd=-1e7;hm=-9899; - for ND = 1:2 - MN = ML + ND - 1; % SAVEKfL_h(ML,ND)=KfL_h(ML,ND); - if SWCC == 1 - if SFCC == 1 - if abs(hh(MN)) >= abs(hd) - Gama_hh(MN) = 0; - elseif abs(hh(MN)) >= abs(hm) - Gama_hh(MN) = log(abs(hd) / abs(hh(MN))) / log(abs(hd) / abs(hm)); - else - Gama_hh(MN) = 1; - end - Theta_m(ML) = Gama_hh(MN) * Theta_r(J) + (Theta_s(J) - Gama_hh(MN) * Theta_r(J)) * (1 + abs(Alpha(J) * (-1))^n(J))^m(J); % Theta_UU==>Theta_LL Theta_LL==>Theta_UU - if Theta_m(ML) >= POR(J) - Theta_m(ML) = POR(J); - elseif Theta_m(ML) <= Theta_s(J) - Theta_m(ML) = Theta_s(J); - end - - if hh(MN) >= -1 - Theta_UU(ML, ND) = Theta_s(J); - DTheta_LLh(ML, ND) = 0; % hh_frez(MN)=h_frez(MN); - if (hh_frez(MN)) >= 0 - Theta_LL(ML, ND) = Theta_s(J); - DTheta_UUh(ML, ND) = 0; % Se(ML,ND)=Theta_LL(ML,ND)/POR(J); - Se(ML, ND) = 1; - else - if Thmrlefc - if (hh(MN) + hh_frez(MN)) <= hd - Theta_LL(ML, ND) = 0; - DTheta_UUh(ML, ND) = 0; % (Theta_s(J)-Gama_hh(MN)*Theta_a(J))*Alpha(J)*n(J)*abs(Alpha(J)*hh(MN))^(n(J)-1)*(-m(J))*(1+abs(Alpha(J)*hh(MN))^n(J))^(-m(J)-1); - Se(ML, ND) = 0; - else - Theta_LL(ML, ND) = Gama_hh(MN) * Theta_r(J) + (Theta_m(ML) - Gama_hh(MN) * Theta_r(J)) / (1 + abs(Alpha(J) * (hh(MN) + hh_frez(MN)))^n(J))^m(J); % Theta_UU==>Theta_LL Theta_LL==>Theta_UU - DTheta_UUh(ML, ND) = (-Theta_r(J)) / (abs(hh(MN) + hh_frez(MN)) * log(abs(hd / hm))) * (1 - (1 + abs(Alpha(J) * (hh(MN) + hh_frez(MN)))^n(J))^(-m(J))) - Alpha(J) * n(J) * m(J) * (Theta_m(ML) - Gama_hh(MN) * Theta_r(J)) * ((1 + abs(Alpha(J) * (hh(MN) + hh_frez(MN)))^n(J))^(-m(J) - 1)) * (abs(Alpha(J) * (hh(MN) + hh_frez(MN)))^(n(J) - 1)); % (Theta_s(J)-Gama_hh(MN)*Theta_a(J))*Alpha(J)*n(J)*abs(Alpha(J)*hh(MN))^(n(J)-1)*(-m(J))*(1+abs(Alpha(J)*hh(MN))^n(J))^(-m(J)-1); - Se(ML, ND) = Theta_LL(ML, ND) / POR(J); - end - - else - if abs(hh(MN) - h(MN)) < 1e-3 - DTheta_UUh(ML, ND) = (Theta_m(ML) - Theta_r(J)) * Alpha(J) * n(J) * abs(Alpha(J) * (hh(MN) + hh_frez(MN)))^(n(J) - 1) * (-m(J)) * (1 + abs(Alpha(J) * (hh(MN) + hh_frez(MN)))^n(J))^(-m(J) - 1); - else - DTheta_UUh(ML, ND) = (Theta_LL(ML, ND) - Theta_L(ML, ND)) / (hh(MN) + hh_frez(MN) - h(MN) - h_frez(MN)); - end - end - end - else - if Thmrlefc - if Gama_hh(MN) == 0 - Theta_UU(ML, ND) = 0; - DTheta_LLh(ML, ND) = 0; % (Theta_s(J)-Gama_hh(MN)*Theta_a(J))*Alpha(J)*n(J)*abs(Alpha(J)*hh(MN))^(n(J)-1)*(-m(J))*(1+abs(Alpha(J)*hh(MN))^n(J))^(-m(J)-1); - if (hh(MN) + hh_frez(MN)) <= hd - Theta_LL(ML, ND) = 0; - DTheta_UUh(ML, ND) = 0; % (Theta_s(J)-Gama_hh(MN)*Theta_a(J))*Alpha(J)*n(J)*abs(Alpha(J)*hh(MN))^(n(J)-1)*(-m(J))*(1+abs(Alpha(J)*hh(MN))^n(J))^(-m(J)-1); - Se(ML, ND) = 0; - else - Theta_LL(ML, ND) = Gama_hh(MN) * Theta_r(J) + (Theta_m(ML) - Gama_hh(MN) * Theta_r(J)) / (1 + abs(Alpha(J) * (hh(MN) + hh_frez(MN)))^n(J))^m(J); % Theta_UU==>Theta_LL Theta_LL==>Theta_UU - DTheta_UUh(ML, ND) = (-Theta_r(J)) / (abs(hh(MN) + hh_frez(MN)) * log(abs(hd / hm))) * (1 - (1 + abs(Alpha(J) * (hh(MN) + hh_frez(MN)))^n(J))^(-m(J))) - Alpha(J) * n(J) * m(J) * (Theta_m(ML) - Gama_hh(MN) * Theta_r(J)) * ((1 + abs(Alpha(J) * (hh(MN) + hh_frez(MN)))^n(J))^(-m(J) - 1)) * (abs(Alpha(J) * (hh(MN) + hh_frez(MN)))^(n(J) - 1)); % (Theta_s(J)-Gama_hh(MN)*Theta_a(J))*Alpha(J)*n(J)*abs(Alpha(J)*hh(MN))^(n(J)-1)*(-m(J))*(1+abs(Alpha(J)*hh(MN))^n(J))^(-m(J)-1); - Se(ML, ND) = Theta_LL(ML, ND) / POR(J); - end - else - Theta_UU(ML, ND) = Gama_hh(MN) * Theta_r(J) + (Theta_m(ML) - Gama_hh(MN) * Theta_r(J)) / (1 + abs(Alpha(J) * hh(MN))^n(J))^m(J); - DTheta_LLh(ML, ND) = (-Theta_r(J)) / (abs(hh(MN)) * log(abs(hd / hm))) * (1 - (1 + abs(Alpha(J) * hh(MN))^n(J))^(-m(J))) - Alpha(J) * n(J) * m(J) * (Theta_m(ML) - Gama_hh(MN) * Theta_r(J)) * ((1 + abs(Alpha(J) * hh(MN))^n(J))^(-m(J) - 1)) * (abs(Alpha(J) * hh(MN))^(n(J) - 1)); % (Theta_s(J)-Gama_hh(MN)*Theta_a(J))*Alpha(J)*n(J)*abs(Alpha(J)*hh(MN))^(n(J)-1)*(-m(J))*(1+abs(Alpha(J)*hh(MN))^n(J))^(-m(J)-1); - if (hh(MN) + hh_frez(MN)) <= hd - Theta_LL(ML, ND) = 0; - DTheta_UUh(ML, ND) = 0; % (Theta_s(J)-Gama_hh(MN)*Theta_a(J))*Alpha(J)*n(J)*abs(Alpha(J)*hh(MN))^(n(J)-1)*(-m(J))*(1+abs(Alpha(J)*hh(MN))^n(J))^(-m(J)-1); - Se(ML, ND) = 0; - else - Theta_LL(ML, ND) = Gama_hh(MN) * Theta_r(J) + (Theta_m(ML) - Gama_hh(MN) * Theta_r(J)) / (1 + abs(Alpha(J) * (hh(MN) + hh_frez(MN)))^n(J))^m(J); % Theta_UU==>Theta_LL Theta_LL==>Theta_UU - DTheta_UUh(ML, ND) = (-Theta_r(J)) / (abs(hh(MN) + hh_frez(MN)) * log(abs(hd / hm))) * (1 - (1 + abs(Alpha(J) * (hh(MN) + hh_frez(MN)))^n(J))^(-m(J))) - Alpha(J) * n(J) * m(J) * (Theta_m(ML) - Gama_hh(MN) * Theta_r(J)) * ((1 + abs(Alpha(J) * (hh(MN) + hh_frez(MN)))^n(J))^(-m(J) - 1)) * (abs(Alpha(J) * (hh(MN) + hh_frez(MN)))^(n(J) - 1)); % (Theta_s(J)-Gama_hh(MN)*Theta_a(J))*Alpha(J)*n(J)*abs(Alpha(J)*hh(MN))^(n(J)-1)*(-m(J))*(1+abs(Alpha(J)*hh(MN))^n(J))^(-m(J)-1); - Se(ML, ND) = Theta_LL(ML, ND) / POR(J); - end - end - else - if abs(hh(MN) - h(MN)) < 1e-3 - DTheta_LLh(ML, ND) = (Theta_m(ML) - Theta_r(J)) * Alpha(J) * n(J) * abs(Alpha(J) * (hh(MN)))^(n(J) - 1) * (-m(J)) * (1 + abs(Alpha(J) * (hh(MN)))^n(J))^(-m(J) - 1); - else - DTheta_LLh(ML, ND) = (Theta_UU(ML, ND) - Theta_U(ML, ND)) / (hh(MN) - h(MN)); - end - end - end - %%%%%%%%%%%%%%%%%%%% Sin function for ice calculation %%%%%%%%%%%%%%%%%%%%% - else - Tf1 = 273.15 + 1; - Tf2 = 273.15 - 3; % XCAP(J)=0.23; - if hh(MN) >= -1e-6 - Theta_UU(ML, ND) = Theta_s(J); - hh(MN) = -1e-6; - DTheta_UUh(ML, ND) = 0; - if TT(MN) + 273.15 > Tf1 - Theta_II(ML, ND) = 0; - Theta_LL(ML, ND) = Theta_s(J); - - elseif TT(MN) + 273.15 >= Tf2 - Theta_II(ML, ND) = 0.5 * (1 - sin(pi() * (TT(MN) + 273.15 - 0.5 * Tf1 - 0.5 * Tf2) / (Tf1 - Tf2))) * XCAP(J); - Theta_LL(ML, ND) = Theta_UU(ML, ND) - Theta_II(ML, ND) * RHOI / RHOL; - - else - Theta_II(ML, ND) = XCAP(J); - Theta_LL(ML, ND) = Theta_UU(ML, ND) - Theta_II(ML, ND) * RHOI / RHOL; - - end - if Theta_LL(ML, ND) <= 0.06 - Theta_LL(ML, ND) = 0.06; - DTheta_LLh(ML, ND) = 0; - Se(ML, ND) = 0; - else - Theta_LL(ML, ND) = Theta_LL(ML, ND); - Se(ML, ND) = Theta_LL(ML, ND) / POR(J); - DTheta_LLh(ML, ND) = DTheta_UUh(ML, ND); - end - elseif hh(MN) <= -1e7 - Theta_UU(ML, ND) = Theta_r(J); - hh(MN) = -1e7; - DTheta_UUh(ML, ND) = 0; - Theta_II(ML, ND) = 0; - Theta_LL(ML, ND) = Theta_r(J); - Se(ML, ND) = 0; - DTheta_LLh(ML, ND) = 0; - else - Theta_UU(ML, ND) = Theta_r(J) + (Theta_s(J) - Theta_r(J)) / (1 + abs(Alpha(J) * hh(MN))^n(J))^m(J); - - if Thmrlefc - DTheta_UUh(ML, ND) = (Theta_s(J) - Theta_r(J)) * Alpha(J) * n(J) * abs(Alpha(J) * hh(MN))^(n(J) - 1) * (-m(J)) * (1 + abs(Alpha(J) * hh(MN))^n(J))^(-m(J) - 1); - else - if abs(hh(MN) - h(MN)) < 1e-3 - DTheta_UUh(ML, ND) = (Theta_s(J) - Theta_r(J)) * Alpha(J) * n(J) * abs(Alpha(J) * hh(MN))^(n(J) - 1) * (-m(J)) * (1 + abs(Alpha(J) * hh(MN))^n(J))^(-m(J) - 1); - else - DTheta_UUh(ML, ND) = (Theta_UU(ML, ND) - Theta_U(ML, ND)) / (hh(MN) - h(MN)); - end - end - if TT(MN) + 273.15 > Tf1 - Theta_II(ML, ND) = 0; - Theta_LL(ML, ND) = Theta_r(J) + (Theta_s(J) - Theta_r(J)) / (1 + abs(Alpha(J) * (hh(MN) + hh_frez(MN)))^n(J))^m(J); % Theta_UU==>Theta_LL Theta_LL==>Theta_UU - - elseif TT(MN) + 273.15 >= Tf2 - Theta_II(ML, ND) = 0.5 * (1 - sin(pi() * (TT(MN) + 273.15 - 0.5 * Tf1 - 0.5 * Tf2) / (Tf1 - Tf2))) * XCAP(J); - Theta_LL(ML, ND) = Theta_UU(ML, ND) - Theta_II(ML, ND) * RHOI / RHOL; % Theta_UU(ML,ND) - - else - Theta_II(ML, ND) = XCAP(J); - Theta_LL(ML, ND) = Theta_UU(ML, ND) - Theta_II(ML, ND) * RHOI / RHOL; % Theta_UU(ML,ND) - - end - if Theta_LL(ML, ND) <= 0.06 - Theta_LL(ML, ND) = 0.06; - DTheta_LLh(ML, ND) = 0; - Se(ML, ND) = 0; - else - Theta_LL(ML, ND) = Theta_LL(ML, ND); - Se(ML, ND) = Theta_LL(ML, ND) / POR(J); - DTheta_LLh(ML, ND) = DTheta_UUh(ML, ND); - end - - end - end - else - if hh(MN) >= Phi_s(J) - Theta_UU(ML, ND) = Theta_s(J); - hh(MN) = Phi_s(J); - DTheta_UUh(ML, ND) = 0; - if hh(MN) + hh_frez(MN) <= -1e7 - Theta_LL(ML, ND) = Theta_r(J); - DTheta_LLh(ML, ND) = 0; - Se(ML, ND) = 0; - elseif hh(MN) + hh_frez(MN) >= Phi_s(J) - Theta_LL(ML, ND) = Theta_s(J); - DTheta_LLh(ML, ND) = 0; - Se(ML, ND) = 1; - else - Theta_LL(ML, ND) = Theta_s(J) * ((hh(MN) + hh_frez(MN)) / Phi_s(J))^(-1 * Lamda(J)); % Theta_r(J)+(Theta_s(J)-Theta_r(J))/(1+abs(Alpha(J)*(hh(MN)+hh_frez(MN)))^n(J))^m(J); %Theta_UU==>Theta_LL Theta_LL==>Theta_UU - if Thmrlefc - DTheta_LLh(ML, ND) = Theta_s(J) / Phi_s(J) * ((hh(MN) + hh_frez(MN)) / Phi_s(J))^(-1 * Lamda(J) - 1); % (Theta_s(J)-Theta_r(J))*Alpha(J)*n(J)*abs(Alpha(J)*(hh(MN)+hh_frez(MN)))^(n(J)-1)*(-m(J))*(1+abs(Alpha(J)*(hh(MN)+hh_frez(MN)))^n(J))^(-m(J)-1); - else - if abs(hh(MN) - h(MN)) < 1e-3 - DTheta_LLh(ML, ND) = (Theta_s(J) - Theta_r(J)) * Alpha(J) * n(J) * abs(Alpha(J) * (hh(MN) + hh_frez(MN)))^(n(J) - 1) * (-m(J)) * (1 + abs(Alpha(J) * (hh(MN) + hh_frez(MN)))^n(J))^(-m(J) - 1); - else - DTheta_LLh(ML, ND) = (Theta_LL(ML, ND) - Theta_L(ML, ND)) / (hh(MN) + hh_frez(MN) - h(MN) - h_frez(MN)); - end - end - Se(ML, ND) = Theta_LL(ML, ND) / POR(J); - end - elseif hh(MN) <= -1e7 - Theta_LL(ML, ND) = Theta_r(J); - Theta_UU(ML, ND) = Theta_r(J); - Theta_II(ML, ND) = 0; - hh(MN) = -1e7; - hh_frez(MN) = -1e-6; - DTheta_UUh(ML, ND) = 0; - DTheta_LLh(ML, ND) = 0; - Se(ML, ND) = 0; - else - Theta_UU(ML, ND) = Theta_s(J) * ((hh(MN)) / Phi_s(J))^(-1 * Lamda(J)); % Theta_r(J)+(Theta_s(J)-Theta_r(J))/(1+abs(Alpha(J)*(hh(MN)+hh_frez(MN)))^n(J))^m(J); %Theta_UU==>Theta_LL Theta_LL==>Theta_UU - if Thmrlefc - DTheta_UUh(ML, ND) = Theta_s(J) / Phi_s(J) * ((hh(MN)) / Phi_s(J))^(-1 * Lamda(J) - 1); % (Theta_s(J)-Theta_r(J))*Alpha(J)*n(J)*abs(Alpha(J)*(hh(MN)+hh_frez(MN)))^(n(J)-1)*(-m(J))*(1+abs(Alpha(J)*(hh(MN)+hh_frez(MN)))^n(J))^(-m(J)-1); - else - if abs(hh(MN) - h(MN)) < 1e-3 - DTheta_UUh(ML, ND) = (Theta_s(J) - Theta_r(J)) * Alpha(J) * n(J) * abs(Alpha(J) * (hh(MN)))^(n(J) - 1) * (-m(J)) * (1 + abs(Alpha(J) * (hh(MN)))^n(J))^(-m(J) - 1); - else - DTheta_UUh(ML, ND) = (Theta_UU(ML, ND) - Theta_U(ML, ND)) / (hh(MN) - h(MN)); - end - end - - if hh(MN) + hh_frez(MN) <= -1e7 - Theta_LL(ML, ND) = Theta_r(J); - DTheta_LLh(ML, ND) = 0; - Se(ML, ND) = 0; - elseif hh(MN) + hh_frez(MN) >= Phi_s(J) - Theta_LL(ML, ND) = Theta_s(J); - DTheta_LLh(ML, ND) = 0; - Se(ML, ND) = 1; - else - Theta_LL(ML, ND) = Theta_s(J) * ((hh(MN) + hh_frez(MN)) / Phi_s(J))^(-1 * Lamda(J)); % Theta_r(J)+(Theta_s(J)-Theta_r(J))/(1+abs(Alpha(J)*(hh(MN)+hh_frez(MN)))^n(J))^m(J); %Theta_UU==>Theta_LL Theta_LL==>Theta_UU - if Thmrlefc - DTheta_LLh(ML, ND) = Theta_s(J) / Phi_s(J) * ((hh(MN) + hh_frez(MN)) / Phi_s(J))^(-1 * Lamda(J) - 1); % (Theta_s(J)-Theta_r(J))*Alpha(J)*n(J)*abs(Alpha(J)*(hh(MN)+hh_frez(MN)))^(n(J)-1)*(-m(J))*(1+abs(Alpha(J)*(hh(MN)+hh_frez(MN)))^n(J))^(-m(J)-1); - else - if abs(hh(MN) - h(MN)) < 1e-3 - DTheta_LLh(ML, ND) = (Theta_s(J) - Theta_r(J)) * Alpha(J) * n(J) * abs(Alpha(J) * (hh(MN) + hh_frez(MN)))^(n(J) - 1) * (-m(J)) * (1 + abs(Alpha(J) * (hh(MN) + hh_frez(MN)))^n(J))^(-m(J) - 1); - else - DTheta_LLh(ML, ND) = (Theta_LL(ML, ND) - Theta_L(ML, ND)) / (hh(MN) + hh_frez(MN) - h(MN) - h_frez(MN)); - end - end - Se(ML, ND) = Theta_LL(ML, ND) / POR(J); - - end - end - end - if Se(ML, ND) >= 1 - Se(ML, ND) = 1; - elseif Se(ML, ND) <= 0 - Se(ML, ND) = 0; - end - if isnan(Se(ML, ND)) == 1 - warning('\n case "isnan(Se(ML, ND)) == 1" happens. Dont know what to do! \r'); - end - Theta_II(ML, ND) = (Theta_UU(ML, ND) - Theta_LL(ML, ND)) * RHOL / RHOI; % ice water content - if Theta_UU(ML, ND) ~= 0 - Ratio_ice(ML, ND) = RHOI * Theta_II(ML, ND) / (RHOL * Theta_UU(ML, ND)); % ice ratio - else - Ratio_ice(ML, ND) = 0; - end - if KIT - MU_WN = Constants.MU_W0 * exp(Constants.MU1 / (8.31441 * (20 + 133.3))); - if TT(MN) < -20 - MU_W(ML, ND) = 3.71e-2; % CKT(MN)=0.2688; - elseif TT(MN) > 150 - MU_W(ML, ND) = 1.81e-3; - % CKT(MN)=5.5151; % kgm^-1s^-1 --> 10 g.cm^-1.s^-1; J.cm^-2---> kg.m^2.s^-2.cm^-2--> 1e7g.cm^2.s^-2.cm^-2 - else - MU_W(ML, ND) = Constants.MU_W0 * exp(Constants.MU1 / (8.31441 * (TT(MN) + 133.3))); - - end - % CKT(MN)=CKTN/(50+2.575*TT(MN)); - CKT(MN) = MU_WN / MU_W(ML, ND); - if Se(ML, ND) == 0 - KL_h(ML, ND) = 0; - else - KL_h(ML, ND) = CKT(MN) * Ks(J) * (Se(ML, ND)^(0.5)) * (1 - (1 - Se(ML, ND)^(1 / m(J)))^m(J))^2; - end - - CORF = 1; - FILM = 1; % indicator for film flow parameterization; =1, Zhang (2010); =2, Lebeau and Konrad (2010) - if FILM == 1 - % %%%%%%%%% see Zhang (2010) - AGR(J) = 0.00035; - - % % %%%%%%%%% see Zeng (2011) and Zhang (2010) - RHOW0 = 1e3; - GVA = 9.81; - uw0 = 2.4152e-5; % Pa s - u1 = 4.7428; % kJ mol-1 - R = 8.314472; % J mol-1 K-1 - e = 78.54; - e0 = 8.85e-12; - B_sig = 235.8E-3; - Tc = 647.15; - e_sig = 1.256; - c_sig = -0.625; - kb = 1.381e-23; % Boltzmann constant - Bi = 1; - Ba = 1.602e-19 * Bi; - uw = uw0 * exp(u1 / R / (TT(MN) + 133.3)); - sigma = B_sig * ((Tc - (TT(MN) + 273.15)) / Tc)^e_sig * (1 + c_sig * (Tc - (TT(MN) + 273.15)) / Tc); - B(J) = 2^0.5 * pi()^2 * RHOW0 * GVA / uw * (e * e0 / 2 / sigma)^1.5 * (kb * (TT(MN) + 273.15) / (Bi * Ba))^3; - - Coef_f = 0.0373 * (2 * AGR(J))^3.109; - Ks_flm(ML, ND) = B(J) * (1 - POR(J)) * (2 * AGR(J))^0.5; % m2 - if hh(MN) < -1 - Kr(ML, ND) = Coef_f * (1 + 2 * AGR(J) * RHOW0 * GVA * abs(hh(MN) / 100) / 2 / sigma)^(-1.5); - else - Kr(ML, ND) = 1; - end - KL_h_flm(ML, ND) = Ks_flm(ML, ND) * Kr(ML, ND) * 1e4; % m2 --> cm2 - else - %%%%%%%%% see sutraset --> perfilm.f based on Lebeau and Konrad (2010) - % EFFECTIVE DIAMETER - ASVL = -6e-20; - GVA = 9.8; - PERMVAC = 8.854e-12; % (VAC)CUM (PERM)EABILITY [PERMFS] [C2J-1M-1 OR F/M OR S4A2/KG/M3] - ELECTRC = 1.602e-19; - BOTZC = 1.381e-23; - RHOW0 = 1000; - PSICM = 1000; - SWM = 0.1; % THE SATURATION WHEN SURFACE ROUGHNESS OF THE SOLID GRAIN BECOMES NEGNIGIBLE (-) IT IS EQUIVALENT TO - % SATURATION BEEN 1000M FROM TULLER AND OR (2005) - RELPW = 78.54; - % SWG=1; - ED = 6.0 * (1 - POR(J)) * (-ASVL / (6.0 * pi() * RHOW0 * GVA * PSICM))^(1.0 / 3.0) / POR(J) / SWM; - % FILM THICKNESS (WARNING) THIS IS NOT THE FULL PART - DEL = (RELPW * PERMVAC / (2 * RHOW0 * GVA))^0.50 * (pi() * BOTZC * 298.15 / ELECTRC); - Ks_flm(ML, ND) = CORF * 4.0 * (1 - POR(J)) * DEL^3.0 / pi() / ED; - - if hh(MN) <= -1 - Kr(ML, ND) = (1 - Se(ML, ND)) * abs(hh(MN) / 100)^(-1.50); - else - Kr(ML, ND) = 1; - end - KL_h_flm(ML, ND) = Ks_flm(ML, ND) * Kr(ML, ND) * 1e4; % m2 --> cm2 - end - if KL_h_flm(ML, ND) <= 0 - KL_h_flm(ML, ND) = 0; - elseif KL_h_flm(ML, ND) >= 1e-6 - KL_h_flm(ML, ND) = 1e-6; - end - if KL_h(ML, ND) <= 1E-20 - KL_h(ML, ND) = 1E-20; - end - if Gama_hh(MN) ~= 1 - KfL_h(ML, ND) = KL_h(ML, ND) * 10^(-1 * Imped(MN) * Ratio_ice(ML, ND)); % +KL_h_flm(ML,ND); % hydraulic conductivity for freezing soil - % KL_h(ML,ND)=KL_h(ML,ND)+KL_h_flm(ML,ND); - else - KfL_h(ML, ND) = KL_h(ML, ND) * 10^(-1 * Imped(MN) * Ratio_ice(ML, ND)); % hydraulic conductivity for freezing soil - end - if isnan(KL_h(ML, ND)) == 1 - KL_h(ML, ND) = 0; - warning('\n case "isnan(KL_h(ML, ND)) == 1", set "KL_h(ML, ND) = 0" \r'); - end - if ~isreal(KL_h(ML, ND)) - warning('\n case "~isreal(KL_h(ML, ND))", dont know what to do! \r'); - end - KfL_T(ML, ND) = helpers.heaviside1(TT_CRIT(MN) - (TT(MN) + T0)) * L_f * 1e4 / (g * (T0)); % thermal consider for freezing soil - else - KL_h(ML, ND) = 0; - KfL_h(ML, ND) = 0; - KfL_T(ML, ND) = 0; - end - end - end - - %%%%%%%%% Unit of KL_h is determined by Ks, which would be given at the%%%% - %%%%%%%%% beginning.Import thing is to keep the unit of matric head hh(MN) - %%%%%%%%% as 'cm'.%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% diff --git a/src/CondT_coeff.m b/src/CondT_coeff.m deleted file mode 100644 index a3faa50b..00000000 --- a/src/CondT_coeff.m +++ /dev/null @@ -1,44 +0,0 @@ -function [c_unsat, Lambda_eff, ZETA, ETCON, EHCAP, TETCON, EfTCON] = CondT_coeff(Theta_LL, Lambda1, Lambda2, Lambda3, RHO_bulk, Theta_g, RHODA, RHOV, c_a, c_V, c_L, NL, nD, ThmrlCondCap, ThermCond, HCAP, SF, TCA, GA1, GA2, GB1, GB2, HCD, ZETA0, CON0, PS1, PS2, XWILT, XK, TT, POR, DRHOVT, L, D_A, Theta_V, Theta_II, TCON_dry, Theta_s, XSOC, TPS1, TPS2, TCON0, TCON_s, FEHCAP, RHOI, RHOL, c_unsat, Lambda_eff, ETCON, EHCAP, TETCON, EfTCON, ZETA) - - if ThmrlCondCap == 1 - % run EfeCapCond; - [ETCON, EHCAP, TETCON, EfTCON, ZETA] = EfeCapCond(HCAP, SF, TCA, GA1, GA2, GB1, GB2, HCD, ZETA0, CON0, PS1, PS2, XWILT, XK, TT, NL, POR, Theta_LL, DRHOVT, L, D_A, RHOV, Theta_V, Theta_II, TCON_dry, Theta_s, XSOC, ThermCond, TPS1, TPS2, TCON0, TCON_s, FEHCAP, RHOI, RHOL, ETCON, EHCAP, TETCON, EfTCON, ZETA); - for ML = 1:NL - for ND = 1:nD - - if ThermCond == 1 - Lambda_eff(ML, ND) = ETCON(ML, ND); - elseif ThermCond == 2 - Lambda_eff(ML, ND) = EfTCON(ML, ND) / 100; - elseif ThermCond == 3 - Lambda_eff(ML, ND) = TETCON(ML, ND); - elseif ThermCond == 4 - Lambda_eff(ML, ND) = EfTCON(ML, ND) / 100; - end - - if Lambda_eff(ML, ND) <= 0 - Lambda_eff(ML, ND) = 0.0008; - elseif Lambda_eff(ML, ND) >= 0.02 - Lambda_eff(ML, ND) = 0.02; - else - Lambda_eff(ML, ND) = Lambda_eff(ML, ND); - end - c_unsat(ML, ND) = EHCAP(ML, ND); - end - end - else - MN = 0; - for ML = 1:NL - for ND = 1:nD - MN = ML + ND - 1; - Lambda_eff(ML, ND) = Lambda1 + Lambda2 * Theta_LL(ML, ND) + Lambda3 * Theta_LL(ML, ND)^0.5; % 3.6*0.001*4.182; % It is possible to add the dispersion effect here to consider the heat dispersion. - - c_unsat(ML, ND) = 837 * RHO_bulk / 1000 + Theta_LL(ML, ND) * c_L + Theta_g(ML, ND) * (RHODA(MN) * c_a + RHOV(MN) * c_V); % 9.79*0.1*4.182;% - end - end - end - - %%%%% Unit of Lambda_eff is (J.m^-1.s^-1.Cels^-1), While c_unsat is (J.m^-3.Cels^-1) - %%%%% UnitC needs to be used here to convert 'm' to 'cm' . 837 in J.kg^-1.Cels^-1. RHO_bulk in kg.m^-3 %%%%% - %%%%% c_a, c_v,would be in J.g^-1.Cels^-1 as showed in - %%%%% globalization. RHOV and RHODA would be in g.cm^-3 diff --git a/src/CondV_DE.m b/src/CondV_DE.m deleted file mode 100644 index 9f83f6bf..00000000 --- a/src/CondV_DE.m +++ /dev/null @@ -1,34 +0,0 @@ -function [D_V, Eta, D_A] = CondV_DE(Theta_LL, TT, fc, Theta_s, NL, nD, Theta_g, POR, ThmrlCondCap, ZETA, XK, DVT_Switch, Theta_UU) - MN = 0; - for ML = 1:NL - J = ML; - for ND = 1:nD - MN = ML + ND - 1; - - if ThmrlCondCap - if Theta_UU(ML, ND) < XK(J) - EnhnLiqIsland(ML, ND) = POR(J); % -Theta_II(ML,ND)*RHOI/RHOL - else - EnhnLiqIsland(ML, ND) = Theta_g(ML, ND) * (1 + Theta_UU(ML, ND) / (POR(J) - XK(J))); % -Theta_II(ML,ND)*RHOI/RHOL - end - - f0(ML, ND) = Theta_g(ML, ND)^(7 / 3) / Theta_s(J)^2; % Theta_g(ML,ND)^0.67; - - D_A(MN) = 0.229 * (1 + TT(MN) / 273)^1.75; % cm2/s---------5.8*10^(-7)*(273.15+TT(MN))^2.3/(UnitC^2);% - - if DVT_Switch == 1 - Eta(ML, ND) = ZETA(ML, ND) * EnhnLiqIsland(ML, ND) / (f0(ML, ND) * Theta_g(ML, ND)); % 0; % - else - Eta(ML, ND) = 0; - end - else - f0(ML, ND) = Theta_g(ML, ND)^0.67; % - D_A(MN) = 1e4 * 2.12 * 10^(-5) * (1 + TT(MN) / 273.15)^2; % 2.12e-5*(1+T/273.15)^2 m2/s---> 1e4*cm2/s - Eta(ML, ND) = 8 + 3 * Theta_LL(ML, ND) / Theta_s(J) - 7 * exp(-((1 + 2.6 / sqrt(fc)) * Theta_LL(ML, ND) / Theta_s(J))^3); - end - - D_V(ML, ND) = f0(ML, ND) * Theta_g(ML, ND) * D_A(MN); - - end - end - %%%%%%%%%%%%% With UnitC^2, m^2.s^-1 would be converted as cm^2.s^-1 %%%%%%%%%%%%% diff --git a/src/CondV_DVg.m b/src/CondV_DVg.m deleted file mode 100644 index b08eef46..00000000 --- a/src/CondV_DVg.m +++ /dev/null @@ -1,33 +0,0 @@ -function [D_Vg, V_A, Beta_g, DPgDZ, Beta_gBAR, Alpha_LgBAR] = CondV_DVg(P_gg, Theta_g, Sa, V_A, k_g, MU_a, DeltZ, Alpha_Lg, KaT_Switch, Theta_s, Se, NL, DPgDZ, Beta_gBAR, Alpha_LgBAR, Beta_g) - - MN = 0; - for ML = 1:NL - J = ML; - for ND = 1:2 - MN = ML + ND - 1; - Sa(ML, ND) = 1 - Se(ML, ND); - f0(ML, ND) = Theta_g(ML, ND)^(7 / 3) / Theta_s(J)^2; % Theta_g(ML,ND)^0.67; - Beta_g(ML, ND) = (k_g(ML, ND) / MU_a); - Alpha_Lg(ML, ND) = 0.078 * (13.6 - 16 * Sa(ML, ND) + 3.4 * Sa(ML, ND)^5) * 100; % - end - end - - for ML = 1:NL - Beta_gBAR(ML) = (Beta_g(ML, 1) + Beta_g(ML, 2)) / 2; - DPgDZ(ML) = (P_gg(ML + 1) - P_gg(ML)) / DeltZ(ML); - Alpha_LgBAR(ML) = (Alpha_Lg(ML, 1) + Alpha_Lg(ML, 2)) / 2; - end - - for ML = 1:NL - V_A(ML) = -Beta_gBAR(ML) * DPgDZ(ML); % 0; % - if KaT_Switch == 1 - D_Vg(ML) = Alpha_LgBAR(ML) * abs(V_A(ML)); % 0; %0; % - else - D_Vg(ML) = 0; - end - end - - %%%%%%%%%%%% Unit of kg is cm^2, MU_a is g.cm^-1.s^-1,V_A is cm.s^-1, D_Vg - %%%%%%%%%%%% is cm^2.s^-1, The unit of soil air pressure should be Pa=kg.m^-1.s^-2=10g.cm^-1.s^-2 %%%%%%%%%%%%% - %%%%%%%%%%%% Notice that '10'in V_A is because Pa needs to be converted as10g.cm^-1.s^-2; has been done in the StarInit subroutine %%%%%%%%%%%%% - %%%%%%%%%%%% MU_a's unit has been changed %%%%%%%%% diff --git a/src/Condg_k_g.m b/src/Condg_k_g.m deleted file mode 100644 index 645be131..00000000 --- a/src/Condg_k_g.m +++ /dev/null @@ -1,21 +0,0 @@ -function [k_g] = Condg_k_g(POR, NL, m, Theta_g, g, MU_W, Ks, RHOL, SWCC, Imped, Ratio_ice, Soilairefc, MN) - - for ML = 1:NL - J = ML; - for ND = 1:2 - Sa(ML, ND) = Theta_g(ML, ND) / POR(J); - if SWCC == 1 - k_g(ML, ND) = Ks(J) * MU_W(ML, ND) * (1 - Sa(ML, ND)^0.5) * (1 - (1 - (1 - Sa(ML, ND)^(1 / m(J))))^m(J))^2 / (g * RHOL); - else - k_g(ML, ND) = 0; - end - if Soilairefc == 1 - k_g(ML, ND) = (Ks(J) * MU_W(ML, ND) * (1 - Sa(ML, ND)^0.5) * (1 - (1 - (1 - Sa(ML, ND)^(1 / m(J))))^m(J))^2 / (g * RHOL)) * 10^(-1 * Imped(MN) * Ratio_ice(ML, ND)); - else - k_g(ML, ND) = 0; - end - end - end - - %%% Unit of k_g is m^2 , with UnitC^2, unit has been converted as cm^2 % - %%% 10^(-12)is used to convert the micrometer to meter % diff --git a/src/EfeCapCond.m b/src/EfeCapCond.m deleted file mode 100644 index fa4b06df..00000000 --- a/src/EfeCapCond.m +++ /dev/null @@ -1,107 +0,0 @@ -function [ETCON, EHCAP, TETCON, EfTCON, ZETA] = EfeCapCond(HCAP, SF, TCA, GA1, GA2, GB1, GB2, HCD, ZETA0, CON0, PS1, PS2, XWILT, XK, TT, NL, POR, Theta_LL, DRHOVT, L, D_A, RHOV, Theta_V, Theta_II, TCON_dry, Theta_s, XSOC, ThermCond, TPS1, TPS2, TCON0, TCON_s, FEHCAP, RHOI, RHOL, ETCON, EHCAP, TETCON, EfTCON, ZETA) - global TCON - MN = 0; - TCON(1) = 1.37e-3 * 4.182; - for ML = 1:NL - for ND = 1:2 - MN = ML + ND - 1; - J = ML; % J=IS(ML); - XXX = Theta_LL(ML, ND); - XII = Theta_II(ML, ND) * RHOI / RHOL; % MODIFIED 201905 - if Theta_LL(ML, ND) < XK(J) - XXX = XK(J); - end - - if XXX < XWILT(J) - SF(2) = GA2 + GB2(J) * XXX; - else - SF(2) = GA1 + GB1(J) * XXX; - end - D_A(MN) = 0.229 * (1 + TT(MN) / 273)^1.75; - TCON(2) = TCA + D_A(MN) * L(MN) * DRHOVT(MN); % TCA+D_A(MN)*L(MN)*DRHOVT(MN); % Revised from ""(D_V(ML,ND)*Eta(ML,ND)+D_Vg(ML))*L(MN)*DRHOVT(MN) - TCON(6) = 5.2e-3 * 4.182; % ice thermal conductivity - SF(6) = 0.125; - TARG(2) = TCON(2) / TCON(1) - 1; - TARG(6) = TCON(6) / TCON(1) - 1; - GRAT(2) = 0.667 / (1 + TARG(2) * SF(2)) + 0.333 / (1 + TARG(2) * (1 - 2 * SF(2))); - GRAT(6) = 0.667 / (1 + TARG(6) * SF(6)) + 0.333 / (1 + TARG(6) * (1 - 2 * SF(6))); - ETCON(ML, ND) = (PS1(J) + XXX * TCON(1) + (POR(J) - XXX - XII) * GRAT(2) * TCON(2) + XII * GRAT(6) * TCON(6)) / (PS2(J) + XXX + (POR(J) - XXX - XII) * GRAT(2) + XII * GRAT(6)); - ZETA(ML, ND) = GRAT(2) / (GRAT(2) * (POR(J) - XXX) + XXX + PS2(J)); % ita_T enhancement factor - ZETA(ML, ND) = GRAT(2) / (GRAT(2) * (POR(J) - XXX - XII) + XXX + PS2(J) + XII * GRAT(6)); - if Theta_LL(ML, ND) == XXX - EHCAP(ML, ND) = HCD(J) + HCAP(1) * Theta_LL(ML, ND); - EHCAP(ML, ND) = EHCAP(ML, ND) + (0.448 * RHOV(MN) * 4.182 + HCAP(2)) * Theta_V(ML, ND) + HCAP(6) * XII; % The Calorie should be converted as J - else - ZETA(ML, ND) = ZETA0(J) + (ZETA(ML, ND) - ZETA0(J)) * Theta_LL(ML, ND) / XXX; - ETCON(ML, ND) = CON0(J) + (ETCON(ML, ND) - CON0(J)) * Theta_LL(ML, ND) / XXX; - EHCAP(ML, ND) = HCD(J) + HCAP(1) * Theta_LL(ML, ND); - EHCAP(ML, ND) = EHCAP(ML, ND) + (0.448 * RHOV(MN) * 4.182 + HCAP(2)) * Theta_V(ML, ND) + HCAP(6) * XII; % The Calorie should be converted as J - end - %%%%%% simpliefied Vries method by Tian 2016 - TSF(2) = 0.333 * (1 - (POR(J) - XXX - XII) / POR(J)); - TSF(6) = 0.333 * (1 - XII / POR(J)); - TSF(7) = 0.5; % shape factor for organic matter - D_A(MN) = 0.229 * (1 + TT(MN) / 273)^1.75; - TCON(2) = TCA + D_A(MN) * L(MN) * DRHOVT(MN); % TCA+D_A(MN)*L(MN)*DRHOVT(MN); % Revised from ""(D_V(ML,ND)*Eta(ML,ND)+D_Vg(ML))*L(MN)*DRHOVT(MN) - TCON(6) = 5.2e-3 * 4.182; % ice conductivity - TCON(7) = 0.25e-2; % soil organic matter conductivity - TARG(2) = TCON(2) / TCON(1) - 1; - TARG(6) = TCON(6) / TCON(1) - 1; - TARG(7) = TCON(7) / TCON(1) - 1; - TTGRAT(2) = 0.667 / (1 + TARG(2) * TSF(2)) + 0.333 / (1 + TARG(2) * (1 - 2 * TSF(2))); - TTGRAT(6) = 0.667 / (1 + TARG(6) * TSF(6)) + 0.333 / (1 + TARG(6) * (1 - 2 * TSF(6))); - TTGRAT(7) = 0.667 / (1 + TARG(7) * TSF(7)) + 0.333 / (1 + TARG(7) * (1 - 2 * TSF(7))); - if XSOC(J) == 0 - TETCON(ML, ND) = (TPS1(J) + XXX * TCON(1) + (POR(J) - XXX - XII) * TTGRAT(2) * TCON(2) + XII * TTGRAT(6) * TCON(6)) / (TPS2(J) + XXX + (POR(J) - XXX - XII) * TTGRAT(2) + XII * TTGRAT(6)); - else - TETCON(ML, ND) = (TPS1(J) + XXX * TCON(1) + (POR(J) - XXX - XII) * TTGRAT(2) * TCON(2) + XII * TTGRAT(6) * TCON(6) + XSOC(J) * TTGRAT(7) * TCON(7)) / (TPS2(J) + XXX + (POR(J) - XXX - XII) * TTGRAT(2) + XII * TTGRAT(6) + XSOC(J) * TTGRAT(7)); - end - TZETA(ML, ND) = TTGRAT(2) / (TTGRAT(2) * (POR(J) - XXX) + XXX + TPS2(J)); % ita_T enhancement factor - if Theta_LL(ML, ND) == XXX - EHCAP(ML, ND) = HCD(J) + HCAP(1) * Theta_LL(ML, ND); - EHCAP(ML, ND) = EHCAP(ML, ND) + (0.448 * RHOV(MN) * 4.182 + HCAP(2)) * Theta_V(ML, ND) + HCAP(6) * XII; % The Calorie should be converted as J - else - ZETA(ML, ND) = ZETA0(J) + (ZETA(ML, ND) - ZETA0(J)) * Theta_LL(ML, ND) / XXX; - TETCON(ML, ND) = TCON0(J) + (TETCON(ML, ND) - TCON0(J)) * Theta_LL(ML, ND) / XXX; - EHCAP(ML, ND) = HCD(J) + HCAP(1) * Theta_LL(ML, ND); - EHCAP(ML, ND) = EHCAP(ML, ND) + (0.448 * RHOV(MN) * 4.182 + HCAP(2)) * Theta_V(ML, ND) + HCAP(6) * XII; % The Calorie should be converted as J - end - %%%%%%%%%%%%%%%%%% Farouki thermal parameter method %%%%%%%%%%% - %%%%%%%%%%%%%%%%% change 21 July 2017 switch of thermal conductivity%%%%%% - if ThermCond == 4 %|| ThermCond==2 - % FEHCAP(ML,ND)=(2.128*Theta_sa(J)+2.385*Theta_cl(J))/(Theta_sa(J)+Theta_cl(J))*1e6; %J m-3 K-1 - % TCON_s(J)=(8.8*Theta_sa(J)+2.92*Theta_cl(J))/(Theta_sa(J)+Theta_cl(J)); % W m-1 K-1 - FFEHCAP(ML, ND) = (FEHCAP(J) * (1 - POR(J)) + 1000 * 4217.7 * Theta_LL(ML, ND) + 917 * 2117.27 * XII) * 1e-6; % The Calorie should be converted as J - EHCAP(ML, ND) = FFEHCAP(ML, ND); - end - %%%%%%%%%%%%%%%%%% Johansen thermal conductivity method %%%%%%%%%%% - TCON_qtz = 7.7; - TCON_o = 2.0; - TCON_L = 0.57; % Theta_qtz(J)=0.47; % thermal conductivities of soil quartz, other soil particles and water; unit W m-1 K-1 - if ThermCond == 2 - TCON_s(J) = TCON_qtz^(Theta_qtz(J)) * TCON_o^(1 - Theta_qtz(J)); % solid soil thermal conductivity Unit W m-1 K-1 - end - TCON_i(J) = 2.2; - %%%%% Peters-Lidard et al. 1998 frozen soil thermal conductivity method %%%%%%% - Satr(ML, ND) = (Theta_LL(ML, ND) + XII) / Theta_s(J); - if Theta_II(ML, ND) <= 0 - if Theta_LL(ML, ND) / POR(J) > 0.1 - K_e(ML, ND) = log10(Theta_LL(ML, ND) / POR(J)) + 1; % Kersten coefficient, weighting dry and wet thermal conductivity - elseif Theta_LL(ML, ND) / POR(J) > 0.05 - K_e(ML, ND) = 0.7 * log10(Theta_LL(ML, ND) / POR(J)) + 1; % Kersten coefficient, weighting dry and wet thermal conductivity - else - K_e(ML, ND) = 0; - end - Coef_k = 1.9; % [4.6 3.55 1.9; 1.7 0.95 0.85]; - K_e1(ML, ND) = Coef_k * (Theta_LL(ML, ND) / POR(J)) / (1 + (Coef_k - 1) * (Theta_LL(ML, ND) / POR(J))); % Kersten coefficient, weighting dry and wet thermal conductivity - ALPHA = 0.9; % [1.05,0.9,0.58]; - K_e2(ML, ND) = exp(ALPHA * (1 - (Theta_LL(ML, ND) / POR(J))^(ALPHA - 1.33))); % Lu and Ren 2007; ALPHA=1.05,0.9,0.58 - TCON_sat(ML, ND) = TCON_s(J)^(1 - Theta_s(J)) * TCON_L^(Theta_s(J)); % saturated soil thermal conductivity Unit W m-1 K-1 - EfTCON(ML, ND) = K_e(ML, ND) * (TCON_sat(ML, ND) - TCON_dry(J)) + TCON_dry(J); - else - K_e(ML, ND) = (Theta_LL(ML, ND) + Theta_II(ML, ND)) / POR(J); % +Theta_II(ML,ND) - TCON_sat(ML, ND) = TCON_s(J)^(1 - Theta_s(J)) * TCON_i(J)^(Theta_s(J) - Theta_LL(ML, ND)) * TCON_L^(Theta_LL(ML, ND)); % saturated soil thermal conductivity Unit W m-1 K-1 - EfTCON(ML, ND) = K_e(ML, ND) * (TCON_sat(ML, ND) - TCON_dry(J)) + TCON_dry(J); - end - end - end diff --git a/src/STEMMUS_SCOPE.m b/src/STEMMUS_SCOPE.m index 14adf912..6c6709f7 100644 --- a/src/STEMMUS_SCOPE.m +++ b/src/STEMMUS_SCOPE.m @@ -46,7 +46,7 @@ ModelSettings = io.getModelSettings(); global J rwuef SWCC Thmrlefc Soilairefc hThmrl KT TIME Delt_t NN ML nD -global W_Chg ThmrlCondCap ThermCond SSUR fc T0 rroot SAVE NL DeltZ +global fc T0 rroot SAVE NL DeltZ NL = ModelSettings.NL; DeltZ = ModelSettings.DeltZ; DeltZ_R = ModelSettings.DeltZ_R; @@ -55,10 +55,6 @@ Thmrlefc = ModelSettings.Thmrlefc; Soilairefc = ModelSettings.Soilairefc; hThmrl = ModelSettings.hThmrl; -W_Chg = ModelSettings.W_Chg; -ThmrlCondCap = ModelSettings.ThmrlCondCap; -ThermCond = ModelSettings.ThermCond; -SSUR = ModelSettings.SSUR; fc = ModelSettings.fc; T0 = ModelSettings.T0; rwuef = ModelSettings.rwuef; @@ -83,15 +79,15 @@ global MN ND hOLD TOLD h hh T TT P_g P_gg Evap QMT hN Trap RWU EVAP theta_s0 Ks0 global Precip frac SUMTIME TTT Theta_LLL CHK Theta_LL Theta_L Theta_UUU Theta_UU -global Theta_U Theta_III Theta_II AVAIL0 TIMEOLD SRT ALPHA alpha_h bx Srt CTT_PH +global Theta_U Theta_III Theta_II AVAIL0 SRT alpha_h bx Srt CTT_PH global CTT_LT CTT_g CTT_Lg c_unsat DhDZ DTDZ DRHOVZ QL QL_h QL_T QV Qa KL_h Chh ChT global Khh KhT Resis_a KfL_h KfL_T TT_CRIT h_frez L_f CTT EPCT DTheta_LLh DTheta_LLT -global DTheta_UUh CKT Lambda_eff EfTCON TETCON DDhDZ DEhBAR DRHOVhDz EtaBAR D_Vg -global DRHOVTDz KLhBAR KLTBAR DTDBAR SAVEDTheta_LLh SAVEDTheta_UUh QVT QVH Sa HR QVa -global QLH QLT DVH DVT Se QL_a DPgDZ k_g V_A Theta_V W WW D_Ta Ratio_ice thermal Xaa -global XaT Xah KL_T DRHOVT DRHOVh DRHODAt DRHODAz Theta_g Alpha_Lg Beta_g D_V D_A Eta -global ZETA MU_W Ks RHODA RHOV ETCON EHCAP L Evapo Beta_gBAR Alpha_LgBAR Gvc -global sfactortot sfactor fluxes lEstot lEctot NoTime Tss +global DTheta_UUh Lambda_eff DDhDZ DEhBAR DRHOVhDz EtaBAR D_Vg +global DRHOVTDz KLhBAR KLTBAR DTDBAR SAVEDTheta_LLh SAVEDTheta_UUh QVT QVH HR QVa +global QLH QLT DVH DVT Se QL_a DPgDZ V_A Theta_V W WW D_Ta thermal Xaa +global XaT Xah KL_T DRHOVT DRHOVh DRHODAt DRHODAz Theta_g Beta_g D_V Eta +global Ks RHODA RHOV L Evapo Gvc +global sfactortot sfactor fluxes lEstot lEctot Tss % Get initial values InitialValues = init.defineInitialValues(TimeProperties.Dur_tot); @@ -100,17 +96,9 @@ Srt = InitialValues.Srt; SAVEDTheta_UUh = InitialValues.SAVEDTheta_UUh; SAVEDTheta_LLh = InitialValues.SAVEDTheta_LLh; -Ratio_ice = InitialValues.Ratio_ice; -KL_T = InitialValues.KL_T; Lambda_eff = InitialValues.Lambda_eff; -W = InitialValues.W; -WW = InitialValues.WW; -MU_W = InitialValues.MU_W; -D_Ta = InitialValues.D_Ta; D_V = InitialValues.D_V; Eta = InitialValues.Eta; -D_A = InitialValues.D_A; -EHCAP = InitialValues.EHCAP; Chh = InitialValues.Chh; ChT = InitialValues.ChT; Khh = InitialValues.Khh; @@ -118,10 +106,7 @@ QL = InitialValues.QL; QL_h = InitialValues.QL_h; QL_T = InitialValues.QL_T; -k_g = InitialValues.k_g; -Sa = InitialValues.Sa; V_A = InitialValues.V_A; -Alpha_Lg = InitialValues.Alpha_Lg; Beta_g = InitialValues.Beta_g; c_unsat = InitialValues.c_unsat; CTT_PH = InitialValues.CTT_PH; @@ -176,11 +161,9 @@ hOLD = InitialValues.hOLD; TOLD = InitialValues.TOLD; -global f0 L_WT Kha Vvh VvT Chg C1 C2 C3 C4 C5 C6 Cah CaT Caa Kah KaT Kaa Vah VaT Vaa Cag CTh CTa KTh KTT KTa -global VTT VTh VTa CTg Kcva Kcah KcaT Kcaa Ccah CcaT Ccaa Ksoil SMC bbx wfrac Ta Ts U HR_a Rns Rnl Rn +global Kha Vvh VvT Chg C1 C2 C3 C4 C5 C6 Cah CaT Caa Kah KaT Kaa Vah VaT Vaa Cag CTh CTa KTh KTT KTa +global VTT VTh VTa CTg Kcva Kcah KcaT Kcaa Ccah CcaT Ccaa SMC bbx Ta Ts U HR_a Rns Rn global RHOV_s DRHOV_sT Tbtm r_a_SOIL Rn_SOIL SH MO Zeta_MO TopPg Tp_t RHS C7 C9 -f0 = InitialValues.f0; -L_WT = InitialValues.L_WT; Kha = InitialValues.Kha; Vvh = InitialValues.Vvh; VvT = InitialValues.VvT; @@ -217,16 +200,13 @@ Ccah = InitialValues.Ccah; CcaT = InitialValues.CcaT; Ccaa = InitialValues.Ccaa; -Ksoil = InitialValues.Ksoil; SMC = InitialValues.SMC; bbx = InitialValues.bbx; -wfrac = InitialValues.wfrac; Ta = InitialValues.Ta; Ts = InitialValues.Ts; U = InitialValues.U; HR_a = InitialValues.HR_a; Rns = InitialValues.Rns; -Rnl = InitialValues.Rnl; Rn = InitialValues.Rn; SH = InitialValues.SH; MO = InitialValues.MO; @@ -245,16 +225,12 @@ %% 1. define Constants Constants = io.define_constants(); -global g RHOL RHOI Rv RDA MU_a Lambda1 Lambda2 Lambda3 c_a c_V c_L Hc c_i Gamma0 Gamma_w RHO_bulk Rl +global g RHOL RHOI Rv RDA c_a c_V c_L Hc c_i Gamma0 Gamma_w Rl g = Constants.g; RHOL = Constants.RHOL; RHOI = Constants.RHOI; Rv = Constants.Rv; RDA = Constants.RDA; -MU_a = Constants.MU_a; -Lambda1 = Constants.Lambda1; -Lambda2 = Constants.Lambda2; -Lambda3 = Constants.Lambda3; c_L = Constants.c_L; c_V = Constants.c_V; c_a = Constants.c_a; @@ -440,32 +416,19 @@ [SoilVariables, VanGenuchten, ThermalConductivity] = StartInit(SoilVariables, SoilProperties, VanGenuchten); %% get variables that are defined global and are used by other scripts -global hm hd hh_frez XWRE POR IH IS XK XWILT KLT_Switch DVT_Switch KaT_Switch ISFT Imped XSOC -global Lamda Phi_s XCAP Gama_hh Gama_h SAVEhh COR CORh m n Alpha TCON_dry TPS1 TPS2 TCON0 TCON_s -global Theta_s Theta_r Theta_f HCAP SF TCA GA1 GA2 GB1 GB2 HCD ZETA0 CON0 PS1 PS2 FEHCAP +global hd hh_frez POR KaT_Switch XSOC +global XCAP SAVEhh COR CORh m n Alpha +global Theta_s Theta_r Theta_f -% get soil constants for StartInit +% get soil constants SoilConstants = io.getSoilConstants(); -hm = SoilConstants.hm; hd = SoilConstants.hd; -XWRE = SoilVariables.XWRE; POR = SoilVariables.POR; -IH = SoilVariables.IH; -IS = SoilVariables.IS; XK = SoilVariables.XK; -XWILT = SoilVariables.XWILT; -KLT_Switch = SoilVariables.KLT_Switch; -DVT_Switch = SoilVariables.DVT_Switch; KaT_Switch = SoilVariables.KaT_Switch; -ISFT = SoilVariables.ISFT; -Imped = SoilVariables.Imped; XSOC = SoilVariables.XSOC; -Lamda = SoilVariables.Lamda; -Phi_s = SoilVariables.Phi_s; XCAP = SoilVariables.XCAP; -Gama_hh = SoilVariables.Gama_hh; -Gama_h = SoilVariables.Gama_h; SAVEhh = SoilVariables.SAVEhh; Theta_s = VanGenuchten.Theta_s; Theta_r = VanGenuchten.Theta_r; @@ -473,25 +436,6 @@ Alpha = VanGenuchten.Alpha; n = VanGenuchten.n; m = VanGenuchten.m; -HCAP = ThermalConductivity.HCAP; -SF = ThermalConductivity.SF; -TCA = ThermalConductivity.TCA; -GA1 = ThermalConductivity.GA1; -GA2 = ThermalConductivity.GA2; -GB1 = ThermalConductivity.GB1; -GB2 = ThermalConductivity.GB2; -HCD = ThermalConductivity.HCD; -ZETA0 = ThermalConductivity.ZETA0; -CON0 = ThermalConductivity.CON0; -PS1 = ThermalConductivity.PS1; -PS2 = ThermalConductivity.PS2; -TCON_s = ThermalConductivity.TCON_s; -TCON_dry = ThermalConductivity.TCON_dry; -RHO_bulk = ThermalConductivity.RHO_bulk; -TPS1 = ThermalConductivity.TPS1; -TPS2 = ThermalConductivity.TPS2; -FEHCAP = ThermalConductivity.FEHCAP; -TCON0 = ThermalConductivity.TCON0; %% these vars are defined as global at the begining of this script %% because they are both input and output of StartInit @@ -508,7 +452,7 @@ BoundaryCondition = init.setBoundaryCondition(SoilVariables, ForcingData, landcoverClass(1)); %% get global vars -global NBCh NBCT NBChB NBCTB BCh DSTOR DSTOR0 RS NBChh DSTMAX IRPT1 IRPT2 NBCP BChB BCTB BCPB BCT BCP BtmPg +global NBCh NBCT NBChB NBCTB BCh DSTOR0 NBChh NBCP BChB BCTB BCPB BCT BCP BtmPg NBCh = BoundaryCondition.NBCh; NBCT = BoundaryCondition.NBCT; NBChB = BoundaryCondition.NBChB; @@ -580,7 +524,7 @@ L_f = 0; % ignore Freeze/Thaw, see issue 139 TT_CRIT(NN) = T0; % unit K hOLD_frez = []; - if IRPT1 == 0 && IRPT2 == 0 && ISFT == 0 + if IRPT1 == 0 && IRPT2 == 0 && SoilVariables.ISFT == 0 for MN = 1:NN hOLD_frez(MN) = h_frez(MN); h_frez(MN) = hh_frez(MN); @@ -837,15 +781,38 @@ DTheta_UUh = SoilVariables.DTheta_UUh; Theta_II = SoilVariables.Theta_II; - [KL_T] = CondL_T(NL); + % Reset KL_T here. CondL_T script is replaced by this line + % see issue 181, item 4 + KL_T = InitialValues.KL_T; + [RHOV, DRHOVh, DRHOVT] = Density_V(TT, hh, g, Rv, NN); - [W, WW, MU_W, D_Ta] = CondL_Tdisp(InitialValues, POR, Theta_LL, Theta_L, SSUR, RHOL, TT, Theta_s, h, hh, W_Chg, NL, nD, Delt_t, Theta_g, KLT_Switch); + + % update inputs + SoilVariables.Theta_L = Theta_L; + TransportCoefficient = conductivity.calculateTransportCoefficient(InitialValues, SoilVariables, VanGenuchten, Delt_t); + W = TransportCoefficient.W; + WW = TransportCoefficient.WW; + D_Ta = TransportCoefficient.D_Ta; + [L] = Latent(TT, NN); [Xaa, XaT, Xah, DRHODAt, DRHODAz, RHODA] = Density_DA(T, RDA, P_g, Rv, DeltZ, h, hh, TT, P_gg, Delt_t, NL, NN, DRHOVT, DRHOVh, RHOV); - [c_unsat, Lambda_eff, ZETA, ETCON, EHCAP, TETCON, EfTCON] = CondT_coeff(Theta_LL, Lambda1, Lambda2, Lambda3, RHO_bulk, Theta_g, RHODA, RHOV, c_a, c_V, c_L, NL, nD, ThmrlCondCap, ThermCond, HCAP, SF, TCA, GA1, GA2, GB1, GB2, HCD, ZETA0, CON0, PS1, PS2, XWILT, XK, TT, POR, DRHOVT, L, D_A, Theta_V, Theta_II, TCON_dry, Theta_s, XSOC, TPS1, TPS2, TCON0, TCON_s, FEHCAP, RHOI, RHOL, c_unsat, Lambda_eff, ETCON, EHCAP, TETCON, EfTCON, ZETA); - [k_g] = Condg_k_g(POR, NL, m, Theta_g, g, MU_W, Ks, RHOL, SWCC, Imped, Ratio_ice, Soilairefc, MN); - [D_V, Eta, D_A] = CondV_DE(Theta_LL, TT, fc, Theta_s, NL, nD, Theta_g, POR, ThmrlCondCap, ZETA, XK, DVT_Switch, Theta_UU); - [D_Vg, V_A, Beta_g, DPgDZ, Beta_gBAR, Alpha_LgBAR] = CondV_DVg(P_gg, Theta_g, Sa, V_A, k_g, MU_a, DeltZ, Alpha_Lg, KaT_Switch, Theta_s, Se, NL, DPgDZ, Beta_gBAR, Alpha_LgBAR, Beta_g); + + Theta_LL = SoilVariables.Theta_LL; + ThermalConductivityCapacity = conductivity.calculateThermalConductivityCapacity(InitialValues, ThermalConductivity, SoilVariables, VanGenuchten, DRHOVT, L, RHOV); + c_unsat = ThermalConductivityCapacity.c_unsat; + Lambda_eff = ThermalConductivityCapacity.Lambda_eff; + + k_g = conductivity.calculateGasConductivity(InitialValues, TransportCoefficient, VanGenuchten, SoilVariables); + + VaporVariables = conductivity.calculateVaporVariables(InitialValues, SoilVariables, VanGenuchten, ThermalConductivityCapacity, TT); + D_V = VaporVariables.D_V; + Eta = VaporVariables.Eta; + + GasDispersivity = conductivity.calculateGasDispersivity(InitialValues, SoilVariables, P_gg, k_g); + D_Vg = GasDispersivity.D_Vg; + V_A = GasDispersivity.V_A; + Beta_g = GasDispersivity.Beta_g; + DPgDZ = GasDispersivity.DPgDZ; run h_sub; if NBCh == 1 diff --git a/src/StartInit.m b/src/StartInit.m index 6e32517f..7c6ec780 100644 --- a/src/StartInit.m +++ b/src/StartInit.m @@ -24,9 +24,7 @@ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ThermalConductivity = init.calculateInitialThermal(SoilVariables, VanGenuchten); - % these will be updated in UpdateSoilWaterContent - SoilVariables.COR = []; - SoilVariables.CORh = []; + % SoilVariables will be updated in UpdateSoilWaterContent SoilVariables = UpdateSoilWaterContent(KIT, LatentHeatOfFreezing, SoilVariables, VanGenuchten); Theta_L = SoilVariables.Theta_L; diff --git a/src/UpdateSoilWaterContent.m b/src/UpdateSoilWaterContent.m index 61b19ba9..813a6907 100644 --- a/src/UpdateSoilWaterContent.m +++ b/src/UpdateSoilWaterContent.m @@ -3,92 +3,48 @@ Update SoilWaterContent i.e. Theta_LL. %} - % these lines can be removed after issue 181 % get model settings ModelSettings = io.getModelSettings(); - NN = ModelSettings.NN; - NL = ModelSettings.NL; - SWCC = ModelSettings.SWCC; - Thmrlefc = ModelSettings.Thmrlefc; - T0 = ModelSettings.T0; - Tr = ModelSettings.Tr; - hThmrl = ModelSettings.hThmrl; - Hystrs = ModelSettings.Hystrs; - Theta_s = VanGenuchten.Theta_s; - Theta_r = VanGenuchten.Theta_r; - Alpha = VanGenuchten.Alpha; - n = VanGenuchten.n; - m = VanGenuchten.m; - - POR = SoilVariables.POR; - h = SoilVariables.h; - Lamda = SoilVariables.Lamda; - Phi_s = SoilVariables.Phi_s; - - h_frez = SoilVariables.h_frez; - hh_frez = SoilVariables.hh_frez; - TT = SoilVariables.TT; - hh = SoilVariables.hh; - COR = SoilVariables.COR; - CORh = SoilVariables.CORh; - XWRE = SoilVariables.XWRE; - IH = SoilVariables.IH; - Ks = SoilVariables.Ks; - Imped = SoilVariables.Imped; - XCAP = SoilVariables.XCAP; - - Theta_L = SoilVariables.Theta_L; - Theta_LL = SoilVariables.Theta_LL; - Theta_V = SoilVariables.Theta_V; - Theta_g = SoilVariables.Theta_g; - Se = SoilVariables.Se; - KL_h = SoilVariables.KL_h; - DTheta_LLh = SoilVariables.DTheta_LLh; - KfL_T = SoilVariables.KfL_T; - Theta_II = SoilVariables.Theta_II; - Theta_I = SoilVariables.Theta_I; - Theta_UU = SoilVariables.Theta_UU; - Theta_U = SoilVariables.Theta_U; - TT_CRIT = SoilVariables.TT_CRIT; - KfL_h = SoilVariables.KfL_h; - DTheta_UUh = SoilVariables.DTheta_UUh; - % get Constants - Constants = io.define_constants(); - g = Constants.g; - RHOI = Constants.RHOI; - RHOL = Constants.RHOL; + COR = []; + CORh = []; - if hThmrl == 1 - for MN = 1:NN + if ModelSettings.hThmrl == 1 + for MN = 1:ModelSettings.NN CORh(MN) = 0.0068; - COR(MN) = exp(-1 * CORh(MN) * (TT(MN) - Tr)); % *COR21(MN) + COR(MN) = exp(-1 * CORh(MN) * (SoilVariables.TT(MN) - ModelSettings.Tr)); % *COR21(MN) if COR(MN) == 0 COR(MN) = 1; end end else - for MN = 1:NN + for MN = 1:ModelSettings.NN COR(MN) = 1; end end - for MN = 1:NN - hhU(MN) = COR(MN) * hh(MN); - hh(MN) = hhU(MN); + for MN = 1:ModelSettings.NN + hhU(MN) = COR(MN) * SoilVariables.hh(MN); + SoilVariables.hh(MN) = hhU(MN); end - [Theta_LL, Se, KfL_h, KfL_T, DTheta_LLh, hh, hh_frez, Theta_UU, DTheta_UUh, Theta_II, KL_h] = CondL_h(SoilVariables, Theta_r, Theta_s, Alpha, hh, hh_frez, h_frez, n, m, Ks, NL, Theta_L, h, KIT, TT, Thmrlefc, POR, SWCC, Theta_U, XCAP, Phi_s, RHOI, RHOL, Lamda, Imped, L_f, g, T0, TT_CRIT, Theta_II, KfL_h, KfL_T, KL_h, Theta_UU, Theta_LL, DTheta_LLh, DTheta_UUh, Se); - for MN = 1:NN - hhU(MN) = hh(MN); - hh(MN) = hhU(MN) / COR(MN); + + SoilVariables = conductivity.calculateHydraulicConductivity(SoilVariables, VanGenuchten, KIT, L_f); + Theta_LL = SoilVariables.Theta_LL; + DTheta_LLh = SoilVariables.DTheta_LLh; + Theta_UU = SoilVariables.Theta_UU; + Theta_II = SoilVariables.Theta_II; + + for MN = 1:ModelSettings.NN + hhU(MN) = SoilVariables.hh(MN); + SoilVariables.hh(MN) = hhU(MN) / COR(MN); end - if Hystrs == 0 - for i = 1:NL + if ModelSettings.Hystrs == 0 + for i = 1:ModelSettings.NL for ND = 1:2 - Theta_V(i, ND) = POR(i) - Theta_UU(i, ND); % -Theta_II(i,ND); % Theta_LL==>Theta_UU + Theta_V(i, ND) = SoilVariables.POR(i) - Theta_UU(i, ND); % -Theta_II(i,ND); % Theta_LL==>Theta_UU if Theta_V(i, ND) <= 1e-14 Theta_V(i, ND) = 1e-14; end @@ -96,34 +52,34 @@ end end else - for i = 1:NL + for i = 1:ModelSettings.NL for ND = 1:2 - if IH(i) == 2 - if XWRE(i, ND) < Theta_LL(i, ND) - Theta_V(i, ND) = POR(i) - Theta_LL(i, ND) - Theta_II(i, ND); + if SoilVariables.IH(i) == 2 + if SoilVariables.XWRE(i, ND) < Theta_LL(i, ND) + Theta_V(i, ND) = SoilVariables.POR(i) - Theta_LL(i, ND) - Theta_II(i, ND); else XSAVE = Theta_LL(i, ND); - Theta_LL(i, ND) = XSAVE * (1 + (XWRE(i, ND) - Theta_LL(i, ND)) / Theta_s(i)); + Theta_LL(i, ND) = XSAVE * (1 + (SoilVariables.XWRE(i, ND) - Theta_LL(i, ND)) / Theta_s(i)); if KIT > 0 DTheta_LLh(i, ND) = DTheta_LLh(i, ND) * (Theta_LL(i, ND) / XSAVE - XSAVE / Theta_s(i)); end - Theta_V(i, ND) = POR(i) - Theta_LL(i, ND) - Theta_II(i, ND); + Theta_V(i, ND) = SoilVariables.POR(i) - Theta_LL(i, ND) - Theta_II(i, ND); end end - if IH(i) == 1 - if XWRE(i, ND) > Theta_LL(i, ND) + if SoilVariables.IH(i) == 1 + if SoilVariables.XWRE(i, ND) > Theta_LL(i, ND) XSAVE = Theta_LL(i, ND); Theta_LL(i, ND) = (2 - XSAVE / Theta_s(i)) * XSAVE; if KIT > 0 DTheta_LLh(i, ND) = 2 * DTheta_LLh(i, ND) * (1 - XSAVE / Theta_s(i)); end - Theta_V(i, ND) = POR(i) - Theta_LL(i, ND) - Theta_II(i, ND); + Theta_V(i, ND) = SoilVariables.POR(i) - Theta_LL(i, ND) - Theta_II(i, ND); else - Theta_LL(i, ND) = Theta_LL(i, ND) + XWRE(i, ND) * (1 - Theta_LL(i, ND) / Theta_s(i)); + Theta_LL(i, ND) = Theta_LL(i, ND) + SoilVariables.XWRE(i, ND) * (1 - Theta_LL(i, ND) / Theta_s(i)); if KIT > 0 - DTheta_LLh(i, ND) = DTheta_LLh(i, ND) * (1 - XWRE(i, ND) / Theta_s(i)); + DTheta_LLh(i, ND) = DTheta_LLh(i, ND) * (1 - SoilVariables.XWRE(i, ND) / Theta_s(i)); end - Theta_V(i, ND) = POR(i) - Theta_LL(i, ND) - Theta_II(i, ND); + Theta_V(i, ND) = SoilVariables.POR(i) - Theta_LL(i, ND) - Theta_II(i, ND); end end if Theta_V(i, ND) <= 1e-14 % consider the negative conditions? @@ -133,20 +89,11 @@ end end end - SoilVariables.KfL_T = KfL_T; - SoilVariables.Theta_II = Theta_II; - SoilVariables.Theta_UU = Theta_UU; - - SoilVariables.hh_frez = hh_frez; - SoilVariables.hh = hh; SoilVariables.COR = COR; SoilVariables.CORh = CORh; - SoilVariables.Se = Se; - SoilVariables.KL_h = KL_h; - SoilVariables.Theta_LL = Theta_LL; - SoilVariables.DTheta_LLh = DTheta_LLh; - SoilVariables.KfL_h = KfL_h; SoilVariables.Theta_V = Theta_V; SoilVariables.Theta_g = Theta_g; - SoilVariables.DTheta_UUh = DTheta_UUh; + SoilVariables.Theta_LL = Theta_LL; + SoilVariables.DTheta_LLh = DTheta_LLh; + end