Simulation of the ice thickness of the Heilongjiang River and application of SD models to a river ice model

SD

System dynamic

STNT

Steadily turn negative temperature

CNT

Cumulative negative temperature

RMSE

Root-mean-square error

R²

Coefficient of determination

NSE

Nash–Sutcliffe efficiency

CMDC

China Meteorology Data Service Center

LHS

Latin hypercube sampling

AT

Average temperature

MAT

Maximum temperature

MIT

Minimum temperature

RH

Relative humidity

AP

Air pressure

SVP

Saturation vapor pressure

ITMR

Ice thickness in the middle of river

ITRS

Ice thickness at the riverside

VP

Vapor pressure

The cryosphere is a part of the climate system and stores nearly 80% of all freshwater (Lunt et al. 2004; Serreze & Rigor 2007). River ice is an important part of the cryosphere and has considerable impacts on climate change (Beltaos & Prowse 2009). Global warming enhances the melting of snow and ice, ultimately leading to spring floods caused by ice jam (Blackburn & Hicks 2003; Pagneux et al. 2011). River ice is a complex system that is affected by hydrodynamic, mechanical, and thermal processes under the influence of meteorological and hydrological conditions (Jin et al. 2019). It has significant meaning for understanding ice growth and melting, conducting research on ice thickness as a system, and calculating river ice thickness. The Heilongjiang Province is located in Northeast China, and ice is present in nearly the entire province for months. The upstream of the Heilongjiang River often experiences spring floods due to the ice dam, which causes substantial economic damage and massive casualties. However, the Heilongjiang River lacks monitoring equipment and data because of its geography, and no river ice thickness models are available for the river. An applicable river ice thickness model in the Heilongjiang River should be established for forecasting spring flood disasters.

Considerable efforts have been made to simulate ice thickness. Numerical simulation has been developed and matured internationally. For example, Stefan (1891) established the heat balance equation of ice thickness on the assumption that the surface temperature of ice is equal to the atmospheric temperature and deduced a model of freezing days simply by simulating ice thickness. Stefan's model shows that only at the temperature after the obtainment of ice can the thickness be reduced; in fact, the measured data showed that when the temperature was lower than 0 °C, ice thickness began to decrease. To simulate the melting of ice thickness, Shen & Yapa (1985) added the day factor. The accumulation of temperature began at the moment ice formed in Shen's degree–day model, but river ice is actually not formed at the moment temperature decreased.

On the basis of the simple numerical model, semiempirical methods are presently developed into the most important numerical simulation and prediction methods for simulating river ice thickness (Prowse & Beltaos 2002; Shen 2002; Beltaos & Prowse 2009; Shen 2010). They are based on the physical mechanism of river ice, such as growth and melting, and physical equations, such as hydraulics, thermodynamics, and hydrology (Shen 2003a, 2016; Shen & Chiang 1984). Scientists have performed mathematical simulation of ice growth and melting in accordance with the thermal mechanism and established semiempirical ice thickness simulation models, which consider the interaction between ice thickness and heat. For instance, Ashton (1979), Greene (1981), Shen & Yapa (1984), Yapa & Shen (1986), Sarraf (1990), Sarraf & Plouffe (1991), Zhang (1992), and Sarraf & Zhang (1996) developed river ice models by using thermal effluent, and the convection–diffusion of thermal energy in the river water was simulated numerically. Shen & Ho (1986) and Hirayama (1986) developed one-dimensional extension models, which provided the description for two-dimensional models. Svensson et al. (1989) developed a two-dimensional river ice model of the thermal growth of border ice, which was based on the conservation equations for mass, momentum, and thermal energy. Although the river ice models are mature and applied in many regions, only a few studies have regarded river ice as a whole system to investigate. A river ice system has dynamic changes. In a thermal-dominant river ice system, every element has an influence on heat flux and ice thickness. In the melting stage, ice thickness is influenced by ice melting capacity, and the ice melting capacity of the next time step is influenced by the ice thickness of this time step. The interaction between river ice and elements constitutes a river ice system. The complex interrelationships among different elements in the river ice system cannot be shown.

Most river ice models were built and applied specifically for the Yellow River in China. For example, Fu et al. (2014) developed the YR1DM river ice model on the basis of the ice condition in the Ning–Meng reach of the Yellow River in China. Numerical ice flood models were developed by Mao et al. (2014) and Wang (2018) for the Ning–Meng reach, leading to enhanced forecasting of the ice regime and improved decision support for upstream reservoir regulation and taking appropriate measures for disaster risk reduction. The Heilongjiang River is different from the Yellow River because of its special location and high latitude. The models built for the Yellow River cannot be used for the Heilongjiang River. Most studies on the Heilongjiang River ice concentrated on ice condition description (Yu 2006; Guo et al. 2009). Although some studies (Li et al. 2009) based on experiments have been performed, data monitoring, river ice simulation, and ice disaster forecasting for the Heilongjiang River are lacking because of its severe environment and backward economy. Few river ice models have been built and applied to the Heilongjiang River because of its specific location and the inadequate data.

