Influence of Variations of Hydrothermal Conditions on Normalized Difference Vegetation Index in Typical Temperature Zones Along the East Coast of China

The eastern coastal areas of China span multiple climatic zones, and the impacts of climate warming on their ecological environment show regional differences. In this research, the Normalized Difference Vegetation Index (NDVI) was used as the indicator to characterize the ecological environment, and selected Guangdong, Jiangsu, and Liaoning as its typical research areas. The authors selected the NDVI, average temperature, and precipitation data of the yearly growth season, respectively, from 1982 to 2016. This study adopted the copula functions model based on Markov Chain Monte Carlo to carry out the research of bivariate joint distribution so as to calculate the joint probability, the joint exceedance probability, the joint return period and the co-occurrence return period. The results showed that 1) the temperature and precipitation in the three regions were respectively related to the NDVI sequence showing the characteristic that was correlated at the upper tail and asymptotically independent at the lower tail, which demonstrated that the temperature and precipitation had little effect on NDVI when they reached their minimum values, and the temperature and precipitation had obvious effect on NDVI when they reached their maximum values. 2) The shorter the return period was, the wider the ranges of the climate factor and the NDVI were, showing that when the climate factor was constant, the probability of the NDVI having a shorter return period was higher. The greater the climate factor was, the longer the return period was, indicating that the probability of plant growth inhibition was higher when the climate factor exceeded a certain threshold. 3) The suitable temperature and precipitation for vegetation growth in the three regions gradually decreased from south to north. These results provide some theoretical guidance and scientific foundation for the protection of regional ecological environment and enhance the understanding of the impact of climate change on the ecosystem.

The eastern coastal areas of China span multiple climatic zones, and the impacts of climate warming on their ecological environment show regional differences. In this research, the Normalized Difference Vegetation Index (NDVI) was used as the indicator to characterize the ecological environment, and selected Guangdong, Jiangsu, and Liaoning as its typical research areas. The authors selected the NDVI, average temperature, and precipitation data of the yearly growth season, respectively, from 1982 to 2016. This study adopted the copula functions model based on Markov Chain Monte Carlo to carry out the research of bivariate joint distribution so as to calculate the joint probability, the joint exceedance probability, the joint return period and the co-occurrence return period. The results showed that 1) the temperature and precipitation in the three regions were respectively related to the NDVI sequence showing the characteristic that was correlated at the upper tail and asymptotically independent at the lower tail, which demonstrated that the temperature and precipitation had little effect on NDVI when they reached their minimum values, and the temperature and precipitation had obvious effect on NDVI when they reached their maximum values. 2) The shorter the return period was, the wider the ranges of the climate factor and the NDVI were, showing that when the climate factor was constant, the probability of the NDVI having a shorter return period was higher. The greater the climate factor was, the longer the return period was, indicating that the probability of plant growth inhibition was higher when the climate factor exceeded a certain threshold. 3) The suitable temperature and precipitation for vegetation growth in the three regions gradually decreased from south to north. These results provide some theoretical guidance and scientific foundation for the protection of regional ecological environment and enhance the understanding of the impact of climate change on the ecosystem.

