Modeling of flow in porousmedia is a classical problem in both civil and petrolium engineering, and itplays as a critical role in analysis and design of civil engineering structures1, 2. Because of its mathematical challenges, it is also interesting inmathematical points of view. Thus, a considerable amount of literature has beenpublished on flow in porous media considering the notion of the phenomenon,numerical modeling and analytical solution 3,4,5,6,7,8. In traditionalstudies, Darcy equation (Darcy’s law) is used which represents a linearrelation between gradient of pressure and the velocity vectors 9,10,11,12.This assumption is valid as long as the Reynolds number is very low, while it is not areasonable assumption in many other engineering boundary value problems, wherethe grain sizes of porous media are not fine and the flow velocity isconsiderable.

A number of cross sectional studies have been made to modifyDarcy’s law 13, 14. The idea of including velocity by the power of two (ormore) in Darcy equation was introduced by Joceph et al. 15 known as Forchimerequation converting the linear Darcy’s equation to a nonlinear one. Brinkmansuggested another method to modify Darcy equation being more relayable in porousmedia with large porosity 16. The logic of the mentioned modification was theforce due to viscosity of fluid (water) to be taken in account in largeporosities.

This force has been considered accurately by the Brinkmanadditional term. Several attempts have been made to use some combination ofForchimer and Brinkman terms to form the momentum equation 17, 18, 19, 20, 21,22. Recently investigators have used Navier-Stockes type equations in modelingflow in porous media 23, 24. One of the challenges of modelling of flow can be seen in theproblems where the area of concerned consists of both porous and non-porousmedia. 23, 24, 25, 26.

These two parts are traditionally solved in twodifferent procedures, and the interaction of two domain are imposed by usingthe result of the model in previous time step on each part for setting theboundary condition of the other part and vice versa. Recently, a fewresearcheres have tried to solve the coupled equations for surface water andsubsurface flow to obtain more accurate results 23, 25, 29