The recharge and discharge of groundwater are mainly vertical surface infiltration recharge, evaporation discharge (evaporation can be shown as negative value of infiltration) and lateral surface water recharge and discharge. One of the most basic problems in hydrogeology is the evaluation criteria of sustainable exploitation of groundwater, which involves the increment of recharge and the decrease of discharge. Therefore, the prediction model of groundwater exploitation must include the above two kinds of recharge and discharge factors, otherwise the evaluation of groundwater cannot meet the requirements. However, the classic Theis unsteady well flow model (1935) only involves the recharge and discharge of the side boundary, and does not considers the surface infiltration recharge of the upper boundary, even when groundwater is pumped near the river. In this way, the analytical model cannot be basically used for prediction, and can only be used for well flow test to obtain the aquifer parameters in dry season. Therefore, the goal of this paper is to develop the Theis unsteady well flow model with surface infiltration recharge. For the problem of phreatic flow, the equation of groundwater flow cannot be established by using head as the dependent variable in a confined aquifer. We use the potential function of the second linearization method to establish groundwater flow equation in the phreatic aquifer. For the solution of the generalized mathematical model of the complex hydrogeological problem with rainfall infiltration and a pumping well, we adopt the method of decomposing it into several simple sub-models and synthesizing them into the solution of the original complex mathematical model. Based on the principle of mass conservation and assuming that the seepage obeys the Darcy’s law and satisfies the Dupuit's assumptions, the differential equation of groundwater flow is established. Then, for the well flow problem with uniform and stable infiltration recharge in two parallel rivers and two kinds of strip regions formed by a river parallel to an impermeable boundary, the general equation of groundwater flow and several kinds of water table equations under specific conditions and their flux equations are obtained. In addition, the “boundary to boundary reflection method” is proposed and applied to solve the same problem in a strip of region between a river boundary in parallel to an impermeable boundary, which reduces many derivation processes. Finally, as a preliminary application of the above theoretical results, it is also an important application, that is, there is a pumping well near the river whose water quality cannot meet the requirements, and the critical flux equation of the pumping well is calculated on the premise of not absorbing the river water. The important and concise relation equation is obtained. The equation can also be used to calculate the critical pumping rate of pumping wells in coastal areas without seawater intrusion. In this paper, the groundwater flow network diagram at a certain time in the process of unstable well flow under the above conditions is given. Compared with the flow network of well flow near the river, which is commonly seen in the literature, the flow network has obvious characteristics.