Flow PKs#

Flow PKs describe the conservation of mass of water as it flows both above and below-ground. Subsurface flow PKs are based on 3D Richards equation, which describes variably saturated flow in porous media. Variations on this include the incorporation of freeze-thaw processes. Surface flow PKs are based on a diffusion wave equation and Manning’s model for sheet flow. Variations on this also include the incorporation of freeze-thaw processes. Finally we include in flow a “snow distribution” algorithm which takes as input precipitation and applies it based on the existing surface level (elevation + water + snowpack), thereby “filling in” low-lying areas preferentially. This makes for more accurate snowpacks at fine scales.

Richards PK#

src/physics/ats/src/pks/flow/richards.hh

Two-phase, variable density Richards equation.

Solves Richards equation:

\[\frac{\partial \Theta}{\partial t} - \nabla \cdot \frac{k_r n_l}{\mu} K ( \nabla p + \rho g \hat{z} ) = Q_w\]

“PK type” = “richards flow”

pk-richards-flow-spec

  • “domain” [string] “domain” Defaults to the subsurface mesh.

  • “primary variable key” [string] DOMAIN-pressure The primary variable associated with this PK.

  • “boundary conditions” [list] Defaults to Neuman, 0 normal flux. See Flow-specific Boundary Conditions

IF

  • “source term” [bool] false Is there a source term?

THEN

  • “source key” [string] DOMAIN-water_source Typically not set, as the default is good. \([mol \, s^{-1}]\)

  • “explicit source term” [bool] false Apply the source term from the previous timestep (source is lagged).

END

  • “permeability rescaling” [double] 1e7 Typically 1e7 or order \(sqrt(K)\) is about right. This rescales both the relative and absolute permeabilities to avoid multiplying small numbers (absolute permeability) by a large number (\(\rho / \mu\)).

Math and solver algorithm options:

  • “absolute error tolerance” [double] 2750.0 in units of \([mol]\).

  • “inverse” [inverse-spec] optional The inverse used for preconditioning in a non-globally coupled problem. See Inverse.

  • “diffusion” [pde-diffusion-typedinline-spec] The (forward) diffusion operator, see Diffusion.

  • “diffusion preconditioner” [pde-diffusion-spec] optional The inverse of the diffusion operator. See Diffusion. Typically this is only needed to set Jacobian options, as all others probably should match those in “diffusion”, and default to those values.

  • “compute boundary values” [bool] false Used to include boundary face unknowns on diffusion discretizations that are cell-only (e.g. FV). This can be useful for surface flow or other wierd boundary conditions. Usually provided by MPCs that need them.

  • “accumulation preconditioner” [pde-accumulation-spec] optional The inverse of the accumulation operator. See Accumulation. Typically not provided by users, as defaults are sufficient.

  • “surface rel perm strategy” [string] none Approach for specifying the relative permeabiilty on the surface face. “clobber” is frequently used for cases where a surface rel perm will be provided. TODO: this should be absorbed into an evaluator. One of:

    • “none” : use the upwind direction to determine whether to use the boundary face or internal cell

    • “clobber” : always use the boundary face rel perm

    • “max” : use the max of the boundary face and internal cell values

    • “unsaturated” : Uses the boundary face when the internal cell is not saturated.

  • “relative permeability method” [string] upwind with Darcy flux Relative permeability is defined on cells, but must be calculated on faces to multiply a flux. There are several methods commonly used. Note these can significantly change answers – you don’t want to change these unless you know what they mean. One of:

    • “upwind with Darcy flux” First-order upwind method that is most common

    • “upwind with gravity” Upwinds according to the gravitational flux direction

    • “cell centered” This corresponds to the harmonic mean, and is most accurate if the problem is always wet, but has issues when it is dry.

    • “arithmetic mean” Face value is the mean of the neighboring cells. Not a good method.

