PorousFlowAquiferBC

Applies a Robin (aquifer) boundary condition: flux = conductance * (P_model - P_aquifer(z)), where P_aquifer is the far-field aquifer pressure. Using aquifer_head automatically yields zero flux on any hydrostatic boundary whose model head equals the aquifer head, even when the boundary has vertical extent.

PorousFlowAquiferBC applies a Robin (Cauchy) boundary condition that couples the model boundary to a far-field aquifer. The mass flux leaving the domain is proportional to the difference between the model pore pressure and the aquifer pressure evaluated at the same elevation: (1) where is the conductance (kg.m.Pa.s), is the nodal pore pressure, and is the elevation at the current quadrature point. A positive value of means fluid is leaving the domain (the boundary acts as a sink); a negative value means fluid is entering (a source).

The elevation is computed per quadrature point from the position vector and the user-supplied gravity vector : where is the unit vector in the direction of gravity (e.g. for gravity acting in : , so , which is the -coordinate). Because is re-evaluated at every quadrature point, the hydrostatic correction is correct on boundaries of any orientation, including vertical faces.

The fluid density used in the hydrostatic correction, , is taken from PorousFlow's own nodal material property. This makes the aquifer pressure self-consistent with the equation of state used in the rest of the model.

Reference-pressure formulations

Exactly one of the following two parameters must be supplied to define the far-field aquifer pressure.

Hydraulic head formulation (aquifer_head)

where is the far-field hydraulic head (m). This form has the important property that the flux is identically zero at every point on a boundary whose model pressure is at hydrostatic equilibrium with the aquifer, even when the boundary has vertical extent. This is because the hydraulic head is constant throughout a hydrostatic column, so everywhere.

The conductance must be supplied directly via aquifer_conductance (kg.m.Pa.s). A physically reasonable estimate is , where is the permeability of the intervening material (m), is the fluid viscosity (Pa.s), and is the distance to the far-field aquifer (m).

Pressure-at-datum formulation (aquifer_pressure_at_datum)

where is the aquifer pressure at reference elevation . This form is convenient when the aquifer pressure at a specific depth is known from field data.

For this formulation the conductance is computed automatically from the permeability and viscosity material properties: where is the permeability of the material between the boundary and the far-field aquifer, is the nodal fluid viscosity, and is the distance to the far-field aquifer supplied via aquifer_distance (m).

By default is the boundary permeability material property projected onto the outward boundary normal. This is appropriate in the common case where the mesh truncates a laterally extensive formation and the aquifer is a continuation of the boundary material; heterogeneity along the boundary and anisotropy are then handled automatically. When the intervening material differs from the boundary material (for example a low-permeability reservoir edge connected to a regional aquifer through higher-permeability strata), the optional aquifer_permeability parameter (m) overrides the boundary value.

Effect of aquifer distance

For the pressure-at-datum formulation the single parameter controls the strength of the hydraulic coupling. For the head formulation an equivalent distance can be recovered from the supplied conductance as . In either case:

  • : . An arbitrarily large flux is generated by any nonzero pressure difference, driving . The boundary approaches a Dirichlet condition pinned to the aquifer pressure.

  • : . The flux approaches zero regardless of the pressure difference. The boundary approaches an impermeable (Neumann) condition.

The equilibration time constant for a compressible domain of length , porosity , fluid bulk modulus , permeability , and viscosity is approximately so can be chosen to match a desired response timescale.

Comparison with PorousFlowPiecewiseLinearSink

PorousFlowPiecewiseLinearSink with pt_vals/multipliers encoding , flux_function = C, and PT_shift = Pe applies where is a single scalar. This is equivalent to PorousFlowAquiferBC only when the boundary is horizontal (constant elevation) and the density is approximately constant, so that can be pre-computed at a single elevation .

On a vertical or inclined boundary varies with elevation and a single scalar PT_shift matches the aquifer pressure at only one height. PorousFlowAquiferBC evaluates and per quadrature point, so it is correct on boundaries of any orientation.

A further difference is that PorousFlowAquiferBC uses from PorousFlow's own equation of state in the hydrostatic correction, so the aquifer pressure is consistent with the model fluid properties. A PorousFlowPiecewiseLinearSink with a pre-computed scalar PT_shift uses a fixed reference density, which introduces a small error for compressible fluids.

See boundary conditions for more details and examples.

