Optimal Operation Model of Drainage Works for Minimizing Waterlogging Loss in Paddy Fields

: The risk of ﬂood or waterlogging in irrigation districts has increased due to global climate change and intensive human activities. A Model of Optimal Operation of Drainage Works (MOODW) for ﬂat irrigation district was established by incorporating the hydrological model of waterlogging process and waterlogging loss estimation, which was solved by an optimization method of genetic algorithm. The model of waterlogging process was built based on a modiﬁed Tank model and hydrodynamic model for the ditch-river system. The waterlogging loss is calculated under the condition of inconstant inundated depth by linear interpolation. The adaptive genetic algorithm with the global optimization function was selected to solve the model. With an extreme rainfall events in Gaoyou irrigation district as cases, results showed that operation time and numbers of pumps increased; thus, operating costs were 1.4 times higher than before, but the yield loss of rice decreased by 35.4% observably. Finally, the total waterlogging loss was reduced by 33.8% compared with the traditional operation of waterlogging work. The most signiﬁcant improvement was found in units with high waterlogging vulnerability. The MOODW can provide the waterlogging information visually and assist the district manager in making a reasonable decision.


Introduction
Flood is one of the main natural disasters because of its high frequency, wide distribution, heavy disaster degree, and significant economic loss. The agricultural system is more sensitive to flood disasters than urban or other systems and is essential to disaster research [1][2][3]. More than one-third of the world's irrigated land suffers from waterlogging worldwide [4,5], reporting that waterlogging affected the crop in many ways, leading to the decline in crop yield. Waterlogging can cause different degrees of harm to different kinds of crops or plants in different duration periods. According to statistics of the Ministry of Water Resources, R. P. China, floods in China affected 47.666 million people and 6,684,000 hectares of crops, of which 1,3215,000 hectares failed to harvest in 2019. The direct economic loss was 192.27 billion yuan, accounting for 0.19% of that year's GDP [6].
Irrigated districts with flat terrain, large population, and fertile land are essential grain production bases in China. However, in recent years, with the impact of global climate change, the frequency of extreme rainfall events has been increasing [7]. At the same time, affected by human activities of urbanization, the irrigated districts are prone to flood disasters, which result in colossal damage to agricultural production [8][9][10].
Understanding the occurrence and evolution of waterlogging in the flat irrigated district is the basis for disaster control. Although the plain irrigated district is characterized by a complex drainage system and variable river direction, which makes the waterlogging in the irrigated district a complicated process, this situation becomes more complex due to the impact of human activity [11].
Simulating or describing the waterlogging process accurately is the basis of optimal operation of waterlogging mitigation works in irrigated districts. Researchers have developed hydrological models to simulate the characteristics of the runoff process over the different underlying surfaces and under different scenarios. Yet, most of the existing hydrological models were developed over natural watersheds, which cannot directly describe the water movement in irrigated districts under the intervention of human activities [12][13][14]. Traditional hydrological and hydraulic models such as SWAT, WEP, and Stanford are mainly developed for natural watersheds suitable for large-and medium-sized basins. These models always provide runoff processes at the outlet of the basin. Nonetheless, they are unable to give results about waterlogging, such as the depth of the hydrological response unit, especially in farmland.
The Tank model was proposed by Sugawara in 1972 [15]. This model has been widely used in various terrain and climatic conditions because of its simple principle and flexible structure. Studies have shown that the Tank model can accurately simulate the rainfall and runoff process in the paddy field [16,17]. An improved tank model is also applicable to flat, low-lying agricultural areas' drainage and inundation analysis [18]. Thus, the tank model is a potential tool in the drainage analysis of paddy fields.
Generally, there is a flood control system with multiple drainage works in a specific irrigated district. Optimizing the operation of multiple drainage works is essential for disaster control. The optimal operation aims to give full play to the function of works in the river system, improve waterlogging control, and reduce the loss caused by waterlogging. The application of intelligent optimization algorithms to the optimal dispatch of flood control systems has become a hot topic in recent years. Genetic algorithm (GA) is a typical intelligent optimization algorithm applied for optimization research, especially for complex, multistage decision-making problems [19]. Oliveira and Loucks firstly applied GA to a reservoir swarm optimization dispatching system and achieved good results [20]. Then, researchers widely used GA in real-time reservoir flood control dispatching [21,22] and multidatabase joint dispatching systems [23,24]. As a result, the optimal operation of flood control systems has become mature in theory and practice. However, the research work mainly focuses on reservoir flood control, and there is still little research on the optimal drainage system operation. Optimal scheduling of drainage works is a complex, multistage decision-making problem. GA will be helpful for solving this problem.
The objectives of this study were to develop an optimal operation model that can describe the waterlogging process, assess the waterlogging losses in paddy, and provide an optimal dispatching scheme of multiple drainage works to minimize the loss caused by waterlogging. It was tested with data collected in Gaoyou Irrigation Area. It will offer an effective tool for optimizing the operation of multiple drainage works for reducing the loss caused by waterlogging in flat irrigated districts.

