Simulating the hydrologic cycle in coal mining subsidence areas with a distributed hydrologic model

Large-scale ground subsidence caused by coal mining and subsequent water-filling leads to serious environmental problems and economic losses, especially in plains with a high phreatic water level. Clarifying the hydrologic cycle in subsidence areas has important practical value for environmental remediation, and provides a scientific basis for water resource development and utilisation of the subsidence areas. Here we present a simulation approach to describe interactions between subsidence area water (SW) and several hydrologic factors from the River-Subsidence-Groundwater Model (RSGM), which is developed based on the distributed hydrologic model. Analysis of water balance shows that the recharge of SW from groundwater only accounts for a small fraction of the total water source, due to weak groundwater flow in the plain. The interaction between SW and groundwater has an obvious annual cycle. The SW basically performs as a net source of groundwater in the wet season, and a net sink for groundwater in the dry season. The results show there is an average 905.34 million m3 per year of water available through the Huainan coal mining subsidence areas (HCMSs). If these subsidence areas can be integrated into water resource planning, the increasingly precarious water supply infrastructure will be strengthened.

. Diagram showing formation of ground surface subsidence caused by coal mining and surface subsidence isoclines to the face 4. Numbers 1, 2, 3 and 4 represent the locations of working faces. W1, W2, W3 and W4 are the subsidence basins corresponding to the working faces. H0 is the depth to the coal seams. Yinghe, Cihuaixin and Huaihe Rivers can be used as head boundary as they are perennial rivers. In sum, it is reasonable and feasible to take the area surrounded by Yinghe, Cihuaixin and Huaihe Rivers as the study area (Fig. S3).    sediments are 30-400m, 40-100m and 25-45m respectively. There is a less permeable layer composed of loam and clay on the top of Q3 formation, and its thickness is about 10-35m. Other aquifers consist of fine sand, coarse sand, and silty sand, and contain several discontinuous weak permeable layers.
The landscape in this region is flat with a topographic slope of 1/15000-1/8000. The groundwater flow direction controlled by the topography is generally along the surface water flow. The small hydraulic gradient (approximately 1/10000) indicates that groundwater flow is slow and the lateral recharge is possibly limited.
Supplementary Information 2: Process scheme and conceptual model of MODCYCLE and its module "river-subsidence-groundwater" simulation module (RSGM)

The Modularized Model for Basin Scale Hydrologic Cycle (MODCYCLE)
The general model structure and hydrologic cycle simulation principle of the MODCYCLE code have been introduced in numerous Chinese monographs, journal papers and patents 2-4 , but not yet reported in international journals. Here we present the framework and part principle of the model, and hope to help readers to understand the role of the model in the hydrologic cycle simulation of the subsidence.

Model structure and hydrologic cycle path in the MODCYCLE
The MODCYCLE is a distributed hydrologic simulation model based on physical mechanisms. In the planar structure, the model firstly divides the river basin into several subbasins using the digital elevation model (DEM). Relationship among subbasins is built by the topology of main channels, which are generated following the corresponding subbasins. Then each subbasin is further divided into a plurality of basic simulation units (BSUs) according to the space distribution of land use & land cover (LUC), soil and land management operation (MGT). The BSU in fact is an aggregate that has the same LUC, soil and MGT. The BSUs may be scattered in the subbasin, and are not connected each other. In addition to the BSUs, the subbasins can include open waters such as swamps, wetlands, ponds, and lakes. The aquifers below the vadose zone are divided into two layers, the shallow layer and the deep layer, in each subbasin. To simplify the vertical flow simulation, the model only simulates the interaction between the shallow and deep groundwater system. The groundwater in subbasin aquifer can be independent or interconnected depending on the specific hydraulic conditions. Figures