Input Parameters

  • PorousFlowDictatorThe UserObject that holds the list of PorousFlow variable names

    C++ Type:UserObjectName

    Controllable:No

    Description:The UserObject that holds the list of PorousFlow variable names

  • boundaryThe list of boundary IDs from the mesh where this object applies

    C++ Type:std::vector<BoundaryName>

    Controllable:No

    Description:The list of boundary IDs from the mesh where this object applies

  • gravityGravitational acceleration vector (m/s^2), e.g. '0 0 -9.81'. The elevation at each quadrature point is computed as the component of the position vector in the direction opposite to gravity.

    C++ Type:libMesh::VectorValue<Real>

    Unit:(no unit assumed)

    Controllable:No

    Description:Gravitational acceleration vector (m/s^2), e.g. '0 0 -9.81'. The elevation at each quadrature point is computed as the component of the position vector in the direction opposite to gravity.

  • variableThe name of the variable that this residual object operates on

    C++ Type:NonlinearVariableName

    Unit:(no unit assumed)

    Controllable:No

    Description:The name of the variable that this residual object operates on

Required Parameters

  • aquifer_conductanceConductance per unit boundary area (kg/(m^2*Pa*s)). The mass flux leaving the domain is conductance * (P_model - P_aquifer). Required when using the aquifer_head formulation. Can be estimated as rho * k / (mu * L) where k is aquifer permeability (m^2), mu is fluid viscosity (Pa.s), and L is the distance to the far-field (m).

    C++ Type:Real

    Unit:(no unit assumed)

    Controllable:No

    Description:Conductance per unit boundary area (kg/(m^2*Pa*s)). The mass flux leaving the domain is conductance * (P_model - P_aquifer). Required when using the aquifer_head formulation. Can be estimated as rho * k / (mu * L) where k is aquifer permeability (m^2), mu is fluid viscosity (Pa.s), and L is the distance to the far-field (m).

  • aquifer_distanceDistance from the boundary to the far-field aquifer (m). Required when using the aquifer_pressure_at_datum formulation. The conductance is computed internally as rho * k_nn / (mu * aquifer_distance), where k_nn is the permeability projected onto the boundary normal, rho is the nodal fluid density, and mu is the nodal fluid viscosity.

    C++ Type:Real

    Unit:(no unit assumed)

    Range:aquifer_distance > 0

    Controllable:No

    Description:Distance from the boundary to the far-field aquifer (m). Required when using the aquifer_pressure_at_datum formulation. The conductance is computed internally as rho * k_nn / (mu * aquifer_distance), where k_nn is the permeability projected onto the boundary normal, rho is the nodal fluid density, and mu is the nodal fluid viscosity.

  • aquifer_headFar-field hydraulic head of the aquifer (m above model datum). P_aq(z) = rho_nodal * |g| * (aquifer_head - z), where rho_nodal is the PorousFlow nodal fluid density at the boundary node. Mutually exclusive with aquifer_pressure_at_datum.

    C++ Type:Real

    Unit:(no unit assumed)

    Controllable:No

    Description:Far-field hydraulic head of the aquifer (m above model datum). P_aq(z) = rho_nodal * |g| * (aquifer_head - z), where rho_nodal is the PorousFlow nodal fluid density at the boundary node. Mutually exclusive with aquifer_pressure_at_datum.

  • aquifer_permeabilityPermeability of the material between the boundary and the far-field aquifer (m^2), used as k_nn in the conductance formula of the aquifer_pressure_at_datum formulation. If not supplied, the boundary permeability projected onto the boundary normal is used, which is appropriate when the aquifer is a continuation of the boundary material.

    C++ Type:Real

    Unit:(no unit assumed)

    Range:aquifer_permeability > 0

    Controllable:No

    Description:Permeability of the material between the boundary and the far-field aquifer (m^2), used as k_nn in the conductance formula of the aquifer_pressure_at_datum formulation. If not supplied, the boundary permeability projected onto the boundary normal is used, which is appropriate when the aquifer is a continuation of the boundary material.

  • aquifer_pressure_at_datumFar-field aquifer pressure at datum_elevation (Pa). P_aq(z) = aquifer_pressure_at_datum + rho_nodal * |g| * (datum_elevation - z), where rho_nodal is the PorousFlow nodal fluid density at the boundary node. Mutually exclusive with aquifer_head.

    C++ Type:Real

    Unit:(no unit assumed)

    Controllable:No

    Description:Far-field aquifer pressure at datum_elevation (Pa). P_aq(z) = aquifer_pressure_at_datum + rho_nodal * |g| * (datum_elevation - z), where rho_nodal is the PorousFlow nodal fluid density at the boundary node. Mutually exclusive with aquifer_head.

  • datum_elevation0Elevation of the reference point for aquifer_pressure_at_datum (m).

    Default:0

    C++ Type:Real

    Unit:(no unit assumed)

    Controllable:No

    Description:Elevation of the reference point for aquifer_pressure_at_datum (m).

  • fluid_phaseIf supplied, then this BC will potentially be a function of fluid pressure, and you can use mass_fraction_component, use_mobility, use_relperm, use_enthalpy and use_energy. If not supplied, then this BC can only be a function of temperature

    C++ Type:unsigned int

    Controllable:No

    Description:If supplied, then this BC will potentially be a function of fluid pressure, and you can use mass_fraction_component, use_mobility, use_relperm, use_enthalpy and use_energy. If not supplied, then this BC can only be a function of temperature

  • flux_function1The flux. The flux is OUT of the medium: hence positive values of this function means this BC will act as a SINK, while negative values indicate this flux will be a SOURCE. The functional form is useful for spatially or temporally varying sinks. Without any use_*, this function is measured in kg.m^-2.s^-1 (or J.m^-2.s^-1 for the case with only heat and no fluids)

    Default:1

    C++ Type:FunctionName

    Unit:(no unit assumed)

    Controllable:No

    Description:The flux. The flux is OUT of the medium: hence positive values of this function means this BC will act as a SINK, while negative values indicate this flux will be a SOURCE. The functional form is useful for spatially or temporally varying sinks. Without any use_*, this function is measured in kg.m^-2.s^-1 (or J.m^-2.s^-1 for the case with only heat and no fluids)

  • mass_fraction_componentThe index corresponding to a fluid component. If supplied, the flux will be multiplied by the nodal mass fraction for the component

    C++ Type:unsigned int

    Controllable:No

    Description:The index corresponding to a fluid component. If supplied, the flux will be multiplied by the nodal mass fraction for the component

  • use_enthalpyFalseIf true, then fluxes are multiplied by enthalpy. In this case bare_flux is measured in kg.m^-2.s^-1 / (J.kg). This can be used in conjunction with other use_*

    Default:False

    C++ Type:bool

    Controllable:No

    Description:If true, then fluxes are multiplied by enthalpy. In this case bare_flux is measured in kg.m^-2.s^-1 / (J.kg). This can be used in conjunction with other use_*

  • use_internal_energyFalseIf true, then fluxes are multiplied by fluid internal energy. In this case bare_flux is measured in kg.m^-2.s^-1 / (J.kg). This can be used in conjunction with other use_*

    Default:False

    C++ Type:bool

    Controllable:No

    Description:If true, then fluxes are multiplied by fluid internal energy. In this case bare_flux is measured in kg.m^-2.s^-1 / (J.kg). This can be used in conjunction with other use_*

  • use_mobilityFalseIf true, then fluxes are multiplied by (density*permeability_nn/viscosity), where the '_nn' indicates the component normal to the boundary. In this case bare_flux is measured in Pa.m^-1. This can be used in conjunction with other use_*

    Default:False

    C++ Type:bool

    Controllable:No

    Description:If true, then fluxes are multiplied by (density*permeability_nn/viscosity), where the '_nn' indicates the component normal to the boundary. In this case bare_flux is measured in Pa.m^-1. This can be used in conjunction with other use_*

  • use_relpermFalseIf true, then fluxes are multiplied by relative permeability. This can be used in conjunction with other use_*

    Default:False

    C++ Type:bool

    Controllable:No

    Description:If true, then fluxes are multiplied by relative permeability. This can be used in conjunction with other use_*

  • use_thermal_conductivityFalseIf true, then fluxes are multiplied by thermal conductivity projected onto the normal direction. This can be used in conjunction with other use_*

    Default:False

    C++ Type:bool

    Controllable:No

    Description:If true, then fluxes are multiplied by thermal conductivity projected onto the normal direction. This can be used in conjunction with other use_*