INTRODUCTION
As the main body of the terrestrial ecosystem, vegetation is an important part connecting the atmosphere, soil water, and energy cycle (Huryna and Pokorný, 2016;Jin et al., 2017), and thus has great significance to regional climate regulation, surface energy balance and soil and water conservation, etc. Duveiller et al., 2018;He et al., 2018). Climate resources such as light, heat, and water are the basis of vegetation growth, and a vegetation variation is a key indicator to measure a regional environment and a climate change (Li et al., 2017;Zewdie et al., 2017;Wan et al., 2018). The interaction between the climate change and the vegetation has aroused widespread interest of scholars and become a research hotspot in recent years (Huete, 2016;Seidl et al., 2017).
There is a bidirectional feedback mechanism between the climate and the terrestrial ecosystem: a variation of a hydrothermal condition causes a change in a temporal and spatial distribution pattern of the terrestrial ecosystem (Huang et al., 2016;Huang et al., 2019); and the change of the ecosystem patterns affects a regional carbon cycle and water cycle process, which in turn reacts on the regional climate (Sippel et al., 2018;Quan et al., 2019). Under the background of climate warming, temperature, precipitation, radiation, and other factors in different regions are obviously different, and the influence of the climate change on the vegetation in different regions varies. The study held that the changes in factors such as temperature, water, and light affect the vegetation by participating in photosynthesis, respiration, and transpiration of the vegetation (Zandalinas et al., 2018;Dusenge et al., 2019). Meanwhile, the vegetation is sensitive and adaptable to the climate, resulting in differences in response rules of the regional vegetation to the climate change. Therefore, the interaction between the regional climate and the vegetation belongs to a nonlinear category (Wen et al., 2017).
Normalized difference vegetation index (NDVI), as an important representative indicator of the growth of surface vegetation, is widely used in the study on dynamic changes of the vegetation (Xu et al., 2016). China's climate resources show obvious regional differences, vegetation types are diverse, and a relationship between the vegetation changes and the regional climate is heterogeneous. Existing researches on the relationship between the climate changes and the NDVI mainly focused on the overall correlation between various elements and the NDVI on the regional scale, or the interrelationship between them on the multitime scale; while a relationship between different action subintervals of climate elements and their corresponding NDVI has not been excavated. Thus, the dynamic correlation analysis of climate factors and the NDVI cannot be achieved, nor can the optimal hydrothermal condition for the NDVI-based growth under the conditions of the regional climate and climate variability differences be described. Besides, no zoning discussion has been done according to the differences in the climatic zones.
Traditional correlation coefficients have been adopted as a popular method to study the correlation between variables, but cannot distinguish the dynamic correlation between the variables in different intervals or exclude the possibility of "pseudo correlation" between the variables (Ramos-Cordoba et al., 2016). In addition, although linear regression fitting is also commonly used in the modeling of multi-element interaction, the application prerequisite of such method is to assume that there is a linear relationship between the variables, which requires the data distribution to conform to the stationarity. In the context of climate warming, the temperature, precipitation, and NDVI all have the changing trend, and the data no longer meet the precondition of the hypothesis of the stationarity. These problems can be effectively solved by Copula function theory (Abdi et al., 2017). The Copula function can reflect a correlation structure of random variables independently from the marginal distribution of random variables. Any marginal distribution can be constructed into a joint distribution through the Copula function. There will be no information distortion and loss as the information of a single variable is included in the marginal distribution (Lee, 2018;Papalexiou and Serinaldi, 2020). The multi-dimensional joint distribution was built on the basis of the marginal distribution and correlation of the variables, and a probability distribution function was established through the Copula function theory to depict in detail the dependence of temperature-NDVI and precipitation-NDVI in different value ranges, and the probability of was obtained. This research has studied the influence of hydrothermal condition variations on an NDVI, explored the spatial differentiation rule of temperature, precipitation, and NDVI changes in different temperature zones and quantified the effect of hydrothermal combination to the vegetation growth, which can provide a theoretical basis for a deep understanding of the interaction between regional scale climate and its ecosystem. And this will be conducive to an active response to the climate changes, and serving environmental protection and decision-making plan of coping with the climate warming better.
This research consisted of five sections: Materials and Methods presented an outline introduction to the research areas and the used data and method. Through the establishment of joint probability distribution and exceedance probability distribution, Results Analysis analyzed the interaction relationship among the NDVI, the average temperature and the precipitation, and calculate the corresponding joint return period. Discussion section focused on the causes of NDVI dynamic change, and analyzed the response of NDVI to temperature and precipitation in detail. Finally, the conclusions were presented in Conclusion.

MATERIALS AND METHODS
An Outline Introduction to the Research Areas Figure 1 shows the geographical location of the study area. Guangdong is a southern province of China with a tropical and subtropical monsoon climate. Therefore, it has long summers and Frontiers in Earth Science | www.frontiersin.org September 2020 | Volume 8 | Article 574101 2 warm winters. The average annual temperature is from 20.4 to 23.1°C. Jiangsu is an eastern-central province of China and is situated in a transitional area between the temperate and subtropical zone. Generally, toward the south of the Huaihe River and sub-northern irrigation canal, humid subtropical climate zone is experienced, whereas, warm temperate climate zone is experienced toward the north. The annual mean temperature is from 13 to 16°C, increasing from the northeast to southwest area. Liaoning is a northern coastal province in Northeast of China and is located in temperate monsoon climate zone. The temperature in this area is characterized by uneven spatial distributions, decreasing from the plain to the mountain area. The annual average temperature in Liaoning is from 7 to 11°C (Figure 1).

