1D Transient Flow (Pressure), BCs of 1st Kind

The Problem Description

The PFLOTRAN Input File (GENERAL Mode)

The PFLOTRAN Input File (TH Mode)

The PFLOTRAN Input File (RICHARDS Mode)

The Python Script

The Problem Description

This problem is adapted from Kolditz, et al. (2015), Thermo-Hydro-Mechanical-Chemical Processes in Fractured Porous Media: Modelling and Benchmarking, Closed Form Solutions, Springer International Publishing, Switzerland. Section 2.2.7, pg.32, “A Transient 1D Pressure Distribution, Time-Dependent Boundary Conditions of 1st Kind.”

The domain is a 20x1x1 meter rectangular column extending along the positive and negative x-axis and is made up of 20x1x1 cubic grid cells with dimensions 1x1x1 meters. The domain material is assigned a porosity :math”phi = 0.25, and permeability \(k\) = 1e-14 m2. The fluid has a viscosity \(\mu\) = 1.728 mPa-sec, and a constant compressibility \(K\) = 1e-8 1/Pa.

The pressure is initially uniform at p(t=0) = 0.25 MPa. The boundary pressures are transient:

\[ \begin{align}\begin{aligned}p(-L,t) = p_b t + p_{t=0}\\p(L,t) = p_b t + p_{t=0}\end{aligned}\end{align} \]

where L = 10 m and \(p_b\) = 2 MPa/day. The transient pressure distribution is governed by,

\[\phi K {{\partial p} \over {\partial t}} = {k \over \mu} {{\partial^{2} p} \over {\partial x^{2}}}\]

With the given boundary conditions, the solution is defined by,

\[ \begin{align}\begin{aligned}\chi = {k \over {\phi \mu K}}\\p(x,t) = p_b t + {{p_b(x^2-L^2)}\over{2\chi}} + {{16 p_b L^2}\over{\chi\pi^3}} \sum_{n=0}^{\infty}{{(-1)^n}\over{(2n+1)^3}} cos{\left({{(2n+1)x\pi}\over{2L}}\right)} e^{\left({-\chi(2n+1)^2\pi^2{{t}\over{4L^2}}}\right)}\end{aligned}\end{align} \]
_images/visit_figure6.png

The PFLOTRAN domain set-up.

If you do not see this image, you must run the QA test suite to generate this figure.

Comparison of the PFLOTRAN vs. analytical solution for GENERAL mode.

If you do not see this image, you must run the QA test suite to generate this figure.

Comparison of the PFLOTRAN vs. analytical solution for TH mode.

If you do not see this image, you must run the QA test suite to generate this figure.

Comparison of the PFLOTRAN vs. analytical solution for RICHARDS mode.

The PFLOTRAN Input File (GENERAL Mode)

The GENERAL Mode PFLOTRAN input file can be downloaded here.

 
SIMULATION
  SIMULATION_TYPE SUBSURFACE
  PROCESS_MODELS
    SUBSURFACE_FLOW flow
      MODE GENERAL
    /
  /
END

SUBSURFACE

GRID
  TYPE structured
  NXYZ 10 1 1
  DXYZ
   2.0d0
   1.0d0
   1.0d0
  END
END

REGION all
  COORDINATES
    0.0d0 0.0d0 0.0d0
    20.0d0 1.0d0 1.0d0
  /
END

REGION left_end
  FACE WEST
  COORDINATES
    0.d0 0.d0 0.d0
    0.d0 1.0d0 1.0d0
  /
END

REGION right_end
  FACE EAST
  COORDINATES
    20.0d0 0.d0 0.d0
    20.0d0 1.0d0 1.0d0
  /
END

MATERIAL_PROPERTY beam
  ID 1
  POROSITY 0.25
  TORTUOSITY 1.d0
  ROCK_DENSITY 2.5d3 kg/m^3
  HEAT_CAPACITY 1.5 J/kg-C
  THERMAL_CONDUCTIVITY_DRY 2.0 W/m-C
  THERMAL_CONDUCTIVITY_WET 2.0 W/m-C
  CHARACTERISTIC_CURVES cc1
  PERMEABILITY
    PERM_X 1.d-14
    PERM_Y 1.d-14
    PERM_Z 1.d-14
  /
END

CHARACTERISTIC_CURVES cc1
  SATURATION_FUNCTION VAN_GENUCHTEN
    LIQUID_RESIDUAL_SATURATION 0.5d-1
    M 0.75
    ALPHA 1.d-3
  /
  PERMEABILITY_FUNCTION MUALEM
    LIQUID_RESIDUAL_SATURATION 0.1d0
    M 0.5d0
  /
  PERMEABILITY_FUNCTION MUALEM_VG_GAS
    LIQUID_RESIDUAL_SATURATION 0.1d0
    GAS_RESIDUAL_SATURATION 0.1d0
    M 0.5d0
  /
END

EOS WATER
  #ref. density [kg/m3], ref. pressure [Pa], water compressibility [1/Pa?]
  DENSITY EXPONENTIAL 1000. 0. 1.d-8
  VISCOSITY CONSTANT 1.728d-3 Pa-s
END

FLUID_PROPERTY
  DIFFUSION_COEFFICIENT 1.d-9
END

STRATA
  REGION all
  MATERIAL beam
END

TIME
  FINAL_TIME 0.75 day
  INITIAL_TIMESTEP_SIZE 1.d-4 day
  MAXIMUM_TIMESTEP_SIZE 0.0001 day at 0.d0 day
END

OUTPUT
  SNAPSHOT_FILE
    TIMES day 0.10 0.25 0.50
    #NO_PRINT_INITIAL
    FORMAT HDF5
  /
END

FLOW_CONDITION initial
  TYPE
    TEMPERATURE dirichlet
    LIQUID_PRESSURE dirichlet
    MOLE_FRACTION dirichlet
  /
  TEMPERATURE 25.0d0 C
  LIQUID_PRESSURE 0.25 MPa
  MOLE_FRACTION 1.d-10
END

FLOW_CONDITION left_end
  TYPE
    TEMPERATURE dirichlet
    LIQUID_PRESSURE dirichlet
    MOLE_FRACTION dirichlet
  /
  TEMPERATURE 25.0 C
  LIQUID_PRESSURE LIST
    # p = p_b*t + offset; p_b=2e6 Pa/day; offset = 0.25 MPa
    TIME_UNITS day
    DATA_UNITS Pa
    INTERPOLATION LINEAR
    #time  #pressure
    0.00d0 0.25d6
    0.25d0 0.75d6
    0.50d0 1.25d6
    1.00d0 2.25d6
  /
  MOLE_FRACTION 1.d-10
END

FLOW_CONDITION right_end
  TYPE
    TEMPERATURE dirichlet
    LIQUID_PRESSURE dirichlet
    MOLE_FRACTION dirichlet
  /
  TEMPERATURE 25.0 C
  LIQUID_PRESSURE LIST
    # p = p_b*t + offset; p_b=2e6 Pa/day; offset = 0.25 MPa
    TIME_UNITS day
    DATA_UNITS Pa
    INTERPOLATION LINEAR
    #time  #pressure
    0.00d0 0.25d6
    0.25d0 0.75d6
    0.50d0 1.25d6
    1.00d0 2.25d6
  /
  MOLE_FRACTION 1.d-10
END

INITIAL_CONDITION initial
  REGION all
  FLOW_CONDITION initial
END

BOUNDARY_CONDITION left_end
  REGION left_end
  FLOW_CONDITION left_end
END

BOUNDARY_CONDITION right_end
  FLOW_CONDITION right_end
  REGION right_end
END

END_SUBSURFACE

The PFLOTRAN Input File (TH Mode)

The TH Mode PFLOTRAN input file can be downloaded here.

 
SIMULATION
  SIMULATION_TYPE SUBSURFACE
  PROCESS_MODELS
    SUBSURFACE_FLOW flow
      MODE TH
    /
  /
END

SUBSURFACE

GRID
  TYPE structured
  NXYZ 10 1 1
  DXYZ
   2.0d0
   1.0d0
   1.0d0
  END
END

REGION all
  COORDINATES
    0.0d0 0.0d0 0.0d0
    20.0d0 1.0d0 1.0d0
  /
END

REGION left_end
  FACE WEST
  COORDINATES
    0.d0 0.d0 0.d0
    0.d0 1.0d0 1.0d0
  /
END

REGION right_end
  FACE EAST
  COORDINATES
    20.0d0 0.d0 0.d0
    20.0d0 1.0d0 1.0d0
  /
END

MATERIAL_PROPERTY beam
  ID 1
  POROSITY 0.25
  TORTUOSITY 1.d0
  ROCK_DENSITY 2.5d3 kg/m^3
  SPECIFIC_HEAT 1.5 J/kg-C
  THERMAL_CONDUCTIVITY_DRY 2.0 W/m-C
  THERMAL_CONDUCTIVITY_WET 2.0 W/m-C
  SATURATION_FUNCTION default
  PERMEABILITY
    PERM_X 1.d-14
    PERM_Y 1.d-14
    PERM_Z 1.d-14
  /
END

SATURATION_FUNCTION default
  SATURATION_FUNCTION_TYPE VAN_GENUCHTEN
  RESIDUAL_SATURATION 0.5d-1
  LAMBDA 0.75
  ALPHA 1.d-3
END

EOS WATER
  #ref. density [kg/m3], ref. pressure [Pa], water compressibility [1/Pa?]
  DENSITY EXPONENTIAL 1000. 0. 1.d-8
  VISCOSITY CONSTANT 1.728d-3 Pa-s
END

FLUID_PROPERTY
  DIFFUSION_COEFFICIENT 1.d-9
END

STRATA
  REGION all
  MATERIAL beam
END

TIME
  FINAL_TIME 0.75 day
  INITIAL_TIMESTEP_SIZE 1.d-4 day
  MAXIMUM_TIMESTEP_SIZE 0.0001 day at 0.d0 day
END

OUTPUT
  SNAPSHOT_FILE
    TIMES day 0.10 0.25 0.50
    #NO_PRINT_INITIAL
    FORMAT HDF5
  /
END

FLOW_CONDITION initial
  TYPE
    TEMPERATURE dirichlet
    PRESSURE dirichlet
  /
  TEMPERATURE 25.0d0 C
  PRESSURE 0.25 MPa
END

FLOW_CONDITION left_end
  TYPE
    TEMPERATURE dirichlet
    PRESSURE dirichlet
  /
  TEMPERATURE 25.0 C
  PRESSURE LIST
    # p = p_b*t + offset; p_b=2e6 Pa/day; offset = 0.25 MPa
    TIME_UNITS day
    DATA_UNITS Pa
    INTERPOLATION LINEAR
    #time  #pressure
    0.00d0 0.25d6
    0.25d0 0.75d6
    0.50d0 1.25d6
    1.00d0 2.25d6
  /
END

FLOW_CONDITION right_end
  TYPE
    TEMPERATURE dirichlet
    PRESSURE dirichlet
  /
  TEMPERATURE 25.0 C
  PRESSURE LIST
    # p = p_b*t + offset; p_b=2e6 Pa/day; offset = 0.25 MPa
    TIME_UNITS day
    DATA_UNITS Pa
    INTERPOLATION LINEAR
    #time  #pressure
    0.00d0 0.25d6
    0.25d0 0.75d6
    0.50d0 1.25d6
    1.00d0 2.25d6
  /
END

INITIAL_CONDITION initial
  REGION all
  FLOW_CONDITION initial
END

BOUNDARY_CONDITION left_end
  REGION left_end
  FLOW_CONDITION left_end
END

BOUNDARY_CONDITION right_end
  FLOW_CONDITION right_end
  REGION right_end
END

END_SUBSURFACE

The PFLOTRAN Input File (RICHARDS Mode)

The RICHARDS Mode PFLOTRAN input file can be downloaded here.

 
SIMULATION
  SIMULATION_TYPE SUBSURFACE
  PROCESS_MODELS
    SUBSURFACE_FLOW flow
      MODE RICHARDS
    /
  /
END

SUBSURFACE