Optional Parameters

  • absolute_value_vector_tagsThe tags for the vectors this residual object should fill with the absolute value of the residual contribution

    C++ Type:std::vector<TagName>

    Controllable:No

    Description:The tags for the vectors this residual object should fill with the absolute value of the residual contribution

  • extra_matrix_tagsThe extra tags for the matrices this Kernel should fill

    C++ Type:std::vector<TagName>

    Controllable:No

    Description:The extra tags for the matrices this Kernel should fill

  • extra_vector_tagsThe extra tags for the vectors this Kernel should fill

    C++ Type:std::vector<TagName>

    Controllable:No

    Description:The extra tags for the vectors this Kernel should fill

  • matrix_onlyFalseWhether this object is only doing assembly to matrices (no vectors)

    Default:False

    C++ Type:bool

    Controllable:No

    Description:Whether this object is only doing assembly to matrices (no vectors)

  • matrix_tagssystemThe tag for the matrices this Kernel should fill

    Default:system

    C++ Type:MultiMooseEnum

    Options:nontime, system

    Controllable:No

    Description:The tag for the matrices this Kernel should fill

  • vector_tagsnontimeThe tag for the vectors this Kernel should fill

    Default:nontime

    C++ Type:MultiMooseEnum

    Options:nontime, time

    Controllable:No

    Description:The tag for the vectors this Kernel should fill

