Philip Function

public function Philip(rain, psic, ksat, theta, thetas, thetar, psdi, itheta, dt, cuminf) result(inf)

calculates the actual rate of infiltration (m/s) according to Philip's equation. Use time compression approximation.

References:

Philip, J. R. 1954. An infiltration equation with physical significance. Soil Science, 77: 153-157.

Milly, P. C. D., An event-based simulation model of moisture and energy fluxes at a bare soil surface, Water Resour. Res., 22, 1680-! 692, 1986

Sivapalan, M., Beven, K., and Wood, E. F. (1987), On hydrologic similarity: 2 . A scaled model of storm runoff production, Water Resour. Res., 23( 12), 2266– 2278, doi:10.1029/WR023i012p02266.

Arguments

Type IntentOptional Attributes Name
real(kind=float), intent(in) :: rain

rainfall rate (m/s)
real(kind=float), intent(in) :: psic

bubbling pressure (m)
real(kind=float), intent(in) :: ksat

saturated hydraulic conductivity (m/s)
real(kind=float), intent(in) :: theta

volumetric water content (m3/m3)
real(kind=float), intent(in) :: thetas

saturated volumetric water content (m3/m3)
real(kind=float), intent(in) :: thetar

residual volumetric water content (m3/m3)
real(kind=float), intent(in) :: psdi

Brooks & Corey pore size distribution index (-)
real(kind=float), intent(in) :: itheta

initial soil moisture
integer(kind=short), intent(in) :: dt

time step (s)
real(kind=float), intent(inout) :: cuminf

cumulative infiltration (m)

Return Value real(kind=float)

infiltration rate (m/s)


Variables

Type Visibility Attributes Name Initial
real(kind=float), public :: bcgamm

2.0 + 3.0 * psdi
real(kind=float), public :: cc = 0.

gravity term in Philip's infiltration equation
real(kind=float), public :: deltrz

difference between soil moisture at beginning of storm event and residual soil moisture
real(kind=float), public :: infcap

infiltration capacity (m/s)
real(kind=float), public :: sorp

sorptivity

Source Code

FUNCTION Philip &
!
(rain, psic, ksat, theta, thetas, thetar, psdi, itheta, dt, cuminf) &
!
RESULT (inf)


IMPLICIT NONE

!Arguments with intent in
REAL (KIND = float), INTENT(IN) :: rain !!rainfall rate (m/s)
REAL (KIND = float), INTENT(IN) :: psic !!bubbling pressure (m)
REAL (KIND = float), INTENT(IN) :: ksat !!saturated hydraulic conductivity (m/s)
REAL (KIND = float), INTENT(IN) :: theta !!volumetric water content (m3/m3)
REAL (KIND = float), INTENT(IN) :: thetas !!saturated volumetric water content (m3/m3)
REAL (KIND = float), INTENT(IN) :: thetar !!residual volumetric water content (m3/m3)
REAL (KIND = float), INTENT(IN) :: psdi !!Brooks & Corey pore size distribution index (-)
INTEGER (KIND = short), INTENT(IN) :: dt !!time step (s)
REAL (KIND = float), INTENT(IN)  :: itheta !!initial soil moisture

!Arguments with intent inout:
REAL (KIND = float), INTENT(INOUT) :: cuminf !!cumulative infiltration (m)

!local declarations:
REAL (KIND = float)  :: inf !!infiltration rate (m/s)
REAL (KIND = float)  :: sorp  !!sorptivity
REAL (KIND = float)  :: cc = 0. !!gravity term in Philip's infiltration equation
REAL (KIND = float)  :: deltrz  !! difference between soil moisture at beginning
                                !! of storm event and residual soil moisture
REAL (KIND = float) :: bcgamm !!2.0 + 3.0 * psdi
REAL (KIND = float) :: infcap !!infiltration capacity (m/s)


!-------------------------end of declarations----------------------------------

  ! Sivapalan et al. (1987)
  bcgamm = 2.0 + 3.0 * psdi
  sorp = (((2. * ksat *( (thetas-itheta)**2.) * psic ) / &
         (thetas-thetar))*((1./(bcgamm + 0.5 * psdi - 1.)) + &
         ((thetas-thetar)/ (thetas-itheta) ) ) ) ** 0.5

  deltrz = itheta - thetar
  IF (deltrz <= 0.) deltrz = 0.
  cc = 0.5 * (1. + ( (deltrz / (thetas - thetar))**(bcgamm / psdi) ) )


!If soil is saturated then set infiltration to zero
IF (theta >= thetas) THEN
  inf = 0.
  
!If surface is unsaturated then calculate infiltration  
ELSE
  !If first step with infiltration then all precipitation is
  !infiltrated to avoid division by zero
  IF (cuminf == 0.) THEN
    inf = rain 

    
  !calculate infiltration capacity as a function of cumulative
  ! infiltration for all other time steps of storm (Milly, 1986)  
  ELSE
     infcap = cc * ksat * (1. + (1./ ((( 1. + (( 4. * cc * ksat * cuminf ) / &
              (sorp**2.)))**0.5) - 1.))) 
     !Take the actual infiltration rate as the minimum of the
     ! precipitation rate or the infiltration capacity.
     inf = MIN (infcap,rain)
  END IF
  
  
  !update cumulative infiltration
  cuminf = cuminf + inf * dt

END IF

RETURN

END FUNCTION Philip