As stated, there are some problems in current research. First, the lag days between river ice and temperature are not considered in the degree–day model; second, the river ice is not regarded as a whole system to investigate; and lastly, there is no applicable river ice model in the Heilongjiang River. The purpose of this study is to improve the degree–day model and establish a river ice del as a system in the Heilongjiang River. To solve the above-mentioned problems in previous models, the correlation between ice thickness conditions and meteorological elements is analyzed using statistical methods. A degree–day model considering the lag days between temperature and river ice formation is built for the ice growth model. A Russian river ice model is applied to the Heilongjiang River directly in the melting phase. A river ice simulation model coupling with the system dynamic (SD) simulation approach is built, and sensitivity analysis and uncertainty analysis are conducted in SD models.

In general, the river ice process includes the growth and melting of ice shown in Figure 1. At the beginning of winter, the water surface cools to 0 °C some days after the air temperature is below 0 °C. Future heat loss will lead to a below zero water temperature, then the river ice begins to freeze. Future cooling will lead to rapid thermal growth of river ice. During this time, air temperature plays an important role in river ice growth. When the air temperature reaches the minimum value, the river ice grows slowly and even stops growing because the loss of energy reaches its maximum. At the beginning of spring, the river ice begins to melt. The thermal erosion of the river ice cover and the deterioration of river ice cover through internal melting by solar radiation lead to the reduction in cover strength. In the melting process, the ice cover is destroyed, while the upper and bottom surfaces of the ice cover are melted. Solar radiation is absorbed by the ice, and heat is converted into a surface of crystals and heterogeneous inclusions, which lead to a large interlayer between pores and liquid. When the ice melts from inside, the ice begins to change from solid phase to liquid phase. The presence of liquid phase in the melting ice is the main reason for the strength reduction of the ice sheet. When the liquid interlayer among the crystals reaches a thickness in which the adhesion is completely lost, the strength of the ice is completely lost, that is, the ice is separated into individual crystals. From the preceding analysis, temperature and heat flux play important roles in the river ice growth and melting (Morse & Hicks 2010; Shen 2010; Ashton 2011; Prowse & Marsh 2011; Wang et al. 2013; Beltaos 2015).

The lag days between river ice thickness and temperature were confirmed by calculating the person coefficient beginning on day 1 and moving toward day 15 and choosing the most significant correlation between river ice thickness and temperature. The steadily turn negative temperature (STNT) date was after the day the temperature became below 0 °C, and lag days served as the first day of CNT. The relationship between river ice formation date and special temperature date was concluded using a statistical method. From the conclusions, the last day of CNT and when to melt were confirmed. After the freezing and melting dates were determined, the coefficient of the freezing model was determined using the least square method.

The degree–day model shows that ice thickness is reduced when the temperature is above zero. However, the measured data show that when the temperature is lower than 0 °C, ice thickness begins to decrease. Hence, in this study, the degree–day model is only used in the freezing model.

The turbulent exchange between ice surface and air is the empirical formula of sensible heat flux, which is calculated using the temperature difference between water and air, wind speed, and other factors.

We evaluate the quality of the river ice model by using root-mean-square error (RMSE), coefficient of determination (R²), and Nash–Sutcliffe efficiency (NSE).

The SD approach is a well-established methodology for understanding, visualizing, and analyzing complex dynamic feedback systems by using stocks, flows, internal feedback loops, and time delays in SD models (Forrester 1958). The SD approach is commonly applied in water resource (Simonovic 2002) and hydropower safety (King et al. 2017). It is especially useful for the physical relationship between river ice and meteorological factors because it has a detailed representation of how a river ice model operates. A professional SD software package, Vensim PLE (www.vensim.com), is used to simulate the integrated system.

On the basis of the river ice process, a river ice dynamic system is built (Figure 2), which represents the status of the river ice model structure. The design of the river ice model includes two subsystems: river ice growth and river ice melting. It consists of 4 level stocks, 22 auxiliary variables, and 11 constant parameters. The key variable in this model is ice thickness. In accordance with the system stock-flow diagram and collected data, the rate, state, and calculation formulas of auxiliary variables involved in the model are shown in Table 1.