Assumptions
In the Model of Optimal Operation of Drainage Works (MOODW), drainage sluices, ditches, and pumping stations were taken as the scheduling object. The following assumptions were made according to the practice of drainage works management.

1.
Ignore the duration it takes to open and close a specific drainage work, and all drainage works are scheduled once only for one waterlogging event.

2.
All drainage works only have two states-fully open and fully closed-the sluice is opened with opening height above the water surface or the pumping station is in operation at the rated flow, respectively. For the pump station with gate, the pump was launched only when the gate was closed. A pump station with multiple pumps is considered as one work. 3.
The flow rate and power consumption of a specific pump station are assumed to be constant and are calculated according to the rated value.

Objective Functions and Calculation Formula
As shown in Figure 1, according to the requirements of the operation of drainage works, the starting time (t 0 ) and closing time (t 1 ) of the drainage work were selected as decision variables in the model. When the number of drainage works is greater than 1, the decision variable is a multidimensional vector, denoted as t 0 = t 1 0 , t 2 0 , . . . , t N 0 , t 1 = t 1 1 , t 2 1 , . . . , t N 1 . N is the total number of drainage works. was launched only when the gate was closed. A pump station with multiple pumps is considered as one work. 3. The flow rate and power consumption of a specific pump station are assumed to be constant and are calculated according to the rated value.

Objective Functions and Calculation Formula
As shown in Figure 1, according to the requirements of the operation of drainage works, the starting time ( ) and closing time ( ) of the drainage work were selected as decision variables in the model. When the number of drainage works is greater than 1, the decision variable is a multidimensional vector, denoted as is the total number of drainage works. The optimization objectives were the minimum loss of waterlogging and cost of drainage. For the cost of drainage, only the energy cost of pumping station for pumping the excessive water out was considered.
In Equation (1), is the minimum summation of work operation cost and flood loss (ten thousand Yuan); is the number of hydrological units; ( ) is yield loss rate of rice as a function of the flooding depth of , %; is the price of rice, yuan/kg; is the average yield of rice in normal years, kg/m 2 ; is the area of the paddy field in unit , hm 2 ; is the number of drainage pumping stations; ( ) is rated power of the drainage pumping station , kWh; ( ) is working efficiency coefficient of drainage pumping station , %; is price of electricity, yuan/kW • h; ( ) is the running time of the pumping station , h.

Waterlogging Loss Estimation
Model for waterlogging loss estimation was calculated as a function of inundated water depth varied along with waterlogging process in every grid of flood. According to the planting structure of crops in the study area, the total loss of waterlogging of crops in the disaster area is calculated by multiplying the flooding loss rate with the total economic output of each crop in average years. The mathematical expression of the waterlogging loss in the unit is as follows: The optimization objectives were the minimum loss of waterlogging and cost of drainage. For the cost of drainage, only the energy cost of pumping station for pumping the excessive water out was considered.
In Equation (1), Z is the minimum summation of work operation cost and flood loss (ten thousand Yuan); N is the number of hydrological units; Y(H i ) is yield loss rate of rice as a function of the flooding depth of H i , %; β is the price of rice, yuan/kg; V is the average yield of rice in normal years, kg/m 2 ; A i is the area of the paddy field in unit i, hm 2 ; K is the number of drainage pumping stations; P(j) is rated power of the drainage pumping station j, kWh; η(j) is working efficiency coefficient of drainage pumping station j, %; δ is price of electricity, yuan/kW·h; T(j) is the running time of the pumping station j, h.