Contribution To Tagged Field Data Parameters

  • control_tagsAdds user-defined labels for accessing object parameters via control logic.

    C++ Type:std::vector<std::string>

    Controllable:No

    Description:Adds user-defined labels for accessing object parameters via control logic.

  • diag_save_inThe name of auxiliary variables to save this BC's diagonal jacobian contributions to. Everything about that variable must match everything about this variable (the type, what blocks it's on, etc.)

    C++ Type:std::vector<AuxVariableName>

    Unit:(no unit assumed)

    Controllable:No

    Description:The name of auxiliary variables to save this BC's diagonal jacobian contributions to. Everything about that variable must match everything about this variable (the type, what blocks it's on, etc.)

  • enableTrueSet the enabled status of the MooseObject.

    Default:True

    C++ Type:bool

    Controllable:Yes

    Description:Set the enabled status of the MooseObject.

  • implicitTrueDetermines whether this object is calculated using an implicit or explicit form

    Default:True

    C++ Type:bool

    Controllable:No

    Description:Determines whether this object is calculated using an implicit or explicit form

  • save_inThe name of auxiliary variables to save this BC's residual contributions to. Everything about that variable must match everything about this variable (the type, what blocks it's on, etc.)

    C++ Type:std::vector<AuxVariableName>

    Unit:(no unit assumed)

    Controllable:No

    Description:The name of auxiliary variables to save this BC's residual contributions to. Everything about that variable must match everything about this variable (the type, what blocks it's on, etc.)

  • search_methodnearest_node_connected_sidesChoice of search algorithm. All options begin by finding the nearest node in the primary boundary to a query point in the secondary boundary. In the default nearest_node_connected_sides algorithm, primary boundary elements are searched iff that nearest node is one of their nodes. This is fast to determine via a pregenerated node-to-elem map and is robust on conforming meshes. In the optional all_proximate_sides algorithm, primary boundary elements are searched iff they touch that nearest node, even if they are not topologically connected to it. This is more CPU-intensive but is necessary for robustness on any boundary surfaces which has disconnections (such as Flex IGA meshes) or non-conformity (such as hanging nodes in adaptively h-refined meshes).

    Default:nearest_node_connected_sides

    C++ Type:MooseEnum

    Options:nearest_node_connected_sides, all_proximate_sides

    Controllable:No

    Description:Choice of search algorithm. All options begin by finding the nearest node in the primary boundary to a query point in the secondary boundary. In the default nearest_node_connected_sides algorithm, primary boundary elements are searched iff that nearest node is one of their nodes. This is fast to determine via a pregenerated node-to-elem map and is robust on conforming meshes. In the optional all_proximate_sides algorithm, primary boundary elements are searched iff they touch that nearest node, even if they are not topologically connected to it. This is more CPU-intensive but is necessary for robustness on any boundary surfaces which has disconnections (such as Flex IGA meshes) or non-conformity (such as hanging nodes in adaptively h-refined meshes).

  • seed0The seed for the master random number generator

    Default:0

    C++ Type:unsigned int

    Controllable:No

    Description:The seed for the master random number generator

  • skip_execution_outside_variable_domainFalseWhether to skip execution of this boundary condition when the variable it applies to is not defined on the boundary. This can facilitate setups with moving variable domains and fixed boundaries. Note that the FEProblem boundary-restricted integrity checks will also need to be turned off if using this option

    Default:False

    C++ Type:bool

    Controllable:No

    Description:Whether to skip execution of this boundary condition when the variable it applies to is not defined on the boundary. This can facilitate setups with moving variable domains and fixed boundaries. Note that the FEProblem boundary-restricted integrity checks will also need to be turned off if using this option

  • use_displaced_meshFalseWhether or not this object should use the displaced mesh for computation. Note that in the case this is true but no displacements are provided in the Mesh block the undisplaced mesh will still be used.

    Default:False

    C++ Type:bool

    Controllable:No

    Description:Whether or not this object should use the displaced mesh for computation. Note that in the case this is true but no displacements are provided in the Mesh block the undisplaced mesh will still be used.

Advanced Parameters

  • prop_getter_suffixAn optional suffix parameter that can be appended to any attempt to retrieve/get material properties. The suffix will be prepended with a '_' character.

    C++ Type:MaterialPropertyName

    Unit:(no unit assumed)

    Controllable:No

    Description:An optional suffix parameter that can be appended to any attempt to retrieve/get material properties. The suffix will be prepended with a '_' character.

  • use_interpolated_stateFalseFor the old and older state use projected material properties interpolated at the quadrature points. To set up projection use the ProjectedStatefulMaterialStorageAction.

    Default:False

    C++ Type:bool

    Controllable:No

    Description:For the old and older state use projected material properties interpolated at the quadrature points. To set up projection use the ProjectedStatefulMaterialStorageAction.

Material Property Retrieval Parameters

Input Files