Infiltration and drainage test descriptions

The 2-phase analytic infiltration solution

The physical setup studied in this section is a 1D column that is initially unsaturated, and which is subject to a constant injection of fluid from its top. This is of physical importance because it is a model of constant rainfall recharge to an initially dry groundwater system. The top surface becomes saturated, and this saturated zone moves downwards into the column, diffusing as it goes. The problem is of computational interest because under certain conditions an analytic solution is available for the saturation profile as a function of depth and time.

The Richards' equation for an incompressible fluid in one spatial dimension () reads where and Here which is the capillary pressure, and recall that .

The analytic solution of this nonlinear diffusion-advection relevant to constant infiltration to groundwater has been derived by Broadbridge and White (Broadbridge and White, 1988) for certain functions and . Broadbridge and White assume the hydraulic conductivity is (1) where and the parameters obey , , and . The diffusivity is of the form . This leads to very complicated relationships between the capillary pressure, , and the saturation, except in the case where is small, when they are related through with being the final parameter introduced by Broadbridge and White.

Broadbridge and White derive time-dependent solutions for constant recharge to one end of a semi-infinite line. Their solutions are quite lengthy, will not be written here. To compare with MOOSE, the following parameters are used — the hydraulic parameters are those used in Figure 3 of Broadbridge and White:

Table 1: Parameter values used in the 2-phase tests

ParameterValue
Bar length20m
Bar porosity0.25
Bar permeability1
Gravity0.1m.s
Fluid density10kg.m
Fluid viscosity4Pa.s
0m.s
1m.s
0m.s
1m.s
1.5
2Pa
Recharge rate 0.5

Broadbridge and white consider the case where the initial condition is , but this yields , which is impossible to use in a MOOSE model. Therefore the initial condition Pa is used which avoids any underflow problems. The recharge rate of corresponds in the MOOSE model to a recharge rate of kg.m.s. Note that m.s, so that the and may be encoded as and in the relative permeability function Eq. (1) in a straightforward way.

Figure 1 shows good agreement between the analytic solution of Broadbridge and White and the MOOSE implementation. There are minor discrepancies for small values of saturation: these get smaller as the temporal and spatial resolution is increased, but never totally disappear due to the initial condition of Pa.

Figure 1: Comparison of the Broadbridge and White analytical solution with the MOOSE solution for 3 times. This figure is shown in the standard format used in the Broadbridge-White paper: the constant recharge is applied to the top (where depth is zero) and gravity acts downwards in this figure.

Two tests are part of the automatic test suite (one is marked "heavy" because it is a high-resolution version).

[Mesh]
  type = GeneratedMesh
  dim = 2
  nx = 400
  ny = 1
  xmin = -10
  xmax = 10
  ymin = 0
  ymax = 0.05
[]

[GlobalParams]
  PorousFlowDictator = dictator
[]

[Functions]
  [dts]
    type = PiecewiseLinear
    y = '1E-5 1E-2 1E-2 1E-1'
    x = '0 1E-5 1 10'
  []
[]