Meteorological Data and Preprocessing
The daily temperature records from January 1, 1982 to December 31, 2016 were collected for the three provinces from 90 Frontiers in Earth Science | www.frontiersin.org September 2020 | Volume 8 | Article 574101 3 meteorological stations provided by the China Meteorological Administration (https://data.cma.cn/site/index.html). The stations were selected on the basis of the length of the time series, data completeness (missing values less than 5%), and spatial coverage. A series of quality control tests were applied to identify outliers at all stations and were marked with a quality control flag. The average daily temperatures were then extracted. It should be noted that the level of data completeness in the study exceeded the minimum requirement of 95%. At this level, scaling indices and scaling behaviors of the time series were not affected. The monthly regional average temperatures and precipitations of Guangdong, Jiangsu, and Liaoning provinces were calculated respectively according to the daily average air temperature and precipitation.

Normalized Difference Vegetation Index Data and Preprocessing
Two kinds of NDVI datasets including the latest updated version of the third generation Global Inventory Monitoring and Modeling System (GIMMS, available at https://ecocast.arc.nasa.gov/data/pub/ gimms/3g.v1/ as nc4 files) and MODIS NDVI were both used in the study. The GIMMS NDVI 3g.v1 is generated from National Oceanic and Atmospheric Administration's Advanced Very High Resolution Radiometer data, and the spatial resolution is 1/12°. Its temporal resolution is 15-days intervals with span from 1982 to 2013. This data set eliminates the effects of volcanic eruption, solar altitude angle, and sensor sensitivity changes over time, making the quality of the GIMMS data set superior to other NDVI data sets. The GIMMS data is the longest time series of NDVI data at present. It has a better correlation with other high-resolution data sets and has been extensively used in the world. The MODIS NDVI data (available at https://ladsweb.modaps.eosdis.nasa.gov/) is derived from the MODIS vegetation index product developed by NASA MODIS Land Product Team according to the unified algorithm. The MODIS NDVI used in the study is MOD13A2, i.e., the vegetation index with a resolution of 1 km and temporal resolution of 16 days.
The data sets of GIMMS NDVI and MODIS NDVI in Guangdong, Jiangsu, and Liaoning were obtained respectively after NDVI data were preprocessed through quality inspection, image mosaic, subset extraction, cropping, format, and projection conversion and so on. In order to get rid of the differences in temporal resolution there between, the Maximum Value Composite (MVC) was used to obtain monthly scale NDVI data, so as to further remove the impact of clouds and aerosols (Liu, 2017;Seong et al., 2020), and reduce the effect of the phenological cycle in the month (Fensholt and Proud, 2012). The averaged NDVI for growing season was calculated for analysis. The pixel values with an average of growing season NDVI <0.1 were masked as non-vegetated areas.
Due to the different spatial resolutions of the GIMMS data and the MODIS data, this research has adopted the GIMMS data with the time span from 1982 to 2013, and the MODIS data with the time series from 2001 to 2016. The correlation analysis of the monthly NDVI data was carried out in terms of data in overlapping periods. The correlation coefficient of the two data is 0.927, and the linear regression equation is NDVI GIMMS 0.781 * NDVI MODIS − 0.0469(r 2 0.9016, p-value < 0.05). Using the regression equation and combining MODIS-NDVI data from 2014 to 2016 to interpolate monthly GIMMS-NDVI data, the time span of GIMMS NDVI data set was extended to 1982-2016.

Maximum Value Composites
The MVC was proposed by Holben (1986). The specific formula of MVC (Kundu et al., 2018) is as follows: where NDVI i refers to the NDVI in the ith month or the ith year; and NDVI ij refers to the NDVI data on the jth 15-day in the ith month or on the jth month in the ith year.

Copula Function Theory
The Copula function raised by Sklar (Li et al., 2018) can use the marginal distribution and correlation framework to build a multidimensional joint distribution Copula function model (Salvadori and De Michele, 2007). The study selected eight Copula function clusters (Clayton, 1978;Genest and Favre, 2007;Li et al., 2013;Sraj et al., 2015), including 1) BB1, 2) Clayton, 3) Frank, 4) Gaussian, 5) Gumbel, 6) Joe, 7) t, and 8) Tawn. These eight forms of the Copula functions have always been common choices for related models due to their performances. The selected Copula functions were used to establish a two-dimensional joint distribution of climate elements and NDVI.

Parameter Estimation
The parameters of the Copula function were calculated by the non-parametric estimation method (Genest and Rivest, 1993). This technique is mainly related to the parameter θ of the Copula. See Table 1 for various Copula function forms. Formula (2) shows the relationship between θ and τ (Kendall correlation coefficient). By calculating τ from the measured data, the corresponding joint distribution parameters can be obtained (2)

Verification and Evaluation
To quantitatively evaluate the fitting error and select the appropriate Copula function, Akaike information criterion (AIC), Bayesian information criterion (BIC) (Pho et al., 2019),

Copula function name
Mathematical description September 2020 | Volume 8 | Article 574101 and root mean square error (Dodangeh et al., 2017) were employed as the criteria for selecting the Copula function clusters.
K was the number of estimated parameters in the model including the intercept and ι( θ |y) was a log-likelihood at its maximum point of the estimated model; n was a sample size. The rule of selection was that the smaller the value of AIC was, the better the model was, and so did the BIC.
where n is the number of observations, X c is the theoretical probability of copula, and X o is the empirical probabilities of observations.