Waterlogging Loss Estimation
Model for waterlogging loss estimation was calculated as a function of inundated water depth varied along with waterlogging process in every grid of flood. According to the planting structure of crops in the study area, the total loss of waterlogging of crops in the disaster area is calculated by multiplying the flooding loss rate with the total economic output of each crop in average years. The mathematical expression of the waterlogging loss in the unit i is as follows: where AD is the total loss of farmland in irrigated district, ten thousand Yuan; n is number of crop species in the study area; D m is economic losses of crop k per unit area in this flood, yuan/m 2 ; CRP a -the planting area of crop k in cell i, hm 2 ; mm is loss coefficient of crop k at different growth stages; CP k is unit price of crop k, yuan/kg; F k is yield of crop k per unit area in average years, kg/m 2 ; DC k is yield reduction rate of crop k per unit area (as a percentage of average annual yield), %. Flood disasters in the flat irrigated district mainly occurred in summer in the Yangtze River Basin, synchronous with rice growing period. Therefore, rice was taken as a typical crop to analyze the vulnerability characteristics of crops and build a vulnerability model.
In this paper, yield loss was selected as the evaluation index to study the inundation loss of rice. The yield reduction was mainly affected by the rice growth stage, duration of inundation, and depth of inundation [25]. Generally, the crop waterlogging production function (WPF) is a mathematical description of the flooded depth H, the flooded duration T, and the rice yield reduction rate DC rice . In addition, researchers have fitted some functional relations through statistical analysis, among which the widely used relations are in exponential form, and the fundamental equation of the mathematical model is as follows: where H is the percentage of flooded water depth to plant height, %; T is flooded duration, day; a, b, and c are parameters of the model. The rice yield reduction rate model was calibrated using the data collected from Jiangsu Province. The model was verified using the experimental field data, except for the booting stage, by [26]. The calibration results are shown in Table 1.

Waterlogging Process Simulation
The tank model simulated the rainfall-runoff process in the paddy field. According to the hydrological and geological characteristics of the irrigated district, the rice field between two channels was regarded as a hydrological unit, and a two-layer tank model developed by [27] was run on each unit. In the rainfall drainage process, the excess water from the paddy field was drained into the river channel through the outlet of the drainage ditch. Figure 2 showed the structure of the two-layer tank model.
In this model, the rice field drainage was uniformly drained into the canal or river along with the water flow and used as the sidestream of river channel. In the two-layer tank model, two side holes and a bottom hole were set in the first layer, the lower hole was used to simulate the side seepage of the ridge of the paddy field, and the upper hole was used to simulate the overflow runoff from the paddy field. The second layer of the Tank model has two side holes and a bottom hole for simulating soil flow. According to the study of [28], the broad-crested weir overflow formula can be applied to calculate the discharge from the paddy field. Therefore, the submersed discharge formula and the free discharge formula of the broad-crested weir overflow formula were used to calculate the discharge from the side holes of the first layer of the Tank model.
The 1D Saint-Venant equations [29] were used as the basic confluence equations of plain river network. This paper uses Preissmann's implicit scheme to solve the Saint-Venant equations, including dynamics and continuity equations. In this model, the rice field drainage was uniformly drained into the canal or river along with the water flow and used as the sidestream of river channel. In the two-layer tank model, two side holes and a bottom hole were set in the first layer, the lower hole was used to simulate the side seepage of the ridge of the paddy field, and the upper hole was used to simulate the overflow runoff from the paddy field. The second layer of the Tank model has two side holes and a bottom hole for simulating soil flow. According to the study of [28], the broad-crested weir overflow formula can be applied to calculate the discharge from the paddy field. Therefore, the submersed discharge formula and the free discharge formula of the broad-crested weir overflow formula were used to calculate the discharge from the side holes of the first layer of the Tank model.
The 1D Saint-Venant equations [29] were used as the basic confluence equations of plain river network. This paper uses Preissmann's implicit scheme to solve the Saint-Venant equations, including dynamics and continuity equations.

Constraint Conditions
Constraints for optimal operation of drainage system in plain river network district should be considered from the hydraulic calculation of river network, the boundary flow capacity of sluice station, and the mode of dispatch of sluice station. The primary constraints are listed as follows: Water level constraints: For nodes without the capacity of water storage and regulation, the flow connection and water level connection must be satisfied in the process of hydraulic calculation of the river network as shown below: Figure 2. Structure of the two-layer tank model. R represents the outflow of side hole, p represents the outflow of bottom hole, S represents water depth of the tank, and E means evaporation.