Simulating principle of surface watergroundwater interaction
In the current version of the MODCYCLE, subbasins can be divided into mountain area subbasins and plain area subbasins to simulate the groundwater respectively.
where xx K , yy K , and zz K are the components of hydraulic conductivity along coordinate X, Y, and Z, which are assumed to be parallel to the major axes of hydraulic conductivity (LT -1 ) , h is the potentiometric head (L), W is a volumetric flux per unit volume representing sources and/or sinks of water (T -1 ), Ss is the specific storage of the porous material (L -1 ) , and t is time (T).

Spatial superposition between subbasins and groundwater numerical simulation grid cells
The superposition of subbasins and grid cells is diagramed in Fig  By using the above method, the spatial correlation between subbasins for hydrologic model and grid cells for groundwater flow simulation is established, which provides the basis for the two-way information exchange between the two models.  In a daily time step, the above "synchronous coupling" includes the following sub steps:

Information exchange between hydrologic model and groundwater flow model
Step 1: The hydrologic model implements hydrologic cycle simulation in the time step firstly, and calculate every vertical flux of groundwater in each subbasin: Step 2: The vertical fluxes of each subbasin simulated by the hydrologic model are transferred to the groundwater numerical model, which is based on the spatial relation between the subbasins and the grid cells and the area proportion of grid cell to subbasin.
Water sources and sinks of each grid cell can be calculated by: where iCell W is the total net water source (+) or sink (-) of the grid cell with number iCell(L 3 T -1 ), sum is the total number of subbasins in the whole simulation space, and the meaning of other symbols can be found in the above.
Step 3: After the above steps, the groundwater flow model has got the sources or sinks in the time step (W in equation S1). Groundwater head and depth of all grid cells are updated after converge of the groundwater flow simulation.
Step 4: The groundwater tables in grid cells are transferred to the hydrologic model.
The average groundwater heads and depths in subbasins will be updated by the area weighted average method, which also based on the spatial relation between the subbasins and the grid cells and the area proportion of grid cell to subbasin. The MODCYCLE can partition subbasins and grid cells flexibly according to the different conditions of landscape and aquifer, reduce unnecessary data redundancy and therefore can be applied to large basins or regions. In addition, the MODCYCLE can be divided into hilly areas and plain areas, and conduct water balance calculations or numerical simulation respectively. Hydraulic information is generally limited for aquifers in hill areas, therefore, it is more appropriate to use water balance calculations, which can reduce the amount of computation.

RSGM module of the MODCYCLE for the simulation of hydrologic cycle in coal mining subsidence areas
Governing equation, spatial discretization, and simulation principle of the conversion relationship between the subsidence area water and the groundwater for the RSGM have been introduced in the main text. The calculation principle of the recharge and discharge for the subsidence areas are presented in this report.

Precipitation and evaporation on the water areas
Precipitation P (L 3 T -1 ) and evaporation E (L 3 T -1 ) on the water areas are calculated using meteorological data and water area:

Flow from upstream
Flow from upstream up Q is the total water of all rivers into the subsidence area which derived from the hydrologic simulation of the MODCYCLE: where i up Q is the flow into the subsidence from the i-th upstream river (L 3 T -1 ).

No-water area runoff
No-water area runoff nw R (L 3 T -1 ) is calculated by the modified SCS curve number where S is the retention parameter, which varies spatially due to changes in soils, land use, management and slope and temporally due to changes in soil water content(LT -1 ), CN is the curve number for a daily time step, nw A is the no-water area(L 2 ), and the meaning of other symbol can be found above.

Water withdrawal
Water withdrawal W (L 3 T -1 ) from the subsidence area is use for irrigation, industry, livelihood, and ecology: where irri w , indu w , live w , and eco w are water withdrawal from the subsidence area for irrigation, industry, livelihood, and ecology, respectively.

Drainage through outlet
Drainage trough the outlet of the subsidence where m is the dimensionless discharge coefficient (-), B0 is the channel breadth (L), g is the gravity acceleration (LT -2 ), and H0 is the upstream total head above the crest (L).