Correlation Analysis and Establishment of Marginal Distribution Function
In order to determine whether there was a correlation between the NDVI and monthly mean temperature and monthly precipitation, this paper established a joint distribution function and adopted Kendall, Pearson, and Spearman rank correlation coefficients to analyze the correlation between the climate elements and the NDVI. Multivariate Copula Analysis Toolbox, as a general software toolbox, uses Markov Chain Monte Carlo simulations to estimate Copula parameters (Sadegh et al., 2017). It was adopted to study the dependency structure between variables and select the optimal marginal distribution function for each variable. Distribution functions include 1) Beta, 2) Birnbaum-Saunders, 3) exponential, 4) extreme value, 5) Gamma, 6) generalized extreme value, 7) generalized Pareto, 8) inverse Gaussian, 9) logistic, 10) log-logistic, 11) lognormal, 12) Nakagami, 13) normal, 14) Rayleigh, 15) Rician, 16) T location scale, and 17) Weibull distributions. The parameters are estimated by the maximum likelihood approach. Detailed descriptions about these distributions refer to Sadegh et al. (2018).

Joint Probability Distribution
For the sake of the study of the joint probability of Tavg-NDVI and Pre-NDVI, the marginal distributions of Tavg, Pre, and NDVI were calculated respectively, and the parameters of this function were obtained. The comparison between a fitting result of the function and the actual data was evaluated by Quantile-Quantile plot. Based on the univariate marginal function, a twodimensional Copula function was constructed, and three goodness-of-fit evaluation indexes of AIC, BIC, and root mean square error were used to select the optimal Copula function type from the Copula function clusters (Chen and Sun, 2015).

The Return Period of Normalized Difference Vegetation Index and Mean Temperature/Precipitation
The return period refers to the time when the value of the random variable appears in a longer period (Singh and Zhang, 2018).
Calculating the return periods of the NDVI under different temperature and precipitation conditions can provide valuable information for a more meticulous study of how heat and water affect the NDVI. This paper calculates a bivariate joint return period and conditional return period because the univariate recurrence interval or return period often leads to overestimation or underestimation of the risk rate of an event (Shiau, 2006). It is defined that in the joint return period, X ≥ x and Y ≥ y, and in the conditional return period, X ≥ x or Y ≥ y.
In the above formulas, T joint represented the joint return period for X ≥ x and Y ≥ y; T conditional indicated the conditional return period for X ≥ x or Y ≥ y; and E(L) showed an expected value of the time interval at which continuous events start. Detailed discussions on the relationships between univariate, bivariate, and conditional return periods could be found in Shiau (2006).

RESULTS ANALYSIS Trend Analysis of Mean Temperature, Precipitation, and Normalized Difference Vegetation Index
From Guangdong to Liaoning, it spanned tropical, subtropical, and temperate zones, and climate and NDVI change rates in different regions vary markedly. Figure 2 showed the variation trends of mean temperature (Tavg), precipitation (Pre), and NDVI in Guangdong, Jiangsu, and Liaoning. It can be seen from Figure 2A that T avg in Guangdong, Jiangsu, and Liaoning provinces showed an obvious rising trend from 1982 to 2016 and displayed a larger regional difference in a descending order as follows: GuangDong > JiangSu > Liaoning. However, the mean temperature change rate was the smallest in Guangzhou and the largest in Jiangsu, which may be caused by the fact that climate regulation of sea surface temperature enhanced the stability of a climate system of Guangdong greatly influenced by marine factors. Located at the junction of subtropical monsoon climate zone and temperate continental monsoon climate zone, serving as the junction of China's zero isotherm, Jiangsu was affected by alternating cold and warm air masses. Factors capable of changing the Jiangsu's climate were complex and changeable, resulting in a large temperature change rate. In Figure 2B, the descending order of precipitations in the three regions was the same as the mean temperature, but the precipitation change rates were as follows: GuangDong > JiangSu > Liaoning. Here, the precipitation in Liaoning Province showed a significant downward trend, which was due to the law that the annual precipitation in Liaoning Province decreases from southeast to northwest at almost equal intervals, and the significant decrease of the Frontiers in Earth Science | www.frontiersin.org September 2020 | Volume 8 | Article 574101 5 precipitation mainly arose from dramatically reduced summer precipitation in southeast and winter precipitation in southern Liaoning. Figure 2C showed the variation trends of the NDVI in the three provinces. The NDVI of the three regions showed a highly consistent increase trend during the growing season. The NDVI change rate in Guangdong was notably higher than that in the other two provinces, and the NDVI of each of the three regions had a huge transition around 1998.