Constraint Conditions
Constraints for optimal operation of drainage system in plain river network district should be considered from the hydraulic calculation of river network, the boundary flow capacity of sluice station, and the mode of dispatch of sluice station. The primary constraints are listed as follows: Water level constraints: For nodes without the capacity of water storage and regulation, the flow connection and water level connection must be satisfied in the process of hydraulic calculation of the river network as shown below: In Equations (5) and (6), Q j i is the flow of river j into node i, m 3 /s; A i is water surface area of i, m 2 -for nodes without storage capacity, its value can be considered as 0; Z j i is the water level of node i in river j, m; Z i is average water level of node i, m.
Hydraulic constraint: The common connecting structures are controlling gates and culverts. Therefore, the hydraulic connection between channels can be expressed by the overflow formula: In Equation (7), Q is overflow discharge of structure, m 3 /s; Z u is upstream water level of overflow structure, m; Z d is downstream water level of overflow structure, m; Z 0 is floor elevation of overflow structure, m; B is overcurrent width of overflow structure, m.
Time constraint: Start time of specific drainage works should be greater than or equal to 0. The time at the end minus the time at the beginning of the drainage work is equal to the operation time of the drainage work. In the drainage period, the process of the drainage work is uninterrupted.
In Equations (8) and (9), Constraint of river network level: Considering all kinds of water demand such as water supply, irrigation, and ecological water use, the flood discharge capacity of the river must be fully exerted on the premise of maintaining a certain water level at the end of the flood season. According to the dispatching requirements, the water level should recover to near the preflood level at the end of the flood season. In addition, the water level of each calculated section must be higher than the ecological water level of the section at any stage to prevent ecological damage.
In Equation (10), Z i is water level of channel section i, m; Z e i is ecological water level of channel section i, m.
The Water level constraint, Hydraulic constraint, and Time constraint are automatically satisfied in the design of a one-dimensional unsteady flow model of the river network. As for the river network level at the end of flood season, it can be satisfied when the water level falls below the maximum waterlogging level. Therefore, the only constraints that need to be solved in the optimization model are drainage continuity and river level constraints.
The operation rationality of the pumping station should be considered in the drainage continuity constraint: where GCloseT, GOpenT are closing time and opening time of gate, respectively; BOpenT, BCloseT are closing time and opening time of pump station.