Globalization and other process-based hacks:

  • “modify predictor with consistent faces” [bool] false In a face+cell diffusion discretization, this modifies the predictor to make sure that faces, which are a DAE, are consistent with the predicted cells (i.e. face fluxes from each sides match).

  • “modify predictor for flux BCs” [bool] false Infiltration into dry ground can be hard on solvers – this tries to do the local nonlinear problem to ensure that face pressures are consistent with the prescribed flux in a predictor.

  • “modify predictor via water content” [bool] false Modifies the predictor using the method of Krabbenhoft [??] paper. Effectively does a change of variables, extrapolating not in pressure but in water content, then takes the smaller of the two extrapolants.

  • “max valid change in saturation in a timestep [-]” [double] -1 Rejects timesteps whose max saturation change is greater than this value. This can be useful to ensure temporally resolved solutions. Usually a good value is 0.1 or 0.2.

  • “max valid change in ice saturation in a timestep [-]” [double] -1 Rejects timesteps whose max ice saturation change is greater than this value. This can be useful to ensure temporally resolved solutions. Usually a good value is 0.1 or 0.2.

  • “limit correction to pressure change [Pa]” [double] -1 If > 0, this limits an iterate’s max pressure change to this value. Not usually helpful.

  • “limit correction to pressure change when crossing atmospheric [Pa]” [double] -1 If > 0, this limits an iterate’s max pressure change to this value when they cross atmospheric pressure. Not usually helpful.

Algorithmic physics control:

IF

  • “coupled to surface via flux” [bool] false If true, apply surface boundary conditions from an exchange flux. Note, if this is a coupled problem, it is probably set by the MPC. No need for a user to set it.

THEN

  • “surface-subsurface flux key” [string] DOMAIN-surface_subsurface_flux

END

  • “coupled to surface via head” [bool] false If true, apply surface boundary conditions from the surface pressure (Dirichlet).

INCLUDES:

EVALUATORS: - “conserved quantity” - “mass density” - “molar density” - “permeability” - “conductivity” - “saturation”

KEYS:

  • “conserved quantity” DOMAIN-water_content

  • “mass density” DOMAIN-mass_density_liquid

  • “molar density” DOMAIN-molar_density_liquid

  • “permeability key” DOMAIN-permeability

  • “conductivity key” DOMAIN-relative_permeability

  • “upwind conductivity key” [string] DOMAIN-upwind_relative_permeability upwinded (face-based) scalar coefficient of the permeability. Note the units of this are strange, but this represents \(\frac{n_l k_r}{\mu}\) [mol kg^-1 s^1 m^-2]

  • “darcy flux key” DOMAIN-water_flux water flux across a face [mol s^-1]

  • “darcy flux direction key” DOMAIN-water_flux_direction direction of the darcy flux (used in upwinding \(k_r\)) [??]

  • “darcy velocity key” DOMAIN-darcy_velocity darcy velocity vector, interpolated from faces to cells [m s^-1]

  • “saturation key” DOMAIN-saturation_liquid [-]

Permafrost Flow PK#

src/physics/ats/src/pks/flow/permafrost.hh

A three-phase, thermal Richard’s equation with water, water vapor, and ice for permafrost applications.

Note that the only difference between permafrost and richards is in constitutive relations – the WRM changes to provide three saturations, while the water content changes to account for water in ice phase. As these are now drop-in field evaluators, there is very little to change in the PK.

In the future, this will not need a different PK.

“PK type” = “permafrost flow”

pk-permafrost-flow-spec

  • “saturation ice key” [string] “DOMAIN-saturation_ice” volume fraction of the ice phase (only when relevant) [-] Typically the default is correct.

INCLUDES:

Overland Flow PK#

src/physics/ats/src/pks/flow/overland_pressure.hh

The overland flow PK solves the diffusion wave equation for overland flow with pressure as a primary variable:

\[\frac{\partial \Theta}{\partial t} - \nabla n_l k \nabla h(p) = Q_w\]

“PK type” = “overland flow, pressure basis”

pk-overland-flow-pressure-basis-spec

  • “domain” [string] “surface” Defaults to the extracted surface mesh.

  • “primary variable” [string] DOMAIN-pressure The primary variable associated with this PK.

  • “boundary conditions” [list] Defaults to Neuman, 0 normal flux. See Flow-specific Boundary Conditions