Joint Probability Distribution of Mean Temperature, Precipitation, and Normalized Difference Vegetation Index
The mean temperature and precipitation provide necessary heat and water conditions for vegetation growth. There are regional differences in the influence of internal hydrothermal conditions on NDVI-based growth in different temperature zones. Correlation  analysis was conducted on T avg , Pre, and NDVI during the growing season in Guangdong, Jiangsu, and Liaoning. Kendall rank correlation coefficients, Pearson correlation coefficients, and Spearman rank correlation coefficients were used to measure the correlation of two-dimensional variables. The calculation results are shown in Table 2.
According to Table 2, the T avg and the NDVI had a relatively obvious positive correlation during the growing season in the three regions, and all had passed the hypothesis test (p < 0.01). The correlation intensity was Liaoning > Jiangsu > Guangdong. There was a comparatively pronounced positive correlation between the Pre and the NDVI in Jiangsu and Liaoning, and the correlation in Liaoning was greater than that in Jiangsu. The monthly Pre and NDVI in Guangdong showed a relatively remarkable negative correlation, and both had passed the hypothesis test (p < 0.05). The correlation between the mean temperature and the precipitation decreased from north to south with the decrease of latitude, while the absolute value of a correlation coefficient corresponding to the precipitation was less than the value of the mean temperature in the same region, which indicated that the influence of the mean temperature on vegetation growth was notably higher than that of the precipitation in various temperature zones along the east coast of China.
As it could be seen from Tables 3 and 4, the optimal Copula functions for the Tavg-NDVI Copulas in Guangdong, Jiangsu, and Liaoning were BB1, Tawn, and BB1, respectively, and their corresponding three evaluation index values were all less than those of the other seven Copula functions. Similarly, the optimal Copula functions for the Pre-NDVI Copulas in the three provinces were Gaussian, Tawn, and Gaussian, respectively. It indicated that these Copula functions had the best fitting effect and were more suitable to describe the joint distribution characteristics of Tavg-NDVI and Pre-NDVI. Therefore, this paper chose BB1 Copula, Tawn Copula, and BB1 Copula as well as Gaussian Copula, Tawn Copula, Gaussian Copula functions to establish two-dimensional joint probability distribution models of the mean temperature, precipitation, and NDVI in the three regions, respectively. Figure 3 showed the joint probability distribution relationship between the monthly mean temperatures and the NDVI in the growing season from 1982 to 2016 in Guangdong ( Figure 3A), Jiangsu ( Figure 3B) and Liaoning ( Figure 3C). The joint Frontiers in Earth Science | www.frontiersin.org September 2020 | Volume 8 | Article 574101 8 probability of the mean temperature and the NDVI at any point could be obtained from the figure. The joint probability distribution function could clearly reflect the correlation between the mean temperatures and the NDVI in different ranges. The most prominent feature in the figure was an asymmetric and inclined dependency structure of monthly scale data. Guangdong probability isoline had the smallest inclination, and Liaoning's probability isoline had the largest inclination, indicating that the correlation between the mean temperature and the NDVI increased from south to north. The limit change of the tail intervals of the commonly used twodimensional Copula distribution was very small, and the tail interval mainly presented a characteristic of approximately independent distribution or relatively intensive asymptotic correlation. The probability distribution function of the mean temperature and the NDVI in the figure showed the characteristics that the upper tail was higher than the lower tail. That is, the mean temperature was correlated with the upper tail on the NDVI sequence, and the lower tail was gradually independent, demonstrating that the minimum mean temperature had little influence on the NDVI, while the maximum mean temperature had significant influence on the NDVI. According to the isogram, the interval distribution of the mean temperature and the NDVI could be obtained when the joint probability is 0.1-0.9.

