UNESCO-IHP Water Programme for Environmental Sustainability - WPA II
The CATchment HYdrology (CATHY) Model Working Group: Claudio Paniconi and Mauro Sulis INRS-ETE, University of Quebec, Canada Mario Putti and Matteo Camporese University of Padova, Italy Stefano Orlandini, Giovanni Moretti, Maurizio Cingi and Alice Cusi University of Modena & Reggio Emilia, Italy First “Climate Change...” Project Meeting, May 29 – June 4, 2009, Brazil
Context
Surface and subsurface waters are not isolated components of the hydrologic system, but interact in response to topographic, soil, geologic, and climatic factors.
Groundwater seepage can take the dual role of source and sink: • Groundwater may seep out from shallow aquifer and feed networks of irregular channels. • Surface water may seep back into the ground, possibly depleting streams until they run dry.
Context
The interaction of groundwater and surface flow is a key focus of interest for: • Water resources management (under the effects of climate and land use and demographic changes) • Water quality (e.g., role hyporheic fluxes in aquatic habitats) • Geomorphology (e.g., action of groundwater seepage on channel initiation)
A number of modeling tools and approaches have been developed for studying coupled surface water-groundwater systems: Three-dimensional equation for variably saturated subsurface flow, i.e., Richards equation, coupled with a one- or two-dimensional approximation of the Saint-Venant equations for overland and channel flow, represent the current state-of-the-art in catchment-aquifer models.
Outline 1. Model description: •
Mathematical formulation
•
Numerical discretization
•
Surface–subsurface interactions
•
Surface flow conceptualization
2. Test cases and applications
3. Data Requirements
Model description: Mathematical formulation
∂ψ σ ( Sw ) = ∇ ⋅ [ K s K rw ( S w )(∇ψ + η z ) ] + qs ( h ) ∂t
∂Q ∂Q ∂ 2Q + ck = Dh 2 + ck qL ( h,ψ ) ∂t ∂s ∂s σ Sw
θ θs
Ss
φ ψ
t Ks Krw
ηz
general storage term [1/L]: σ = SwSs + φ(dSw/dψ) water saturation = θ/θs [/] volumetric moisture content [L3/L3] saturated moisture content [L3/L3] specific storage [1/L] porosity (= θs if no swelling/shrinking) pressure head [L] time [T] saturated conductivity tensor [L/T] relative hydraulic conductivity [/] zero in x and y and 1 in z direction
z qs h s Q ck Dh qL
(1)
(2)
vertical coordinate +ve upward [L] subsurface equation coupling term (more generally, source/sink term) [L3/L3T] ponding head (depth of water on surface of each cell) [L] hillslope/channel link coordinate [L] discharge along s [L3/T] kinematic wave celerity [L/T] hydraulic diffusivity [L2/T] surface equation coupling term (overland flow rate) [L3/LT]
(1) Paniconi & Wood, Water Resour. Res., 29(6), 1993 ; Paniconi & Putti, Water Resour. Res., 30(12), 1994 (2) Orlandini & Rosso, J. Hydrologic Engrg., ASCE, 1(3), 1996 ; Orlandini & Rosso, Water Resour. Res., 34(8), 1998 (1)+(2) Bixio et al., CMWR Proceedings, 2000 ; Putti & Paniconi, CMWR Proceedings, 2004
Model description: Numerical discretization Surface: • PDE of the kinematic wave solved by a finite difference (FD) scheme • Numerical dispersion arising from the truncation error of the scheme is used to simulate the physical dispersion • Unconditional stability reached by matching numerical and physical diffusivities through the temporal weighting factor used to discretize the kinematic wave model
Model description: Numerical discretization Subsurface: • PDE solved by a three-dimensional finite element (FE) spatial integrator and by a weighted finite difference (FD), i.e. Euler or Crank-Nicolson, scheme in time • Nonlinearity arising from the storage σ(Sw) and conductivity Krw(Sw) terms are handled via a Picard or Newton linearization scheme • Time varying boundary conditions: prescribed head (Dirichlet type) or flux (Neumann type), atmospheric fluxes, source/sink terms, and seepage faces
Model description: Surface–subsurface interactions The coupling between the land surface and the subsurface is handled by an automatic boundary condition (BC) switching algorithm acting on the source/sink terms qs(h) and qL(h,ψ). The coupling term is computed as the balance between atmospheric forcing (rainfall and potential evaporation) and the amount of water that can actually infiltrate or exfiltrate the soil. The switching check is done surface node by surface node in order to account for soil and topographic variability. The switching check is done at each time the surface equation is solved (according to the values of ponding heads at the surface) and at each subsurface time or iteration.
Model description: Surface–subsurface interactions
Unified Flow Direction Algorithm The mathematical formulation implemented is based on the use of: y
Triangular facets introduced by Tarboton (1997, WRR)
y
Path-based analysis introduced by Orlandini et al. (2003, WRR)
y
Plan curvature computation introduced by Zevenbergen and Thorne (1987, ESPL)
Unified Flow Direction Algorithm
Orlandini et al. (2003, WRR)
Orlandini and Moretti (2009, JGR)
Unified Flow Direction Algorithm
Orlandini and Moretti (2009, WRR)
Validation Using Contour Elevation Data
Moretti and Orlandini (2008, WRR)
Drainage Basin and Drainage Slope
Orlandini and Moretti (2009, WRR)
Prediction of the Drainage Network
Hydraulic Geometry (Leopold and Maddock, 1953)
W = a Qb
Ym = c Q f
U = kQ
m
Sf = tQ
z
kS = r Q y
Parameterization of Stream Channel Geometry (Channels and Hillslope Rivulets) b ′ W =aQ
′
( at-a-station relationship )
W = a′′ Q bf′′ ( downstream relationship ) Q f = u Aw ( given frequency discharge, bankfull discharge )
W = W ( A,1) Q b′ W ( A,1) = W ( As , Q f ) Q f ( As )
Orlandini and Rosso (1998, WRR)
− b′
(A
As )
w( b′′− b′)
Parameterization of Conductance Coefficients (Channels and Hillslope Rivulets)
k S = r′ Q y′ ( at-a-station relationship ) k S = r′′ Q fy′′ ( downstream relationship ) Q f = u Aw ( given frequency discharge, bankfull discharge ) k S = k S ( A,1) Q
y′
k S ( A,1) = k S ( As , Q f ) Q f ( As ) Orlandini (2002, WRR)
− y′
(A
As )
w( y ′′− y ′ )
Diffusion Wave Modeling: Mathematical Model
Kinematic wave model
∂Q ∂Q + ck = ck qL ∂t ∂s ck =
dQ dΩ
Diffusion wave model
∂Q ∂Q ∂ 2Q + ck = Dh 2 + ck qL ∂t ∂s ∂s
∂S f ck = − ∂Ω
∂S f ∂Q
⎛ W ∂S f ⎞ Dh = 1 ⎜ ⎟ β cos ∂ Q ⎝ ⎠
Diffusion Wave Modeling: Parameterization of the Drainage System Gauckler-Manning-Strickler Equation
Q = k S W −2 3 Ω5 3 S f 1 2 Sf = Incorporating the variation of stream channel geometry
Sf =
Q2 kS W (Q ) 2
−4 3
Ω
10 3
( kS
= 1 n ,W ≈ P )
Q2 k S 2 W −4 3 Ω10 3 Incorporating the variation of conductance coefficient and stream channel geometry
Sf =
Q2 kS (Q ) W (Q ) 2
−4 3
Ω10 3
Diffusion Wave Modeling: Constitutive Equations
Flow Rating Curve
Q = k S ( A,1)
1 1− y ′+ 2 3b′
W ( A,1)
2 − 3 (1− y ′+ 2 3b′)
Ω
5 3 (1− y ′+ 2 3b′)
Sf
1 2 (1− y ′+ 2 3b′)
Kinematic Celerity 2⎛
1 ⎞
3 (1− 2 y ′+ 4 3b′) 2 (1− 2 y ′+ 4 3b′) 3 (1− 2 y ′+ 4 3b′) ⎜ 1+ y ′− b′ ⎟ 5 1 − 5⎝ 2 3 ⎠ ′ ′ ′ ′ ′ ′ y b y b y b − + − + − + 5 1 2 3 5 1 2 3 10 1 2 3 ( ) ( ) ( ) ck = k S ( A,1) W ( A,1) S f ( A,1) Q 3 (1 − y ′ + 2 3 b′)
Hydraulic Diffusivity
Dh =
Q
1− b ′
cos β
2 (1 − y ′ + 2 3 b′ ) W ( A,1) S f
1
Diffusion Wave Modeling: Muskingum-Cunge Method with Variable Parameters Qi +j +11 = C1 Qi j +1 + C2 Qi j + C3 Qi +j 1 + C4 qLij++11 C1 =
ck ( Δt Δs ) − 2 X 2 (1 − X ) + ck ( Δt Δs )
ck ( Δt Δs ) + 2 X C2 = 2 (1 − X ) + ck ( Δt Δs ) 2 (1 − X ) − ck ( Δt Δs ) C3 = 2 (1 − X ) + ck ( Δt Δs ) C4 =
2 ck Δt 2 (1 − X ) + ck ( Δt Δs )
Diffusion Wave Modeling: Muskingum-Cunge Method with Variable Parameters
Dn = ck Δs (1 2 − X )
Dn = Dh
⎛ W ∂S f ⎞ Dh = 1 ⎜ ⎟ ∂ cos Q β ⎝ ⎠
∂S f 1 W X= − 2 ck Δs cos β ∂Q
The Muskingum-Cunge method with variable parameters is: • Unconditionally stable (Dn = Dh). • Accurate for Courant numbers not too far from 1 (∆s ≈ ck ∆t). • Independent of structural parameters ∆s and ∆t. (Cunge, 1969, JHR; Ponce, 1986, JHE; Orlandini and Rosso, 1996, JHE)
Test cases and applications 1D sloping plane with homogeneous subsurface* Sloping plane catchment (400 m x 320 m) S0=0.0005 (slope) Δx=Δy=80 m (spatial discretization) Ts = 300 min (simulation time) Tr = 200 min (rainfall time) Td = 100 min (recession time) qr = 20 mm/h (rainfall intensity) Saturation excess (Dunnian process): Ksat=41.7 mm/h Water table heights: 0.5 and 1.0 m Excess infiltration (Hortonian process): Water table height: 1.0 m Ksat=0.417 and 4.17 mm/h * Results in agreement with Kollet & Maxwell, Adv. Water. Res., 29, 2006.
Test cases and applications: Saturation excess
Test cases and applications: Infiltration excess
Test cases and applications Application of CATHY model in a climate change frameork to the 720 km2 des Anglais catchment located in southern Quebec - Data provided by the Canadian Regional Climate Model (CRCM V4.2.0) - Precipitation (rain, snow) - Temperature (max, min) - Data provided for a 45-km horizontal grid-size mesh with a daily temporal resolution: - Model calibrated to distributed (water table) and aggregated (flow rate) data
Test cases and applications
Test cases and applications
Data Requirements: Catchment Properties Grid-based or contour-based digital elevation models Soil retention characteristics • residual and saturated water content • pore-size distribution index • saturated soil matrix potential • hydraulic conductivity Geologic stratigraphy, bedrock depth Land cover and use Channel network • cartographic blue lines • cross sections • roughness coefficients Man-made hydraulic structures
Data Requirements: Hydrologic Variables Precipitation • rainfall • occult precipitation • snow Evaporative demand (estimated using the Penman-Monteith equation) • solar radiation • air humidity • wind speed • air temperature • barometric pressure Streamflow Soil moisture Water table depth