IF

  • “source term” [bool] false Is there a source term?

THEN

  • “source key” [string] DOMAIN-water_source Typically not set, as the default is good. [m s^-1] or [mol s^-1]

  • “water source in meters” [bool] true Is the source term in [m s^-1]?

END

Math and solver algorithm options:

  • “absolute error tolerance” [double] 550. Defaults to 1 cm of water (in mols). A small, but significant, amount of water.

  • “inverse” [inverse-spec] optional The inverse used for preconditioning in a non-globally coupled problem. See Inverse.

  • “diffusion” [pde-diffusion-typedinline-spec] The (forward) diffusion operator, see Diffusion.

  • “diffusion preconditioner” [pde-diffusion-typedinline-spec] optional The inverse of the diffusion operator. See Diffusion. Typically this is only needed to set Jacobian options, as all others probably should match those in “diffusion”, and default to those values.

  • “accumulation preconditioner” [pde-accumulation-spec] optional The inverse of the accumulation operator. See Accumulation. Typically not provided by users, as defaults are sufficient.

Globalization:

  • “limit correction to pressure change [Pa]” [double] -1 If > 0, this limits an iterate’s max pressure change to this value. Not usually helpful.

  • “limit correction to pressure change when crossing atmospheric [Pa]” [double] -1 If > 0, this limits an iterate’s max pressure change to this value when they cross atmospheric pressure. Not usually helpful.

  • “allow no negative ponded depths” [bool] false Modifies all correction updates to ensure only positive ponded depth is allowed.

  • “min ponded depth for velocity calculation” [double] 1.e-2 For ponded depth below this height, declare the velocity 0.

  • “min ponded depth for tidal bc” [double] 0.02 Control on the tidal boundary condition. TODO: This should live in the BC spec?

Algorithmic physics control:

  • “coupled to subsurface via flux” [bool] false Set by MPC.

  • “coupled to subsurface via head” [bool] false Set by MPC.

INCLUDES:

KEYS:

  • “conserved quantity” DOMAIN-water_content [mol]

  • “elevation” DOMAIN-elevation [m]

  • “slope magnitude” DOMAIN-slope_magnitude [-]

EVALUATORS:

  • “conserved quantity”

  • “water content”

  • “cell volume”

  • “surface_subsurface_flux”

  • “elevation”

  • “slope magnitude”

  • “overland_conductivity”

  • “ponded_depth”

  • “pres_elev”

  • “source”

Overland Flow with Ice#

src/physics/ats/src/pks/flow/icy_overland.hh

This modifies the diffusion wave equation for overland flow that includes freeze-thaw processes.

“PK type” = “overland flow with ice”

pk-overland-flow-with-ice-spec

INCLUDES:

Snow Distribution PK#

src/physics/ats/src/pks/flow/snow_distribution.hh

This PK is a heuristic PK that distributes incoming snow precipitation using a diffusion wave equation. Think of it as an analogue to overland flow – it effectively ensures that new snow “flows downhill,” due to a uniformly distributed random direction wind, and lands on the lowest lying areas.

The Manning coefficient determines with how uniform the snow layer becomes. Most of the parameters are set depending upon the snow precipitation input data interval.

“PK type” = “snow distribution”

pk-snow-distribution-spec

  • “distribution time” [double] 86400. Interval of snow precip input dataset. [s]

  • “precipitation function” [function-spec] Snow precipitation function, see Functions.

  • “diffusion” [pde-diffusion-spec] Diffusion drives the distribution. Typically we use finite volume here. See Diffusion.

  • “diffusion preconditioner” [pde-diffusion-spec] Inverse of the above. Likely only Jacobian term options are needed here, as the others default to the same as the “diffusion” list. See Diffusion.

  • “inverse” [inverse-typed-spec] Inverse method for the solve.

  • “accumulation preconditioner” [pde-accumulation-spec] optional See Accumulation.