Climate Response of Normalized Difference Vegetation Index to Mean Temperature and Precipitation
From the isogram, it could be seen that the greater the NDVI was when the mean temperature took a fixed value or the greater the mean temperature was when the NDVI was fixed, the greater the joint probability was. The joint probabilities of the mean temperatures and the NDVI in different intervals on the same isoline were also obviously different. The joint scenario occurrence probability of the mean temperature and the NDVI in the middle of the isoline was greater than that at either end of the isoline. It could be seen from Figure 3 that the joint probability of three scenarios of the minimum temperature and the minimum NDVI, the minimum temperature and the maximum NDVI, and the maximum temperature and the minimum NDVI was lower. However, the joint probability of the maximum temperature and the maximum NDVI was higher. When the monthly mean temperature in Guangdong was lower than 23°C, Jiangsu was lower than 17°C, and Liaoning is lower than 7°C, the probability of a lower NDVI was higher, which is not conducive to plant growth. However, when the mean temperature was between 27 and 30°C, 23 and 27°C, and 20 and 24°C, the probability of a larger NDVI (greater than 0.6) was higher and was greater than 0.7, which indicated that the vegetation in the three regions grew well. Figure 4 showed isograms of the joint exceedance probabilities of mean temperature-NDVI in Guangdong, Jiangsu, and Liaoning. According to the joint exceedance probability graph, the joint exceedance probabilities when the mean temperature and NDVI were arbitrary values could be obtained. The figure showed different combinations of the mean temperatures and the NDVI which were both greater than or equal to a specific value when the joint exceedance probabilities are 0.1-0.9, respectively. The smaller the values of the mean temperature and the NDVI were, the greater the joint exceedance probability was and vice versa. This indicated that the occurrence probability that both the mean temperature and the NDVI exceeded a smaller value was greater than the occurrence probability that both exceeded a larger value. Figure 5 showed two-dimensional isolines for joint return periods of mean temperatures and NDVI in Guangdong, Jiangsu, and Liaoning. By using the joint return period diagram, the return period corresponding to the mean temperature or NDVI greater than or equal to a specific value could be calculated. Figure 5 reflected the combinations of mean temperatures or NDVI greater than or equal to a specific value when the return periods are 2, 5, 10, 25, and 50 years, respectively. The shorter the joint return period was, the larger the value range of the mean temperature and the NDVI was, which indicated that when the mean temperature was fixed, the probability of the NDVI having a shorter return period was higher. When the mean temperature was in a higher range like 25-30°C (Guangdong), 22-26°C (Jiangsu), and 17-22°C (Liaoning), the probability of the NDVI having the shorter return period is maximum. The values of the NDVI in the above temperature ranges were all greater than 0.5, stating that the occurrence frequency of excellent plant growth was the highest. The difference in heat conditions among the three regions led to a gradual decrease of the suitable growth temperature of plants from south to north. Figure 6 showed the co-occurrence return periods of the mean temperatures and NDVI in Guangdong, Jiangsu, and Liaoning. From the co-occurrence return period diagram, return periods where the mean temperature and NDVI were both greater than or equal to a specific value are known. Figure 6 demonstrated the combinations of the mean temperatures and NDVI which were both greater than or equal to a specific value when the return periods were 2, 5, 10, 25, and 50 years. For the case that the values of the mean temperature and the NDVI were fixed, their cooccurrence return period was obviously longer than the joint return period. For the co-occurrence return period, the longer the given return period was, the greater the values of the mean temperature and the NDVI were. Even when the values of the mean temperature and the NDVI were large enough, the co-occurrence return period exceeds 50 years. Accordingly, under the background of climate warming, the increase of the mean temperature in the historical period spurred the growth of vegetation in the three regions. However, with the further increase of the temperature in the future, it may not necessarily promote the growth of plants, or even the temperature exceeded a certain threshold, the probability of inhibiting the growth of plants was higher. Figure 7 showed a dependency between monthly precipitations and NDVI in the growing season from 1982 to 2016 in Guangdong ( Figure 7A), Jiangsu ( Figure 7B), and Liaoning ( Figure 7C), reflecting values of the precipitations and the NDVI when the joint probability was 0.1-0.9. The most prominent feature in the figure was an asymmetric and inclined dependency structure of monthly scale data. Guangdong's probability isoline had the smallest inclination, and Liaoning's probability isoline had the largest inclination, which reflected that the correlation between the precipitation and the NDVI gradually increases from south to north. The greater the NDVI was when the precipitation took a fixed value or the greater the precipitation was when the NDVI was fixed, the greater the joint probability was. The probability distribution function of the precipitation and the NDVI showed the characteristics that the upper tail was higher than the lower tail. That is, the precipitation was correlated with the upper tail on the NDVI sequence, and the lower tail was gradually independent, demonstrating that the minimum precipitation had little influence on the NDVI, while the maximum precipitation had significant influence on the NDVI. When the monthly precipitation was higher than a threshold: 200 mm (Guangdong), 110 mm (Jiangsu), and 60 mm (Liaoning), NDVI values in the three regions were all greater than 0.5, and each of the joint probabilities of the two was greater than 0.5, indicating that these precipitation intervals had a greater impact on the NDVI. When the monthly precipitations in Guangdong, Jiangsu, and Liaoning were in the ranges of 250-700, 160-400, and 110-350 mm, respectively, the probability of a higher NDVI (>0.6) in each region was more than 0.7. When the precipitations in the three regions were located in the above-mentioned intervals, the precipitations enhanced the plant growth, wherein on the probability isolines at 0.7, when the precipitations in the three regions were located in the intervals of 250-550, 150-170, and 110-280 mm, the occurrence number of well vegetation growth (NDVI> 0.6) was larger than that in other intervals. Figure 8 showed isograms of the joint exceedance probabilities of precipitation-NDVI in Guangdong, Jiangsu, and Liaoning. It showed a variety of combinations of the precipitations and NDVI which were both greater than or equal to a specific value when the joint exceedance probabilities are 0.1-0.9, respectively. The smaller the values of the precipitation and the NDVI were, the greater the joint exceedance probability was and vice versa. This indicated that the occurrence probability that both the precipitation and the NDVI exceeded a smaller value was greater than the occurrence probability that both exceeded a larger value. Figure 9 showed two-dimensional isolines for joint return periods of precipitations and NDVI in Guangdong, Jiangsu, and Liaoning. Figure 9 reflected the combinations of the precipitations or NDVI greater than or equal to a specific value when the return periods were 2, 5, 10, 25 and, 50 years, respectively. The shorter the joint return period was, the larger the value range of the precipitation and the NDVI was, which indicated that when the precipitation was fixed, the probability of the NDVI having a shorter return period was higher. Figure 10 showed the co-occurrence return periods of the precipitations and NDVI in Guangdong, Jiangsu, and Liaoning. Figure 10 demonstrated the combinations of the precipitations and NDVI which are both greater than or equal to a specific value when the return periods were 2, 5, 10, 25, and 50 years. Under the same combination, the co-occurrence return period of the precipitation and the NDVI was obviously longer than the joint return period. For the co-occurrence return period, the longer the given return period was, the greater the values of the precipitation and the NDVI were. Even when the values of the precipitation and the NDVI were large enough, the cooccurrence return period exceeded 50 years. Therefore, if the precipitation continued to increase, it may not necessarily promote the growth of plants. When the precipitation exceeded a certain threshold, it was more likely to inhibit the growth of the plants.