Groundwater discharges from aquifer under the water area
The calculations of recharge and discharge of groundwater in the subsidence area are similar to MODFLOW in simulations of river, reservoir and lake. Groundwater discharges from aquifers into the subsidence area water when the groundwater level is higher than the water level in subsidence area (Fig. S10).

Groundwater discharge from aquifer under the no-water area
When the groundwater table is higher than the water level in the subsidence area, there is a no-water area between the two water levels, where groundwater will discharge into subsidence area (Fig. S11).        Figure S14. (a) Subbasins and (b) main channels in the study area. Maps were created using ArcGIS 9.3 (http://www.esri.com/software/arcgis/arcgis-for-desktop).

Grid for groundwater flow simulation
In our study, the model grid was divided into two layers, representing the shallow and deep aquifers, respectively. Each model layer consists of 132 rows and 276 columns with a uniform size of 500m×500m. Each layer has 36432 cells in which 16166 cells are in the study area and others are inactive (Fig. S15). Figure S15. Grid used in the groundwater flow simulation. Map was created using ArcGIS 9.3 (http://www.esri.com/software/arcgis/arcgis-for-desktop).

Land use and cover
Land use and cover (LUC) types in the study area include agricultural land, forest land, grassland, water area, industrial and mining land, resident land, etc., in which the agricultural land has the largest area, accounting for approximately 80% of the total area, followed by the residential land (about 14%), and other types (less than 5%). Spatial distribution of the LUC types is shown in Fig. S16. Figure S16. Land use and cover (LUC) types in the study area. Map was created using ArcGIS 9.3 (http://www.esri.com/software/arcgis/arcgis-for-desktop).

Soil types
Spatial distribution of soil types in study area is presented in Fig. S17

Weather stations in subbasins
Spatial distribution of weather data measured at stations is processed using a method similar to the Thiessen polygons. Weather information at the station in a subbasin or the nearest one is used as the weather input for the subbasin, as shown in Fig. S18 and

Coal mining subsidence data
Coal mining subsidence data was provided by the China University of Mining and Technology (CUMT). The data are presented in the form contour maps of subsidence depth (Fig. S20). Figure S20. Contour map of the coal mining subsidence depth (in unit of meter). Map was created using ArcGIS 9.3 (http://www.esri.com/software/arcgis/arcgis-for-desktop).

Hydrogeological parameters
Main hydrogeological information and parameters used in the groundwater flow

Weather information
Weather data used in the MODCYCLE, including daily precipitation, maximum/minimum air temperature, solar radiation, wind speed and relative humidity during 2001 to 2010, come from 17 precipitation stations and 3 national basic weather stations mentioned above (Fig. S22).

Boundary conditions of the groundwater flow model
The three rivers surrounding the study area are represented as specified head boundaries in the groundwater flow simulation, and the specified head can be

Crop growth period and parameters
Crop is an important participant in the hydrologic cycle of the study area. The MODCYCLE can simulate the field water moving process associated with crop growth.
The growth periods and parameters of main crops are shown in Table S1.  Table S1. Parameters of main crops for MODCYCLE in the study area.

Water usage
According to the stability of water usage, it can be divided into two types of artificial water: agricultural irrigation water and other water use. Other water usage includes urban/rural domestic water, industrial / thermal power, and ecological/environment water usage, etc.
The applied agricultural irrigation is generally not measured or recorded. In MODCYCLE, the automatic irrigation method is used to estimate the applied agricultural irrigation. The method takes the soil moisture condition as a judgment condition, when the soil moisture in a BSU is lower than a given threshold, the model will automatically pump water from reservoir, channel, or aquifer, etc. for farmland irrigation, then complete the temporal and spatial distribution of the irrigation. Other water usage data from local water resources bulletin is distributed into subbasins according to the area of urban/rural land use distributions.

Supplementary Information 4: Other results of model calibration
The model was calibrated using river stages data from a temporary gage station and groundwater table data from 6 observation wells. In the main text, the comparison between the measured and calculated river stages and groundwater levels at an observation well is given. The comparisons at other 5 wells are presented in Fig. S24.