Sensitivity analysis is a critical tool to test a model if changes in parameters or inputs lead to important changes in simulation. It may provide information on which of the model parameters or inputs are important for the simulation output. The analysis includes numerical sensitivity, policy sensitivity, and behavior sensitivity. Numerical sensitivity evaluates the sensitivity of output values to the change in model assumptions, whereas behavior mode sensitivity determines the sensitivity of output behavior to the alterations in the model. Policy sensitivity is defined as the change in desirability or suitability of an existing policy when a change in model parameters or structure occurs. Behavior pattern sensitivity is to cover the effect of changing inputs on the output patterns of the model. The step of sensitivity analysis of SD models starts with the determination of parameter ranges and distributions. Then, the estimation of behavior measures and regression analysis should be performed sequentially.

The distribution function and range of each parameter are determined first. Usually, of parameter values are used as distribution ranges. To determine the parameters better, of values are selected in this research. A sampling strategy is selected after determining the parameter ranges. Sampling strategies include random sampling, stratified sampling, Taguchi method, and Latin hypercube sampling (LHS). Among them, LHS is an effective technique. In this study, LHS is used in sensitivity analysis. Then, the estimation of behavior measures should be performed (Hekimoglu & Barlas 2016).

In sum, input parameters are determined, and 50% of parameter values are selected first. The changing input model behavior curve is compared with the ice thickness model behavior curve, and model sensitivity is estimated through changing different inputs. Then, of values and LHS are selected. Lastly, the estimation of behavior is performed. The larger the bandwidth is, the more sensitive it is; the smaller it is, the less sensitive it is.

Uncertainty analysis aims at quantifying the variability of the output that is due to the variability of the input. In the research, input variables, such as ice surface temperature, wind speed, AP, VP, temperature, ice surface VP, longwave radiation, and shortwave radiation, are increased 25, 50, and 75%. The ice thickness range is performed based on different changed input variables, and the uncertainty elements are confirmed.

The SD simulation model has been built for the Mohe stream (E122°21′50″, N53°28′48″), a tributary of the Heilongjiang River (Amur River), which is the largest international boundary river in China (Figure 3). The basin has a cold temperate continental monsoon climate, characterized by uneven seasonal distribution, large annual variation in temperature, and different humidity in winter and summer. Given that the Siberian cold air activity is active, the monsoon climate is highly characterized in this region. The winter is long, cold, and dry. The temperature generally drops to 0 °C from mid-to-late October.

The Mohe experiment station consists of ice monitoring station and automatic meteorological monitoring station. When river ice is frozen completely, automatic ice monitoring equipment is installed at two sites on the surface of river ice, one of which is close to the river bank and the other one to the river center. Automatic solar radiation monitoring equipment is installed at the height of 2 m in 10 m anemometer mast at the Mohe hydrological station. Automatic meteorological monitoring equipment is installed at the Mohe hydrological station. An ultrasonic wind sensor is installed at the height of 10 m, and a rotary wind sensor is installed at the height of 2 m in 10 m anemometer mast. Temperature, humidity, and pressure sensors are installed at the height of 2 m in 10 m anemometer mast. Longwave and shortwave radiation sensors are installed on the ground.

Elements of the Mohe experiment station, including hourly and shortwave and longwave radiation from 27 November 2016 to 31 October 2017, are recorded. Hourly temperature, elements of hourly 2 and 10 m wind speed, hourly relative humidity (RH), hourly dew point temperature, and hourly local AP from 1 January 2014 to 30 October 2016 are measured. Hourly ice temperatures from 1 January 2014 to 2 April 2014, 9 December 2014 to 8 April 2015, and 1 January 2016 to 2 April 2016 are also measured. The ice condition data, such as ice frozen date and ice melting date, from 2013 to 2016 are obtained from the Hydrology Bureau of Heilongjiang Province. The meteorological data of the China Meteorology Data Service Center (CMDC) from 2013 to 2016 are acquired online at http://data.cma.cn/en. These data are real and exact.

In this experiment, meteorological data are insufficient because of the environment, monitoring equipment, and other factors, and they should be interpolated. Figure 4 shows the temperature monitoring data, including average temperature (AT), maximum temperature (MAT), and minimum temperature (MIT), during 2013–2016.

Figure 5 shows the RH and AP during 2013–2016. Missing data are determined among the monitoring data. In addition, some input data are not monitored directly; for example, VP, sensible flux, and latent flux are calculated using other elements.

AT, MAT, and MIT are interpolated through the temperature of the experiment station based on the correlation between the data in the experiment station and the data of the CMDC from October to April in 2014–2016.

The interpolation of RH and AP is based on the correlation between the data of the experiment station and the data of the CMDC from October to April in 2014–2016.