Dynamic Change of Normalized Difference Vegetation Index
The monthly NDVI in Guangdong, Jiangsu, and Liaoning showed an overall upward trend during the growing season from 1982 to 2016. This result was consistent with the research results of the NDVI change trends on different spatial scales such as global , Eurasia (Piao et al., 2011), China (Peng et al., 2011), and eastern China (Jiang et al., 2014). The change rates of the NDVI in the three regions showed a law of increasing from north to south. However, the NDVI experienced a huge jump around 1998, which may be related to the intensified warming and strong El Nino around 1998. Climate warming leads to an increase in the daily temperature range (Ma et al., 2019). Photosynthesis of plants increased with the rise of temperature in daytime, and respiration (consumption) of the plants decreased with the reduction of temperature at night, resulting in a higher net photosynthetic accumulation of the plants, which was conducive to the growth and development of vegetation and vice versa (Wen et al., 2018). Although the vegetation in Liaoning was generally on the rise, the rising rate was lower, which could be related to the gradual warming and drying of Liaoning. Guangdong and Jiangsu were relatively rich in water resources, and the evapotranspiration caused by temperature rise would take away excess water, which in turn promoted the growth of vegetation.

Response of Normalized Difference Vegetation Index to Climate Change
Compared with the correlation study of NDVI with mean temperature and precipitation on the annual scale, the response of the NDVI to the temperature and precipitation on the monthly scale could reveal the influence of hydrothermal changes on the NDVI more deeply (Pei et al., 2019;Guo et al., 2020). On the one hand, the annual-scale dependency between the NDVI and climate factors reflects the long-term change trend relationship between the two, which failed to "peel off" other factors, such as urbanization, land utilization/change, and solar radiation in the course of the year, which affect vegetation. On the other hand, the monthly scale time series removed seasonal signals, and to some extent, eliminated the interference of signals having seasonal variation characteristics like the temperature, precipitation, solar radiation and agricultural production, etc., on the research results. It was more scientific to explore the response mechanism of the NDVIs to the climate factors (Van Gelder et al., 2008). Therefore, this paper discussed the response characteristics of monthly scale NDVI to the temperature and precipitation.

Response of Normalized Difference Vegetation Index to Mean Temperature
The difference in latitudes among Guangdong, Jiangsu, and Liaoning caused unbalanced heat. The suitable growth temperature of plants in the three regions decreases from south to north. When the mean temperature was in an appropriate range, the sufficient heat and accelerated physiological activities of plants were bringing the plants on nicely, making the NDVI higher. When the mean temperature  was in a lower range, insufficient heat would slow down the physiological activities of the plants, resulting in slower growth even stop and wither of the plants. When the temperature was in a higher range, evapotranspiration of trees was intensified by the increased temperature, which would lead to the excessive consumption of soil water and eventually affect the normal growth of vegetation (Kong et al., 2017;Zhang et al., 2019). It could be seen from Table 2 that there was a relatively obvious positive correlation between the mean temperatures and NDVI in the three regions; and the correlation intensity was Liaoning > Jiangsu > Guangdong; while the heat condition in the three regions was Liaoning < Jiangsu < Guangdong; and the demand of plant growth for heat: LiaoNing > JiangSu > Guangdong. The effect of the mean temperature on the plant growth in Liaoning was greater than that in either of Jiangsu and Guangdong, showing that the positive correlation between the mean temperature and the NDVI in Liaoning was higher than that in either of the other two regions.

