Abstract
Subsurface temperatures depend on climate and groundwater flow. A lack of observations of subsurface temperature collected over decades limits interpretation of the combined influences of surface warming and groundwater flow on subsurface thermal regimes. Subsurface temperature-depth profile data acquired for Kumamoto Plain, Japan, between 1987 and 2012 were collected and analyzed to elucidate regional groundwater and heat flows. The observed and simulated temperature-depth profiles showed the following: subsurface water flows from northeast to southwest in the study area; the combined influence of surface warming and water flow perturbation produces different temporal changes in thermal profiles in recharge, intermediate, and discharge areas; and aquifer thermal properties contribute more than hydraulic parameters to the perturbation of temperature-depth profiles. Spatial and temporal evolution features of subsurface thermal regimes may be utilized to investigate the influence of surface warming events on subsurface water and heat flows at the basin scale.
1. Introduction
Temperature-depth profiles may be disturbed by both surface warming events and subsurface water flow [1, 2]. Consequently, analyses of observed temperature-depth profiles yield information useful for past climate reconstruction and regional groundwater flow detection. Changes in ground surface temperature propagate downward and impart a thermal signature to the subsurface environment that can be easily measured. Moreover, this information can be used for the reconstruction of past climate [3]. The disturbance of subsurface temperature by surface warming is recorded as a deviation of the temperature-depth profile from the steady state. The observed thermal profiles can be interpreted inversely to reconstruct the climate records of the last five centuries [3–7], extending the limit of knowledge much further into the past compared to that of meteorological observation records.
The studies cited above emphasize the analyses of heat conduction in soils and rocks, ignoring heat convection caused by subsurface water flow. Stallman [8] first presented the basic governing equations of simultaneous transfer of heat and water in a porous medium. Subsequently, this general equation was solved analytically [1, 9–11] and numerically [12, 13] to interpret groundwater flow patterns from the observed temperature-depth profiles. Studies on the interpretation of groundwater flow using heat as a groundwater tracer [1, 9, 10] indicate that downward movement of groundwater flow produces concave-upward thermal profiles, while upward movement produces convex-upward profiles. The vertical distribution of the thermal gradient provides information about the vertical component of groundwater flux. Thus, heat carried by groundwater acts as a tracer for identifying surface water infiltration, flow through fractures, and flow patterns in groundwater basins [14].
Due to severe surface warming caused by urbanization and climate change, researchers during the last two decades have considered the cohesive effects of surface warming and convective groundwater flow on subsurface temperature [1, 2, 15–18]. Due to the influence of surface warming, the temperature-depth profiles of surficial zones shift rightward on a temperature-depth diagram [1]. Although important responses to surface warming commonly occur in shallow aquifers [19], surface warming effects may reach more than 100 m underground, where they may influence groundwater flow systems [18] in urban areas such as Tokyo and Osaka. However, the lack of repeated temperature-depth profile measurements over decades limits interpretation of the temporal and spatial variations in subsurface temperature, which are influenced by both surface temperature and groundwater flow. Subsurface temperature is commonly measured to reconstruct ground surface temperature from the observed temperature-depth profiles, and the signals of climate change or urbanization might be isolated in these analyses [3, 5, 6]. Repeated measurements in different seasons are also needed to detect seasonal fluctuations of subsurface temperature. Little attention has been paid to the observation and analyses of subsurface temperature changes over a scale of years, except studies by Taniguchi et al. [17], Bense et al. [20], and Bense and Kurylyk [21]. The repeated temperature-depth profiles were analyzed to trace the groundwater flows in these studies. This study extends analyses of observed subsurface temperatures in different decades to theoretical analyses using numerical modeling and investigates temporal changes in subsurface temperature under the combined influence of surface warming and regional groundwater flow.
This study was designed to detect temporal and spatial changes in subsurface temperature caused by surface warming and groundwater flow in Kumamoto Plain, Japan. The main objectives were to evaluate changes in groundwater flow systems between 1987 and 2012 based on repeated measurements of subsurface temperature, to elucidate the temporal and spatial evolution of the subsurface temperature regime at the basin scale using a numerical model and to quantify the propagation of surface warming events in aquifers with different groundwater velocities.
2. Data and Method
2.1. Study Area and Data
The investigated aquifer is located in Kumamoto Plain, on the western coast of the island of Kyushu, Japan (Figure 1). Three rivers (Kikuchi, Shira, and Midori rivers) flow through the plain, discharging into the Ariake Sea. The eastern part of this plain is covered with pyroclastic flow terraces and alluvial deposits, while the western area consists of alluvial delta clay [17]. Artesian springs, which form lakes and ponds, are present in the middle portion of the plain [22]. The geological section along line AA’ is presented in Figure 2. The bedrock is granite formed in the Mesozoic era. The aquifers mainly consist of pyroclastic rock from four large eruptions of the volcano Aso 90,000 to 260,000 years ago. These deposits were designated Aso-1, Aso-2, Aso-3, and Aso-4 in order from the oldest volcanic eruption event to the most recent. Among the pyroclastic rock strata is a distribution of deposits from periods of volcanic inactivity. The groundwater system in this region is a two-layer system. The aquifer above the lacustrine clay stratum (red color in Figure 2) is an unconfined aquifer consisting of Aso-4. Below this stratum lies a confined aquifer consisting of pyroclastic flow deposits from Aso-1, Aso-2, and Aso-3 and densely crannied Togawa lava. The primary underground water resources of the Kumamoto region are supplied from this semiconfined aquifer [23].
The annual means of air temperature and precipitation from 1890 to 2012 in the Kumamoto Meteorological Station were 15.8°C and 2011 mm, respectively (Figure 3). Figure 3(a) shows that the air temperature in the study area is almost constant from the 1890s to 1950s, while a significant increase occurs from the 1950s to 2010s. This significant increase in air temperature might be associated with the gradual increase in the urban area on Kumamoto Plain [24]. The growth of the urban population resulted in redevelopment of farmland for urban use, and the urban area within Kumamoto Plain expanded nearly threefold from the 1960s to present [23]. No significant trend was detected in the time series of annual precipitation (Figure 3(b)). However, the annual precipitation amounts in 1987, 2001, and 2012 are anomalous, which might disturb the subsurface temperature by changing the groundwater recharge/discharge rates.
(a)
(b)
Temperature-depth profiles were measured in 28 boreholes in Kumamoto Plain (Figure 1) by the hydrology laboratory of Kumamoto University around July in 2012 to detect local subsurface water flow features. The boreholes are managed by the Kumamoto City government for the purpose of groundwater flow system monitoring. The borehole temperatures were measured with a thermistor thermometer (Technol Seven Ltd., Tokyo, Japan), which can be read to 0.01°C. The thermometer was connected to a monitor with a 500 m cable, and the borehole temperatures were measured by moving the thermometer at 2 m depth intervals. Three sets of measurements were repeated to reduce the error. The wells were drilled between 1977 and 1978, and at least one year had passed between the drilling of the boreholes and measurement. Therefore, the water temperatures measured in the boreholes should represent the temperature of groundwater surrounding the boreholes. However, vertical temperature profiles were observed in some wells (such as well numbers 4, 8, 11, 13, 14, and 16), which might be associated with the vertical convection effects caused by density-driven groundwater flow.
To evaluate the effects of climate change on the subsurface temperature distribution in the study area, we compared repeated measurements of thermal profiles collected over the last 25 years. The measured thermal profiles from 1987 and 2001 were obtained from the studies of Shimano et al. [25] and Taniguchi et al. [17], respectively. A summary of the boreholes sampled for temperature measurements in 1987, 2001, and 2012 is presented in Table 1. Three thermal profiles (well numbers 8, 14, and 25) were observed in both 1987 and 2012. Thermal profiles of 10 boreholes (well numbers 5, 7, 12, 15, 17, 18, 20, 21, 22, and 26) were observed in both 2001 and 2012. The thermal profiles of nine other boreholes (well numbers 1, 2, 3, 4, 6, 10, 11, 19, and 28) were measured all three times.
2.2. Conceptual Model
Taniguchi et al. [1] developed an analytical solution for one-dimensional (1D) heat transfer in homogeneous porous media that accounts for the effects of linear surface warming. Domenico and Palciauskas [10] derived an analytical solution for steady 2D heat flow disturbed by groundwater flow advection. The conceptual model used in this study was constructed based on those theories, as shown in Figure 4. The left and right boundaries were rendered as impermeable hydraulic and insulated thermal boundaries. The bottom boundary was set as an impermeable and constant thermal gradient, while the top was set as a cosine hydraulic head and temporally variable Dirichlet temperature boundary. The governing equation for the subsurface temperature response to 2D thermal flow in a moving medium is expressed as
In this equation, signifies the temperature, denotes the depth from the surface (positive downward), and represents the distance in the horizontal direction (positive rightward). In addition, stands for time, is the thermal diffusivity of the aquifer, and denote the components of fluid velocity in the and directions, respectively, and and represent the volumetric heat capacities of water and the aquifer, respectively ( is the mass-specific heat capacity and is the bulk density). The water table distribution equals (Figure 4). and are constants, and is the length of a cross-section in the horizontal direction. To determine the influence of surface warming, the upper thermal boundary condition was substituted with the observed surface temperature data . The lower boundary condition was expressed as a constant temperature, as shown below.
In this equation, is the surface temperature and is the thermal gradient. Because our focus was on the effects of groundwater flow perturbation on the thermal regime, the influence of subsurface temperature changes on groundwater flow was ignored. Determination of the hydraulic head () distribution was conducted prior to heat flow simulation. The classic hydraulic head assumption [10] was used in this paper, providing a good explanation of the geometrical controls on the spatial distribution of the potential. Under this assumption, the hydraulic head equals [10]. Thus, the hydraulic head in the 2D diagram could be calculated and regional groundwater flow lines were depicted (lines with arrows in Figure 4). Furthermore, hydraulic gradients in the vertical and horizontal directions could be determined as the derivatives of the hydraulic head, which are expressed as
The velocities of the vertical and horizontal groundwater fluxes in each simulation grid can be determined according to Darcy’s law. The analytical solution [10] for the 2D forced convective heat transfer problem in a steady state (with constant surface temperature) was used to represent the initial thermal distribution. This solution is expressed as
Therein, is the hydraulic conductivity. Under the assumptions of these boundary and initial conditions, the differential equation (1) can be resolved numerically using the Crank–Nicolson scheme [26]. Calculations were conducted in a FORTRAN.90 environment. And the results were employed as the initial condition in the subsequent model simulation of 2D heat flow disturbed by both groundwater flowing and surface warming.
3. Results and Discussion
3.1. Classification of Thermal Profiles
Heat can be an effective tracer for the detection of subsurface water flow features [14]. In the 1D model, downward vertical movement of groundwater produces concave–upward thermal profiles, while upward movement results in convex–upward profiles [9]. This basic theory has been used to characterize groundwater recharge or discharge in numerous field studies [1, 2, 12, 27]. Although the scale was extended into two dimensions, similar features to those described for the 1D model were examined in recharge and discharge areas [10]. Furthermore, convection of water flux produces a high subsurface temperature at the same depth in the discharge area. These features are used to depict groundwater flow in Kumamoto Plain, where 28 temperature-depth profiles were collected in July 2012 to detect local subsurface water flow features.
Most observed near-surface borehole temperatures ranged from 18°C to 20°C, while the average air temperature from 1890 to 2012 was 15.8°C. This difference might be due to the following: the ground surface temperature is often 1 or 2 degrees higher than the air temperature [28, 29]; and most observations were collected later than 1987, and the average air temperature exceeded 17°C between 1987 and 2012 in the study area. The observed thermal profiles were classified into recharge, intermediate, and discharge types (Figure 5) according to the observed temperature values. The intermediate type is defined as the thermal profile shows no disturbances from the vertical groundwater flow. It can be observed that abrupt shifts in the measured temperature profile gradients were detected in most boreholes, which might be associated with the changes of soil thermal properties caused by vertical heterogeneity. Statistical information related to the observed temperature-depth profiles is presented in Table 2. Relatively high subsurface temperatures were detected in the discharge profiles. As a result of surface warming events and downward water fluxes, negative thermal gradients were observed in most of the recharge profiles.
(a) Recharge type
(b) Intermediate type
(c) Discharge type
To show the relationship between the observed temperature profiles and the groundwater flow system, each type of thermal profile is shown with different symbols in Figure 6. Visual examination of Figure 6 reveals that most of the recharge and intermediate profiles are distributed in the northeastern part of Kumamoto Plain, whereas the discharge types are located in the southwestern coastal area. The thermal profile distribution implies that groundwater flows from the northeast to the southwest in Kumamoto Plain, as indicated by the hydraulic head measurements in these wells. This groundwater flow system interpretation corresponds well with previous investigations in the study area [17, 30]. Most recharge areas occur near the middle stream of Shira River, where paddy fields predominate. Furthermore, some recharge-type wells are located in the northern upland area. Groundwater discharge in the study area commonly occurs near Ezu Lake in the form of springs, artesian wells, and open lakes.
3.2. Temporal Changes of Temperature-Depth Profiles
To evaluate the effects of climate change on subsurface temperature distributions in the study area, comparison of repeated thermal profile measurements taken during the last 25 years was conducted. Temporal changes in twice-measured thermal profiles are presented in Table 3. Thermal profiles observed in all three years are presented in Figure 7. In Table 3, the subsurface temperature changes in shallow (above 40 m) and deep (below 40 m) aquifers are shown separately. The average absolute temperature change in all measured boreholes below 10 m is 0.2°C (the maximum change is 2.1°C and minimum change is 0°C). Therefore, thermal profiles with absolute temperature changes not exceeding 0.2°C were identified as being in a steady state, whereas those exceeding 0.2°C were identified as increased or decreased. A significant increase was detected in thermal profiles in the shallow aquifer, which was thought to be associated with surface warming events. In the deep aquifer, temperatures in the recharge-type wells generally show an increasing trend. Temperatures in intermediate and discharge wells often remain stable for 25 years. It should also be noted that the subsurface temperature increase between 2000 and 2012 was about 1°C greater than that between 1988 and 2000 and the identical temporal variation can be observed in the air temperature.
(a) Well number 1
(b) Well number 2
(c) Well number 3
(d) Well number 4
(e) Well number 6
(f) Well number 10
(g) Well number 11
(h) Well number 19
(i) Well number 28
These temporal variation features were examined in Figure 7. The triple measurements of temperature-depth profiles in well numbers 1, 2, 3, 4, 6, 10, 11, 19, and 28 are shown in Figures 7(a)–7(i). Well numbers 1, 2, 3, 4, 11, and 19 were identified as recharge wells. Concave–upward thermal profiles were found in these boreholes. The subsurface temperature increments between 2001 and 2012 were greater than those between 1987 and 2001, in accordance with the variation in surface temperature. The thermal profiles in both shallow and deep aquifers for recharge profiles increased continuously, which indicates that the surface warming event propagates deeply with vertical groundwater infiltration. Borehole number 10 was recognized as an intermediate-type well. The thermal profiles remained steady between 1987 and 2012. In the discharge boreholes (numbers 6 and 28), increasing subsurface temperatures were observed in the shallow aquifer and the increasing tendencies diminish as the depth increase. The subsurface thermal profiles above 40 m in borehole number 6 showed increasing tendency between 1987 and 2012, whereas the temporal change in thermal profiles in borehole number 28 showed increasing tendency below 70 m. The differences in temporal changes between thermal profiles in recharge and discharge areas reflect the influence of vertical groundwater flow on propagation of surface warming. The surface temperature signal propagates deeper in the recharge area due to the heat convection effect. Anomalous temperature changes dominated the surficial zone (above 10 m), which were primarily caused by seasonal surface temperature fluctuations.
3.3. Model Simulation
Instead of attempting to replicate field-measured temperature data, our model focused on theoretical analyses of temporal changes in thermal profiles in the study area. A numerical simulation was used to evaluate the influence of surface warming events on the subsurface thermal regime in the groundwater flow region, and the parameters included are listed in Table 4. The investigated aquifer was assumed to be homogeneous to maintain accordance with the initial thermal distribution. The topography parameters, including and , were estimated according to the length and well depth along the cross-section, marked AA’ in Figure 2; thermodynamic parameters, including , , , and , were calibrated according to previous investigations in this area [17, 25, 31]; the surface temperature variation parameters, including and , were calibrated to the observed air temperature records from the Kumamoto Meteorological Station (Figure 8(a)), and these two parameters were utilized to generate the top temperature boundary in the next 100 years; the hydraulic parameters, including and , were estimated from the measured hydraulic head along cross-section AA’ (Figure 8(b)). Then the thermal regimes over the next 100 years were simulated. The contours of the subsurface temperature distributions at t = 0 and 100 years are shown in Figure 9(a) and 9(b), respectively.
(a)
(b)
The groundwater flow regime in the study area was determined based on the annual average water levels in wells distributed along cross-section AA’. Under the classic hydraulic head assumption [10], parameter B was calibrated using the hydraulic head equation of the upper boundary, as shown in Figure 2. Then the hydraulic gradients in the vertical and horizontal directions were determined according to (3) and (4), respectively. Finally, the Darcy fluxes in each simulation grid were determined according to Darcy’s law. The calculated maximum vertical groundwater velocities in recharge and discharge areas were 0.296 and 0.292 m/year, respectively. The calculated groundwater flow pattern in cross-section AA’ showed great similarity to that in Figure 2. Visual examination of Figure 9 indicates that the temperature-depth profiles in Figures 9(a) and 9(b) were both formed by convective water flow. Furthermore, Figure 9(b) showed a warmer shallow aquifer than Figure 9(a), which was the result of surface warming events. As indicated in previous studies [1, 9, 10, 14], downward groundwater fluxes result in concave-upward thermal profiles, while upward fluxes result in convex-upward profiles. These features are shown in Figure 8; when disturbed by groundwater flow, higher temperatures were observed in a discharge area rather than in a recharge area at the same depth. The main differences between the thermal regimes in Figures 9(a) and 9(b) occurred in the shallow aquifer. For clearer visualization, the difference in temperature between Figures 9(a) and 9(b) is presented in Figure 9(c), which shows that surface warming propagates downward to warm the subsurface environment. Using the estimated parameters, surface warming effects reach a depth greater than 100 m in the recharge area, but only 50 m at the discharge site. The deeper extension of surface warm effects in the recharge area might result from convection of vertical groundwater fluxes. Furthermore, the combined effects of surface warming and convective groundwater flow produce a significant negative thermal gradient in the recharge area, as discussed by Taniguchi et al. [1]. Temperature at 100 m depth gradually increases from 13°C to 25°C as the groundwater flows from the recharge area to the discharge area. Most of the observed temperatures in the discharge area range from 20°C to 24°C, which indicates that the seaward region near A’ might also represent a discharge area, and higher temperatures might be observed in this region.
To evaluate the temporal changes in subsurface thermal regime influent due to convective water flow, temperature-depth profiles at the sites of x = 1000, 15,000, and 29,000 m at various times (t = 0, 25, 50, 75, and 100 years) are depicted in Figure 10. For comparative purposes, the numerical solution of a lateral heat convection problem with an identical surface warming rate was used to indicate temporal changes in thermal profiles without the vertical convective force of water flow. Figure 10(a) emphasizes the temporal variation in the subsurface thermal regime as affected by surface warming events. Figures 10(b)–10(d) portray the temporal evolution of thermal profiles in recharge, intermediate, and discharge areas, respectively. Without the vertical convective force of water flow, surface warming effects propagate downward as thermal waves. The amplitude diminishes exponentially with increasing depth below the surface [3]. Temporal changes in the temperature-depth profile were disturbed by convective water flow. In the recharge area (x = 1000 m), warm water in the shallow aquifer moved downward and transferred surface heat to the deeper aquifer through convection. The amplitude of thermal waves damped more slowly than those in the intermediate and discharge profiles during propagation. The thermal profiles in the intermediate area (x = 15,000 m) show identical temporal variation features to those shown in Figure 10(a) in the shallow aquifer. In the discharge area (x = 29,000 m), thermal profiles were influenced by surface warming, showing an increasing trend in the shallow aquifer. Affected by convective water flow, different responses in thermal profiles to surface warming events along the groundwater flow were detected in the basin investigated. The temporal evolution of thermal profiles estimated by the numerical solutions was confirmed by the observational data.
(a)
(b)
(c)
(d)
Using the numerical solution described above, the influences of surface warming and convective groundwater flow on temperature-depth profiles could be analyzed for various differences including in surface warming rates, aquifer thermal diffusivities, and hydraulic conductivities. The surface warming rate represents the long-term ground temperature increase. The subsurface temperature shift driven by different surface warming rates was analyzed to evaluate the propagation of surface thermal waves in the aquifer. In this model, the aquifer thermal diffusivity () and hydraulic conductivity () incorporate the critical determination of thermodynamic and hydraulic parameters for subsurface heat transfer. Furthermore, sensitivity analyses were used to identify the factors influencing model results among the parameters used.
Sensitivity analyses were conducted to address the contributions of various parameters to the model results. The sensitivity coefficient, defined as the ratio of the change in output to the change in input while other parameters remained constant [32, 33], was used to quantify sensitivity. It should be noted that hydraulic conductivity varies by 12 orders of magnitude more than thermal diffusivity but in the present study, hydraulic conductivity variation was confined to 10−5~10−4 m/s because of the geologic composition. Thus, in natural geologic media, it is possible that natural hydraulic conductivity variations exert more thermal influence than natural thermal diffusivity variations. The investigated parameter was shifted by ±66.7%, while other parameters remained constant. Consequently, temperature shifts with the parameter = 9.5 × 10−5/2.17 × 10−7, 9.5 × 10−5/6.5 × 10−7, 9.5 × 10−5/1.08 × 10−6, 3.17 × 10−5/6.5 × 10−7, and 1.58 × 10−4/6.5 × 10−7 were produced. Other pertinent parameters are listed in Table 4. The sensitivity coefficients of these two parameters at x = 1000 m, 15,000 m, and 29,000 m were calculated globally (Table 5). Table 5 shows that aquifer thermal diffusivity is more influential than hydraulic conductivity for simulation of temperature-depth profiles disturbed by surface warming and convective groundwater flow. Furthermore, a comparison of sensitivity coefficients indicates that convective flow and aquifer thermal properties gradually increase in impact on subsurface thermal regimes along the direction of groundwater flow. The global sensitivity coefficient of aquifer thermal diffusivity (0.944) is more than twice the value of hydraulic conductivity (0.439), which indicates that aquifer thermal properties contribute more than hydraulic parameters to the perturbation of temperature-depth profiles.
4. Conclusions
Subsurface temperature-depth profile data in Kumamoto Plain, Japan, between 1987 and 2012 was collected and analyzed using a 2D numerical model to evaluate the regional groundwater circulation, elucidate the temporal and spatial evolution of the subsurface temperature regime at the basin scale, and quantify the propagation of surface warming events in aquifers with different groundwater velocities. The temperature-depth profiles were obtained using thermometers in 28 boreholes located in Kumamoto Plain in July 2012 to investigate regional groundwater circulation features. The observed profiles were classified into three types: recharge, intermediate, and discharge, according to the observed temperature values and shapes of the thermal profiles. Most recharge and intermediate profiles are distributed in the middle stream of the Shira River and the northern upland area, while discharge in the study area commonly occurs near Ezu Lake in the form of springs, artesian wells, and open lakes. Subsurface water flows from the northeast to southwest in the study area.
Observed subsurface thermal profiles collected over the last 25 years were compared to evaluate the effects of climate change on subsurface temperature distributions in the study area. Potentially associated with surface warming events, a significant increase can be detected in thermal profiles in the shallow aquifer. In the deep aquifer, temperatures in recharge-type wells generally showed an increasing trend. Temperatures in intermediate and discharge-type wells commonly remained stable over the 25 years analyzed. The subsurface temperature increments between 2001 and 2012 are larger than those between 1987 and 2001, which agrees with the change in surface temperature. Anomalous temperature changes dominated the surficial zone, which may be caused by diurnal and seasonal surface temperature fluctuations. The observed features of spatial and temporal evolution were explained using a 2D numerical model. However, the heterogeneity was not considered during the model simulation because of the lack of precise geological information, which might influence the precise estimation of groundwater fluxes using heat tracer. Parameter sensitivity analyses indicated that aquifer thermal properties contribute more than hydraulic parameters to the perturbation of temperature-depth profiles.
Data Availability
The datasets generated during and/or analyzed during the current study are available from the corresponding author on reasonable request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
Special thanks are extended to Professor Taniguchi of the Research Institute for Humanity and Nature and Professor Shimano of the Bunsei University of Art for the provision of the preterit subsurface temperature data used for this research. The authors also thank the crew of the hydrology laboratory of Kumamoto University for their constructive comments and suggestions. The quality of this manuscript was improved considerably through their help. This study was funded by the National Key R&D Program of China (nos. 2017YFC1502506 and 2017YFC1502502) and National Natural Science Foundation of China (no. 41501037).