GRID
  TYPE structured
  NXYZ 10 1 1
  DXYZ
   2.0d0
   1.0d0
   1.0d0
  END
END

REGION all
  COORDINATES
    0.0d0 0.0d0 0.0d0
    20.0d0 1.0d0 1.0d0
  /
END

REGION left_end
  FACE WEST
  COORDINATES
    0.d0 0.d0 0.d0
    0.d0 1.0d0 1.0d0
  /
END

REGION right_end
  FACE EAST
  COORDINATES
    20.0d0 0.d0 0.d0
    20.0d0 1.0d0 1.0d0
  /
END

MATERIAL_PROPERTY beam
  ID 1
  POROSITY 0.25
  TORTUOSITY 1.d0
  ROCK_DENSITY 2.5d3 kg/m^3
  CHARACTERISTIC_CURVES cc1
  PERMEABILITY
    PERM_X 1.d-14
    PERM_Y 1.d-14
    PERM_Z 1.d-14
  /
END

CHARACTERISTIC_CURVES cc1
  SATURATION_FUNCTION VAN_GENUCHTEN
    LIQUID_RESIDUAL_SATURATION 0.5d-1
    M 0.75
    ALPHA 1.d-3
  /
  PERMEABILITY_FUNCTION MUALEM
    LIQUID_RESIDUAL_SATURATION 0.1d0
    M 0.5d0
  /
  PERMEABILITY_FUNCTION MUALEM_VG_GAS
    LIQUID_RESIDUAL_SATURATION 0.1d0
    GAS_RESIDUAL_SATURATION 0.1d0
    M 0.5d0
  /
END

EOS WATER
  #ref. density [kg/m3], ref. pressure [Pa], water compressibility [1/Pa?]
  DENSITY EXPONENTIAL 1000. 0. 1.d-8
  VISCOSITY CONSTANT 1.728d-3 Pa-s
END

FLUID_PROPERTY
  DIFFUSION_COEFFICIENT 1.d-9
END

STRATA
  REGION all
  MATERIAL beam
END

TIME
  FINAL_TIME 0.75 day
  INITIAL_TIMESTEP_SIZE 1.d-4 day
  MAXIMUM_TIMESTEP_SIZE 0.0001 day at 0.d0 day
END

OUTPUT
  SNAPSHOT_FILE
    TIMES day 0.10 0.25 0.50
    #NO_PRINT_INITIAL
    FORMAT HDF5
  /
END

FLOW_CONDITION initial
  TYPE
    PRESSURE dirichlet
  /
  PRESSURE 0.25 MPa
END

FLOW_CONDITION left_end
  TYPE
    PRESSURE dirichlet
  /
  PRESSURE LIST
    # p = p_b*t + offset; p_b=2e6 Pa/day; offset = 0.25 MPa
    TIME_UNITS day
    DATA_UNITS Pa
    INTERPOLATION LINEAR
    #time  #pressure
    0.00d0 0.25d6
    0.25d0 0.75d6
    0.50d0 1.25d6
    1.00d0 2.25d6
  /
END

FLOW_CONDITION right_end
  TYPE
    PRESSURE dirichlet
  /
  PRESSURE LIST
    # p = p_b*t + offset; p_b=2e6 Pa/day; offset = 0.25 MPa
    TIME_UNITS day
    DATA_UNITS Pa
    INTERPOLATION LINEAR
    #time  #pressure
    0.00d0 0.25d6
    0.25d0 0.75d6
    0.50d0 1.25d6
    1.00d0 2.25d6
  /
END

INITIAL_CONDITION initial
  REGION all
  FLOW_CONDITION initial
END

BOUNDARY_CONDITION left_end
  REGION left_end
  FLOW_CONDITION left_end
END

BOUNDARY_CONDITION right_end
  FLOW_CONDITION right_end
  REGION right_end
END

END_SUBSURFACE

The Python Script

def flow_transient_1D_BC1stkind(path,input_prefix,remove,screen_on,pf_exe,
				                             mpi_exe,num_tries):
#==============================================================================#
# Based on:
# Kolditz, et al. (2015) Thermo-Hydro-Mechanical-Chemical Processes in
# Fractured Porous Media: Modelling and Benchmarking, Closed Form Solutions,
# Springer International Publishing, Switzerland.
# Section 2.2.7, pg.32
# "A Transient 1D Pressure Distribution, Time-Dependent Boundary Conditions
# of 1st Kind"
# With some modification of the parameter values given in Kolditz (2015)
#
# Author: Jenn Frederick
# Date: 08/23/2016
# ******************************************************************************
  nxyz = np.zeros(3) + 1
  dxyz = np.zeros(3) + 1.
  lxyz = np.zeros(3) + 1.
  error_analysis = np.zeros(num_tries)
  dxyz_record = np.zeros((num_tries,3))
  test_pass = False
  try_count = 0
  
  # initial discretization values
  lxyz[0] = 20.                # [m] lx
  nxyz[0] = 10                 # [m] nx 
  dxyz = lxyz/nxyz             # [m]
  k = 1.0e-14                  # [m^2]
  mu = 1.728e-3                # [Pa*sec]
  por = 0.25                   # [-]
  K = 1.0e-8                   # [1/Pa]
  pb_day = 2.0e6               # [Pa/day]
  pb_sec = 2.0e6/(24.0*3600.0) # [Pa/sec]
  p_offset = 0.25e6            # [Pa]
  chi = k/(por*K*mu)
  ierr = 0
  
  while (not test_pass) and (try_count < num_tries): 
    print_discretization(lxyz,nxyz,dxyz)
    nx = int(nxyz[0]); ny = int(nxyz[1]); nz = int(nxyz[2])
    dx = dxyz[0]; dy = dxyz[1]; dz = dxyz[2]
    Lx = lxyz[0]; Ly = lxyz[1]; Lz = lxyz[2]
    try_count = try_count + 1

    x_soln1 = np.linspace(0.+(dx/2.),(Lx/2.)-(dx/2.),nx/2.)        # [m]
    x_soln2 = np.linspace(-((Lx/2.)-(dx/2.)),-(0.+(dx/2.)),nx/2.)  # [m]
    x_soln = np.concatenate((x_soln2,x_soln1),axis=0)              # [m]
    # Add boundary values to analytical solution
    x_soln = np.concatenate(([-(Lx/2.)],x_soln),axis=0)            # [m]
    x_soln = np.concatenate((x_soln,[Lx/2.]),axis=0)               # [m]
    t_soln = np.array([0.10,0.25,0.50,0.75])                       # [day]
    p_soln = np.zeros((4,(nx+2)))                                  # [Pa]
    x_pflotran1 = np.linspace(0.+(dx/2.),(Lx/2.)-(dx/2.),nx/2.)    # [m]
    x_pflotran2 = np.linspace(-((Lx/2.)-(dx/2.)),-(0.+(dx/2.)),nx/2.)  # [m]
    x_pflotran = np.concatenate((x_pflotran2,x_pflotran1),axis=0)  # [m]
    p_pflotran = np.zeros((4,nx))                                  # [Pa]

    # create the analytical solution
    for time in range(4):
      t = t_soln[time]*(24.0*3600.0)   # [sec]
      p_soln[time,:] = pb_sec*t + ((pb_sec*(pow(x_soln,2)-pow((Lx/2.),2)))/(2.*chi))
      sum_term = np.zeros(nx+2)
      sum_term_old = np.zeros(nx+2)
      n = 0
      epsilon = 1.0
      while epsilon > epsilon_value:
	sum_term_old = sum_term
	sum_term = sum_term_old + (((pow(-1.,n))/(pow(((2*n)+1),3)))*np.cos((math.pi*x_soln*((2*n)+1))/(2*(Lx/2.)))*np.exp(-chi*pow((2*n)+1,2)*pow(math.pi,2)*(t/(4*pow((Lx/2.),2)))))
	epsilon = np.max(np.abs(sum_term_old-sum_term))
	n = n + 1
      p_soln[time,:] = p_offset + p_soln[time,:] + ((16.*pb_sec*pow((Lx/2.),2))/(chi*pow(math.pi,3)))*sum_term

    # run PFLOTRAN simulation
    run_pflotran(input_prefix,nxyz,dxyz,lxyz,remove,screen_on,pf_exe,mpi_exe)
    
    # load data from HDF5
    index_string = 'Time:  1.00000E-01 d/Liquid_Pressure [Pa]'
    p_pflotran[0,:] = read_pflotran_output_1D(path+
		      '/1D_transient_pressure_BC_1st_kind.h5',index_string,False)
    index_string = 'Time:  2.50000E-01 d/Liquid_Pressure [Pa]'
    p_pflotran[1,:] = read_pflotran_output_1D(path+
		      '/1D_transient_pressure_BC_1st_kind.h5',index_string,False)
    index_string = 'Time:  5.00000E-01 d/Liquid_Pressure [Pa]'
    p_pflotran[2,:] = read_pflotran_output_1D(path+
		      '/1D_transient_pressure_BC_1st_kind.h5',index_string,False)
    index_string = 'Time:  7.50000E-01 d/Liquid_Pressure [Pa]'
    p_pflotran[3,:] = read_pflotran_output_1D(path+
		      '/1D_transient_pressure_BC_1st_kind.h5',index_string,remove)
    ierr = check(p_pflotran)
  
    # Convert units [Pa -> MPa]
    p_pflotran = p_pflotran/1.e6
    p_soln = p_soln/1.e6

    max_percent_error = calc_relative_error(p_soln[:,1:nx+1],p_pflotran,ierr)
    record_error(error_analysis,dxyz_record,max_percent_error,dxyz,try_count)
    test_pass = does_pass(max_percent_error,try_count,num_tries)
    nxyz[0] = nxyz[0]*2.
    dxyz = lxyz/nxyz

  # Plot the PFLOTRAN and analytical solutions
  plot_1D_transient(path,t_soln,'day',x_soln,p_soln,x_pflotran,p_pflotran,
                    'Distance [m]','Pressure [MPa]',"{0:.2f}".format(max_percent_error))
  
  # Plot error analysis
  plot_error(error_analysis,dxyz_record,path,try_count)
  
  # Add test result to report card
  add_to_report(path,test_pass,max_percent_error,ierr)

  return;

Refer to section Python Helper Functions for documentation on the qa_tests_helper module, which defines the helper functions used in the Python script above.

1D Transient Flow (Pressure), BCs of 2nd Kind

The Problem Description

The PFLOTRAN Input File (GENERAL Mode)

The PFLOTRAN Input File (TH Mode)

The PFLOTRAN Input File (RICHARDS Mode)

The Python Script

The Problem Description

This problem is adapted from Kolditz, et al. (2015), Thermo-Hydro-Mechanical-Chemical Processes in Fractured Porous Media: Modelling and Benchmarking, Closed Form Solutions, Springer International Publishing, Switzerland. Section 2.2.8, pg.33, “A Transient 1D Pressure Distribution, Time-Dependent Boundary Conditions of 2nd Kind.”

The domain is a 25x1x1 meter rectangular column extending along the positive x-axis and is made up of 50x1x1 cubic grid cells with dimensions 0.5x1x1 meters. The domain material is assigned the following properties: porosity \(\phi\) = 0.25; \(k\) = 1.0e-14 m^2. The fluid viscosity \(\mu\) = 0.864 mPa/sec, and the fluid compressibility \(K\) = 1.0e-8 1/Pa. The accuracy of the numerical solution is very sensitive to the time step size because the fluid flux value is held constant at the value at the beginning of the timestep. Reducing the timestep will make the numerical solution more accurate relative to the analytical solution, but note this is not an issue with time truncation error.

The pressure is initially uniform at p(t=0) = 0 MPa. At the left end, a no fluid flux boundary condition is applied, and at the right end, a transient fluid flux boundary condition is applied (both in units of m/sec):

\[ \begin{align}\begin{aligned}q(0,t) = 0\\q(L,t) = 9.0x10^{-6} t\end{aligned}\end{align} \]

where L = 25 m. The transient pressure distribution is governed by,