Model Solving by Using Genetic Algorithm Genetic Algorithm
The optimal operation of drainage work is a complex, multistage decision-making problem. The model in this paper contains many decision variables and complex constraints. The influence of decision variables on the objective function needs to be calculated through tricky runoff generation and affluxion, as well as waterlogging loss. Therefore, genetic algorithm [30] with the global optimization function was selected to solve the model. The calculation process of the optimal scheduling mode of waterlogging removal work deduced by the genetic algorithm is shown in Figure 3.
The decision variable of the optimal scheduling problem of drainage works is the start and end time of each works, which is a multidimensional decision variable. Therefore, this paper selected the floating-point encoding method to encode decision variables, and its representation is t 0 : t 0 1 , t 0 2 , . . . , t 0 n and t 1 : t 1 1 , t 1 2 , . . . , t 1 n . t 0 is the start operation time of drainage works and t 1 is the end time of drainage work operation. n is the total number of drainage works.
Water 2021, 13, x FOR PEER REVIEW 7 of 18 model. The calculation process of the optimal scheduling mode of waterlogging removal work deduced by the genetic algorithm is shown in Figure 3. The decision variable of the optimal scheduling problem of drainage works is the start and end time of each works, which is a multidimensional decision variable. Therefore, this paper selected the floating-point encoding method to encode decision variables, and its representation is : { , ,⋅⋅⋅, } and : { , ,⋅⋅⋅, }.
is the start operation time of drainage works and is the end time of drainage work operation. is the total number of drainage works.
To prevent interference by the artificial selection, the initial population of decision variables randomly generates according to the working conditions of drainage work, which are as follows: In Equation (12), , , , are, respectively, the upper limit and lower limit of the start time of drainage work , h; , , , are, respectively, the upper limit and lower limit of the end time of drainage work , h; both and are random numbers greater than 0 and less than 1.
The fitness measures the adaptability of a species to the living environment. The objective function of the optimal operation model is to minimize the cost of drainage and the To prevent interference by the artificial selection, the initial population of decision variables randomly generates according to the working conditions of drainage work, which are as follows: In Equation (12), t 0 i,min , t 0 i,max are, respectively, the upper limit and lower limit of the start time of drainage work i, h; t 1 i,min , t 1 i,max are, respectively, the upper limit and lower limit of the end time of drainage work i, h; both α and β are random numbers greater than 0 and less than 1.
The fitness measures the adaptability of a species to the living environment. The objective function of the optimal operation model is to minimize the cost of drainage and the yield loss of waterlogging, which does not meet the requirements of maximizing the fitness evaluation function. Therefore, the objective function is converted into a fitness function as follows: In Equation (13), F is fitness function; Z is objective function; Z max is the maximum possible loss of waterlogging. Since the operating cost of drainage works takes a tiny proportion to the total loss, the total output value of the paddy fields in average years in the study area is taken as the maximum possible loss of waterlogging.
This model used the roulette selection method recommended by Goldberg [31] to select individuals with higher fitness to be inherited to the next generation. The method determined the selection probability or survival probability of the individual according to the proportion of fitness value of each chromosome. This selection process simulates the rotations of a roulette wheel and the frequency is equal to the population size; for each trial, one individual was selected for the new population. The characteristic of the roulette selection method was the random adoption process. First, the arithmetic crossover is used to generate individuals. That is, a linear combination of two-parent bodies produces two new individuals. The next step adopts uniform variation.
Finally, genetic iteration algebra and the number of repeats of the optimal individual were used as termination constraints. In this paper, the maximum number of iterations is 500, while the optimal number of individual repetitions is 50.
In this study, the penalty function is adopted to deal with the river level constraint. The optimization process must meet the requirements of the river's ecological water level. Therefore, when the water level of any section of the river is lower than the ecological water level, a minimal fitness value will be returned, which eliminates the operation scheme.

Study Area
Gaoyou Irrigation District is a typical plain irrigation area in the east of China, located in Yangzhou, Jiangsu Province. The irrigated district covers 650 km 2 , composed of independent polders such that there is no water exchange between them. Therefore, this study only selects a specific agricultural field, Longben Polder, as the object to study the operation mode of drainage work.
Located in the middle of Gaoyou Irrigation Area, Longben polder is a typical agricultural polder area. The total study area is 26.6 km 2 , of which the cultivated area accounts for 90%. It is high in the south and low in the west, and the flow direction of the whole area is from south to north. The drainage system is composed of embankment, drainage ditch, barrier pond, river course, and gate station. There are eight pumping stations and four culverts in the research area, which are combined into seven stations according to the location of the pumping stations, including one inner drainage station and six outer drainage stations. The map of the study area and drainage system is shown in Figure 4.
The generalization of the drainage system is the basis of the hydrodynamic calculation of the river network in the flat irrigation district, and its rationality directly affects the accuracy of the complete analysis. According to the specific characteristics of the drainage system, the backbone rivers are regarded as the skeleton of the whole drainage system for waterlogging.
In irrigation districts, under artificial reconstruction, the engineering layout of alternate irrigation and drainage was basically formed. Among the four boundaries of the generalized field, three of them are channels or villages, and the other side is adjacent to the river or ditch. Therefore, there is no direct water exchange between fields, and field drainage flows directly into rivers or ditches. Each field corresponds to a river, and a unique number connects the rice fields and river channels.
As shown in Figure 5, the study area was divided into 23 fields. There are a total of 32 nodes and 31 riverways (including hydraulic structures) in the study area. Riverways relate to each other by nodes. Due to the small scope of the study area, each channel is divided into five equally spaced sections in the hydrodynamic calculation of the river network. As a result, the flow direction of the river in the plain river network area is uncertain. Therefore, for the sake of convenience, the assumption of the flow direction of the river is made: For the outer river channel, the flow direction is positive from the external node to the internal node-that is, when the calculated flow is positive, the flow direction of the river channel is from the external node to the internal node; otherwise, it is from the internal node to the external node. For the internal channel, the flow direction is assumed to be positive from east to west and from south to north. In this paper, all pumping stations and culverts are generalized as river channels with a length of 0 and are numbered uniformly.
. Kouhe, (2). Nijia, (3). Hongqi, (4). Jiangmahe, (5). Zhongshihe, (6). Lvyanghe. river is made: For the outer river channel, the flow direction is positive from the external node to the internal node-that is, when the calculated flow is positive, the flow direction of the river channel is from the external node to the internal node; otherwise, it is from the internal node to the external node. For the internal channel, the flow direction is assumed to be positive from east to west and from south to north. In this paper, all pumping stations and culverts are generalized as river channels with a length of 0 and are numbered uniformly.

