 Welcome to this lecture number 36. In this particular lecture, we will cover this modeling and management of groundwater. And under this, so our main topic is modeling and management of groundwater and under this topics to be covered are groundwater management model that is for confined and unconfined aquifer. And this linked simulation optimization also will cover that meta model based in the last lecture or lecture number 35. We have talked about groundwater management model for confined aquifers that is basically confined aquifers with one dimensional flow situation under steady state condition. In that one, we have one linear objective function and we have discretized our governing equation using finite difference method and we have utilized those equations as constraints for our optimization model. And finally, the total model or total formulation can be solved using linear programming because all the constraints and objectives are linear in nature. So, today we will start with the unconfined aquifer first and we will try to solve that unconfined confined aquifer flow aquifer flow management model. This is also applicable to one dimensional flow situation. So, steady state groundwater equation for unconfined aquifer can be written as, so this is for steady state flow situation where x is the flow direction, then Tx and this is del H by del x equals to w. So, in this one Tx can be written as Kh, K is the hydraulic conductivity and H is the hydraulic head. So, with this expression, if we replace it in our original equation, we can write the whole equation as Kh del H by del x equals to w or we can write it as del H by del x half K and this is basically H square x square. So, for homogeneous type of aquifer, we can take this K out of this derivative thing, homogeneous uniform condition, we can write it as K and we can transfer the two on the other side. So, del 2 H 2 divided by this is del x 2 and we can write it as 2w. So, from this one final equation will be our del 2 H 2 del x 2 or 2w K or simply we can write this as, because H is only a function of x. So, we can write it as our ordinary differential equation. Now, we can apply our finite difference discretization to this left hand portion, but the problem is that it has got this H square term. So, what we can do, we can substitute another secondary variable with let us say this is small w, we can write it as H square. So, with this replacement our equation will become a linear equation. So, with this substitution our final equation that will look like d 2 w, this is dx 2 and 2w by K and finally, if we apply the finite difference method. So, we can write it as w i plus 1 minus 2 w i plus w i minus 1 divided by del x square and on the right hand side it will be w i divided by K. We are using a single K value for the whole domain, because we have considered that aquifer is homogeneous in nature. So, now, this particular equation can be used as a constraint for our optimization problem. So, our optimization problem then it will become maximize our z, which is summation of all small w i's, where i belongs to that set capital i which is a set of all wells. Now, this is subject to our total pumping for all pumping wells that should be greater than equal to minimum pumping value and other constraints that we need are small w should be greater than equals to 0, head cannot be negative and also our pumping cannot be negative. So, these are the constraints which need to be satisfied for our objective or optimization problem and we have our this particular equation that will act as binding constraint for the optimization problem or equality constraint. This approach is basically our embedding technique approach, where we directly utilize the governing equation within the optimization problem as constraints. So, with this h i or head at a particular location that can be determined from our this root over w small w i thing. And we can also satisfy additional constraints, because for this particular location we can have constraints like w 5 that should be less than equal to w 4 small w. Similarly, our w 4 should be less than w 3 and w 3 should be less than w 2. Similarly, we have w 2 less than equals to w 1 and w 1 is less than equals to our w naught. So, like our previous confined aquifer problem, we can set another problem where this is your impervious datum and we have 4 wells. So, this is the 4 well situation both the left hand and right hand boundaries are defined in terms of hydraulic head that is basically h 5 and h naught. So, these are basically specified boundary conditions or specified head boundary condition. And we have q 1, q 2, q 3 and q 4 for this situation. Then this corresponding to this we have h 1, this is h 2, this is h 3, this is h 4. So, we can see that this particular equation, this particular set of constraints can be directly utilized within the optimization model to get the solution. Because the head here at this downstream or down gradient point should be should have a less value compared to one up gradient point. So, in this particular problem, this is the aquifer portion, we have the problem as maximize the z where small w 1, w 2, w 3, w 4. This is our objective function subject to finite difference equations that is 2 w 1 plus w 2 minus del x square by k w 1 minus w naught and w 1 minus 2 w 2 plus w 3 minus 2 del x square by k w 2, this is 0. Again we have w 2 minus 2 w 3 plus w 4 minus 2 del x square divided by k and w 3, this is 0. And the final equation, we can write it as w 3 to w 4 minus 2 del x square by k w 4 equal to minus w 5. All these are small w's, these two are small w's in this side, these are basically small w's and these are basically capital W values or pumping values. These corresponds to our hydraulic head. Now, this is this set is the binding constraint, this set is the binding constraint and this particular set which we have already written that is valid for the optimization problem that will give a proper result because any down gradient hydraulic head cannot be greater than the up gradient or cannot be greater than compared to a up gradient value. So, what are the other constraints that will be required for this one? Those are related to production rate, production rate. In case of production rate we have w 1, w 2, w 3, w 4, this should be greater than w minimum and small w i's should be greater than 0 for i equals to 1 to 4, capital W i greater than equals to 0 for i 1 to 4. Again all the objectives and constraints are linear in nature. So, we can directly solve it using linear programming algorithm. Now, we have already covered this confined and unconfined one dimensional flow management problem for ground waters. Now, we will try to see what is there in two dimensional thing. So, for any steady state homogeneous confined aquifer system the equation can be written as this is this T is having single value because we are considering constant hydraulic conductivity over the two dimensional aquifer. So, if we discretize this. So, discretization should be i j minus 2 i j plus h i minus j divided by del x square plus h i j plus 1 to h i j plus h i j minus 1 divided by del y square equals to W i j and divided by T. So, if we consider that del x equal to our del y then we can rewrite this equation as h i plus 1 j minus 4 h i j plus h i minus 1 j plus h i j plus 1 plus h i j minus 1 equal to del x square W i j divided by T. Again if we want to write a two dimensional management problem 2D management problem for this kind of aquifers we can write it as maximize z. In this case we should have h i j where i comma j this belongs to the set i for welds and subject to our binding constraints that we have already derived for 2D aquifer and i j all i W i j greater than W minimum and h i j greater than 0 W i j greater than equals to 0. So, this can be solved using again by LP method or linear programming method. So, if you have a transient problem let us consider transient problem with confined aquifer system. So, transient problem with confined aquifer. So, in this case we can have this del this is for one dimensional case and this is transient 1D thing. So, the equation del 2 h by del x 2 that can be expressed as Crankt-Ekelson thing. In Crankt-Ekelson half of the derivatives are evaluated at the present time step and half of the derivatives are evaluated at the future time step. Let us n denotes the time level and i plus 1 denotes our space level. So, this is basically 2 h i n plus h i minus 1 n divided by del x square plus our h i plus 1 this is n plus 1 level to h i n plus 1 level plus h i minus 1 that is also at n plus 1 level. So, for a single the second order derivative we are using this Crankt-Ekelson scheme. So, in this scheme we can discretize the whole thing like this otherwise we can also discretize this thing using our any particular time step. So, Crankt-Ekelson scheme is basically bounded 1. So, we have this is nth n plus 1 level this is nth level and this is spatially i th derivative or i th location this is i minus 1 this is i plus 1 this is again i i minus 1 i plus 1. So, it involves all these points within our calculation and that way it is it will give advantage during the solution process. And if you discretize our this time dependent or transient term like this this is i n plus 1 minus h i n th level and del t also we can discretize our w that is average between i n th level and i n plus 1 level. So, if we substitute these terms in our original equation we can get one linear equation and again if we employ our transient problem our transient values within our optimization framework then it will be a time dependent optimization problem. So, we can write this thing as maximization of z and our i all i and this is i tau tau is basically our last time step. So, our subject to we can use the constraints as i w i n or we can directly use it here t level and this is greater than equal to w minimum and which is at t level and this t varies from 1 to that tau. So, at the last time step our head values should be maximum or we should maximize our head values at the last time step. And we have this constraint that for each time period or management period we can have the total pumping which will be greater than our minimum value of pumping. Finally, if we have any 2D unconfined aquifer situation then we have the equation like this and with our usual substitution for unconfined aquifers we can rewrite the equation as del x square del 2 w y 2 2 w k and we can discretize this particular equation with our finite difference formulation. Now, this is time independent or steady state problem for unconfined aquifer situation. Now, let us talk about the simulation optimization problem. In case of simulation optimization problem we directly use the simulation model as binding constraint within our optimization model, but we link it externally using any ready made simulation model or we can write our code, but we link it externally with the optimization model. For any flow and transport related problem we can have one extra constraint that is your concentration value should be below the specified concentration limit. So, in this case this has been taken from our salt water inclusion management problem. So, we have this maximization of pumping from production wells then minimization of pumping from extraction or barrier wells because these pumped waters cannot be directly used for our water supply systems. We need to have some kind of reverse osmosis or any other kind of plant to treat these water. So, this capital Q is for any x i location and T k is the time level. So, time period so, at any specific spatial location at any particular time period. So, summation over all spatial locations and all spatial all time periods we can maximize this production function or production well pumping and minimize the pumping from extraction barrier well. So, that is for location j and for time period k and this is subject to this externally linked simulation model C we can directly get from this external g function which is not exactly one proper expression, but it is some kind of linked thing and this is concentration related constraint that it should not be greater than any specified concentration limit these are the limit for our production well and barrier well. So, we can have this kind of situation where this is our optimization model. So, during evaluation of objective and constraints we can go to f m water f m water is density dependent flow and transport simulation model. So, during this objective function and constraint calculation we need the information about concentration. So, what we are passing through this? So, we are passing the information related to Q and small Q that is pumping values for production well, pumping values for barrier well and what we are getting out of this f m water model we are getting the concentration values that way it is externally linked. Now, with this algorithm we can have this final solutions that is this is this is the final solution set one solution set that is for a number of objective function values. So, one objective is maximization of f 1, f 1 is the production one is minimize the f 2. So, that way there is conflict. So, if we maximize one thing there will be increase in other variables value. So, there is tradeoff between first objective and second objective. So, let us consider one hypothetical example this particular phase let us say this is our ocean phase or sea phase these three wells are extraction or barrier wells and these wells are our production wells and these two boundaries are no flow boundaries and this is our inland phase. So, this is the size 1800 meters 100 meter thickness of the aquifer and 1400 meters and this has been modeled using f m water. So, these are the well screen levels. So, these are the values for the hydraulic conductivity 25 meters per day for x and y direction and 0.25 meters per day in z direction longitudinal dispersivity alpha l is 50 meters alpha t that is 20 meters molecular diffusion coefficient that is 0.69 meter square per day. Soil porosity that is 0.2 and density reference ratio interestingly previously we have used this alpha c divided by C s. Now, this alpha by C s is basically this particular epsilon. So, this is density reference ratio this is our vertical recharge although it is having that small q notation, but it is different from or pumping from extraction barrier thing. So, different parameters let us say we have that C max value is 500 mg per liter. So, this is our vertical recharge so this is standard for secondary maximum contaminant level. So, this is 500 mg for any multiple objective optimization because we have two objectives. So, we can designate it as multi objective optimization. So, this is called as Pareto front. So, one objective is spreading in these two directions spreading of final solutions and second and this is the second goal and the first goal is movement of solutions and the final front should be almost coinciding with our Pareto front or it is nearing to our Pareto front. So, with different number of iterations or generations let us say this is generation number one there is first objective is can be realized because this is 100 generation 200 300. So, it is moving towards 300 number in 300 number generation it is moving towards Pareto front and with the increase in generation number it is also spreading in both the directions. So, let us say this is our final front. So, in this particular front 11 and 14 are two points and these two points are showing two different results two different means 11 and 14 will give you two sets of capital Q and small Q values and with that you can explain the process. One thing is that for point number 11 we have F 1 value which is less than our F 1 value compared to 14, but it is having a better value for F 2 compared to solution number 14. So, this we can see that there is some kind of tradeoff between the solutions. So, it took around 24 into 800 simulations it took 30 days of running time. So, the problem is that running time is a problem for linked simulation optimization. So, what we can do we can use our metamodel based approach for management purpose. First thing is that response matrix approach, other thing response surface methodology or other metamodel based approaches artificial network radial basis function, support vector machine, relevance vector machine and creaking model and GP. So, what is the difference between our original simulation model and metamodel? So, original simulation model will give you the exact value of H or C, but in case of metamodel there will be some amount of error involved with it. So, what we can do is although we can gain in terms of solution time or gain in terms of solving a particular problem in management related issues. So, first thing is this response matrix thing, response matrix approach and this response matrix approach let us say draw down for any location s draw down s for any location k n, k is the special location or cell and n at the end of nth time period can be written in terms of this expression. So, in this case this unit response function, unit response function it tells something that is in the kth cell it tells something that kth cell at the end of nth time period due to unit. So, it is basically denoting the change in draw down for nth cell the kth cell at nth time period due to due to unit pumpage from jth cell and pth time period. So, this is some kind of unit response for kth cell at the end of nth time period due to the unit pumpage or injection in the jth cell during pth time period. So, for kth cell at the end of nth time period we need to consider all time periods which are less than nth time period. So, starting from 1 to n and j we need to consider all possible cell locations. So, this is all possible cell locations. So, we can replace our binding constraints with this particular expression these unit response functions can be obtained by simulation of the original simulation model. Now, we can use ANN, ANN has got architecture of something. So, we have some input and we can get some output and we can evaluate the performance using total error. Total error is desired and actual one and again desired and actual one average absolute relative error can be computed using this. So, this is basically number of, let us say layers and we can also find out this correlation coefficient for this one to check the accuracy of the ANN models. So, what is there in the link simulation optimization with meta model? We have physical simulation model here and we have meta model here and we have optimization model here. So, from simulation model we can generate input patterns using that this is Latin hyper cube sampling, this is one sampling strategy with which we can generate our input pumping values and we can input that in FEM water model and we can generate training and testing data set. We can train that ANN model or any SVM or GP model, then we have that train model and train model will pass that capital Q and small q again will get C value part of this and we can use this C value for optimization and again we can get this multi objective solutions. So, these many iterations are required for the testing part and the 3000 training data set was used and 600 testing data set and this was the architecture 33 inputs, these are two hidden layers and this is the output thing. So, we can see that meta model based approach is giving almost equal or better result compared to our direct link simulation optimization model. There can be another approach where we can use this screening model thing. So, we have this one. This physical simulation model from physical simulation model we can have our meta model from meta model. We can run our optimization model and we can get some intermediate Pareto front or intermediate solutions and we can pass that particular solution to the final objective or final optimization model where we have original optimization, original simulation model is linked with the optimization model. So, that way we can get a good accuracy in terms of results and also we can reduce the number of iterations that will be required for simulation of original FU water model during management optimization simulation process. So, partially trained meta model linked simulation optimization and with 100 generation with 1000 generation these are giving this is giving screening model based approach. This is with link simulation based approach direct link simulation based approach it is giving always the better result because these are the correct values and it is more near to our Pareto front because as it moves in the right hand direction it is more near to our Pareto front. There can be link simulation optimization with operation uncertainty we can run this simulation model for multiple realizations. We can give Q and small Q values and we can generate Q plus del Q combinations del Q is small variation because in reality there will be variation of these values in field situations. So, with this combinations we can run the model and we can get the average objective function value and standard deviation that we can utilize here and we can get some kind of C bar with some seed standard deviation C bar means some mean value for the concentration we can pass it and we can get a robust optimal solution for the management problem. So, integrated planning mechanism so, we have robust optimal solution from this multiple realization approach we can get final front or final Pareto front out of that we will select one strategy and we will implement that single strategy in the field there will be monitoring and we will collect some information again we will with those information we can update the simulation model that is the actual process for any integrated planning approach. Thank you.