\[ \begin{align}\begin{aligned}\chi = {k \over {\phi \mu K}}\\{\phi K} {{\partial p} \over {\partial t}} = {k \over \mu} {{\partial^{2} p} \over {\partial x^{2}}}\end{aligned}\end{align} \]

With the given boundary conditions, the solution is defined by,

\[p(x,t) = {{8q\sqrt{\chi t^3}} \over {k \over \mu}} \sum_{n=0}^{\infty} \left[{ i^3erfc{{(2n+1)L-x}\over{2\sqrt{\chi t}}} + i^3erfc{{(2n+1)L+x}\over{2\sqrt{\chi t}}} }\right]\]

where \(i^3erfc(g)\) represents the third repeated integral of the complimentary error function, given by,

\[i^3erfc(g) = {2 \over \pi} \int^{\infty}_{g} {{(s-g)^3}\over{3!}} e^{-s^2} ds\]
_images/visit_figure7.png

The PFLOTRAN domain set-up.

If you do not see this image, you must run the QA test suite to generate this figure.

Comparison of the PFLOTRAN vs. analytical solution for GENERAL mode.

If you do not see this image, you must run the QA test suite to generate this figure.

Comparison of the PFLOTRAN vs. analytical solution for TH mode.

If you do not see this image, you must run the QA test suite to generate this figure.

Comparison of the PFLOTRAN vs. analytical solution for RICHARDS mode.

The PFLOTRAN Input File (GENERAL Mode)

The General Mode PFLOTRAN input file can be downloaded here.

 
SIMULATION
  SIMULATION_TYPE SUBSURFACE
  PROCESS_MODELS
    SUBSURFACE_FLOW flow
      MODE GENERAL
    /
  /
END

SUBSURFACE

GRID
  TYPE structured
  NXYZ 10 1 1
  DXYZ
   2.5d0
   1.0d0
   1.0d0
  END
END

REGION all
  COORDINATES
    0.0d0 0.0d0 0.0d0
    25.0d0 1.0d0 1.0d0
  /
END

REGION left_end
  FACE WEST
  COORDINATES
    0.d0 0.d0 0.d0
    0.d0 1.0d0 1.0d0
  /
END

REGION right_end
  FACE EAST
  COORDINATES
    25.0d0 0.d0 0.d0
    25.0d0 1.0d0 1.0d0
  /
END

MATERIAL_PROPERTY beam
  ID 1
  POROSITY 0.20
  TORTUOSITY 1.d0
  ROCK_DENSITY 2.0d3 kg/m^3
  HEAT_CAPACITY 100. J/kg-C
  THERMAL_CONDUCTIVITY_DRY 1.16 W/m-C
  THERMAL_CONDUCTIVITY_WET 1.16 W/m-C
  CHARACTERISTIC_CURVES cc1
  PERMEABILITY
    PERM_X 1.d-14
    PERM_Y 1.d-14
    PERM_Z 1.d-14
  /
END

CHARACTERISTIC_CURVES cc1
  SATURATION_FUNCTION VAN_GENUCHTEN
    LIQUID_RESIDUAL_SATURATION 0.5d-1
    M 0.75
    ALPHA 1.d-3
  /
  PERMEABILITY_FUNCTION MUALEM
    LIQUID_RESIDUAL_SATURATION 0.1d0
    M 0.5d0
  /
  PERMEABILITY_FUNCTION MUALEM_VG_GAS
    LIQUID_RESIDUAL_SATURATION 0.1d0
    GAS_RESIDUAL_SATURATION 0.1d0
    M 0.5d0
  /
END

EOS WATER
  # ref. density, ref. pressure, water compressibility
  DENSITY EXPONENTIAL 1000. 101325. 1.d-8
  VISCOSITY CONSTANT 0.864d-3 Pa-s
END

FLUID_PROPERTY
  DIFFUSION_COEFFICIENT 1.d-9
END

LINEAR_SOLVER FLOW
  SOLVER DIRECT
  ZERO_PIVOT_TOL 1.d-20
END

STRATA
  REGION all
  MATERIAL beam
END

TIME
  FINAL_TIME 0.20 day
  INITIAL_TIMESTEP_SIZE 1.d-4 day
  MAXIMUM_TIMESTEP_SIZE 0.0001 day at 0.d0 day
END

OUTPUT
  SNAPSHOT_FILE
    TIMES day 0.01 0.04 0.09 0.12
    NO_PRINT_INITIAL
    FORMAT HDF5
  /
END

FLOW_CONDITION initial
  TYPE
    TEMPERATURE dirichlet
    LIQUID_PRESSURE dirichlet
    MOLE_FRACTION dirichlet
  /
  TEMPERATURE 25.0d0 C
  LIQUID_PRESSURE 101325 Pa
  MOLE_FRACTION 1.d-10
END

FLOW_CONDITION left_end
  TYPE
    LIQUID_FLUX neumann
    TEMPERATURE dirichlet
    GAS_FLUX neumann
  /
  LIQUID_FLUX 0.d0 m/sec
  TEMPERATURE 25.0d0 C
  GAS_FLUX 0.d0 m/sec
END

FLOW_CONDITION right_end
  TYPE
    LIQUID_FLUX neumann
    TEMPERATURE dirichlet
    GAS_FLUX neumann
  /
  LIQUID_FLUX LIST
    # liquid_flux = (9.0d-6)*t [m/sec]
    TIME_UNITS day
    DATA_UNITS m/sec
    INTERPOLATION LINEAR
    #time  #liquid_flux
    0.00d0 0.0d0
    1.00d0 9.0d-6
  /
  TEMPERATURE 25.0d0 C
  GAS_FLUX 0.d0 m/sec
END

INITIAL_CONDITION initial
  REGION all
  FLOW_CONDITION initial
END

BOUNDARY_CONDITION left_end
  REGION left_end
  FLOW_CONDITION left_end
END

BOUNDARY_CONDITION right_end
  FLOW_CONDITION right_end
  REGION right_end
END

END_SUBSURFACE

The PFLOTRAN Input File (TH Mode)

The TH Mode PFLOTRAN input file can be downloaded here.

 
SIMULATION
  SIMULATION_TYPE SUBSURFACE
  PROCESS_MODELS
    SUBSURFACE_FLOW flow
      MODE TH
    /
  /
END

SUBSURFACE

GRID
  TYPE structured
  NXYZ 10 1 1
  DXYZ
   2.5d0
   1.0d0
   1.0d0
  END
END

REGION all
  COORDINATES
    0.0d0 0.0d0 0.0d0
    25.0d0 1.0d0 1.0d0
  /
END

REGION left_end
  FACE WEST
  COORDINATES
    0.d0 0.d0 0.d0
    0.d0 1.0d0 1.0d0
  /
END

REGION right_end
  FACE EAST
  COORDINATES
    25.0d0 0.d0 0.d0
    25.0d0 1.0d0 1.0d0
  /
END

MATERIAL_PROPERTY beam
  ID 1
  POROSITY 0.20
  TORTUOSITY 1.d0
  ROCK_DENSITY 2.0d3 kg/m^3
  SPECIFIC_HEAT 100.0 J/kg-C
  THERMAL_CONDUCTIVITY_DRY 1.16 W/m-C
  THERMAL_CONDUCTIVITY_WET 1.16 W/m-C
  SATURATION_FUNCTION default
  PERMEABILITY
    PERM_X 1.d-14
    PERM_Y 1.d-14
    PERM_Z 1.d-14
  /
END

SATURATION_FUNCTION default
  SATURATION_FUNCTION_TYPE VAN_GENUCHTEN
  RESIDUAL_SATURATION 0.5d-1
  LAMBDA 0.75
  ALPHA 1.d-3
END

EOS WATER
  # ref. density, ref. pressure, water compressibility
  DENSITY EXPONENTIAL 1000. 101325. 1.d-8
  VISCOSITY CONSTANT 0.864d-3 Pa-s
END

FLUID_PROPERTY
  DIFFUSION_COEFFICIENT 1.d-9
END

LINEAR_SOLVER FLOW
  SOLVER DIRECT
  ZERO_PIVOT_TOL 1.d-20
END

STRATA
  REGION all
  MATERIAL beam
END

TIME
  FINAL_TIME 0.20 day
  INITIAL_TIMESTEP_SIZE 1.d-4 day
  MAXIMUM_TIMESTEP_SIZE 0.0001 day at 0.d0 day
END

OUTPUT
  SNAPSHOT_FILE
    TIMES day 0.01 0.04 0.09 0.12
    NO_PRINT_INITIAL
    FORMAT HDF5
  /
END

FLOW_CONDITION initial
  TYPE
    TEMPERATURE dirichlet
    PRESSURE dirichlet
  /
  TEMPERATURE 25.0d0 C
  PRESSURE 101325 Pa
END

FLOW_CONDITION left_end
  TYPE
    FLUX neumann
    TEMPERATURE dirichlet
  /
  TEMPERATURE 25.0 C
  FLUX 0.d0 m/sec
END

FLOW_CONDITION right_end
  TYPE
    TEMPERATURE dirichlet
    FLUX neumann
  /
  TEMPERATURE 25.0d0 C
  FLUX LIST
    # liquid_flux = (9.0d-6)*t [m/sec]
    TIME_UNITS day
    DATA_UNITS m/sec
    INTERPOLATION LINEAR
    #time  #liquid_flux
    0.00d0 0.0d0
    1.00d0 9.0d-6
  /
END

INITIAL_CONDITION initial
  REGION all
  FLOW_CONDITION initial
END

BOUNDARY_CONDITION left_end
  REGION left_end
  FLOW_CONDITION left_end
END

BOUNDARY_CONDITION right_end
  FLOW_CONDITION right_end
  REGION right_end
END

END_SUBSURFACE

The PFLOTRAN Input File (RICHARDS Mode)

The RICHARDS Mode PFLOTRAN input file can be downloaded here.

 
SIMULATION
  SIMULATION_TYPE SUBSURFACE
  PROCESS_MODELS
    SUBSURFACE_FLOW flow
      MODE RICHARDS
    /
  /
END

SUBSURFACE

GRID
  TYPE structured
  NXYZ 10 1 1
  DXYZ
   2.5d0
   1.0d0
   1.0d0
  END
END

REGION all
  COORDINATES
    0.0d0 0.0d0 0.0d0
    25.0d0 1.0d0 1.0d0
  /
END

REGION left_end
  FACE WEST
  COORDINATES
    0.d0 0.d0 0.d0
    0.d0 1.0d0 1.0d0
  /
END

REGION right_end
  FACE EAST
  COORDINATES
    25.0d0 0.d0 0.d0
    25.0d0 1.0d0 1.0d0
  /
END

MATERIAL_PROPERTY beam
  ID 1
  POROSITY 0.20
  TORTUOSITY 1.d0
  ROCK_DENSITY 2.0d3 kg/m^3
  CHARACTERISTIC_CURVES cc1
  PERMEABILITY
    PERM_X 1.d-14
    PERM_Y 1.d-14
    PERM_Z 1.d-14
  /
END

CHARACTERISTIC_CURVES cc1
  SATURATION_FUNCTION VAN_GENUCHTEN
    LIQUID_RESIDUAL_SATURATION 0.5d-1
    M 0.75
    ALPHA 1.d-3
  /
  PERMEABILITY_FUNCTION MUALEM
    LIQUID_RESIDUAL_SATURATION 0.1d0
    M 0.5d0
  /
  PERMEABILITY_FUNCTION MUALEM_VG_GAS
    LIQUID_RESIDUAL_SATURATION 0.1d0
    GAS_RESIDUAL_SATURATION 0.1d0
    M 0.5d0
  /
END

EOS WATER
  # ref. density, ref. pressure, water compressibility
  DENSITY EXPONENTIAL 1000. 101325. 1.d-8
  VISCOSITY CONSTANT 0.864d-3 Pa-s
END

FLUID_PROPERTY
  DIFFUSION_COEFFICIENT 1.d-9
END

LINEAR_SOLVER FLOW
  SOLVER DIRECT
  ZERO_PIVOT_TOL 1.d-20
END

STRATA
  REGION all
  MATERIAL beam
END