Calibration and Vertification
Detailed information of the data measurement and model calibration can be found in the paper by Xiong [26] and in Table 2, where the calibration results are listed.

Calibration and Vertification
Detailed information of the data measurement and model calibration can be found in the paper by Xiong [26] and in Table 2, where the calibration results are listed.  The application and effectiveness evaluation of the waterlogging process simulation model can be found in our previous study [32]. Here, we chose ten rainfalls to verify the feasibility of the model with three hydrological variables: paddy field runoff, water level, and flooded water depth in the paddy field, and the result is shown in Table 2. The result confirmed that the overall performance of the waterlogging process model is excellent.

Model Application
The study area has the characteristics of large rainfall and concentration. Through the analysis of rainfall occurrence probability, the 24-h rainfall corresponding to 1%, 2%, and 5% rainfall return period are 250.46 mm, 228.41 mm, and 196.1 mm, respectively (Table 3). The rainfalls were recorded by an automatic weather station (WatchDog2000 series, the U.S.A.) installed in the study area from 2012 to 2015. We selected two typical waterlogging processes for engineering optimization to verify the effectiveness of MOODW.
Rainfall Event 1 occurred from 13 August to 18 August 2014 (Figure 6a). The cumulative rainfall was 133.8 mm, the maximum 24-h rainfall was 108.8 mm, and the peak intensity was 18.4 mm/h. As a result, the lowest water level of the Nanchengzi River was 1.77 m, and the highest water level was 2.54 m. Furthermore, all drainage works were not scheduled during the rainfall, i.e., all of the pumping stations were closed during the deluge.  The rainfalls were recorded by an automatic weather station (WatchDog2000 series, the U.S.A.) installed in the study area from 2012 to 2015. We selected two typical waterlogging processes for engineering optimization to verify the effectiveness of MOODW.
Rainfall Event 1 occurred from 13 August to 18 August 2014 (Figure 6a). The cumulative rainfall was 133.8 mm, the maximum 24-h rainfall was 108.8 mm, and the peak intensity was 18.4 mm/h. As a result, the lowest water level of the Nanchengzi River was 1.77 m, and the highest water level was 2.54 m. Furthermore, all drainage works were not scheduled during the rainfall, i.e., all of the pumping stations were closed during the deluge. Rainfall Event 2 occurred from 9 to 21 August 2015 (Figure 6b). The rainstorm was of a typhoon-type, with the cumulative rainfall reaching 271.6 mm, the cumulative maximum rainfall in 24 h reaching 269.2 mm, and the maximum rainfall intensity reaching 47.4 mm/h. Such rainfall intensity is enormous, equivalent to the standard of once in 100 years in Gaoyou City. The river water level was 1.53 m at the initial rainfall stage, and the peak water level was 2.95 m.
According to the survey, the paddy field water depth did not exceed 20% of the plant in Rainfall Event 1, so there was no waterlogging loss and pump station operation. After Rainfall Event 2 occurred from 9 to 21 August 2015 (Figure 6b). The rainstorm was of a typhoon-type, with the cumulative rainfall reaching 271.6 mm, the cumulative maximum rainfall in 24 h reaching 269.2 mm, and the maximum rainfall intensity reaching 47.4 mm/h. Such rainfall intensity is enormous, equivalent to the standard of once in 100 years in Gaoyou City. The river water level was 1.53 m at the initial rainfall stage, and the peak water level was 2.95 m.
According to the survey, the paddy field water depth did not exceed 20% of the plant in Rainfall Event 1, so there was no waterlogging loss and pump station operation.
After running the optimization model, the results of the optimized operation scheme for pumping stations in Rainfall Event 2 were determined (as listed in Figure 7b). Compared with Figure 7a, the optimized operation scheme changed the time interval and duration of the pumps and gates. After optimization, the pump of P1 was not enabled, but the operation time of gate was increased by 32 h. The operation time of each gate from P2 to P7 was reduced by about 46 h, but the scheduling changes of the pumps were different after optimization. The pumping stations P2, P3, and P7 were added to the drainage scheduling, and the operation time was increased by 60 h, 60 h, and 65 h, respectively. The state of P4 pump became continuous with little change in operation time, while the pumps of P5 and P6, on the contrary, reduced the operation time by 32% and 6%, respectively. Overall, after optimization, the total working hours of the gate decreased by 244 h, accounting for 28.4% of the original. At the same time, the number of enabled pumps increased, the working hours became longer, rising by 49.5%. Therefore, after optimization, the operating electricity cost increased from 55,500 yuan to 76,800 yuan. running the optimization model, the results of the optimized operation scheme for pumping stations in Rainfall Event 2 were determined (as listed in Figure 7b). Compared with Figure 7a, the optimized operation scheme changed the time interval and duration of the pumps and gates. After optimization, the pump of P1 was not enabled, but the operation time of gate was increased by 32 h. The operation time of each gate from P2 to P7 was reduced by about 46 h, but the scheduling changes of the pumps were different after optimization. The pumping stations P2, P3, and P7 were added to the drainage scheduling, and the operation time was increased by 60 h, 60 h, and 65 h, respectively. The state of P4 pump became continuous with little change in operation time, while the pumps of P5 and P6, on the contrary, reduced the operation time by 32% and 6%, respectively. Overall, after optimization, the total working hours of the gate decreased by 244 h, accounting for 28.4% of the original. At the same time, the number of enabled pumps increased, the working hours became longer, rising by 49.5%. Therefore, after optimization, the operating electricity cost increased from 55,500 yuan to 76,800 yuan. As shown in Figure 8a,b, inundated water depth in different hydrological units mostly reached its maximum at the third day of rainfall events. Afterwards, it gradually decreased until there was a slight rise on the eighth day caused by another rainfall. The area with deep water in the paddy field is located east of the study area. For fields with higher topography in the southwest, the flood water was shallow and decreased rapidly with drainage. In comparison, the paddy water was deep and subsided slowly for the lower topography in the northeast. As shown in Figure 8a,b, inundated water depth in different hydrological units mostly reached its maximum at the third day of rainfall events. Afterwards, it gradually decreased until there was a slight rise on the eighth day caused by another rainfall. The area with deep water in the paddy field is located east of the study area. For fields with higher topography in the southwest, the flood water was shallow and decreased rapidly with drainage. In comparison, the paddy water was deep and subsided slowly for the lower topography in the northeast.
The optimized scheme resulted in a reduction of inundated water depth compared with local practice. In Figure 8c, the reduction in inundated depth varied greatly among units. The maximum reduction for most units occurred on the third day since raining. About 30% of fields had no or very slight decrease in inundated depth, primarily located in the southwest and northeast corner of the study area, which is relatively flat and not vulnerable to waterlogging. Waterlogging was eliminated faster for the optimized scheme, which was evident during the 4th to 7th day. As shown in Figure 8b, units with a considerable reduction in flooding water depth were mainly located in the low-lying areas with severe waterlogging. It implies that the optimization procedure tends to improve the practice in drainage work management and plays a more critical role in the part with high waterlogging crisis or vulnerability.  The optimized scheme resulted in a reduction of inundated water depth compared with local practice. In Figure 8c, the reduction in inundated depth varied greatly among units. The maximum reduction for most units occurred on the third day since raining. About 30% of fields had no or very slight decrease in inundated depth, primarily located It shows that the waterlogging vulnerability of the paddy field is greatly affected by both the discharge condition of the paddy field itself and the river level. For the area with higher terrain, the jacking effect caused by river water level is negligible. Hence, the variation of flooded water depth was caused by the combined effect of rainfall intensity and the drainage capacity of ditches. While in low-lying areas, the water depth of the paddy field is prone to being affected by river level. It will form a jacking when the river level is too high, which will affect the drainage speed and even cause backward water flow. The distribution of losses caused by waterlogging for the original and optimized scheme is shown in Figure 9. In the original operation, the heavy losses, 229,000~302,810 yuan/km 2 , occurred in a small part located at the eastern part of the study area. In Figure 8b, the reduction of average loss per km 2 was 13,670~166,530 yuan, which showed different decreases compared with the original operation in other regions. We found the most significant improvement in units with high vulnerability. By adopting the optimized scheme, the total yield loss of rice decreased by 35.4% to 1.589 million yuan, compared with 2.459 million yuan in local practice. The effect was very significant. The total of yield losses and operating costs of pump stations was 2,196,300 yuan, which decreased by about 33.8% compared with the original loss. scheme, which was evident during the 4th to 7th day. As shown in Figure 8b, units with a considerable reduction in flooding water depth were mainly located in the low-lying areas with severe waterlogging. It implies that the optimization procedure tends to improve the practice in drainage work management and plays a more critical role in the part with high waterlogging crisis or vulnerability.
It shows that the waterlogging vulnerability of the paddy field is greatly affected by both the discharge condition of the paddy field itself and the river level. For the area with higher terrain, the jacking effect caused by river water level is negligible. Hence, the variation of flooded water depth was caused by the combined effect of rainfall intensity and the drainage capacity of ditches. While in low-lying areas, the water depth of the paddy field is prone to being affected by river level. It will form a jacking when the river level is too high, which will affect the drainage speed and even cause backward water flow. The distribution of losses caused by waterlogging for the original and optimized scheme is shown in Figure 9. In the original operation, the heavy losses, 229,000~302,810 yuan/km 2 , occurred in a small part located at the eastern part of the study area. In Figure 8b, the reduction of average loss per km 2 was 13,670~166,530 yuan, which showed different decreases compared with the original operation in other regions. We found the most significant improvement in units with high vulnerability. By adopting the optimized scheme, the total yield loss of rice decreased by 35.4% to 1.589 million yuan, compared with 2.459 million yuan in local practice. The effect was very significant. The total of yield losses and operating costs of pump stations was 2,196,300 yuan, which decreased by about 33.8% compared with the original loss.