The interpolation of wind speed data is divided into four processes. (1) Whether the measured sample of the missing data contains dead pixels is determined. The absolute mean method is used as the detection method. For the zero-mean sequence, the mean value of data samples is obtained, multiplied by the coefficient k, then W is generated. When ⁠, is considered the dead pixels. If the data sample is a nonzero mean data sample, it is zero-centered, that is, ⁠. After testing, the formula is ⁠. In general, considering 4–5 for k is reasonable. To retain the original data as much as possible, k is set to 5. (2) To judge whether the sample data before the missing data have stationarity, a model is built and analyzed in accordance with a stationary time series model, and a set of wind speed values of the missing part is interpolated using the forward method. (3) To utilize the sample data after the missing data, the sequence of the sample data after the missing data is reversed, the inverse sequence is modeled and analyzed in accordance with a stationary time series model, and the wind speed value of another group of missing data is interpolated using the inverse method. (4) The corresponding points of the two interpolated wind speed values in steps 2 and 3 are averaged to serve as the measured wind speed interpolation values of each point of the missing part.

Figure 6 shows that the AT, MAT, and MIT of the experiment station and the CMDC have a good liner correlation. The R² are 0.9777, 0.9215, and 0.9348. The calculation equations are y = 0.9805x + 0.3657, y = 0.998x − 7.212, and y = 0.9638x + 6.6936, where y represents the temperature data of the experiment station and x represents the temperature data of the CMDC.

Figure 7 shows that RH and AP have a good linear correlation between the experiment and the CMDC. The R² of RH and AP are, respectively, 0.7853 and 0.9883, and the fitting formulas are, respectively, y = 0.9368x + 7.2371 and y = 1.0643x − 44.139.

Figure 8 shows the river ice thickness process at the river center and riverside from 2013 to 2016. River ice thickness grows rapidly at first because holes should be drilled in the ice when the experimental equipment is installed and the hole at this point needs to be filled. Subsequently, the river ice process shows relatively slow growth before mid-February in 2014, the beginning of February in 2015, and mid-February in 2016. Then, river ice grows slowly and stops growing. Lastly, the river ice begins to melt rapidly.

To compare the relationship between temperature and thickness data, the lag days between daily thickness and daily AT are analyzed and selected. Ice thickness includes ice thickness in the middle of river (ITMR) and ice thickness at the riverside (ITRS). Owing to the inclement working conditions, the equipment did not work during 2015–2016. Only the lag days between daily thickness and daily AT in 2013–2015 are considered. Table 2 shows that coefficients increase and then decrease. From the correlation between thickness and the temperature of 10 lag days, the ITMR of 2014, the ITRS of 2015, and the ITMR of 2015 are the most significant. The coefficient between ITRS thicknesses of 2014 and temperature of 12 delay days reaches the maximum, and the coefficient between ITRS thicknesses of 2014 and temperature of 11 delay days is very close to the coefficient of 12 delay days from the analysis. Thus, the thickness and temperature of 10 delay days have the most remarkable correlation. The 10th day is defined as the STNT date when the temperature is below zero for 10 consecutive days.

As shown in Table 3, maximum ice thickness is determined after AT reaches the minimum. When the highest temperature is greater than 0 °C, ice in the middle of the river begins to melt, and riverside ice begins to melt 1 day in advance. Consequently, ice at the riverside and ice in the middle of the river melt approximately on the same day. In the establishment of the ice thickness growth model, the AT STNT is considered the starting point when the ice thickness starts to grow. When the highest temperature is greater than 0 °C, the ice begins to melt.

Table 4 shows the results of model performance evaluations. Comparison of the simulation value with the measured value indicates that the R² of ITMR and ITRS are 0.67 and 0.69, respectively; the NSEs of ITMR and ITRS are 0.72 and 0.70, respectively; and the RMSEs of ITMR and ITRS are 6.50 and 6.84, respectively. The coefficient shows that the model has good performance.

Sensitivity analysis is conducted on the basis of ITMR during 2013–2014. Temperature, shortwave radiation, longwave radiation, wind speed, atmospheric pressure, and VP at 2 m are selected. One variable is changed, while other variables remain unchanged. Then, river ice thickness is simulated one by one. Figure 10 shows that changing different variables affects ice thickness in varying degrees. Comparison of the changing curves with the model curve indicates a large gap between the simulation curve of changing temperature and the original simulation curve. Gaps appear in the melting process between the simulation curve of changing longwave and shortwave radiation and the original simulation curve, and a small gap exists between the models of other changing elements and the original simulation. The change in temperature has the greatest influence on ice thickness, followed by the change in shortwave radiation, and other inputs have minimal influences on ice thickness.