TIME
  FINAL_TIME 0.20 day
  INITIAL_TIMESTEP_SIZE 1.d-4 day
  MAXIMUM_TIMESTEP_SIZE 0.0001 day at 0.d0 day
END

OUTPUT
  SNAPSHOT_FILE
    TIMES day 0.01 0.04 0.09 0.12
    NO_PRINT_INITIAL
    FORMAT HDF5
  /
END

FLOW_CONDITION initial
  TYPE
    PRESSURE dirichlet
  /
  PRESSURE 101325 Pa
END

FLOW_CONDITION left_end
  TYPE
    FLUX neumann
  /
  FLUX 0.d0 m/sec
END

FLOW_CONDITION right_end
  TYPE
    FLUX neumann
  /
  FLUX LIST
    # fluid_flux = (9.0d-6)*t [m/sec]
    TIME_UNITS day
    DATA_UNITS m/sec
    INTERPOLATION LINEAR
    #time  #fluid_flux
    0.00d0 0.0d0
    1.00d0 9.0d-6
  /
END

INITIAL_CONDITION initial
  REGION all
  FLOW_CONDITION initial
END

BOUNDARY_CONDITION left_end
  REGION left_end
  FLOW_CONDITION left_end
END

BOUNDARY_CONDITION right_end
  FLOW_CONDITION right_end
  REGION right_end
END

END_SUBSURFACE

The Python Script

def flow_transient_1D_BC2ndkind(path,input_prefix,remove,screen_on,pf_exe,
				                             mpi_exe,num_tries):
#==============================================================================#
# Based on:
# Kolditz, et al. (2015) Thermo-Hydro-Mechanical-Chemical Processes in
# Fractured Porous Media: Modelling and Benchmarking, Closed Form Solutions,
# Springer International Publishing, Switzerland.
# Section 2.2.8, pg.20
# "A Transient 1D Pressure Distribution, Time-Dependent Boundary Conditions
# of 2nd Kind"
# With some modification of the parameter values given in Kolditz (2015)
#
# Author: Jenn Frederick
# Date: 08/24/2016
# ******************************************************************************
  nxyz = np.zeros(3) + 1
  dxyz = np.zeros(3) + 1.
  lxyz = np.zeros(3) + 1.
  error_analysis = np.zeros(num_tries)
  dxyz_record = np.zeros((num_tries,3))
  test_pass = False
  try_count = 0
  
  # initial discretization values
  lxyz[0] = 25.         # [m] lx
  nxyz[0] = 10          # [m] nx 
  dxyz = lxyz/nxyz      # [m]
  k = 1.0e-14           # [m2]
  mu = 0.864e-3         # [Pa*sec]
  por = 0.20 
  K = 1.0e-8            # [1/Pa]
  q = 9.0e-6            # [m/s/day]
  p_offset = 101325     # [Pa]
  chi = k/(por*K*mu)
  ierr = 0
  
  while (not test_pass) and (try_count < num_tries): 
    print_discretization(lxyz,nxyz,dxyz)
    nx = int(nxyz[0]); ny = int(nxyz[1]); nz = int(nxyz[2])
    dx = dxyz[0]; dy = dxyz[1]; dz = dxyz[2]
    Lx = lxyz[0]; Ly = lxyz[1]; Lz = lxyz[2]
    try_count = try_count + 1

    x_soln = np.linspace(0.+(dx/2.),Lx-(dx/2.),nx)  # [m]
    x_soln = np.concatenate((x_soln,[Lx]),axis=0)   # [m]
    p_soln = np.zeros((4,nx+1))                     # [Pa]
    x_pflotran = x_soln[0:nx]                       # [m]
    p_pflotran = np.zeros((4,nx))                   # [Pa]
    t_soln = np.array([0.01,0.04,0.09,0.12])        # [day]

    # create the analytical solution
    p1 = np.zeros(nx+1)
    p2 = np.zeros(nx+1)
    i3erfc1 = np.zeros(nx+1)
    i3erfc2 = np.zeros(nx+1)
    for time in range(4):
      t = t_soln[time]*24.0*3600.0  # [sec]
      sum_term = np.zeros(nx+1)
      sum_term_old = np.zeros(nx+1)
      n = 0
      epsilon = 1.0
      while epsilon > epsilon_value:
	p1 = ((((2*n)+1)*Lx)-x_soln)/(2.*math.sqrt(chi*t))
	i3erfc1 = (1./12.)*np.sqrt(math.pi)*p1**3*erf(1.0*p1) - (1./12.)*np.sqrt(math.pi)*p1**3 + (1./3.)*p1**2*np.exp(-1.0*p1**2) + 0.5*p1*(-0.5*p1*np.exp(-1.0*p1**2) + 0.25*np.sqrt(math.pi)*erf(1.0*p1)) - 0.125*np.sqrt(math.pi)*p1 + (1./12.)*np.exp(-1.0*p1**2)
	p2 = ((((2*n)+1)*Lx)+x_soln)/(2.*math.sqrt(chi*t))
	i3erfc2 = (1./12.)*np.sqrt(math.pi)*p2**3*erf(1.0*p2) - (1./12.)*np.sqrt(math.pi)*p2**3 + (1./3.)*p2**2*np.exp(-1.0*p2**2) + 0.5*p2*(-0.5*p2*np.exp(-1.0*p2**2) + 0.25*np.sqrt(math.pi)*erf(1.0*p2)) - 0.125*np.sqrt(math.pi)*p2 + (1./12.)*np.exp(-1.0*p2**2)
	sum_term_old = sum_term
	sum_term = sum_term_old + (2./math.sqrt(math.pi))*i3erfc1 + (2./math.sqrt(math.pi))*i3erfc2
	epsilon = np.max(np.abs(sum_term_old-sum_term))
	n = n + 1
      p_soln[time,:] = p_offset + ((8*q*(1./(24.*3600.))*math.sqrt(chi*pow(t,3)))/(k/mu))*sum_term

    # To calculate the definite integral, use sympy:
    # p1, s = sym.symbols('p1 s')
    # i3erfc1 = sym.integrate(((pow((s-p1),3))/(3.*2.*1.))*sym.exp(-1.*pow(s,2)),(s,p1,sym.oo))
    # p2, s = sym.symbols('p2 s')
    # i3erfc2 = sym.integrate(((pow((s-p2),3))/(3.*2.*1.))*sym.exp(-1.*pow(s,2)),(s,p2,sym.oo))

    # run PFLOTRAN simulation
    run_pflotran(input_prefix,nxyz,dxyz,lxyz,remove,screen_on,pf_exe,mpi_exe)
    
    # load data from HDF5
    index_string = 'Time:  1.00000E-02 d/Liquid_Pressure [Pa]'
    p_pflotran[0,:] = read_pflotran_output_1D(path+
		      '/1D_transient_pressure_BC_2nd_kind.h5',index_string,False)
    index_string = 'Time:  4.00000E-02 d/Liquid_Pressure [Pa]'
    p_pflotran[1,:] = read_pflotran_output_1D(path+
		      '/1D_transient_pressure_BC_2nd_kind.h5',index_string,False)
    index_string = 'Time:  9.00000E-02 d/Liquid_Pressure [Pa]'
    p_pflotran[2,:] = read_pflotran_output_1D(path+
		      '/1D_transient_pressure_BC_2nd_kind.h5',index_string,False)
    index_string = 'Time:  1.20000E-01 d/Liquid_Pressure [Pa]'
    p_pflotran[3,:] = read_pflotran_output_1D(path+
		      '/1D_transient_pressure_BC_2nd_kind.h5',index_string,remove)
    ierr = check(p_pflotran)
    
    # Convert units [Pa -> MPa]
    p_pflotran = p_pflotran/1.e6
    p_soln = p_soln/1.e6

    max_percent_error = calc_relative_error(p_soln[:,0:nx],p_pflotran,ierr)
    record_error(error_analysis,dxyz_record,max_percent_error,dxyz,try_count)
    test_pass = does_pass(max_percent_error,try_count,num_tries)
    nxyz[0] = nxyz[0]*2.
    dxyz = lxyz/nxyz

  # Plot the PFLOTRAN and analytical solutions
  plot_1D_transient(path,t_soln,'day',x_soln,p_soln,x_pflotran,p_pflotran,
                    'Distance [m]','Pressure [MPa]',"{0:.2f}".format(max_percent_error))
  
  # Plot error analysis
  plot_error(error_analysis,dxyz_record,path,try_count)
  
  # Add test result to report card
  add_to_report(path,test_pass,max_percent_error,ierr)

  return;

Refer to section Python Helper Functions for documentation on the qa_tests_helper module, which defines the helper functions used in the Python script above.

1D Transient Flow (Pressure), BCs of 1st & 2nd Kind

The Problem Description

The PFLOTRAN Input File (GENERAL Mode)

The PFLOTRAN Input File (TH Mode)

The PFLOTRAN Input File (RICHARDS Mode)

The Dataset

The Python Script

The Problem Description

This problem is adapted from Kolditz, et al. (2015), Thermo-Hydro-Mechanical-Chemical Processes in Fractured Porous Media: Modelling and Benchmarking, Closed Form Solutions, Springer International Publishing, Switzerland. Section 2.2.9, pg.35, “A Transient 1D Pressure Distribution, Non-Zero Initial Pressure, Boundary Conditions of 1st and 2nd Kind.”

The domain is a 100x1x1 meter rectangular beam extending along the positive x-axis and is made up of 50x1x1 hexahedral grid cells with dimensions 2x1x1 meters. The domain is composed of a single material and is assigned the following properties: porosity \(\phi\) = 0.20; permeability \(k\) = 1.0e-14 m^2; rock density \(\rho\) = 2,000 kg/m^3; fluid compressibility \(K\) = 1.0e-9 1/Pa; fluid viscosity \(\mu\) = 1.728e-3 Pa-sec.

The pressure is initially distributed according to p(x,t=0)=f(x), where f(x) is defined (in units of MPa) as

\[ \begin{align}\begin{aligned}f(x) = 0 \hspace{0.25in} 0 \leq x < {L \over 10}\\f(x) = {{10x} \over {3L}}-{1 \over 3} \hspace{0.25in} {L \over 10} \leq x < {{4L} \over 10}\\f(x) = 1 \hspace{0.25in} {{4L} \over 10} \leq x < {{6L} \over 10}\\f(x) = 3-{{10x} \over {3L}} \hspace{0.25in} {{6L} \over 10} \leq x < {{9L} \over 10}\\f(x) = 0 \hspace{0.25in} {{9L} \over 10} \leq x \leq L\end{aligned}\end{align} \]

At the two boundaries, a no fluid flux condition is applied,

\[ \begin{align}\begin{aligned}q(0,t) = 0\\q(L,t) = 0\end{aligned}\end{align} \]

where L = 100 m. The transient pressure distribution is governed by,

\[{\phi K} {{\partial p} \over {\partial t}} = {k \over \mu} {{\partial^{2} p} \over {\partial x^{2}}}\]

With the initial pressure given, the solution is defined by,

\[ \begin{align}\begin{aligned}p(x,t) = {1 \over 2} + \sum_{n=1}^{\infty} exp\left({-\chi n^2 \pi^2 {t \over L^2}}\right)\left({80 \over {3(n\pi)^2}}\right) cos{{n \pi y} \over L} cos{{n\pi} \over 2} sin{{n\pi} \over 4} sin{{3n\pi} \over 20}\\\chi = {{k} \over {\phi \mu K}}\end{aligned}\end{align} \]
_images/visit_figure8.png

The PFLOTRAN domain set-up.

If you do not see this image, you must run the QA test suite to generate this figure.

Comparison of the PFLOTRAN vs. analytical solution for GENERAL mode.

If you do not see this image, you must run the QA test suite to generate this figure.

Comparison of the PFLOTRAN vs. analytical solution for TH mode.

If you do not see this image, you must run the QA test suite to generate this figure.

Comparison of the PFLOTRAN vs. analytical solution for RICHARDS mode.

The PFLOTRAN Input File (GENERAL Mode)

