Uncoupled saturated/unsaturated porous flow with storage term

[nlosacco]

Hi all,

I am trying to solve an uncoupled, purely hydraulic flow problem, to simulate groundwater flow in a slope. This includes both saturated and unsaturated conditions, and the ability to switch from one to the other. As the simulation is purely hydraulic, I want to include some additional storage terms depending on the compressibility of porewater and the soil skeleton (not solid soil grains, which I can realistically assume incompressible, given the relatively low pressures involved). I read the docs thoroughly but still can't figure out how to solve the problem.
For the unsaturated bit I came up with this modified version of Richards' equation, when porewater pressure $P < 0$, and with the simplifying assumption of constant water density $\rho_\mathrm{w}$ (though it may be slightly compressible), and isotropic hydraulic conductivity :

$$\bigl[C(P) - (c_\mathrm{sk} + c_\mathrm{w})Sr \phi\bigr]\dfrac{\partial{P}}{\partial{t}} - \dfrac{K(P)}{\rho_\mathrm{w}g} \nabla^2P=0$$

I think the equation is correct and includes the additional desired storage terms:

$$(c_\mathrm{sk} + c_\mathrm{w})Sr \phi$$

Here:

• $C(P)$ and $K(P)$ are, respectively, the so-called water retention curve and hydraulic conductivity as functions of $P$, described by the Van-Genuchten and Mualem relations (for C and K respectively) ;
• $c_\mathrm{sk}$ and $c_\mathrm{w}$ are the compressibility of the soil skeleton and of water;
• $\phi$ and $S_\mathrm{r}$ are the porosity and degree of saturation;

Can you suggest how to proceed to include the storage term in the PDE?

Also, I think $S_\mathrm{r}$ is calculated from Van Genuchten's (as stated in the documentation), but $\phi$ as well should be updated according to the total volumetric change, i.e.:

$$ d\phi = c_\mathrm{sk} dP $$

[josebastiase]

Im not 100% sure that I understand the issue. Some thoughts, I hope it helps

If you want to use Richard's, your "storage" is somehow hidden in the porosity term.

So you may want to use PorousFlowPorosity and activate only the hydro (fluid) term and use solid_bulk and biot_coefficient. Im not sure if 1-alpha/K equals the compressibility term that you want, c_sk.

For a quick fix, I would set biot_coefficient and solid_bulk so the term 1-alpha/K equals your c_sk. For example, very small alpha and K = 1/c_sk.

Then, to the initial porosity I would set it as porosity_zero times (c_sk + c_w).

But, Im not sure this is what you want to do.

[nlosacco]

The last equation corresponds to a linear relation between $d\phi$ and $dP$, hence I guess a PorousFlowPorosityLinear with P_coeffequat to $c_\mathrm{sk} + c_\mathrm{w}$ and P_ref equal to the initial $P$ would do, right?
About the latter point, I guess I can store the initial value of porepressure in an AuxVariable... is that ok? How can I do that?

[GiudGiud]

If the storage is “hidden” in the porosity, isn’t there an additional equation there? It might be appropriate to make porosity a nonlinear variable

[josebastiase]

@GiudGiud yes is it nonlinear, but most of time is assumed constant or linear.

@nlosacco I thought you wanted to change porosity over time, if not, then a linear porosity should work. I would check the equations if the ones in porousflow match the ones you want to solve. Wang, 2000 is a good reference for poroelasticity, I guess your system is not very well consolidated so alpha is very close to one. Then, your biot modulus (inverse of storage at constant strain) is close to bulk modulus/porosity. See

modules/porous_flow/tutorial_01
modules/porous_flow/tidal

You can simulate your steady state case and then make a restart and then save the inital value of pp in a AuxVar, like in modules/porous_flow/examples/groundwater/ex02_abstraction.i

[nlosacco]

@GiudGiud I don't know, the various PorousFlowPorosity* are already there, and the relation is a simple linear one. Maybe adding porosity as an additional nonlinear variable is overkill for this simple case, what do you think?
In general, though, if I understand your suggestion, one might want to add his own ODE for $d\phi$ and couple it with the PorousFlowMassTimeDerivative, is that what you mean?
In this case, how should one deal with the original PorousFlowPorosity* that MOOSE is expecting as an input?

@josebastiase yes, the only problem with PorousFlowPorosityLinear is that it is not written in incremental terms, but following your suggestion on the AuxVar it should do the job. I'll let you know.
PS: to us, following Terzaghi's principle of effective stress, it is always $\alpha_\mathrm{B} = 1$, so that Eq. 10 here: https://mooseframework.inl.gov/modules/porous_flow/governing_equations.html becomes:

$$\sigma_{ij}^\mathrm{eff} = \sigma_{ij}^\mathrm{tot} + \delta_{ij}P_\mathrm{f}$$

[GiudGiud]

In this case, how should one deal with the original PorousFlowPorosity* that MOOSE is expecting as an input?

MOOSE is expecting a material property for porosity. If that material property is equal to a variable, we can still compute the value of the material property from the value of the variable. If using AD, we can still have the correct derivative terms for the Jacobian. If not, one would need to manually compute the derivatives of the material wrt the variable (here, 1) and include that in the kernels.

Maybe adding porosity as an additional nonlinear variable is overkill for this simple case, what do you think?

Yes you're probably right. Especially since porosity does not have spatial derivative terms so it can be solved for in every element independently. Keeping it as a material property, with an expression coded in a Material, and depending on the other variables works.

yes, the only problem with PorousFlowPorosityLinear is that it is not written in incremental terms

Please feel free to write a new material class for your application. This is normal when using MOOSE, we can't have code in the modules for every single application.