Conclusions
An optimal operation model of drainage work in the flat irrigated district was established by incorporating the tank model and yield loss estimation model. MOODW aims to minimize the sum of loss caused by waterlogging and energy cost of drainage works by scheduling the worktime of each pumping station. Several conclusions are summarized as follows: (1) Compared with the traditional optimal scheduling model for drainage works, the model takes into account the runoff process and the hydraulic connection constraints among the drainage works to achieve optimal scheduling in time and space.
(2) With extreme rainfall events in the Gaoyou irrigation district as cases, an optimized scheme was compared with the local practice and the result showed that the total waterlogging loss under the optimization was decreased by 33.8%. The most significant improvement occurred in hydrological units located in low-lying areas that are more vulnerable to waterlogging. Thus, the MOODW is an effective nonengineering measure for waterlogging control and disaster reduction, which will be an essential decision support tool for the management department of irrigated districts.
Due to the research limitations, the optimal operation model in this study is still in the exploration stage. There are some shortcomings, which are manifested in the following three aspects. First, the application events are insufficient. The MOODW is designed to reduce waterlogging loss from short-term heavy rainfall that usually exceeds drainage standards. However, during the two-year observation, we met only one super-standard rainfall in 2015. In the follow-up study, we believe that the model can further provide reasonable scheduling strategies by applying it to different simulated scenarios, such as rainfall at different frequencies and various water storage depths caused by irrigation strategies.
Secondly, the MOODW based on multiple waterlogging losses assessment needs further study. In this study, only rice inundation losses are considered in the objective function of the optimal operation model. Still, in the actual waterlogging process, multiple losses are caused, such as factory shutdown and village inundation. Therefore, the evaluation of waterlogging losses of different waterlogging-bearing objects and their influence on the optimal scheduling of waterlogging removal projects need to be studied.
Finally, it is not enough for the model to stay in the simulation of existing waterlogging. If the model can predict future waterlogging, it can send early warnings and make early dispatch responses [33,34]. So, it is necessary to strengthen the simulation study of possible future waterlogging and develop the model's early warning and forecasting functions.