The GENERAL Mode PFLOTRAN input file can be downloaded here.

 
SIMULATION
  SIMULATION_TYPE SUBSURFACE
  PROCESS_MODELS
    SUBSURFACE_FLOW flow
      MODE GENERAL
    /
  /
END

SUBSURFACE

GRID
  TYPE structured
  NXYZ 10 1 1
  DXYZ
   10.0d0
   1.0d0
   1.0d0
  END
END

REGION all
  COORDINATES
    0.0d0 0.0d0 0.0d0
    100.0d0 1.0d0 1.0d0
  /
END

REGION left_end
  FACE WEST
  COORDINATES
    0.d0 0.d0 0.d0
    0.d0 1.0d0 1.0d0
  /
END

REGION right_end
  FACE EAST
  COORDINATES
    100.0d0 0.d0 0.d0
    100.0d0 1.0d0 1.0d0
  /
END

MATERIAL_PROPERTY beam
  ID 1
  POROSITY 0.20
  TORTUOSITY 1.d0
  ROCK_DENSITY 2.0d3 kg/m^3  
  HEAT_CAPACITY 100. J/kg-C
  THERMAL_CONDUCTIVITY_DRY 0.5787037 W/m-C
  THERMAL_CONDUCTIVITY_WET 0.5787037 W/m-C
  CHARACTERISTIC_CURVES cc1
  PERMEABILITY
    PERM_X 1.d-14
    PERM_Y 1.d-14
    PERM_Z 1.d-14
  /
END

CHARACTERISTIC_CURVES cc1
  SATURATION_FUNCTION VAN_GENUCHTEN
    LIQUID_RESIDUAL_SATURATION 0.5d-1
    M 0.75
    ALPHA 1.d-3
  /
  PERMEABILITY_FUNCTION MUALEM
    LIQUID_RESIDUAL_SATURATION 0.1d0
    M 0.5d0
  /
  PERMEABILITY_FUNCTION MUALEM_VG_GAS
    LIQUID_RESIDUAL_SATURATION 0.1d0
    GAS_RESIDUAL_SATURATION 0.1d0
    M 0.5d0
  /
END

EOS WATER
  #ref. density [kg/m3], ref. pressure [Pa], water compressibility [1/Pa]
  DENSITY EXPONENTIAL 1000. 101325. 1.0d-9
  VISCOSITY CONSTANT 1.728d-3 Pa-s
END

FLUID_PROPERTY
  DIFFUSION_COEFFICIENT 1.d-9
END

LINEAR_SOLVER FLOW
  SOLVER DIRECT
  ZERO_PIVOT_TOL 1.d-20
END

STRATA
  REGION all
  MATERIAL beam
END

TIME
  FINAL_TIME 0.50 day
  INITIAL_TIMESTEP_SIZE 1.d-4 day
  MAXIMUM_TIMESTEP_SIZE 0.0001 day at 0.d0 day
END

OUTPUT
  SNAPSHOT_FILE
    TIMES day 0.05 0.1 0.25
    #NO_PRINT_INITIAL
    FORMAT HDF5
  /
END

DATASET pressure_initial
  HDF5_DATASET_NAME initial
  FILENAME ../dataset.h5
END

FLOW_CONDITION initial
  TYPE
    LIQUID_PRESSURE dirichlet
    TEMPERATURE dirichlet
    MOLE_FRACTION dirichlet
  /
  LIQUID_PRESSURE DATASET pressure_initial
  TEMPERATURE 25.0 C
  MOLE_FRACTION 1.d-10
END

FLOW_CONDITION left_end
  TYPE
    LIQUID_FLUX neumann
    ENERGY_FLUX neumann
    GAS_FLUX neumann
  /
  LIQUID_FLUX 0.d0 m/yr
  ENERGY_FLUX 0.d0 W/m^2
  GAS_FLUX 0.d0 m/yr
END

FLOW_CONDITION right_end
  TYPE
    LIQUID_FLUX neumann
    ENERGY_FLUX neumann
    GAS_FLUX neumann
  /
  LIQUID_FLUX 0.d0 m/yr
  ENERGY_FLUX 0.d0 W/m^2
  GAS_FLUX 0.d0 m/yr
END

INITIAL_CONDITION initial
  REGION all
  FLOW_CONDITION initial
END

BOUNDARY_CONDITION left_end
  REGION left_end
  FLOW_CONDITION left_end
END

BOUNDARY_CONDITION right_end
  FLOW_CONDITION right_end
  REGION right_end
END

END_SUBSURFACE

The PFLOTRAN Input File (TH Mode)

The TH Mode PFLOTRAN input file can be downloaded here.

 
SIMULATION
  SIMULATION_TYPE SUBSURFACE
  PROCESS_MODELS
    SUBSURFACE_FLOW flow
      MODE TH
    /
  /
END

SUBSURFACE

GRID
  TYPE structured
  NXYZ 10 1 1
  DXYZ
   10.0d0
   1.0d0
   1.0d0
  END
END

REGION all
  COORDINATES
    0.0d0 0.0d0 0.0d0
    100.0d0 1.0d0 1.0d0
  /
END

REGION left_end
  FACE WEST
  COORDINATES
    0.d0 0.d0 0.d0
    0.d0 1.0d0 1.0d0
  /
END

REGION right_end
  FACE EAST
  COORDINATES
    100.0d0 0.d0 0.d0
    100.0d0 1.0d0 1.0d0
  /
END

MATERIAL_PROPERTY beam
  ID 1
  POROSITY 0.20
  TORTUOSITY 1.d0
  ROCK_DENSITY 2.0d3 kg/m^3
  SPECIFIC_HEAT 100. J/kg-C
  THERMAL_CONDUCTIVITY_DRY 0.5787037 W/m-C
  THERMAL_CONDUCTIVITY_WET 0.5787037 W/m-C
  SATURATION_FUNCTION default
  PERMEABILITY
    PERM_X 1.d-14
    PERM_Y 1.d-14
    PERM_Z 1.d-14
  /
END

SATURATION_FUNCTION default
  SATURATION_FUNCTION_TYPE VAN_GENUCHTEN
  RESIDUAL_SATURATION 0.5d-1
  LAMBDA 0.75
  ALPHA 1.d-3
END

EOS WATER
  #ref. density [kg/m3], ref. pressure [Pa], water compressibility [1/Pa]
  DENSITY EXPONENTIAL 1000. 101325. 1.0d-9
  VISCOSITY CONSTANT 1.728d-3 Pa-s
END

FLUID_PROPERTY
  DIFFUSION_COEFFICIENT 1.d-9
END

LINEAR_SOLVER FLOW
  SOLVER DIRECT
  ZERO_PIVOT_TOL 1.d-20
END

STRATA
  REGION all
  MATERIAL beam
END

TIME
  FINAL_TIME 0.50 day
  INITIAL_TIMESTEP_SIZE 1.d-4 day
  MAXIMUM_TIMESTEP_SIZE 0.0001 day at 0.d0 day
END

OUTPUT
  SNAPSHOT_FILE
    TIMES day 0.05 0.1 0.25
    #NO_PRINT_INITIAL
    FORMAT HDF5
  /
END

DATASET pressure_initial
  HDF5_DATASET_NAME initial
  FILENAME ../dataset.h5
END

FLOW_CONDITION initial
  TYPE
    PRESSURE dirichlet
    TEMPERATURE dirichlet
  /
  PRESSURE DATASET pressure_initial
  TEMPERATURE 25.0 C
END

FLOW_CONDITION left_end
  TYPE
    FLUX neumann
    ENERGY_FLUX neumann
  /
  FLUX 0.d0 m/yr
  ENERGY_FLUX 0.d0 W/m^2
END

FLOW_CONDITION right_end
  TYPE
    FLUX neumann
    ENERGY_FLUX neumann
  /
  FLUX 0.d0 m/yr
  ENERGY_FLUX 0.d0 W/m^2
END

INITIAL_CONDITION initial
  REGION all
  FLOW_CONDITION initial
END

BOUNDARY_CONDITION left_end
  REGION left_end
  FLOW_CONDITION left_end
END

BOUNDARY_CONDITION right_end
  FLOW_CONDITION right_end
  REGION right_end
END

END_SUBSURFACE

The PFLOTRAN Input File (RICHARDS Mode)

The RICHARDS Mode PFLOTRAN input file can be downloaded here.

 
SIMULATION
  SIMULATION_TYPE SUBSURFACE
  PROCESS_MODELS
    SUBSURFACE_FLOW flow
      MODE RICHARDS
    /
  /
END

SUBSURFACE

GRID
  TYPE structured
  NXYZ 10 1 1
  DXYZ
   10.0d0
   1.0d0
   1.0d0
  END
END

REGION all
  COORDINATES
    0.0d0 0.0d0 0.0d0
    100.0d0 1.0d0 1.0d0
  /
END

REGION left_end
  FACE WEST
  COORDINATES
    0.d0 0.d0 0.d0
    0.d0 1.0d0 1.0d0
  /
END

REGION right_end
  FACE EAST
  COORDINATES
    100.0d0 0.d0 0.d0
    100.0d0 1.0d0 1.0d0
  /
END

MATERIAL_PROPERTY beam
  ID 1
  POROSITY 0.20
  TORTUOSITY 1.d0
  ROCK_DENSITY 2.0d3 kg/m^3
  CHARACTERISTIC_CURVES cc1
  PERMEABILITY
    PERM_X 1.d-14
    PERM_Y 1.d-14
    PERM_Z 1.d-14
  /
END

CHARACTERISTIC_CURVES cc1
  SATURATION_FUNCTION VAN_GENUCHTEN
    LIQUID_RESIDUAL_SATURATION 0.5d-1
    M 0.75
    ALPHA 1.d-3
  /
  PERMEABILITY_FUNCTION MUALEM
    LIQUID_RESIDUAL_SATURATION 0.1d0
    M 0.5d0
  /
  PERMEABILITY_FUNCTION MUALEM_VG_GAS
    LIQUID_RESIDUAL_SATURATION 0.1d0
    GAS_RESIDUAL_SATURATION 0.1d0
    M 0.5d0
  /
END

EOS WATER
  #ref. density [kg/m3], ref. pressure [Pa], water compressibility [1/Pa]
  DENSITY EXPONENTIAL 1000. 101325. 1.0d-9
  VISCOSITY CONSTANT 1.728d-3 Pa-s
END

FLUID_PROPERTY
  DIFFUSION_COEFFICIENT 1.d-9
END

LINEAR_SOLVER FLOW
  SOLVER DIRECT
  ZERO_PIVOT_TOL 1.d-20
END

STRATA
  REGION all
  MATERIAL beam
END

TIME
  FINAL_TIME 0.50 day
  INITIAL_TIMESTEP_SIZE 1.d-4 day
  MAXIMUM_TIMESTEP_SIZE 0.0001 day at 0.d0 day
END

OUTPUT
  SNAPSHOT_FILE
    TIMES day 0.05 0.1 0.25
    #NO_PRINT_INITIAL
    FORMAT HDF5
  /
END

DATASET pressure_initial
  HDF5_DATASET_NAME initial
  FILENAME ../dataset.h5
END

FLOW_CONDITION initial
  TYPE
    PRESSURE dirichlet
  /
  PRESSURE DATASET pressure_initial
END

FLOW_CONDITION left_end
  TYPE
    FLUX neumann
  /
  FLUX 0.d0 m/yr
END

FLOW_CONDITION right_end
  TYPE
    FLUX neumann
  /
  FLUX 0.d0 m/yr
END

INITIAL_CONDITION initial
  REGION all
  FLOW_CONDITION initial
END

BOUNDARY_CONDITION left_end
  REGION left_end
  FLOW_CONDITION left_end
END

BOUNDARY_CONDITION right_end
  FLOW_CONDITION right_end
  REGION right_end
END

END_SUBSURFACE

The Dataset

The hdf5 dataset required to define the initial/boundary conditions is created with the following python script called create_dataset.py:

import sys
from h5py import *
import numpy as np

filename = 'dataset.h5'
h5file = File(filename,mode='w')