[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = pressure
    number_fluid_phases = 1
    number_fluid_components = 1
  []
  [pc]
    type = PorousFlowCapillaryPressureBW
    Sn = 0.0
    Ss = 1.0
    C = 1.5
    las = 2
  []
[]

[FluidProperties]
  [simple_fluid]
    type = SimpleFluidProperties
    bulk_modulus = 2e9
    viscosity = 4
    density0 = 10
    thermal_expansion = 0
  []
[]

[Materials]
  [massfrac]
    type = PorousFlowMassFraction
  []
  [temperature]
    type = PorousFlowTemperature
  []
  [simple_fluid]
    type = PorousFlowSingleComponentFluid
    fp = simple_fluid
    phase = 0
  []
  [ppss]
    type = PorousFlow1PhaseP
    porepressure = pressure
    capillary_pressure = pc
  []
  [relperm]
    type = PorousFlowRelativePermeabilityBW
    Sn = 0.0
    Ss = 1.0
    Kn = 0
    Ks = 1
    C = 1.5
    phase = 0
  []
  [porosity]
    type = PorousFlowPorosityConst
    porosity = 0.25
  []
  [permeability]
    type = PorousFlowPermeabilityConst
    permeability = '1 0 0  0 1 0  0 0 1'
  []
[]

[Variables]
  [pressure]
    initial_condition = -9E2
  []
[]

[Kernels]
  [mass0]
    type = PorousFlowMassTimeDerivative
    fluid_component = 0
    variable = pressure
  []
  [flux0]
    type = PorousFlowAdvectiveFlux
    fluid_component = 0
    variable = pressure
    gravity = '-0.1 0 0'
  []
[]

[AuxVariables]
  [SWater]
    family = MONOMIAL
    order = CONSTANT
  []
[]

[AuxKernels]
  [SWater]
    type = MaterialStdVectorAux
    property = PorousFlow_saturation_qp
    index = 0
    variable = SWater
  []
[]

[BCs]
  [recharge]
    type = PorousFlowSink
    variable = pressure
    boundary = right
    flux_function = -1.25 # corresponds to Rstar being 0.5 because i have to multiply by density*porosity
  []
[]

[Preconditioning]
  [andy]
    type = SMP
    full = true
    petsc_options = '-snes_converged_reason -ksp_diagonal_scale -ksp_diagonal_scale_fix -ksp_gmres_modifiedgramschmidt -snes_linesearch_monitor'
    petsc_options_iname = '-ksp_type -pc_type -sub_pc_type -sub_pc_factor_shift_type -pc_asm_overlap -snes_atol -snes_rtol -snes_max_it'
    petsc_options_value = 'gmres      asm      lu           NONZERO                   2               1E-10      1E-10      10000'
  []
[]

[VectorPostprocessors]
  [swater]
    type = LineValueSampler
    warn_discontinuous_face_values = false
    variable = SWater
    start_point = '-10 0 0'
    end_point = '10 0 0'
    sort_by = x
    num_points = 101
    execute_on = timestep_end
  []
[]

[Executioner]
  type = Transient
  solve_type = Newton
  petsc_options = '-snes_converged_reason'
  end_time = 8

  [TimeStepper]
    type = FunctionDT
    function = dts
  []
[]

[Outputs]
  file_base = bw01
  sync_times = '0.5 2 8'
  [exodus]
    type = Exodus
    sync_only = true
  []
  [along_line]
    type = CSV
    sync_only = true
  []
[]
(modules/porous_flow/test/tests/infiltration_and_drainage/bw01.i)

The two-phase analytic drainage solution

Warrick, Lomen and Islas (Warrick et al., 1990) extended the analysis of Broadbridge and White to include the case of drainage from a medium.

The setup is an initially-saturated semi-infinite column of material that drains freely from its lower end. This is simulated by placing a boundary condition of at the lower end. To obtain their analytical solutions, Warrick, Lomen and Islas make the same assumptions as Broadbridge and White concerning the diffusivity and conductivity of the medium. Their solutions are quite lengthy, so are not written here

A MOOSE model with the parameters almost identical to those listed in Table 1 is compared with the analytical solutions. The only differences are that the "bar" length is m (to avoid any interference from the lower Dirichlet boundary condition), and since there is no recharge. The initial condition is Pa: the choice leads to poor convergence since by construction the Broadbridge-White capillary function is only designed to simulate the unsaturated zone and a sensible extension to is discontinuous at .

Figure 2 shows good agreement between the analytic solution and the MOOSE implementation. Any minor discrepancies get smaller as the temporal and spatial resolution increase.

Figure 2: Comparison of the Warrick, Lomen and Islas analytical solution with the MOOSE solution for 3 times. This figure is shown in the standard format used in the literature: the top of the model is at the top of the figure, and gravity acts downwards in this figure, with fluid draining from the infinitely deep point.

Two tests are part of the automatic test suite (one is marked "heavy" because it is a high-resolution version).

[Mesh]
  type = GeneratedMesh
  dim = 2
  nx = 1000
  ny = 1
  xmin = -10000
  xmax = 0
  ymin = 0
  ymax = 0.05
[]

[GlobalParams]
  PorousFlowDictator = dictator
[]

[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = pressure
    number_fluid_phases = 1
    number_fluid_components = 1
  []
  [pc]
    type = PorousFlowCapillaryPressureBW
    Sn = 0.0
    Ss = 1.0
    C = 1.5
    las = 2
  []
[]

[FluidProperties]
  [simple_fluid]
    type = SimpleFluidProperties
    bulk_modulus = 2e9
    viscosity = 4
    density0 = 10
    thermal_expansion = 0
  []
[]

[Materials]
  [massfrac]
    type = PorousFlowMassFraction
  []
  [temperature]
    type = PorousFlowTemperature
  []
  [simple_fluid]
    type = PorousFlowSingleComponentFluid
    fp = simple_fluid
    phase = 0
  []
  [ppss]
    type = PorousFlow1PhaseP
    porepressure = pressure
    capillary_pressure = pc
  []
  [relperm]
    type = PorousFlowRelativePermeabilityBW
    Sn = 0.0
    Ss = 1.0
    Kn = 0
    Ks = 1
    C = 1.5
    phase = 0
  []
  [porosity]
    type = PorousFlowPorosityConst
    porosity = 0.25
  []
  [permeability]
    type = PorousFlowPermeabilityConst
    permeability = '1 0 0  0 1 0  0 0 1'
  []
[]

[Variables]
  [pressure]
    initial_condition = -1E-4
  []
[]

[Kernels]
  [mass0]
    type = PorousFlowMassTimeDerivative
    fluid_component = 0
    variable = pressure
  []
  [flux0]
    type = PorousFlowAdvectiveFlux
    fluid_component = 0
    variable = pressure
    gravity = '-0.1 0 0'
  []
[]

[AuxVariables]
  [SWater]
    family = MONOMIAL
    order = CONSTANT
  []
[]

[AuxKernels]
  [SWater]
    type = MaterialStdVectorAux
    property = PorousFlow_saturation_qp
    index = 0
    variable = SWater
  []
[]

[BCs]
  [base]
    type = DirichletBC
    boundary = 'left'
    value = -1E-4
    variable = pressure
  []
[]

[Preconditioning]
  [andy]
    type = SMP
    full = true
    petsc_options = '-snes_converged_reason -ksp_diagonal_scale -ksp_diagonal_scale_fix -ksp_gmres_modifiedgramschmidt -snes_linesearch_monitor'
    petsc_options_iname = '-ksp_type -pc_type -sub_pc_type -sub_pc_factor_shift_type -pc_asm_overlap -snes_atol -snes_rtol -snes_max_it'
    petsc_options_value = 'gmres      asm      lu           NONZERO                   2               1E-10      1E-10      10000'
  []
[]

[VectorPostprocessors]
  [swater]
    type = LineValueSampler
    warn_discontinuous_face_values = false
    variable = SWater
    start_point = '-5000 0 0'
    end_point = '0 0 0'
    sort_by = x
    num_points = 71
    execute_on = timestep_end
  []
[]

[Executioner]
  type = Transient
  solve_type = Newton
  petsc_options = '-snes_converged_reason'
  end_time = 1000
  dt = 1
[]

[Outputs]
  file_base = wli01
  sync_times = '100 500 1000'
  [exodus]
    type = Exodus
    sync_only = true
  []
  [along_line]
    type = CSV
    sync_only = true
  []
[]
(modules/porous_flow/test/tests/infiltration_and_drainage/wli01.i)

Single-phase infiltration and drainage

Forsyth, Wu and Pruess (Forsyth et al., 1995) describe a HYDRUS simulation of an experiment involving infiltration (experiment 1) and subsequent drainage (experiment 2) in a large caisson. The simulation is effectively one dimensional, and is shown in Figure 3.

Figure 3: Two experimental setups from Forsyth, Wu and Pruess. Experiment 1 involves infiltration of water into an initially unsaturated caisson. Experiment 2 involves drainage of water from an initially saturated caisson.

The properties common to each experiment are listed in Table 2

Table 2: Parameter values used in the single-phase infiltration and drainage tests

ParameterValue
Caisson porosity0.33
Caisson permeabilitym
Gravity10m.s
Water density at STP1000kg.m
Water viscosity0.00101Pa.s
Water bulk modulus20MPa
Water residual saturation0.0
Air residual saturation0.0
Air pressure0.0
van Genuchten Pa
van Genuchten 0.336
van Genuchten turnover0.99

In each experiment 120 finite elements are used along the length of the Caisson. The modified van Genuchten relative permeability curve with a "turnover" (set at ) is employed in order to improve convergence significantly. Hydrus also uses a modified van Genuchten curve, although I couldn't find any details on the modification.

In experiment 1, the caisson is initially at saturation 0.303 (Pa), and water is pumped into the top with a rate 0.002315kg.m.s. This causes a front of water to advance down the caisson. Figure 4 shows the agreement between MOOSE and the published result (this result was obtained by extracting data by hand from online graphics).

In experiment 2, the caisson is initially fully saturated at , and the bottom is held at to cause water to drain via the action of gravity. Figure 4 and Figure 5 show the agreement between MOOSE and the published result.

Figure 4: Saturation profile in the caisson after 4.16 days of infiltration. Note that the HYDRUS results are only approximate they were extracted by hand from online graphics.

Figure 5: Saturation profiles in the caisson after drainage from an initially-saturated simulation (4 days and 100 days profiles). Note that the HYDRUS results are only approximate they were extracted by hand from online graphics.

Experiment 1 and the first 4 simulation days of experiment 2 are marked as "heavy" in the PorousFlow test suite since the simulations take around 3 seconds to complete.

The input file for Experiment 1:

[Mesh]
  type = GeneratedMesh
  dim = 2
  nx = 120
  ny = 1
  xmin = 0
  xmax = 6
  ymin = 0
  ymax = 0.05
[]

[GlobalParams]
  PorousFlowDictator = dictator
[]

[Functions]
  [dts]
    type = PiecewiseLinear
    y = '1E-2 1 10 500 5000 5000'
    x = '0 10 100 1000 10000 100000'
  []
[]

[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = pressure
    number_fluid_phases = 1
    number_fluid_components = 1
  []
  [pc]
    type = PorousFlowCapillaryPressureVG
    m = 0.336
    alpha = 1.43e-4
  []
[]

[FluidProperties]
  [simple_fluid]
    type = SimpleFluidProperties
    bulk_modulus = 2e7
    viscosity = 1.01e-3
    density0 = 1000
    thermal_expansion = 0
  []
[]

[Materials]
  [massfrac]
    type = PorousFlowMassFraction
  []
  [temperature]
    type = PorousFlowTemperature
  []
  [simple_fluid]
    type = PorousFlowSingleComponentFluid
    fp = simple_fluid
    phase = 0
  []
  [ppss]
    type = PorousFlow1PhaseP
    porepressure = pressure
    capillary_pressure = pc
  []
  [relperm]
    type = PorousFlowRelativePermeabilityVG
    m = 0.336
    seff_turnover = 0.99
    phase = 0
  []
  [porosity]
    type = PorousFlowPorosityConst
    porosity = 0.33
  []
  [permeability]
    type = PorousFlowPermeabilityConst
    permeability = '0.295E-12 0 0  0 0.295E-12 0  0 0 0.295E-12'
  []
[]

[Variables]
  [pressure]
    initial_condition = -72620.4
  []
[]

[Kernels]
  [mass0]
    type = PorousFlowMassTimeDerivative
    fluid_component = 0
    variable = pressure
  []
  [flux0]
    type = PorousFlowAdvectiveFlux
    fluid_component = 0
    variable = pressure
    gravity = '-10 0 0'
  []
[]

[AuxVariables]
  [SWater]
    family = MONOMIAL
    order = CONSTANT
  []
[]

[AuxKernels]
  [SWater]
    type = MaterialStdVectorAux
    property = PorousFlow_saturation_qp
    index = 0
    variable = SWater
  []
[]

[BCs]
  [base]
    type = PorousFlowSink
    boundary = right
    flux_function = -2.315E-3
    variable = pressure
  []
[]

[Preconditioning]
  [andy]
    type = SMP
    full = true
    petsc_options = '-snes_converged_reason -ksp_diagonal_scale -ksp_diagonal_scale_fix -ksp_gmres_modifiedgramschmidt -snes_linesearch_monitor'
    petsc_options_iname = '-ksp_type -pc_type -sub_pc_type -sub_pc_factor_shift_type -pc_asm_overlap -snes_atol -snes_rtol -snes_max_it'
    petsc_options_value = 'gmres      asm      lu           NONZERO                   2               1E-10      1E-10      10'
  []
[]

[VectorPostprocessors]
  [swater]
    type = LineValueSampler
    warn_discontinuous_face_values = false
    variable = SWater
    start_point = '0 0 0'
    end_point = '6 0 0'
    sort_by = x
    num_points = 121
    execute_on = timestep_end
  []
[]

[Executioner]
  type = Transient
  solve_type = Newton
  petsc_options = '-snes_converged_reason'
  end_time = 359424

  [TimeStepper]
    type = FunctionDT
    function = dts
  []
[]

[Outputs]
  file_base = rd01
  [exodus]
    type = Exodus
    execute_on = 'initial final'
  []
  [along_line]
    type = CSV
    execute_on = final
  []
[]
(modules/porous_flow/test/tests/infiltration_and_drainage/rd01.i)

The input file for the first 4 simulation days of Experiment 2:

[Mesh]
  type = GeneratedMesh
  dim = 2
  nx = 120
  ny = 1
  xmin = 0
  xmax = 6
  ymin = 0
  ymax = 0.05
[]

[GlobalParams]
  PorousFlowDictator = dictator
[]

[Functions]
  [dts]
    type = PiecewiseLinear
    y = '1E-2 1 10 500 5000 50000'
    x = '0 10 100 1000 10000 500000'
  []
[]

[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = pressure
    number_fluid_phases = 1
    number_fluid_components = 1
  []
  [pc]
    type = PorousFlowCapillaryPressureVG
    m = 0.336
    alpha = 1.43e-4
  []
[]

[FluidProperties]
  [simple_fluid]
    type = SimpleFluidProperties
    bulk_modulus = 2e7
    viscosity = 1.01e-3
    density0 = 1000
    thermal_expansion = 0
  []
[]

[Materials]
  [massfrac]
    type = PorousFlowMassFraction
  []
  [temperature]
    type = PorousFlowTemperature
  []
  [simple_fluid]
    type = PorousFlowSingleComponentFluid
    fp = simple_fluid
    phase = 0
  []
  [ppss]
    type = PorousFlow1PhaseP
    porepressure = pressure
    capillary_pressure = pc
  []
  [relperm]
    type = PorousFlowRelativePermeabilityVG
    m = 0.336
    seff_turnover = 0.99
    phase = 0
  []
  [porosity]
    type = PorousFlowPorosityConst
    porosity = 0.33
  []
  [permeability]
    type = PorousFlowPermeabilityConst
    permeability = '0.295E-12 0 0  0 0.295E-12 0  0 0 0.295E-12'
  []
[]

[Variables]
  [pressure]
    initial_condition = 0.0
  []
[]

[Kernels]
  [mass0]
    type = PorousFlowMassTimeDerivative
    fluid_component = 0
    variable = pressure
  []
  [flux0]
    type = PorousFlowAdvectiveFlux
    fluid_component = 0
    variable = pressure
    gravity = '-10 0 0'
  []
[]

[AuxVariables]
  [SWater]
    family = MONOMIAL
    order = CONSTANT
  []
[]

[AuxKernels]
  [SWater]
    type = MaterialStdVectorAux
    property = PorousFlow_saturation_qp
    index = 0
    variable = SWater
  []
[]

[BCs]
  [base]
    type = DirichletBC
    boundary = left
    value = 0.0
    variable = pressure
  []
[]

[Preconditioning]
  [andy]
    type = SMP
    full = true
    petsc_options = '-snes_converged_reason -ksp_diagonal_scale -ksp_diagonal_scale_fix -ksp_gmres_modifiedgramschmidt -snes_linesearch_monitor'
    petsc_options_iname = '-ksp_type -pc_type -sub_pc_type -sub_pc_factor_shift_type -pc_asm_overlap -snes_atol -snes_rtol -snes_max_it'
    petsc_options_value = 'gmres      asm      lu           NONZERO                   2               1E-10      1E-10      10'
  []
[]

[VectorPostprocessors]
  [swater]
    type = LineValueSampler
    warn_discontinuous_face_values = false
    variable = SWater
    start_point = '0 0 0'
    end_point = '6 0 0'
    sort_by = x
    num_points = 121
    execute_on = timestep_end
  []
[]

[Executioner]
  type = Transient
  solve_type = Newton
  petsc_options = '-snes_converged_reason'
  end_time = 345600

  [TimeStepper]
    type = FunctionDT
    function = dts
  []
[]

[Outputs]
  file_base = rd02
  [exodus]
    type = Exodus
    execute_on = 'initial final'
  []
  [along_line]
    type = CSV
    execute_on = final
  []
[]
(modules/porous_flow/test/tests/infiltration_and_drainage/rd02.i)

The input file for the latter 96 days of Experiment 2:

[Mesh]
  file = gold/rd02.e
[]

[GlobalParams]
  PorousFlowDictator = dictator
[]

[Functions]
  [dts]
    type = PiecewiseLinear
    y = '2E4 1E6'
    x = '0 1E6'
  []
[]

[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = pressure
    number_fluid_phases = 1
    number_fluid_components = 1
  []
  [pc]
    type = PorousFlowCapillaryPressureVG
    m = 0.336
    alpha = 1.43e-4
  []
[]

[FluidProperties]
  [simple_fluid]
    type = SimpleFluidProperties
    bulk_modulus = 2e7
    viscosity = 1.01e-3
    density0 = 1000
    thermal_expansion = 0
  []
[]

[Materials]
  [massfrac]
    type = PorousFlowMassFraction
  []
  [temperature]
    type = PorousFlowTemperature
  []
  [simple_fluid]
    type = PorousFlowSingleComponentFluid
    fp = simple_fluid
    phase = 0
  []
  [ppss]
    type = PorousFlow1PhaseP
    porepressure = pressure
    capillary_pressure = pc
  []
  [relperm]
    type = PorousFlowRelativePermeabilityVG
    m = 0.336
    seff_turnover = 0.99
    phase = 0
  []
  [porosity]
    type = PorousFlowPorosityConst
    porosity = 0.33
  []
  [permeability]
    type = PorousFlowPermeabilityConst
    permeability = '0.295E-12 0 0  0 0.295E-12 0  0 0 0.295E-12'
  []
[]

[Variables]
  [pressure]
    initial_from_file_timestep = LATEST
    initial_from_file_var = pressure
  []
[]

[Kernels]
  [mass0]
    type = PorousFlowMassTimeDerivative
    fluid_component = 0
    variable = pressure
  []
  [flux0]
    type = PorousFlowAdvectiveFlux
    fluid_component = 0
    variable = pressure
    gravity = '-10 0 0'
  []
[]

[AuxVariables]
  [SWater]
    family = MONOMIAL
    order = CONSTANT
  []
[]

[AuxKernels]
  [SWater]
    type = MaterialStdVectorAux
    property = PorousFlow_saturation_qp
    index = 0
    variable = SWater
  []
[]

[BCs]
  [base]
    type = DirichletBC
    boundary = left
    value = 0.0
    variable = pressure
  []
[]

[Preconditioning]
  [andy]
    type = SMP
    full = true
    petsc_options = '-snes_converged_reason -ksp_diagonal_scale -ksp_diagonal_scale_fix -ksp_gmres_modifiedgramschmidt -snes_linesearch_monitor'
    petsc_options_iname = '-ksp_type -pc_type -sub_pc_type -sub_pc_factor_shift_type -pc_asm_overlap -snes_atol -snes_rtol -snes_max_it'
    petsc_options_value = 'gmres      asm      lu           NONZERO                   2               1E-10      1E-10      10'
  []
[]

[VectorPostprocessors]
  [swater]
    type = LineValueSampler
    warn_discontinuous_face_values = false
    variable = SWater
    start_point = '0 0 0'
    end_point = '6 0 0'
    sort_by = x
    num_points = 121
    execute_on = timestep_end
  []
[]

[Executioner]
  type = Transient
  solve_type = Newton
  petsc_options = '-snes_converged_reason'
  end_time = 8.2944E6

  [TimeStepper]
    type = FunctionDT
    function = dts
  []
[]

[Outputs]
  file_base = rd03
  [exodus]
    type = Exodus
    execute_on = 'initial final'
  []
  [along_line]
    type = CSV
    execute_on = final
  []
[]
(modules/porous_flow/test/tests/infiltration_and_drainage/rd03.i)

Water infiltration into a two-phase (oil-water) system

An analytic solution of the two-phase Richards' equations with gravity on a semi-infinite line , with a constant water infiltration flux at has been derived by (Rogers et al., 1983) (Unfortunately there must be a typo in the RSC paper as for nonzero gravity their results are clearly incorrect.). The authors assume incompressible fluids; linear relative permeability relationships; the "oil" (or "gas") viscosity is larger than the water viscosity; and, a certain functional form for the capillary pressure. When the oil viscosity is exactly twice the water viscosity, their effective saturation reads (2) where is the capillary pressure, and and are arbitrary parameters to be defined by the user in the PorousFlow implementation. For other oil/water viscosity ratios is more complicated, and note that their formulation allows , but only the particular form Eq. (2) need be used to validate the MOOSE implementation.

RSC's solutions are quite lengthy, so I will not write them here. To compare with MOOSE, the following parameters are used:

Table 3: Parameter values used in the Rogers-Stallybrass-Clements infiltration tests

ParameterValue
Bar length10 m
Bar porosity0.25
Bar permeabilitym
Gravity0 m.s
Water density10 kg.m
Water viscosityPa.s
Oil density20 kg.m
Oil viscosityPa.s
Capillary 10 Pa
Capillary 1 Pa
Initial water pressure0 Pa
Initial oil pressure15 Pa
Initial water saturation0.08181
Initial oil saturation0.91819
Water injection rate1 kg.s.m

The input file:

# RSC test with high-res time and spatial resolution
[Mesh]
  type = GeneratedMesh
  dim = 2
  nx = 600
  ny = 1
  xmin = 0
  xmax = 10 # x is the depth variable, called zeta in RSC
  ymin = 0
  ymax = 0.05
[]

[GlobalParams]
  PorousFlowDictator = dictator
  gravity = '0 0 0'
[]

[Functions]
  [dts]
    type = PiecewiseLinear
    y = '3E-3 3E-2 0.05'
    x = '0 1 5'
  []
[]

[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'pwater poil'
    number_fluid_phases = 2
    number_fluid_components = 2
  []
  [pc]
    type = PorousFlowCapillaryPressureRSC
    oil_viscosity = 2E-3
    scale_ratio = 2E3
    shift = 10
  []
[]

[FluidProperties]
  [water]
    type = SimpleFluidProperties
    bulk_modulus = 2e9
    density0 = 10
    thermal_expansion = 0
    viscosity = 1e-3
  []
  [oil]
    type = SimpleFluidProperties
    bulk_modulus = 2e9
    density0 = 20
    thermal_expansion = 0
    viscosity = 2e-3
  []
[]

[Materials]
  [temperature]
    type = PorousFlowTemperature
  []
  [ppss]
    type = PorousFlow2PhasePP
    phase0_porepressure = pwater
    phase1_porepressure = poil
    capillary_pressure = pc
  []
  [massfrac]
    type = PorousFlowMassFraction
    mass_fraction_vars = 'massfrac_ph0_sp0 massfrac_ph1_sp0'
  []
  [water]
    type = PorousFlowSingleComponentFluid
    fp = water
    phase = 0
    compute_enthalpy = false
    compute_internal_energy = false
  []
  [oil]
    type = PorousFlowSingleComponentFluid
    fp = oil
    phase = 1
    compute_enthalpy = false
    compute_internal_energy = false
  []
  [relperm_water]
    type = PorousFlowRelativePermeabilityCorey
    n = 1
    phase = 0
  []
  [relperm_oil]
    type = PorousFlowRelativePermeabilityCorey
    n = 1
    phase = 1
  []
  [porosity]
    type = PorousFlowPorosityConst
    porosity = 0.25
  []
  [permeability]
    type = PorousFlowPermeabilityConst
    permeability = '1E-5 0 0  0 1E-5 0  0 0 1E-5'
  []
[]

[Variables]
  [pwater]
  []
  [poil]
  []
[]

[ICs]
  [water_init]
    type = ConstantIC
    variable = pwater
    value = 0
  []
  [oil_init]
    type = ConstantIC
    variable = poil
    value = 15
  []
[]

[Kernels]
  [mass0]
    type = PorousFlowMassTimeDerivative
    fluid_component = 0
    variable = pwater
  []
  [flux0]
    type = PorousFlowAdvectiveFlux
    fluid_component = 0
    variable = pwater
  []
  [mass1]
    type = PorousFlowMassTimeDerivative
    fluid_component = 1
    variable = poil
  []
  [flux1]
    type = PorousFlowAdvectiveFlux
    fluid_component = 1
    variable = poil
  []
[]

[AuxVariables]
  [SWater]
    family = MONOMIAL
    order = CONSTANT
  []
  [SOil]
    family = MONOMIAL
    order = CONSTANT
  []
  [massfrac_ph0_sp0]
    initial_condition = 1
  []
  [massfrac_ph1_sp0]
    initial_condition = 0
  []
[]

[AuxKernels]
  [SWater]
    type = MaterialStdVectorAux
    property = PorousFlow_saturation_qp
    index = 0
    variable = SWater
  []
  [SOil]
    type = MaterialStdVectorAux
    property = PorousFlow_saturation_qp
    index = 1
    variable = SOil
  []
[]

[BCs]
# we are pumping water into a system that has virtually incompressible fluids, hence the pressures rise enormously.  this adversely affects convergence because of almost-overflows and precision-loss problems.  The fixed things help keep pressures low and so prevent these awful behaviours.   the movement of the saturation front is the same regardless of the fixed things.
  active = 'recharge fixedoil fixedwater'
  [recharge]
    type = PorousFlowSink
    variable = pwater
    boundary = 'left'
    flux_function = -1.0
  []
  [fixedwater]
    type = DirichletBC
    variable = pwater
    boundary = 'right'
    value = 0
  []
  [fixedoil]
    type = DirichletBC
    variable = poil
    boundary = 'right'
    value = 15
  []
[]

[Preconditioning]
  [andy]
    type = SMP
    full = true
    petsc_options = '-snes_converged_reason -ksp_diagonal_scale -ksp_diagonal_scale_fix -ksp_gmres_modifiedgramschmidt -snes_linesearch_monitor'
    petsc_options_iname = '-ksp_type -pc_type -sub_pc_type -sub_pc_factor_shift_type -pc_asm_overlap -snes_atol -snes_rtol -snes_max_it'
    petsc_options_value = 'gmres      asm      lu           NONZERO                   2               1E-10      1E-10      10000'
  []
[]

[VectorPostprocessors]
  [swater]
    type = LineValueSampler
    warn_discontinuous_face_values = false
    variable = SWater
    start_point = '0 0 0'
    end_point = '7 0 0'
    sort_by = x
    num_points = 21
    execute_on = timestep_end
  []
[]

[Executioner]
  type = Transient
  solve_type = Newton
  petsc_options = '-snes_converged_reason'
  end_time = 5
  [TimeStepper]
    type = FunctionDT
    function = dts
  []
[]

[Outputs]
  file_base = rsc01
  [along_line]
    type = CSV
    execute_vector_postprocessors_on = final
  []
  [exodus]
    type = Exodus
    execute_on = 'initial final'
  []
[]
(modules/porous_flow/test/tests/infiltration_and_drainage/rsc01.i)

In the RSC theory water is injected into a semi-infinite domain, whereas of course the MOOSE implementation has finite extent ( is chosen). Because of the near incompressibility of the fluids (I choose the bulk modulus to be 2 GPa) this causes the porepressures to rise enormously, and the problem can suffer from precision-loss problems. Therefore, the porepressures are fixed at . This does not affect the progress of the water saturation front.

Figure 6 shows good agreement between the analytic solution and the MOOSE implementation. Any minor discrepancies get smaller as the temporal and spatial resolution increase, as is suggested by the two comparisons in that figure.

The "low-resolution" test has 200 elements in and uses 15 time steps is part the automatic test suite that is run every time the code is updated. The "high-resolution" test has 600 elements and uses 190 time steps, and is marked as "heavy".

Figure 6: Water saturation profile after 5 seconds of injection in the Rogers-Stallybrass-Clements test. The initial water saturation is 0.08181, and water is injected at the top of this figure at a constant rate. This forms a water front which displaces the oil. Black line: RSC's analytic solution. Red squares: high-resolution MOOSE simulation. Green triangles: lower resolution MOOSE simulation.

References

  1. P. Broadbridge and I. White. Constant rate rainfall infiltration: a versatile nonlinear model, 1. analytical solution. Water Resources Research, 24:145–154, 1988.[BibTeX]
  2. P. A. Forsyth, Y. S. Wu, and K. Pruess. Robust numerical methods for saturated-unsaturated flow with dry initial conditions in heterogeneous media. Water Resources Research, 18:25–38, 1995.[BibTeX]
  3. C. Rogers, M. P. Stallybrass, and D. L. Clements. On two phase filtration under gravity and with boundary infiltration: application of a Backlund transformation. Nonlinear Analysis, Theory, Methods and Applications, 7:785–799, 1983.[BibTeX]
  4. A. W. Warrick, D. O. Lomen, and A. Islas. An analytical solution to Richards' Equation for a Draining Soil Profile. Water Resources Research, 26:253–258, 1990.[BibTeX]