 Welcome to lecture number 39 of this ground water hydrology course. So, in last lecture, we have seen that regional scale development of ground water. In that one, one important aspect is conjunctive use of surface water and ground water. In this particular lecture, we will try to cover that aspect. So, our main topic is modeling and management of ground water and under that topic to be covered is conjunctive surface, subsurface modeling of overland flow including flow through vados. So, we already know that if this is our ground surface and this is our water table and this level is that our lower portion or the lower most portion of the phreatic aquifer, then this zone from starting from the water table to water table to this bed level or if it is of rocky nature, then we can say that this is impermeable rock and this impermeable starting from this impermeable rock to this water table level, we can say that this is our zone of saturation and above that there is a thin region which is called as capillary zone and near to surface there exists another region which is called as soil water zone. This is important for plants and in between we have intermediate zone or intermediate or vados zone. Sometimes this total region starting from water table to ground surface, this is called as zone of aeration or unsaturated zone, vados zone. So, in this particular lecture we are concerned about the modeling of this zone of aeration or unsaturated zone or vados zone and the water above that ground level or ground surface. So, what are the processes or what is the interaction mechanism between this overland flow and subsurface flow that we will examine. So, to model this particular problem we need certain governing equations for flow representation. So, for mathematical modeling we need governing equations for surface water. Next we need governing equation for subsurface water or in this particular case this is for vados zone, zone flow. So, the whole problem can be represented as this is let us say we have ground surface which is not always straight, it is always having some slope. So, let us say this is x direction and this is x direction which is one dimensional thing and which is in line with the or in the same direction of the our bed slope and we have vertical direction is z direction. This is for infiltration, infiltration from our overland flow. So, infiltration then we have subsurface flow 2D model and our surface flow is modeled as 1D. So, this is our surface flow this is 1D model and this is our surface water level and we have rainfall above that. So, starting from rainfall we have subsurface flow and this is our ground surface below we have this infiltration and due to this infiltration only there will be subsurface flow that is 2D or 2 dimensional in nature. So, considering prismatic channel we can write our governing equation. So, we can write this surface flow equations. So, surface flow equations first assumption is that the surface flow occurs in prismatic channel of rectangular cross section and flow second is flow is 1D shallow water flow. So, in this case we can write our governing equation for surface water flow as del Q by del T where U is vector plus del F is vector by del X and S is another vector or we can say that these are column matrices. So, in this case U is basically H and Q then we have this F which is Q and Q square by H G H square by 2 and S this is R minus I G H S naught minus SF. So, term by term we need to define different variables. So, H is flow depth in meters or any standard unit Q is discharge per unit width meter square per second G is acceleration due to gravity or is volumetric rate of rainfall per unit surface area this is meter per second then I is volumetric infiltration rate per unit area bottom slope SF is friction slope. So, these are the key parameters. Now, energy slope or friction slope can be calculated using Darcy Weisbach formula. So, friction slope SF can be calculated using the formula SF equals to FD Q square 8 G H cube where FD is the frictional resistance coefficient and this depends on instantaneous state of flow. So, flow is laminar by assumption. So, we can write this FD equals to C L by R E where this R E is Reynolds number C L is the constant which depends on rainfall intensity FD depends on rainfall intensity and Reynolds number. And this Reynolds number can be calculated as Q by nu nu is kinematic viscosity of water. So, this is all about the equations related to surface flow. Now, subsurface flow equations are subsurface flow equations are one is related to the flow that is first assumption is 2D flow, second assumption is transient flow in an isotropic porous medium. So, with the conservation principles we can write this equation without any source or sink as del theta by del t del v x del x and this v z z equals to 0 where this theta is volumetric moisture content v x and v z these are Darcy flow velocities Darcy flow velocity in x direction and this is Darcy flow velocity in y direction. So, in this case x and z are distances along 2 coordinate directions and z is taken as positive downwards. So, for Darcy's law the velocity components can be calculated as this is for Darcy's law this our v x can be calculated as minus k psi del psi del x and v z can be calculated as minus k psi del psi by del z minus 1. This psi is pressure head psi is pressure and k psi this is unsaturated hydraulic conductivity. So, by substituting these two expressions in our original equation for conservation of mass this is from conservation of mass we can get our final form of the equation that is del theta by del t equal to del x k psi del psi x and z this is k psi del psi by del z minus 1. So, this particular equation is called as mixed form mixed form equation because it includes both theta and psi in single equation. So, it is important to have the expression for k psi which is unsaturated hydraulic conductivity and this psi k relationship dictates the flow regime and this moisture content and theta there is a relationship and interestingly this k psi or psi k relationship and psi theta relationship these are not unique not unique. So, it is always problem dependent and it is related to soil type and it varies with problem to problem. To solve these two equations we can start discretizing the governing equations or for our surface flow and subsurface flow. To have this first we need to have some kind of conceptual definition of surface flow and subsurface flow. So, we can have situation like this where our this is x direction, this is z direction although these are not having 90 degree, but approximately these are considered to be 90 degree situations then we have regions like this where this is ith cell, this is i minus 1, this is i plus 1 and from here we have 1D flow, 1D subsurface flow, surface flow and above that we have water level and 1D this is surface flow and this is our ground surface. So, with this picture in our mind we can start discretizing our governing equations. So, first is surface flow. So, in surface flow this is can be solved using some kind of predictor corrector method for surface flow. So, our governing equation is del u del t plus del f del x equals to s. So, for that the predictor step can be written using n and n and star level n is known time level, this is predictor level the star is basically predictor level. So, with this two levels we can write this u i star equal to u i n minus del t by del x and this is f i plus half nth level and this is f i minus half this is again at nth level and this is del t s i this is also calculated at nth level. So, this represents this f at i plus half represents numerical flux through the self face between nodes i plus 1 and i and del x is basically grid spacing all terms in the right hand side of this particular equation are at known time level. Therefore, u star can be computed directly or explicitly. So, we can write it as explicit set equation. So, if f i plus half can be calculated as f i plus half equal to half f r plus f l minus u r minus u l. So, what is this? This is basically alpha is positive coefficient and f r is the flux computed using information from the right hand side of the self face. So, f r is basically flux computed with the information from right hand side and this is with the information from left hand side of the self face and u l and u r are obtained using a particular this u l and u r these are calculated based on a particular expression that is called as m u scl approach m u scl is monotone upwind scheme for conservation law. So, under this this u l is calculated as u i plus del u i and half and u r is calculated as u i plus 1 plus half this is minus half del u i plus 1. There are several ways of calculating these two of these two variables. We can calculate it using a mean mode approach that is called mean mode approach. In that one this del u i is calculated as mean mode of u i plus 1 minus u i and this is u i minus u minus 1. What is this mean mode? Mean mode is defined as a b where a if modulus of a is less than b modulus of b and a b is greater than 0, b if modulus of b is less than modulus of a and a b greater than 0 and 0 if a b less than equals to 0 or negative. This positive coefficient alpha in our previous expression is determined as alpha should be greater than maximum of modulus of this q i divided by h i plus root over g h i and this is for all i this is for all i that means considering all i's that is i starting from 1 to n 1 to n, n is the number of grid points we can calculate the minimum value of alpha. So, alpha should be greater than this level. So, next step is character step. So, in this character step the start values or predicted values are utilized to get the corrected value. So, in this one it involves n plus 1 level and star level star level is already known level from our predicted step and n level is unknown level. So, in this one this u i n plus 1 this can be written as u i nth plus u i star minus del t del x this is f i plus half f i minus half star plus del t into s i this is at star level where this f i plus half is calculated as f r star plus f l star minus alpha u r star minus u l star. So, u l and u r this can be calculated using that mean mode approach we have utilized in our predictor step and we can solve this particular equation again it is explicit equation. So, we can directly solve this. So, in this particular equation only the source term is evaluated using predicted value of h and q. So, these two steps we can calculate using explicit equations and finally, we can get n plus 1 level value for this one. So, what are the initial conditions? At t equals to 0 initial boundary conditions. So, at t equals to 0 we need to consider one h initial depth and q initial although it will have impact on the final solution, but we need to specify to remove the discontinuity in the initial condition. So, although it will have impact on the final results this can eliminate the numerical singularity of the solutions. So, at this level we have i 2 to n minus 2 n minus 1 we have internal nodes and for the first node and the last node we can specify the boundary condition. So, values of the variables at upstream and downstream ends domain are using determined by appropriate values. So, discharge at the upstream end is equal to 0. So, discharge at upstream is equal to 0. However, one initial discharge is specified to remove the singularity, but it should be small sufficiently small compared to the other values. The flow depth at upstream we can say that flow depth at upstream level can be determined using negative characteristics equations or MOC or method of characteristics equation and downstream also this MOC can be utilized to get the levels for h. For stability the current number condition or CFL condition this G by h plus root over gh that should be less than 1 the Cn is called current number and in this particular problem delta T and delta X should be such that this G and h values this should be less than 1. We can also say that this is actually maximum value for all i this is g i or this is actually q q i h i plus root over g into h i less than equals to 1. For subsurface flow again we need to discretize those equations subsurface flow. So, in this one one way is that we can discretize our original equation that is theta ij because it is a two dimensional thing. So, theta ij at nth level this is basically V i plus half j bar represents time averaged values and V i minus half j this is there plus we have ij plus half V i j minus half bar divided by del z equal to 0. So, this V bar is basically w into V n plus 1 1 minus w V n and this is time waiting factor and V i j plus half this is calculated using our expression of Darcy's law this is j and psi is ij plus half psi ij and this is minus del z and divided by del z. So, k ij plus half is the unsaturated hydraulic conductivity and it should be evaluated at the inter block phases between ij plus half and ij. So, this is basically calculated based on ij plus half and ij. So, if you discretize our original equation in terms of this one that our V ij plus half and if we utilize this in this equation and we can replace it here then we will get the final equation as del t del x square minus k i plus half n plus 1 psi i plus half i plus 1 n plus 1 psi ij n plus 1 plus k i minus half j plus n psi ij n plus 1 psi i minus 1 n plus 1 j and next one is related to this z which is z square minus k ij plus half this is i plus half j i minus half j this is ij plus half and this is at psi j plus half n plus 1 minus psi ij at n plus 1 level del z plus k ij minus half both are at n plus 1 level psi ij n plus half minus psi ij minus half n plus 1 minus del z. So, with this another term will be there that is theta ij n plus 1 minus theta ij n 1 minus w is del t by del x this is V i plus half j n V i minus half j n minus 1 minus w this is del z V i j plus half n minus V i j minus half n equal to 0. So, in this one the inter block hydraulic conductivity calculation is important. So, that is calculated as gamma k psi ij plus 1 minus gamma k psi ij plus 1. So, gamma is weight coefficient. So, now we can have boundary conditions. So, we can have flux type or pressure head type boundary condition. So, with this two boundary condition we can model that and there will be interaction at our with our concept that is for each cell there will be interaction between the surface and ground water we need to solve our surface water subsurface flow at time level n. Then surface flow calculation and we need to repeat this process to get the convergence and we need to proceed for future time levels. This subsurface flow equation can be solved using Newton Rapson technique. So, this is the method for solving our surface water ground water interaction, but we have reached up to the zone of aeration or unsaturated zone. So, we need to consider the saturated interaction with the saturated zone. So, maybe in the next class or next lecture we will discuss that aspect. Thank you.