# 1d line in x
# Pressure initial condition [Pa]
L = 100.
h5grp = h5file.create_group('initial')
h5grp.attrs['Dimension'] = np.string_('X')
# Delta length between points [m]
h5grp.attrs['Discretization'] = [1.]
# Location of origin
h5grp.attrs['Origin'] = [0.]
# Load the dataset values
nx = 101
p_offset = 0.101325   # [MPa]
rarray = np.zeros(nx,'=f8')
for i in range(nx):
  if (0. <= i < (L/10.)):
    rarray[i] = ( 0. + p_offset ) * 1.0e6  # [Pa]
  if ((L/10.) <= i < (4.*L/10.)):
    rarray[i] = ( (10./(3.*L))*float(i) - (1./3.) + p_offset ) * 1.0e6  # [Pa]
  if ((4.*L/10.) <= i < (6.*L/10.)):
    rarray[i] = ( 1. + p_offset ) * 1.0e6  # [Pa]
  if ((6.*L/10.) <= i < (9.*L/10.)):
    rarray[i] = ( 3. - (10./(3.*L))*float(i) + p_offset ) * 1.0e6  # [Pa]
  if ((9.*L/10.) <= i < L):
    rarray[i] = ( 0. + p_offset ) * 1.0e6  # [Pa]
h5dset = h5grp.create_dataset('Data', data=rarray)

The Python Script

def flow_transient_1D_BC1st2ndkind(path,input_prefix,remove,screen_on,pf_exe,
				                             mpi_exe,num_tries):
#==============================================================================#
# Based on:
# Kolditz, et al. (2015) Thermo-Hydro-Mechanical-Chemical Processes in
# Fractured Porous Media: Modelling and Benchmarking, Closed Form Solutions,
# Springer International Publishing, Switzerland.
# Section 2.2.9, pg.35
# "A Transient 1D Pressure Distribution, Non-Zero Initial Pressure,
# Boundary Conditions of 1st and 2nd Kind"
# With some modification of the parameter values given in Kolditz (2015)
#
# Author: Jenn Frederick
# Date: 09/09/2016
# ******************************************************************************
  nxyz = np.zeros(3) + 1
  dxyz = np.zeros(3) + 1.
  lxyz = np.zeros(3) + 1.
  error_analysis = np.zeros(num_tries)
  dxyz_record = np.zeros((num_tries,3))
  test_pass = False
  try_count = 0
  
  # initial discretization values
  lxyz[0] = 100.        # [m] lx
  nxyz[0] = 10          # [m] nx 
  dxyz = lxyz/nxyz      # [m]
  k = 1.0e-14           # [m2]
  mu = 1.728e-3         # [Pa-s]
  por = 0.20
  K = 1.0e-9            # [1/Pa]
  p_offset = 0.101325   # [MPa]
  chi = k/(por*K*mu)
  ierr = 0
  
  while (not test_pass) and (try_count < num_tries): 
    print_discretization(lxyz,nxyz,dxyz)
    nx = int(nxyz[0]); ny = int(nxyz[1]); nz = int(nxyz[2])
    dx = dxyz[0]; dy = dxyz[1]; dz = dxyz[2]
    Lx = lxyz[0]; Ly = lxyz[1]; Lz = lxyz[2]
    try_count = try_count + 1

    x_soln = np.linspace(0.+(dx/2.),Lx-(dx/2.),nx)      # [m]
    x_pflotran = x_soln                                 # [m]
    t_soln = np.array([0.05,0.10,0.25,0.50])            # [day]
    p_soln = np.zeros((4,nx))                           # [MPa]
    p_pflotran = np.zeros((4,nx))                       # [Pa]

    # create the analytical solution
    for time in range(4):
      t = t_soln[time]*24.0*3600.0  # [sec]
      sum_term_old = np.zeros(nx)
      sum_term = np.zeros(nx)
      n = 1
      epsilon = 1.0
      while epsilon > epsilon_value:
	sum_term_old = sum_term
	sum_term = sum_term_old + (np.cos(n*math.pi*x_soln/Lx)*np.exp(-chi*pow(n,2)*pow(math.pi,2)*t/pow(Lx,2))*(80./(3.*pow((n*math.pi),2)))*np.cos(n*math.pi/2.)*np.sin(n*math.pi/4.)*np.sin(3.*n*math.pi/20.))
	epsilon = np.max(np.abs(sum_term_old-sum_term))
	n = n + 1
      p_soln[time,:] = 0.50 + sum_term + p_offset
      
    # run PFLOTRAN simulation
    run_pflotran(input_prefix,nxyz,dxyz,lxyz,remove,screen_on,pf_exe,mpi_exe)
    
    # load data from HDF5
    index_string = 'Time:  5.00000E-02 d/Liquid_Pressure [Pa]'
    p_pflotran[0,:] = read_pflotran_output_1D(path+
		 '/1D_transient_pressure_BC_1st_2nd_kind.h5',index_string,False)
    index_string = 'Time:  1.00000E-01 d/Liquid_Pressure [Pa]'
    p_pflotran[1,:] = read_pflotran_output_1D(path+
		 '/1D_transient_pressure_BC_1st_2nd_kind.h5',index_string,False)
    index_string = 'Time:  2.50000E-01 d/Liquid_Pressure [Pa]'
    p_pflotran[2,:] = read_pflotran_output_1D(path+
		 '/1D_transient_pressure_BC_1st_2nd_kind.h5',index_string,False)
    index_string = 'Time:  5.00000E-01 d/Liquid_Pressure [Pa]'
    p_pflotran[3,:] = read_pflotran_output_1D(path+
		'/1D_transient_pressure_BC_1st_2nd_kind.h5',index_string,remove)
    ierr = check(p_pflotran)

    # Convert units [Pa -> MPa]
    p_pflotran = p_pflotran/1.e6

    max_percent_error = calc_relative_error(p_soln,p_pflotran,ierr)
    record_error(error_analysis,dxyz_record,max_percent_error,dxyz,try_count)
    test_pass = does_pass(max_percent_error,try_count,num_tries)
    nxyz[0] = nxyz[0]*2.
    dxyz = lxyz/nxyz

  # Plot the PFLOTRAN and analytical solutions
  plot_1D_transient(path,t_soln,'day',x_soln,p_soln,x_pflotran,p_pflotran,
                    'Distance [m]','Pressure [MPa]',"{0:.2f}".format(max_percent_error))

  # Plot error analysis
  plot_error(error_analysis,dxyz_record,path,try_count)
  
  # Add test result to report card
  add_to_report(path,test_pass,max_percent_error,ierr)

  return;

Refer to section Python Helper Functions for documentation on the qa_tests_helper module, which defines the helper functions used in the Python script above.

2D Transient Flow (Pressure), BCs of 1st & 2nd Kind

The Problem Description

The PFLOTRAN Input File (GENERAL Mode)

The PFLOTRAN Input File (TH Mode)

The PFLOTRAN Input File (RICHARDS Mode)

The Dataset

The Python Script

The Problem Description

This problem is adapted from Kolditz, et al. (2015), Thermo-Hydro-Mechanical-Chemical Processes in Fractured Porous Media: Modelling and Benchmarking, Closed Form Solutions, Springer International Publishing, Switzerland. Section 2.2.10, pg.37, “A Transient 2D Pressure Distribution, Non-Zero Initial Pressure, Boundary Conditions of 1st and 2nd Kind.”

The domain is a 100x100x1 meter rectangular plate extending along the positive x-axis and y-axis and is made up of 50x50x1 hexahedral grid cells with dimensions 2x2x1 meters. The domain is composed of a single material and is assigned the following properties: porosity \(\phi\) = 0.20; permeability \(k\) = 1.0e-14 m^2; rock density \(\rho\) = 2,000 kg/m^3; fluid compressibility \(K\) = 1.0e-9 1/Pa; fluid viscosity \(\mu\) = 1.728e-3 Pa-sec.

The pressure is initially distributed according to p(x,y,t=0)=f(x)*f(y)+ \(p_{offset}\), where f(x) is defined (in units of MPa) as

\[ \begin{align}\begin{aligned}f(x) = 0 \hspace{0.25in} 0 \leq x < {L \over 10}\\f(x) = {{10x} \over {3L}}-{1 \over 3} \hspace{0.25in} {L \over 10} \leq x < {{4L} \over 10}\\f(x) = 1 \hspace{0.25in} {{4L} \over 10} \leq x < {{6L} \over 10}\\f(x) = 3-{{10x} \over {3L}} \hspace{0.25in} {{6L} \over 10} \leq x < {{9L} \over 10}\\f(x) = 0 \hspace{0.25in} {{9L} \over 10} \leq x \leq L\end{aligned}\end{align} \]

and f(y) is defined (in units of MPa) as

\[ \begin{align}\begin{aligned}f(y) = 0 \hspace{0.25in} 0 \leq y < {L \over 10}\\f(y) = {{10y} \over {3L}}-{1 \over 3} \hspace{0.25in} {L \over 10} \leq y < {{4L} \over 10}\\f(y) = 1 \hspace{0.25in} {{4L} \over 10} \leq y < {{6L} \over 10}\\f(y) = 3-{{10y} \over {3L}} \hspace{0.25in} {{6L} \over 10} \leq y < {{9L} \over 10}\\f(y) = 0 \hspace{0.25in} {{9L} \over 10} \leq y \leq L\end{aligned}\end{align} \]

At the north and south boundaries, a no fluid flux condition is applied,

\[ \begin{align}\begin{aligned}q(x,0,t) = 0\\q(x,L,t) = 0\end{aligned}\end{align} \]

and at the east and west boundaries, a pressure boundary condition is applied,

\[ \begin{align}\begin{aligned}p(0,y,t) = 0\\p(L,y,t) = 0\end{aligned}\end{align} \]

where L = 100 m and \(p_{offset}\) = 0.101325 MPa. The transient pressure distribution is governed by,

\[{\phi K} {{\partial p} \over {\partial t}} = {k \over \mu} \left({ {{\partial^{2} p} \over {\partial x^{2}}} + {{\partial^{2} p} \over {\partial y^{2}}} }\right)\]

With the initial pressure given, the solution is defined by,

\[ \begin{align}\begin{aligned}p(x,y,t) = p_1(x,t) p_2(y,t) + p_{offset}\\p_1(x,t) = \sum_{n=1}^{\infty} exp\left({-\chi n^2 \pi^2 {t \over L^2}}\right)\left({80 \over {3(n\pi)^2}}\right) sin{{n \pi x} \over L} sin{{n\pi} \over 2} sin{{n\pi} \over 4} sin{{3n\pi} \over 20}\\p_2(y,t) = {1 \over 2} + \sum_{n=1}^{\infty} exp\left({-\chi n^2 \pi^2 {t \over L^2}}\right)\left({80 \over {3(n\pi)^2}}\right) cos{{n \pi y} \over L} cos{{n\pi} \over 2} sin{{n\pi} \over 4} sin{{3n\pi} \over 20}\\\chi = {{k} \over {\phi \mu K}}\end{aligned}\end{align} \]
_images/visit_figure9.png

The PFLOTRAN domain set-up.

If you do not see this image, you must run the QA test suite to generate this figure.

Comparison of the PFLOTRAN vs. analytical solution for GENERAL mode.

If you do not see this image, you must run the QA test suite to generate this figure.

Comparison of the PFLOTRAN vs. analytical solution for TH mode.

If you do not see this image, you must run the QA test suite to generate this figure.

Comparison of the PFLOTRAN vs. analytical solution for RICHARDS mode.

The PFLOTRAN Input File (GENERAL Mode)

The GENERAL Mode PFLOTRAN input file can be downloaded here.

 
SIMULATION
  SIMULATION_TYPE SUBSURFACE
  PROCESS_MODELS
    SUBSURFACE_FLOW flow
      MODE GENERAL
    /
  /
END

SUBSURFACE

GRID
  TYPE structured
  NXYZ 40 40 1
  DXYZ
   2.5d0
   2.5d0
   1.0d0
  END
END

REGION all
  COORDINATES
    0.0d0 0.0d0 0.0d0
    100.0d0 100.0d0 1.0d0
  /
END