Response of Normalized Difference Vegetation Index to Precipitation
Water conditions in Guangdong, Jiangsu, and Liaoning were also significantly different. The correlation coefficients in Table 2 reveal that monthly precipitation and NDVI in Jiangsu and Liaoning have obvious positive correlation, and the correlation in Liaoning was greater than that in Jiangsu; while there was an obvious negative correlation between the precipitation and the NDVI in Guangdong. Guangdong and Jiangsu were sufficient in water and higher in soil water content, and thus, generally avoids shortage of water during the plant growing season. Even if the precipitation was less, the inhibition effect on plant growth was weak, which was manifested in that the smaller precipitation and the NDVI were gradually independent of each other, and have a lower correlation. When the precipitation was excessive, the plant roots carried out oxygen-free respiration, which directly adversely affected the absorption efficiency of nutrients and water of vegetation, and inhibited the growth and development of vegetation. Meanwhile, the increased precipitation reduced light, bringing an adverse influence on the photosynthesis of vegetation, and stunting the growth of vegetation (Ye et al., 2016).
Here, a negative correlation existed between the precipitations and the NDVI. Liaoning was relatively short in water resources and was lower in the soil water content, so that the precipitation increase was beneficial to the growth of vegetation.
Climate warming was asymmetric in diurnal variations, which had a significant impact on vegetation growth (Peng et al., 2013;Du et al., 2019). The physiological effects of day-and-night temperature on plants were different, and the responses of the plants to round-the-clock asymmetric temperature rise would also be different. Therefore, it is necessary to further study the effect of asymmetric circadian warming on vegetation and the impact of continuous precipitation events on the NDVI. In addition, the precipitation events also include precipitation intensity, precipitation days, precipitation frequency, etc. In order to have a deeper and clearer understanding of the response mechanism of the NDVI to precipitation changes, indepth research on the NDVI and more precipitation events is needed.

CONCLUSION
This paper uses the Copula function theory to discuss the response characteristics of the NDVI in typical climate zones along the east coastal areas of China to monthly mean temperature and monthly precipitation. The main conclusions are as follows.
(1) The mean temperatures in Guangdong, Jiangsu, and Liaoning all showed a rising trend, with a rising rate of Jiangsu > Liaoning > Guangdong. The precipitations in Guangdong and Jiangsu showed an increasing trend, with an increasing rate of Guangzhou > Jiangsu; and the precipitation in Liaoning showed a significant downward trend. The NDVI of the three regions were all increased significantly in an ascending order of Guangdong > Jiangsu > Liaoning.
(2) The joint probability distribution functions of the monthly mean temperatures and the NDVI as well as the monthly mean precipitations and the NDVI in Guangdong, Jiangsu, and Liaoning all showed the characteristics that the upper tail was higher than the lower tail. That is, the temperature and precipitation were correlated with the upper tail of the NDVI, and the lower tail was gradually independent, demonstrating that the minimum temperature and precipitation had a little influence on the NDVI, while the maximum temperature and precipitation had a significant influence on the NDVI. (3) The shorter the return period was, the wider the ranges of the climate factor and the NDVI were, showing that when the climate factor was constant, the probability of the NDVI having a shorter return period was higher. The greater the climate factor was, the longer the return period was, indicating that the probability of plant growth inhibition was higher when the climate factor exceeded a certain threshold. (4) The hydrothermal combination conditions required for plant growth in different regions were significantly different. Plant Frontiers in Earth Science | www.frontiersin.org September 2020 | Volume 8 | Article 574101 15 growth had an optimal temperature interval and an optimal precipitation interval, and would be inhibited above or below these thresholds. (5) The subsystem composed of vegetation and climate elements had the characteristics of complexity and nonlinearity. The Copula function is utilized to quantitatively explain the relationship between different hydrothermal conditions and the NDVI through the probabilities, so that its uncertainty was studied, and the return period was calculated.
Vegetation was affected by both climate change and human activities. De-seasonal changes may eliminate signals with seasonal characteristics such as temperature, precipitation, and regular agricultural production to some extent, but the impact of the human activities on the NDVI had not been separated and quantified. In addition, the low spatial resolution of NDVI data may also affect the response degree of climate elements. Further research will be devoted to the construction of high-resolution vegetation index data sets, and the attribution of regional vegetation cover changes would be further explored, concentrating on quantitative research of influencing factors of the vegetation cover changes. It was worth noting that the rate of global warming was uneven on the spatial scale and asymmetric in time. This asymmetrical diurnal warming trend would have a major impact on the growth of global vegetation. Therefore, further studies are necessary to determine the effects of asymmetric day-and-night warming on natural ecosystems.

DATA AVAILABILITY STATEMENT
The meteorological datasets analyzed for this study can be found in the China Meteorological Administration (https://data.cma. cn/site/index.html). The GIMMS NDVI 3g.v1 (the latest updated version of the third generation Global Inventory Monitoring and Modeling System) datasets for this study can be found at https:// ecocast.arc.nasa.gov/data/pub/gimms/3g.v1/ as nc4 files. The MODIS (Moderate Resolution Imaging Spectroradiometer) NDVI dataset for this study can be found at https://ladsweb. modaps.eosdis.nasa.gov/.