For understanding which element has a considerable influence on the ice thickness behavior, detailed sensitivity analyses are conducted. Table 5 shows the parameter distributions.

Figure 11 shows the results of the sensitivity analyses. The bandwidth of the changing temperature model is the largest, the bandwidth of changing shortwave radiation is larger in the melting process, the curves of changing longwave radiation and wind speed have a minimal color change in the melting process, and no change occurs when AP and VP at 2 m vary. The graphs show that temperature has the most sensitive behavior, shortwave radiation and longwave radiation have high sensitivity in the melting process, and other inputs have minimal sensitivity. The conclusions are the same as previous conclusions.

During ice growth and melting, temperature, as a result of thermal state, influences the ice thickness process directly. Hence, ice thickness is mostly sensitive to temperature. Shortwave radiation and longwave radiation have no influence on the ice growth process. On the contrary, shortwave radiation and longwave radiation, as the main heat sources, have substantial influences on the ice melting process. Sensible flux and latent flux have a minimal influence on the ice melting process. To verify this conclusion, different flux sensitivities are analyzed for ice thickness behavior. First, this study determines which one of upper flux and lower flux affects ice thickness behavior more. Figure 12 shows that the graph of changing upper flux has bandwidth in the melting process and the graph of changing lower flux has no bandwidth, which means that changing upper flux indicates higher sensitivity.

After determining that upper flux influences ice thickness behavior more, this study identifies which flux mainly affects upper flux. Figure 13 shows that the graph of changing shortwave and changing longwave has bandwidth in the melting process and the graph of the two other fluxes has no bandwidth. That is, shortwave and longwave mainly affect the upper flux.

In the ice melting process, the upper flux mainly affects the ice thickness behavior, and shortwave radiation and longwave radiation considerably influence the upper flux. Thus, ice thickness is mostly sensitive to shortwave radiation and longwave radiation. Air temperature characterizing the thermal state also affects ice growth and melting.

Uncertainty analysis is conducted on the basis of ITMR during 2013–2014. As shown in Figure 14, increasing temperature causes the largest change of range and the box plot is much higher than the boxplot of ITMR of 2013–2014; temperature increased by 75% causes the AT to increase approximately 30%, which is the most significant change; temperature increased by 50 and 25% causes the AT to increase about 20 and 15%. Increasing longwave radiation and shortwave radiation cause slight changes of range. Through the analysis, the change of temperature is the most important factor causing the change of model. The conclusion is the same as sensitivity analysis.

The understanding of the processes of ice growth and melting has a considerable influence on building an ice thickness model and combating climate change. In particular, studying the ice thickness process of the Heilongjiang River and establishing a river ice model are necessary, which provides an important theoretical foundation on the spring flood forecast for the Heilongjiang River. In this study, river ice conditions and the relationship between elements and ice thickness are summarized through statistical methods, river ice thickness SD models at the riverside and in the middle of the river are built, and sensitivity analyses for ice thickness behavior are performed through Vensim.

With the statistical methods, this study determines that thickness and temperature of 10 delay days have the most remarkable correlation. The maximum ice thickness occurs after AT reaches the minimum. When the highest temperature is greater than 0 °C, ice in the middle of river begins to melt, and riverside ice begins to melt 1 day in advance. ITMR and ITRS can be considered to melt on the same day.

The SD models present the relationship between elements and ice thickness behavior directly. Comparison of the simulation value with the measured value indicates the good performance of the model.

Sensitivity analyses demonstrate that shortwave radiation and longwave radiation mainly influence upper flux, and upper flux mainly affects the ice thickness behavior. Air temperature characterizing the thermal state also affects ice growth and melting. Uncertainty analyses show the temperature is the most important factor to influence the model and it is the same as sensitivity analysis results.

This study only establishes a simple ice growth and melting model through some statistical methods. The mechanism of the Heilongjiang River ice cannot be understood totally. In the future work, additional efforts should be spent on studying the ice growth physical process deeply, determining the time point of ice growth and melting by using physical methods, and improving simulation accuracy. Further coupling the river ice physical process with SD models can provide a good perspective in understanding climate change.

This work was supported by the National Key R&D Program of China (Grant No. 2018YFC1508001), the Special Fund of State Key Laboratory of Hydrology-Water Resources and Hydraulic Engineering (Grant Nos 20195025612, 20195018812, and 520004412), the program of China Scholarships Council (No. 20180671), and the Fundamental Research Funds for the Central Universities (No. 2017B611X14).

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Data cannot be made publicly available; readers should contact the corresponding author for details.