REGION west_face
  FACE WEST
  COORDINATES
    0.d0 0.d0 0.d0
    0.d0 100.0d0 1.0d0
  /
END

REGION east_face
  FACE EAST
  COORDINATES
    100.0d0 0.d0 0.d0
    100.0d0 100.0d0 1.0d0
  /
END

REGION north_face
  FACE NORTH
  COORDINATES
    0.d0 100.0d0 0.d0
    100.0d0 100.0d0 1.0d0
  /
END

REGION south_face
  FACE SOUTH
  COORDINATES
    0.d0  0.d0  0.d0
    100.0d0 0.d0 1.0d0
  /
END

MATERIAL_PROPERTY plate
  ID 1
  POROSITY 0.20
  TORTUOSITY 1.d0
  ROCK_DENSITY 2.0E3
  HEAT_CAPACITY 100. J/kg-C
  THERMAL_CONDUCTIVITY_DRY 0.5787037 W/m-C
  THERMAL_CONDUCTIVITY_WET 0.5787037 W/m-C
  CHARACTERISTIC_CURVES cc1
  PERMEABILITY
    PERM_X 1.d-14
    PERM_Y 1.d-14
    PERM_Z 1.d-14
  /
END

CHARACTERISTIC_CURVES cc1
  SATURATION_FUNCTION VAN_GENUCHTEN
    LIQUID_RESIDUAL_SATURATION 0.5d-1
    M 0.75
    ALPHA 1.d-3
  /
  PERMEABILITY_FUNCTION MUALEM
    LIQUID_RESIDUAL_SATURATION 0.1d0
    M 0.5d0
  /
  PERMEABILITY_FUNCTION MUALEM_VG_GAS
    LIQUID_RESIDUAL_SATURATION 0.1d0
    GAS_RESIDUAL_SATURATION 0.1d0
    M 0.5d0
  /
END

EOS WATER
  DENSITY EXPONENTIAL 1000. 101325. 1.0d-9
  VISCOSITY CONSTANT 1.728d-3 Pa-s
END

FLUID_PROPERTY
  DIFFUSION_COEFFICIENT 1.d-9
END

STRATA
  REGION all
  MATERIAL plate
END

TIME
  FINAL_TIME 0.10 day
  INITIAL_TIMESTEP_SIZE 1.d-4 day
  MAXIMUM_TIMESTEP_SIZE 0.0001 day at 0.d0 day
END

OUTPUT
  SNAPSHOT_FILE 
    #NO_PRINT_INITIAL
    TIMES day 0.01 0.04 0.06
    #NO_PRINT_FINAL
    FORMAT HDF5
  /
END

DATASET pressure_initial
  HDF5_DATASET_NAME initial
  FILENAME ../dataset.h5
END

FLOW_CONDITION initial
  TYPE
    TEMPERATURE dirichlet
    LIQUID_PRESSURE dirichlet
    MOLE_FRACTION dirichlet
  /
  TEMPERATURE 25.0 C
  LIQUID_PRESSURE DATASET pressure_initial
  MOLE_FRACTION 1.d-20
END

FLOW_CONDITION west_face
  TYPE
    TEMPERATURE dirichlet
    LIQUID_PRESSURE dirichlet
    MOLE_FRACTION dirichlet
  /
  TEMPERATURE 25.0 C
  LIQUID_PRESSURE 101325. Pa
  MOLE_FRACTION 1.d-20
END

FLOW_CONDITION east_face
  TYPE
    TEMPERATURE dirichlet
    LIQUID_PRESSURE dirichlet
    MOLE_FRACTION dirichlet
  /
  TEMPERATURE 25.0 C
  LIQUID_PRESSURE 101325. Pa
  MOLE_FRACTION 1.d-20
END

FLOW_CONDITION north_face
  TYPE
    ENERGY_FLUX neumann
    LIQUID_FLUX neumann
    GAS_FLUX neumann
  /
  ENERGY_FLUX 0.d0 W/m^2
  LIQUID_FLUX 0.d0 m/yr
  GAS_FLUX 0.d0 m/yr
END

FLOW_CONDITION south_face
  TYPE
    ENERGY_FLUX neumann
    LIQUID_FLUX neumann
    GAS_FLUX neumann
  /
  ENERGY_FLUX 0.d0 W/m^2
  LIQUID_FLUX 0.d0 m/yr
  GAS_FLUX 0.d0 m/yr
END

INITIAL_CONDITION initial
  REGION all
  FLOW_CONDITION initial
END

BOUNDARY_CONDITION west_face
  REGION west_face
  FLOW_CONDITION west_face
END

BOUNDARY_CONDITION east_face
  FLOW_CONDITION east_face
  REGION east_face
END

BOUNDARY_CONDITION north_face
  REGION north_face
  FLOW_CONDITION north_face
END

BOUNDARY_CONDITION south_face
  FLOW_CONDITION south_face
  REGION south_face
END

END_SUBSURFACE

The PFLOTRAN Input File (TH Mode)

The TH Mode PFLOTRAN input file can be downloaded here.

 
SIMULATION
  SIMULATION_TYPE SUBSURFACE
  PROCESS_MODELS
    SUBSURFACE_FLOW flow
      MODE TH
    /
  /
END

SUBSURFACE

GRID
  TYPE structured
  NXYZ 40 40 1
  DXYZ
   2.5d0
   2.5d0
   1.0d0
  END
END

REGION all
  COORDINATES
    0.0d0 0.0d0 0.0d0
    100.0d0 100.0d0 1.0d0
  /
END

REGION west_face
  FACE WEST
  COORDINATES
    0.d0 0.d0 0.d0
    0.d0 100.0d0 1.0d0
  /
END

REGION east_face
  FACE EAST
  COORDINATES
    100.0d0 0.d0 0.d0
    100.0d0 100.0d0 1.0d0
  /
END

REGION north_face
  FACE NORTH
  COORDINATES
    0.d0 100.0d0 0.d0
    100.0d0 100.0d0 1.0d0
  /
END

REGION south_face
  FACE SOUTH
  COORDINATES
    0.d0  0.d0  0.d0
    100.0d0 0.d0 1.0d0
  /
END

MATERIAL_PROPERTY plate
  ID 1
  POROSITY 0.20
  TORTUOSITY 1.d0
  ROCK_DENSITY 2.0E3
  HEAT_CAPACITY 100. J/kg-C
  THERMAL_CONDUCTIVITY_DRY 0.5787037 W/m-C
  THERMAL_CONDUCTIVITY_WET 0.5787037 W/m-C
  SATURATION_FUNCTION default
  PERMEABILITY
    PERM_X 1.d-14
    PERM_Y 1.d-14
    PERM_Z 1.d-14
  /
END

SATURATION_FUNCTION default
  SATURATION_FUNCTION_TYPE VAN_GENUCHTEN
  RESIDUAL_SATURATION 0.5d-1
  LAMBDA 0.75
  ALPHA 1.d-3
END

EOS WATER
  DENSITY EXPONENTIAL 1000. 101325. 1.0d-9
  VISCOSITY CONSTANT 1.728d-3 Pa-s
END

FLUID_PROPERTY
  DIFFUSION_COEFFICIENT 1.d-9
END

STRATA
  REGION all
  MATERIAL plate
END

TIME
  FINAL_TIME 0.10 day
  INITIAL_TIMESTEP_SIZE 1.d-4 day
  MAXIMUM_TIMESTEP_SIZE 0.0001 day at 0.d0 day
END

OUTPUT
  SNAPSHOT_FILE 
    #NO_PRINT_INITIAL
    TIMES day 0.01 0.04 0.06
    #NO_PRINT_FINAL
    FORMAT HDF5
  /
END

DATASET pressure_initial
  HDF5_DATASET_NAME initial
  FILENAME ../dataset.h5
END

FLOW_CONDITION initial
  TYPE
    TEMPERATURE dirichlet
    PRESSURE dirichlet
  /
  TEMPERATURE 25.0 C
  PRESSURE DATASET pressure_initial
END

FLOW_CONDITION west_face
  TYPE
    TEMPERATURE dirichlet
    PRESSURE dirichlet
  /
  TEMPERATURE 25.0 C
  PRESSURE 101325. Pa
END

FLOW_CONDITION east_face
  TYPE
    TEMPERATURE dirichlet
    PRESSURE dirichlet
  /
  TEMPERATURE 25.0 C
  PRESSURE 101325. Pa
END

FLOW_CONDITION north_face
  TYPE
    ENERGY_FLUX neumann
    FLUX neumann
  /
  ENERGY_FLUX 0.d0 W/m^2
  FLUX 0.d0 m/yr
END

FLOW_CONDITION south_face
  TYPE
    ENERGY_FLUX neumann
    FLUX neumann
  /
  ENERGY_FLUX 0.d0 W/m^2
  FLUX 0.d0 m/yr
END

INITIAL_CONDITION initial
  REGION all
  FLOW_CONDITION initial
END

BOUNDARY_CONDITION west_face
  REGION west_face
  FLOW_CONDITION west_face
END

BOUNDARY_CONDITION east_face
  FLOW_CONDITION east_face
  REGION east_face
END

BOUNDARY_CONDITION north_face
  REGION north_face
  FLOW_CONDITION north_face
END

BOUNDARY_CONDITION south_face
  FLOW_CONDITION south_face
  REGION south_face
END

END_SUBSURFACE

The PFLOTRAN Input File (RICHARDS Mode)

The RICHARDS Mode PFLOTRAN input file can be downloaded here.

 
SIMULATION
  SIMULATION_TYPE SUBSURFACE
  PROCESS_MODELS
    SUBSURFACE_FLOW flow
      MODE RICHARDS
    /
  /
END

SUBSURFACE

GRID
  TYPE structured
  NXYZ 40 40 1
  DXYZ
   2.5d0
   2.5d0
   1.0d0
  END
END

REGION all
  COORDINATES
    0.0d0 0.0d0 0.0d0
    100.0d0 100.0d0 1.0d0
  /
END

REGION west_face
  FACE WEST
  COORDINATES
    0.d0 0.d0 0.d0
    0.d0 100.0d0 1.0d0
  /
END

REGION east_face
  FACE EAST
  COORDINATES
    100.0d0 0.d0 0.d0
    100.0d0 100.0d0 1.0d0
  /
END

REGION north_face
  FACE NORTH
  COORDINATES
    0.d0 100.0d0 0.d0
    100.0d0 100.0d0 1.0d0
  /
END

REGION south_face
  FACE SOUTH
  COORDINATES
    0.d0  0.d0  0.d0
    100.0d0 0.d0 1.0d0
  /
END

MATERIAL_PROPERTY plate
  ID 1
  POROSITY 0.20
  TORTUOSITY 1.d0
  ROCK_DENSITY 2.0E3
  CHARACTERISTIC_CURVES cc1
  PERMEABILITY
    PERM_X 1.d-14
    PERM_Y 1.d-14
    PERM_Z 1.d-14
  /
END

CHARACTERISTIC_CURVES cc1
  SATURATION_FUNCTION VAN_GENUCHTEN
    LIQUID_RESIDUAL_SATURATION 0.5d-1
    M 0.75
    ALPHA 1.d-3
  /
  PERMEABILITY_FUNCTION MUALEM
    LIQUID_RESIDUAL_SATURATION 0.1d0
    M 0.5d0
  /
  PERMEABILITY_FUNCTION MUALEM_VG_GAS
    LIQUID_RESIDUAL_SATURATION 0.1d0
    GAS_RESIDUAL_SATURATION 0.1d0
    M 0.5d0
  /
END

EOS WATER
  DENSITY EXPONENTIAL 1000. 101325. 1.0d-9
  VISCOSITY CONSTANT 1.728d-3 Pa-s
END

FLUID_PROPERTY
  DIFFUSION_COEFFICIENT 1.d-9
END

STRATA
  REGION all
  MATERIAL plate
END

TIME
  FINAL_TIME 0.10 day
  INITIAL_TIMESTEP_SIZE 1.d-4 day
  MAXIMUM_TIMESTEP_SIZE 0.0001 day at 0.d0 day
END

OUTPUT
  SNAPSHOT_FILE 
    #NO_PRINT_INITIAL
    TIMES day 0.01 0.04 0.06
    #NO_PRINT_FINAL
    FORMAT HDF5
  /
END

DATASET pressure_initial
  HDF5_DATASET_NAME initial
  FILENAME ../dataset.h5
END

FLOW_CONDITION initial
  TYPE
    PRESSURE dirichlet
  /
  PRESSURE DATASET pressure_initial
END

FLOW_CONDITION west_face
  TYPE
    PRESSURE dirichlet
  /
  PRESSURE 101325. Pa
END

FLOW_CONDITION east_face
  TYPE
    PRESSURE dirichlet
  /
  PRESSURE 101325. Pa
END

FLOW_CONDITION north_face
  TYPE
    FLUX neumann
  /
  FLUX 0.d0 m/yr
END

FLOW_CONDITION south_face
  TYPE
    FLUX neumann
  /
  FLUX 0.d0 m/yr
END

INITIAL_CONDITION initial
  REGION all
  FLOW_CONDITION initial
END

BOUNDARY_CONDITION west_face
  REGION west_face
  FLOW_CONDITION west_face
END

BOUNDARY_CONDITION east_face
  FLOW_CONDITION east_face
  REGION east_face
END

BOUNDARY_CONDITION north_face
  REGION north_face
  FLOW_CONDITION north_face
END

BOUNDARY_CONDITION south_face
  FLOW_CONDITION south_face
  REGION south_face
END

END_SUBSURFACE

The Dataset

The hdf5 dataset required to define the initial/boundary conditions is created with the following python script called create_dataset.py:

import sys
from h5py import *
import numpy as np
import math
import matplotlib.pyplot as plt

filename = 'dataset.h5'
h5file = File(filename,mode='w')

p0 = 1.             # [MPa]
p_offset = 0.101325 # [MPa]
L = 100.            # [m]
dx = 1.0            # [m]
dy = 1.0            # [m]
dz = 1.0            # [m]

# 2D Surface:
# -----------------------------------------------------------------------------
# Pressure initial condition z=0 face
h5grp = h5file.create_group('initial')
h5grp.attrs['Dimension'] = np.string_('XY')
# Delta length between points [m]
h5grp.attrs['Discretization'] = [dx,dy]
# Location of origin
h5grp.attrs['Origin'] = [0.,0.]
# Load the dataset values
nx = L*dx + 1
ny = L*dy + 1
rarray = np.zeros((nx,ny),'=f8')

fx = np.zeros(nx,'=f8')
for i in range(int(nx)):
  if (0. <= i < (L/10.)):
    fx[i] = 0.
  if ((L/10.) <= i < (4.*L/10.)):
    fx[i] = (10./(3.*L))*float(i) - (1./3.)
  if ((4.*L/10.) <= i < (6.*L/10.)):
    fx[i] = 1.
  if ((6.*L/10.) <= i < (9.*L/10.)):
    fx[i] = 3. - (10./(3.*L))*float(i)
  if ((9.*L/10.) <= i < L):
    fx[i] = 0.
    
fy = np.zeros(ny,'=f8')
for j in range(int(ny)):
  if (0. <= j < (L/10.)):
    fy[j] = 0.
  if ((L/10.) <= j < (4.*L/10.)):
    fy[j] = (10./(3.*L))*float(j) - (1./3.)
  if ((4.*L/10.) <= j < (6.*L/10.)):
    fy[j] = 1.
  if ((6.*L/10.) <= j < (9.*L/10.)):
    fy[j] = 3. - (10./(3.*L))*float(j)
  if ((9.*L/10.) <= j < L):
    fy[j] = 0.

for i in range(int(nx)):
  for j in range(int(ny)):
    rarray[i][j] = ( p0*fx[i]*fy[j] + p_offset ) * 1.0e6  # [Pa]
    
h5dset = h5grp.create_dataset('Data', data=rarray)



The Python Script

def flow_transient_2D_BC1st2ndkind(path,input_prefix,remove,screen_on,pf_exe,
				                             mpi_exe,num_tries):
#==============================================================================#
# Based on:
# Kolditz, et al. (2015) Thermo-Hydro-Mechanical-Chemical Processes in
# Fractured Porous Media: Modelling and Benchmarking, Closed Form Solutions,
# Springer International Publishing, Switzerland.
# Section 2.2.10, pg.37
# "A Transient 2D Pressure Distribution, Non-Zero Initial Pressure,
# Boundary Conditions of 1st and 2nd Kind"
#
# Author: Jenn Frederick
# Date: 09/09/2016
# ******************************************************************************
  nxyz = np.zeros(3) + 1
  dxyz = np.zeros(3) + 1.
  lxyz = np.zeros(3) + 1.
  error_analysis = np.zeros(num_tries)
  dxyz_record = np.zeros((num_tries,3))
  test_pass = False
  try_count = 0
  
  # initial discretization values
  lxyz[0] = 100.        # [m] lx
  lxyz[1] = 100.        # [m] ly
  nxyz[0] = 10          # [m] nx 
  nxyz[1] = 10          # [m] ny 
  dxyz = lxyz/nxyz      # [m]
  k = 1.0e-14           # [m2]
  mu = 1.728e-3         # [Pa-s]
  por = 0.20
  K = 1.0e-9            # [1/Pa]
  p_offset = 0.101325   # [MPa]
  p0 = 1.               # [MPa]
  chi = k/(por*K*mu)
  ierr = 0
  
  while (not test_pass) and (try_count < num_tries): 
    print_discretization(lxyz,nxyz,dxyz)
    nx = int(nxyz[0]); ny = int(nxyz[1]); nz = int(nxyz[2])
    dx = dxyz[0]; dy = dxyz[1]; dz = dxyz[2]
    Lx = lxyz[0]; Ly = lxyz[1]; Lz = lxyz[2]
    try_count = try_count + 1

    p_soln = np.zeros((4,nx,ny))                     # [MPa]
    p_pflotran = np.zeros((4,nx,ny))                 # [Pa]
    t_soln = np.array([0.0,0.04,0.06,0.10])          # [day]
    x_soln = np.linspace(0.+(dx/2.),Lx-(dx/2.),nx)   # [m]
    y_soln = np.linspace(0.+(dy/2.),Ly-(dy/2.),ny)   # [m]
    p1x = np.zeros(int(nx))
    p2y = np.zeros(int(ny))

    # for time=0, use the actual initial condition as the analytical
    # solution, not the analytical solution with t=0:
    fx = np.zeros(int(nx))
    for i in range(int(nx)):
      x = x_soln[i]
      if (0. <= x < (Lx/10.)):
	fx[i] = 0.
      if ((Lx/10.) <= x < (4.*Lx/10.)):
	fx[i] = (10./(3.*Lx))*(x) - (1./3.)
      if ((4.*Lx/10.) <= x < (6.*Lx/10.)):
	fx[i] = 1.
      if ((6.*Lx/10.) <= x < (9.*Lx/10.)):
	fx[i] = 3. - (10./(3.*Lx))*(x)
      if ((9.*Lx/10.) <= x < Lx):
	fx[i] = 0.
      
    fy = np.zeros(int(ny))
    for j in range(int(ny)):
      y = y_soln[j]
      if (0. <= y < (Ly/10.)):
	fy[j] = 0.
      if ((Ly/10.) <= y < (4.*Ly/10.)):
	fy[j] = (10./(3.*Ly))*(y) - (1./3.)
      if ((4.*Ly/10.) <= y < (6.*Ly/10.)):
	fy[j] = 1.
      if ((6.*Ly/10.) <= y < (9.*Ly/10.)):
	fy[j] = 3. - (10./(3.*Ly))*(y)
      if ((9.*Ly/10.) <= y < Ly):
	fy[j] = 0.

    p_soln_t0 = np.zeros((nx,ny))
    for i in range(int(nx)):
      for j in range(int(ny)):
	p_soln_t0[i][j] = ( p0*fx[i]*fy[j] + p_offset )  # [MPa]

    # create the analytical solution for all other times > 0
    for time in range(4):
      t = t_soln[time]*24.0*3600.0  # [sec]
      # create p1y
      sum_term_y = np.zeros(int(ny))
      sum_term_old_y = np.zeros(int(ny))
      n = 1
      epsilon = 1
      #while n < 5001:
      while epsilon > epsilon_value:
	sum_term_old_y = sum_term_y
	sum_term_y = sum_term_old_y + (np.cos(n*math.pi*y_soln/Ly)*np.exp(-chi*pow(n,2)*pow(math.pi,2)*t/pow(Ly,2))*(80./(3.*pow((n*math.pi),2)))*np.cos(n*math.pi/2.)*np.sin(n*math.pi/4.)*np.sin(3.*n*math.pi/20.))
	epsilon = np.max(np.abs(sum_term_old_y-sum_term_y))
	n = n + 1
      p2y = 0.5 + sum_term_y
      # create p1x
      sum_term_x = np.zeros(int(nx))
      sum_term_old_x = np.zeros(int(nx))
      n = 1
      epsilon = 1
      while epsilon > epsilon_value:
	sum_term_old_x = sum_term_x
	sum_term_x = sum_term_old_x + (np.sin(n*math.pi*x_soln/Lx)*np.exp(-chi*pow(n,2)*pow(math.pi,2)*t/pow(Lx,2))*(80./(3.*pow((n*math.pi),2)))*np.sin(n*math.pi/2.)*np.sin(n*math.pi/4.)*np.sin(3.*n*math.pi/20.))
	epsilon = np.max(np.abs(sum_term_old_x-sum_term_x))
	n = n + 1
      p1x = sum_term_x
      for i in range(int(nx)):
	for j in range(int(ny)):
	  p_soln[time,i,j] = ( p0*p1x[i]*p2y[j] + p_offset ) # [MPa]

    p_soln[0,:,:] = p_soln_t0[:,:]
    
    # run PFLOTRAN simulation
    run_pflotran(input_prefix,nxyz,dxyz,lxyz,remove,screen_on,pf_exe,mpi_exe)
    
    # load data from HDF5
    index_string = 'Time:  0.00000E+00 d/Liquid_Pressure [Pa]'
    p_pflotran[0,:,:] = read_pflotran_output_2D(path+
		 '/2D_transient_pressure_BC_1st_2nd_kind.h5',index_string,False)
    index_string = 'Time:  4.00000E-02 d/Liquid_Pressure [Pa]'
    p_pflotran[1,:,:] = read_pflotran_output_2D(path+
		 '/2D_transient_pressure_BC_1st_2nd_kind.h5',index_string,False)
    index_string = 'Time:  6.00000E-02 d/Liquid_Pressure [Pa]'
    p_pflotran[2,:,:] = read_pflotran_output_2D(path+
		 '/2D_transient_pressure_BC_1st_2nd_kind.h5',index_string,False)
    index_string = 'Time:  1.00000E-01 d/Liquid_Pressure [Pa]'
    p_pflotran[3,:,:] = read_pflotran_output_2D(path+
		'/2D_transient_pressure_BC_1st_2nd_kind.h5',index_string,remove)
    ierr = check(p_pflotran)

    # Convert units [Pa -> MPa]
    p_pflotran = p_pflotran/1.e6

    max_percent_error = calc_relative_error(p_soln,p_pflotran,ierr)
    record_error(error_analysis,dxyz_record,max_percent_error,dxyz,try_count)
    test_pass = does_pass(max_percent_error,try_count,num_tries)
    nxyz[0] = nxyz[0]*2.
    nxyz[1] = nxyz[1]*2.
    dxyz = lxyz/nxyz

  # Plot the PFLOTRAN and analytical solution:
  plot_2D_transient(path,t_soln,'day',x_soln,y_soln,p_soln,p_pflotran,
                    'Distance [m]','Pressure [MPa]',"{0:.2f}".format(max_percent_error))
  
  # Plot error analysis
  plot_error(error_analysis,dxyz_record,path,try_count)
  
  # Add test result to report card
  add_to_report(path,test_pass,max_percent_error,ierr)

  return;

Refer to section Python Helper Functions for documentation on the qa_tests_helper module, which defines the helper functions used in the Python script above.