Grassland Phenology Response to Climate Conditions in Biobio, Chile from 2001 to 2020

: Plant phenology is affected by climate conditions and therefore provides a sensitive indicator to changes in climate. Studying the evolution and change in plant phenology aids in a better understanding of and predicting changes in ecosystems. Vegetation Indices (VIs) have been recognized for their utility in indicating vegetation activity. Understanding climatic variables and their relationship to VI support the knowledge base of how ecosystems are changing under a new climatic scenario. This study evaluates grassland growth phenology in the Biobio, Chile, biweekly with Moderate Resolution Imaging Spectroradiometer (MODIS) Normalized Difference Vegetation Index (NDVI) time series. Four growth parameters for the six agro-climatic regions were analyzed from 2001 to 2020: start and end of the season, time and value of maximum NDVI. For this purpose, the NDVI time series were smoothed using Savitzky–Golay ﬁltering. In addition, by using monthly gridded database climate data, we studied correlations between phenology markers and rainfall, maximum temperature and minimum temperature. The results show that both the start and end of the growing season did not signiﬁcantly change; however, all agro-climatic regions grow faster and more vigorously. Thus, climatic conditions in Biobio have become more conducive to grassland growth over the 2001–2020 period.


Introduction
Vegetation phenology is the study of biological patterns in plant growth over time, such as germination, flowering, fruiting, lead emergence, etc. Moreover, it examines how these patterns relate to environmental factors, such as rainfall and air temperature, thus providing information about vegetation productivity, carbon reserve and carbon dynamics [1], which, in turn, yields data on the responses and adaptations of vegetation to climate change [2]. Therefore, monitoring vegetation with phenological metrics provides data to track changes in vegetation linked to events such as drought, fire, climate fluctuations or directional climate change [3]. The impact of climate change on vegetation causes variations in the phenological patterns of vegetation such as temporal displacement of phenological cycles, changes in plant morphology, colonization by other species, or even extinction of other species. Thus, these processes generate important indicators in ecosystem functions and species composition [4,5].
Grasslands are widely distributed around the world and play an important role in carbon storage [6], which is a crucial factor in the mitigation of climate change [7]. Grasslands also provide essential resources for animal life and humankind as well as rics extraction methods can be divided into two categories, threshold-based and change detection methods.
The objective of this study was to track phenological changes in Biobio grasslands. The study was carried out with the specific aims: (i) Are there significant differences in the timing of phenological markers? (ii) What climatic variables contribute to them? (iii) Are there any trends in these markers in the 2001-2020 period?

Study Area
The Biobio Region (Chile) (central coordinates 37 • 15 S, 72 • 30 W, WGS-84) covers an area of 37,068 km 2 ( Figure 1a) and is divided in six agro-climatic regions: Secano Costero, Secano Interior, Depresión Intermedia, Cordón Isla, Precordillera and Cordillera (Figure 1b). This region is in a transition area ranging from a warm Mediterranean climate to a humid and temperate climate which makes it sensitive to ecological changes, such as human land use and the effects of climate change. Its bioclimatic variables per Hijmans et al. [55] and through monthly rainfall and temperature values show a total annual range of precipitation from 750 mm to 2000 mm, with the precipitation in the driest month being lower than 35 mm while in wettest month ranging from 150 mm to 350 mm. The temperature in the warmest month is around 25 • C while the coldest month ranges from 0 • C to 7 • C, with the annual mean temperature equal to 10 • C. Most of the region is covered by forest, followed by croplands and grasslands.

methods.
The objective of this study was to track phenological changes in Biobio grasslands. The study was carried out with the specific aims: (i) Are there significant differences in the timing of phenological markers? (ii) What climatic variables contribute to them? (iii) Are there any trends in these markers in the 2001-2020 period?

Study Area
The Biobio Region (Chile) (central coordinates 37° 15′ S, 72° 30′ W, WGS-84) covers an area of 37,068 km 2 ( Figure 1a) and is divided in six agro-climatic regions: Secano Costero, Secano Interior, Depresión Intermedia, Cordón Isla, Precordillera and Cordillera (Figure 1b). This region is in a transition area ranging from a warm Mediterranean climate to a humid and temperate climate which makes it sensitive to ecological changes, such as human land use and the effects of climate change. Its bioclimatic variables per Hijmans et al. [55] and through monthly rainfall and temperature values show a total annual range of precipitation from 750 mm to 2000 mm, with the precipitation in the driest month being lower than 35 mm while in wettest month ranging from 150 mm to 350 mm. The temperature in the warmest month is around 25 °C while the coldest month ranges from 0 °C to 7 °C, with the annual mean temperature equal to 10 °C. Most of the region is covered by forest, followed by croplands and grasslands.

Evolution of Climatic Variables
This section explains the evolution of the accumulated rainfall and maximum and minimum temperature variables for the six agro-climatic regions of the study area over the period of 2001-2020. The per-year graph densities are shown in Figure 2, which includes the agro-climatic region and climatic variable, the differences between the monthly mean values and their corresponding monthly mean value over the whole time series. In addition, in Appendix A, Figures A1-A3, represented in the boxplot, graphs the temporal evolution of these variables on a monthly scale for the period studied. The time series of accumulated rainfall, maximum temperature and minimum temperature was decomposed into three components: trend, seasonal and irregular components through an additive model. Each agro-climatic region in Figure 3 shows the trend of the climatic variables analyzed. All of the regions show a downward trend related to accumulated rainfall (Figure 3a). This progressive decrease follows sinusoidal behavior, The variations in accumulated rainfall are represented in Figure 2a. When the density graph shifts to the right of 0 in any given year, it means that that year had greater rainfall within the period analyzed. On the other hand, when shifted to the left, that year was dry. In addition, the density graphs in Figure 2a present a binomial distribution, implying an irregular distribution of rainfall between very rainy and very dry periods and presenting similar behaviors across the agro-climatic regions in overall accumulated rainfall. In general terms, between 2001 and 2010, the annual density graphs show a distribution shifted to the right, and therefore accumulated rainfall was higher than the mean of the total period. On the other hand, from 2010 onwards, the curves shifted to the left and, therefore, it was drier. Thus, the accumulated rainfall has been decreasing over the years. There are, however, years with extreme accumulated rainfall values, both due to excess and deficiency. Thus In regard to temperature, density graphs shifted to the right of 0 mean warmer years, while those shifted to the left mean colder years. The density graphs with variations of the maximum temperatures with respect to their monthly mean values ( Figure 2b) show a tendency to shift to the right with respect to the central value of 0, which indicates an increase in maximum temperatures at the end of the temporal series. The year 2016, in particular, stands out for its increase in maximum temperatures to levels much higher than those in other years within the study period. In addition, it can be seen that the peaks of the density graphs are left at 0 at the beginning of the time series, which means that monthly maximum temperatures were lower than the mean values for 2001-2020. As the time series progresses, the crests of the density graphs move to the right, thus showing increasing maximum temperatures. As with accumulated rainfall, this trend appears across all the agro-climatic regions. Finally, the minimum temperatures, shown in Figure 2c, display the same behavior as the maximum temperatures in such a way that the minimum temperatures at the beginning of the time series were lower than the mean monthly values of the time series analyzed and gradually shift to the right indicating an increase.
The time series of accumulated rainfall, maximum temperature and minimum temperature was decomposed into three components: trend, seasonal and irregular components through an additive model. Each agro-climatic region in Figure 3 shows the trend of the climatic variables analyzed. All of the regions show a downward trend related to accumulated rainfall (Figure 3a). This progressive decrease follows sinusoidal behavior, with very wet years interspersed with drier ones. Moreover, this trend is evident from the beginning of all the agro-climatic series. The Cordón Isla region (Figure 3 Figure 4 summarizes the methodology applied to this research. MODIS MOD13Q1 normalized difference vegetation index (NDVI) [57] was used as the remote sensing dataset, while TERRACLIMATE [58] was used as the climate database. The extraction of remote sensing and climate data for the study area between 2001 and 2020 was carried out with Google Earth Engine (GEE) using Earth Engine Python API as the client server to connect with Earth Engine Services. These extracted data, stored in Google Drive, were downloaded for further processing on a local machine. Information on the areas occupied by grasslands as well as the delimitation of the agro-climatic regions was obtained from the Corporación Nacional Forestal (CONAF) [59]. Grassland areas for each region were used to determine accumulated rainfall, maximum and minimum temperature, and NDVI trends. These data were then analyzed to determine the statistical relationships between climatic variables and NDVI. The phenological markers of the grassland were determined according to each agro-climatic region. an average value equal to 5.35 • C. The standard deviation of maximum temperatures for all regions is equal to ±0.5 • C. Based on the resulting trends obtained for accumulated rainfall, maximum temperatures and minimum temperatures in the 2001-2020 period, the six agro-climatic regions in Biobío are undergoing a reduction in water and an increase in temperature. Figure 4 summarizes the methodology applied to this research. MODIS MOD13Q1 normalized difference vegetation index (NDVI) [57] was used as the remote sensing dataset, while TERRACLIMATE [58] was used as the climate database. The extraction of remote sensing and climate data for the study area between 2001 and 2020 was carried out with Google Earth Engine (GEE) using Earth Engine Python API as the client server to connect with Earth Engine Services. These extracted data, stored in Google Drive, were downloaded for further processing on a local machine. Information on the areas occupied by grasslands as well as the delimitation of the agro-climatic regions was obtained from the Corporación Nacional Forestal (CONAF) [59]. Grassland areas for each region were used to determine accumulated rainfall, maximum and minimum temperature, and NDVI trends. These data were then analyzed to determine the statistical relationships between climatic variables and NDVI. The phenological markers of the grassland were determined according to each agro-climatic region.  MODIS MOD13Q1 [57] data were used to study spatiotemporal changes in the grasslands. This dataset was computed from atmospherically corrected surface reflectance scenes with a temporal resolution of 16 days and a spatial resolution of 250 m. These were masked for water, clouds, cloud-shadow and heavy aerosols. Monthly climate data of maximum and minimum temperatures as well as rainfall were extracted from the TERRACLIMATE database.

Datasets and Image Processing
A total of 480 MOD13Q1 NDVI images, as well as 240 images for each climatic variable, maximum and minimum temperature, and rainfall, from TERRACLIMATE from 2001 to 2021 were extracted. In addition, as TERRACLIMATE data have a monthly reso- MODIS MOD13Q1 [57] data were used to study spatiotemporal changes in the grasslands. This dataset was computed from atmospherically corrected surface reflectance scenes with a temporal resolution of 16 days and a spatial resolution of 250 m. These were masked for water, clouds, cloud-shadow and heavy aerosols. Monthly climate data of maximum and minimum temperatures as well as rainfall were extracted from the TERRACLIMATE database.
A total of 480 MOD13Q1 NDVI images, as well as 240 images for each climatic variable, maximum and minimum temperature, and rainfall, from TERRACLIMATE from 2001 to 2021 were extracted. In addition, as TERRACLIMATE data have a monthly resolution, MOD13Q1 NDVI data were processed using the Maximum Value Composites (MVC) method to acquire monthly NDVI. We used the MVC method to avoid the influence of solar altitude angle, clouds and/or atmospheric effects [60], thus permitting the study of the relationships between the two types of variables.

NDVI Time-Series
The phenological markers SOS, EOS, the peak of maximum NDVI and the grassland values from 2001 to 2020 were extracted from the MOD13Q1 dataset. Firstly, referencing land coverage and land use (LCLU) maps from CONAF [59], pure pixels occupied by grasslands were identified. In order to perform this, the polygons categorized as grasslands were extracted and intersected with the MODIS pixel grid. On this spatial division, only those pixels completely occupied by grassland were selected for each agro-climatic region, thereby removing mixed pixels. For each MOD13Q1 scene and agro-climatic region, the median NDVI was determined, obtaining six NDVI time-series, one for each agroclimatic region, with a temporal resolution of 16 days. R-commander was used for this purpose.
The six NDVI time-series were processed with TIMESAT software v3.3 [53] where annual NDVI values were smoothed using the Savitzky-Golay filter [61]. The phenological markers were determined as the 20% of the season amplitude from the left and right minimum values of NDVI fitted values, respectively [62,63]. Based on previous research [64], grassland phenology for each agro-climatic region was summarized as they are relatively homogeneous both in climate and biological diversity.

Statistical Analysis
First, to obtain an overview of the dynamics of the climate variables, the mean value of these variables for the period of 2001-2020 was determined for each month. Subsequently, for each month and year, their difference with the monthly mean value for the whole period was determined. These differences were represented by annual density plots to describe the behavior of the variables over the period analyzed.
Secondly, climatic variables and NDVI time series were decomposed into trend, cycle and residual. Trend corresponds to a long-term process that occurs over time, and the cycle explains the cyclical process that operates for each cycle once the trend is accounted for and residual, obtained after accounting for trend and cycle, explains local process caused by variability between cycles. Next, via cross-correlation analysis, the relationships between NDVI values and climatic variables for each agro-climatic region were examined in order to evaluate the climatic dependency of the Biobio grasslands.
Finally, as in other research [65,66], temporal trends in the data sets were evaluated with a linear regression model in which time was the independent variable while phenological markers were the dependent variables.

NDVI Time-Series
The evolution of the monthly mean NDVI for the six regions studied is shown in Figure 5. Each region shows a different development curve. Thus, in regions such as the Cordillera (Figure 5f), there is little variation in NDVI throughout the year, while in Cordón Isla (Figure 5d Seasonal parameters, SOS and EOS, extracted from TIMESAT for each agro-climatic region are shown in Table 1, where the term season refers to the annual cycle of NDVI. In addition, Figure A4 in Appendix A represents the NDVI time series and phenological markers SOS and EOS for each agro-climatic region. Based on these data, in the Biobío region, the average SOS between 2001 and 2020 is on the day of the year (DOY) 133, while EOS is on 388, resulting in a mean duration of the grassland cycle equal to 255. However, each agro-climatic region presented different phenological markers. The grasslands in the Secano Costero region are the ones that start the cycle earlier (SOS equals DOY 103). As they move away from the coast towards the interior, SOS increases until reaching the Precordillera region, where SOS decreases. Thus, mean SOS in Secano Interior, Depresión Intermedia, Cordón Isla and Precordillera appears on DOY 111, 113, 118 and 107, respectively. Finally, the grasslands in the Cordillera region present the greatest delayed mean SOS with a DOY equal to 245. On the other hand, EOS markers in grasslands according to the agro-climatic region do not show a pattern similar to SOS. The Cordillera region presents a mean EOS in DOY 506, i.e., the grassland cycle is between two years, while Cordon Isla presents in DOY 349. On the other hand, Secano Costero, Secano Interior, Depresión Intermedia and Precordillera present a mean SOS in DOYs equal to 383, 361, 364, 349 and 368, respectively. With these phenological markers, the average grassland cycle in the Biobio regions ranges from 279 days in Secano Costero to 231 days in Cordon Isla. In the Secano Interior and Depresión Intermedia regions, the average grassland cycle is 250 days, while in the Precordillera and Cordillera, it is 260 days. Cordon Isla (Figure 5d) is the region with the greatest differences in NDVI between months. The lowest NDVI value is March (0.383), and the highest is September (0.641). From March to May, the NDVI increases progressively, the latter being the one with the largest interquartile range (0.08). Until September, the NDVI remains stable at high values and then starts to decrease. For the Precordillera region (Figure 5e), the NDVI has its lowest values in February (0.56). From February onwards, the index increases progressively until November (0.76), where it starts to decrease. Finally, the Cordillera region (Figure 5f) is the one with the most stable NDVI value annually, with respective maximum and minimum values of 0.67 in November and 0.58 in May, registering the highest interquartile range in the season of growth for grasslands (December-March).
Seasonal parameters, SOS and EOS, extracted from TIMESAT for each agro-climatic region are shown in Table 1, where the term season refers to the annual cycle of NDVI. In addition, Figure A4 in Appendix A represents the NDVI time series and phenological markers SOS and EOS for each agro-climatic region. Based on these data, in the Biobío region, the average SOS between 2001 and 2020 is on the day of the year (DOY) 133, while EOS is on 388, resulting in a mean duration of the grassland cycle equal to 255. However, each agro-climatic region presented different phenological markers. The grasslands in the Secano Costero region are the ones that start the cycle earlier (SOS equals DOY 103). As they move away from the coast towards the interior, SOS increases until reaching the Precordillera region, where SOS decreases. Thus, mean SOS in Secano Interior, Depresión Intermedia, Cordón Isla and Precordillera appears on DOY 111, 113, 118 and 107, respectively. Finally, the grasslands in the Cordillera region present the greatest delayed mean SOS with a DOY equal to 245. On the other hand, EOS markers in grasslands according to the agro-climatic region do not show a pattern similar to SOS. The Cordillera region presents a mean EOS in DOY 506, i.e., the grassland cycle is between two years, while Cordon Isla presents in DOY 349. On the other hand, Secano Costero, Secano Interior, Depresión Intermedia and Precordillera present a mean SOS in DOYs equal to 383, 361, 364, 349 and 368, respectively. With these phenological markers, the average grassland cycle in the Biobio regions ranges from 279 days in Secano Costero to 231 days in Cordon Isla. In the Secano Interior and Depresión Intermedia regions, the average grassland cycle is 250 days, while in the Precordillera and Cordillera, it is 260 days.

Analysis of Relations between Climatic Variables and NDVI
Previous research projects have shown that the time lag between NDVI to climatic factors ranges from 0 to 3 months [67][68][69]. Each agro-climatic region was analyzed for the correlation coefficient values between concurrent average monthly NDVI and concurrent average monthly climatic variables (lag 0), the monthly average of the previous month (lag 1) and monthly averages of the two previous months (lag 2) ( Table 2). Positive correlation results between accumulated rainfall and NDVI were obtained in lag 0 and 1 for all agro-climatic regions, while for lag 2, the Cordillera region showed a negative correlation, mainly due to the presence of snow. Accumulated rainfall and NDVI showed a significant positive correlation from lag 0 and, therefore, an immediate positive impact on grasslands. Particularly, in Secano Costero, Precordillera and Cordillera, the response of NDVI to accumulated rainfall showed the highest correlation at lag 0 while in Secano Interior, Depresión Intermedia and Cordón Isla was at lag 1, which means that the highest positive effect appears the following month. In Lag 2, high values are maintained above all in the areas where it rains less, Secano Interior, Depresión Intermedia and Cordón Isla. However, in the areas where rainfall is less necessary, a drop in the correlation in lag 2 can be observed. For maximum and minimum temperature, except in the Cordillera region, the highest correlations were obtained in Lag 0 for both maximum and minimum temperatures. The correlation was negative, indicating that higher temperatures have a negative effect on NDVI, reducing it. Moreover, its effect on grassland appears immediately, with lag equal to 0. The Cordon Isla region showed the highest values, corresponding to the driest region of all the analyzed regions. At the opposite extreme, temperatures in the Cordillera region, both maximum and minimum, have a positive effect on vegetation, increasing the NDVI. The result of an increase in temperature is reflected in the following two months (Lag 2). The results obtained show how climatic variables are related to NDVI values. Thus, an increase in rainfall results in an increase in NDVI values, while they decrease with increasing temperatures.
Having proved that climatic variables influence grassland development in Biobio, with a growing tendency towards reduced rainfall and increased maximum and minimum temperatures, the general trend and phenological markers of NDVI in each region between 2001 and 2020 are presented below. Trends of NDVI values for each region are shown in Figure 6. The trends of the NDVI series of the agro-climatic regions in Biobio show a similar pattern, differentiating an initial period with little change in the NDVI and a second period with an increase in the index. Thus, between 2001 and 2007, NDVI values do not show a clear change in their evolution. In 2008 there was a reduction in the NDVI value in each region, and from that year onwards, there was an increase in NDVI, with greater or lesser intensity depending on the agro-climatic region. Such NDVI trends are aligned with the evolution of the climatic variables presented in Section 3.1. Of all the regions, the Secano Costero region is the one with the highest NDVI values and the smallest increase, presenting stability compared to the other regions. On the other hand, the Precordillera region presents the highest increase in NDVI, reaching similar values to the Secano Costero region in recent years. The Secano Interior and Depresión Intermedia regions behave similarly, with shared increases in NDVI, but more moderately than in the Precordillera region. The Cordillera region has NDVI values in the same range as the two previous regions, although with more accentuated fluctuations over time.
regions, the Secano Costero region is the one with the highest NDVI values and the smallest increase, presenting stability compared to the other regions. On the other hand, the Precordillera region presents the highest increase in NDVI, reaching similar values to the Secano Costero region in recent years. The Secano Interior and Depresión Intermedia regions behave similarly, with shared increases in NDVI, but more moderately than in the Precordillera region. The Cordillera region has NDVI values in the same range as the two previous regions, although with more accentuated fluctuations over time.    In general, the trend changes in grassland phenology at the peak of maximum NDVI and its value in the period 2001-2020 are more pronounced than SOS and EOS. In the Secano Costero region, the temporal relationships between SOS (Figure 7(a.1)) and EOS (Figure 7(b.1)) slope negatively, which means both growth periods occur earlier, although the duration of the cycle stays the same. This same behavior appears in the Depresión Intermedia (Figure 7(a.3,b.3) and Precordillera (Figure 7(a.5,b.5)) regions. On the other hand, the Secano Interior and Cordillera region present with positive slopes in SOS (Figure 7(a.2,a.6)) and a negative slope in EOS (Figure 7(b.2,b.6)). This means a reduction in the phenological cycle of the grasslands in these regions due to the delay of SOS and the advancement of EOS. Finally, the Cordon Isla region has a positive slope in both SOS (Figure 7(a.4)) and EOS (Figure 7(b.4)) so that there is a delay in the entire grassland cycle in this region without major changes in the duration of the cycle. Based on the low values obtained in the slopes of the SOS and EOS trends, there appear to be no important changes in the start or end of the phenological cycle of the grasslands.
When analyzing the results obtained in the trend at the moment of reaching the maximum NDVI (Figure 7c), as well as its value (Figure 7d), there is a general trend to reach maximum NDVI earlier (negative slope) and with a higher value in the index (positive slope). It is in the regions of Secano Costero (Figure 7(c.1)) and Depresión Intermedia (Figure 7(c.3)) where there is a more evident advance in reaching maximum vigor, with a higher negative slope. Secano Interior (Figure 7(c.2)), Precordillera (Figure 7(c.5)) and Cordillera (Figure 7(c.6)) show this advance as well but in a less accentuated manner. Cordón Isla (Figure 7(c.4)) is the only region that shows a delay in this marker. Finally, all regions show an upward trend in NDVI, with Cordon Isla (Figure 7(d.4)) showing the greatest increment over time, while the Secano Costero (Figure 7(d.1)) and Cordillera (Figure 7

Discussion
Phenological observations provide localized information with high temporal resolutions in different biological phases. Phenology varies over geographic gradients according to climate zone and vegetation type. Phenology also varies within communities, and the phenology of individual plants plays a key role in determining how ecosystems are structured and how they function [70]. Temperate forest trends are attributed to warming temperatures [23,71]. In addition, boreal and subalpine forests show similar trends to those observed in temperate forests, which are attributed to a warming climate [72,73]; moreover, there is little evidence that rainfall has a significant influence on boreal forest phenology [74]. In Subalpine meadows and the Arctic and alpine tundras, in high latitudes and high altitude ecosystems, there are two main factors that regulate SOS, the timing of snowmelt and the temperatures that follow [75] while warm might delay EOS [76]; however, it has not been extensively studied. On the other hand, the impacts of climate change on tropical forests vary across regional-to-continental scales. Decreasing precipitation and reduced cloudiness in Amazonia may result in different phenological shifts than, for example, in central Africa, where precipitation is projected to increase [77]. Mediterranean areas include a diverse range of plant types, and phenological responses to environmental cues vary accordingly. Temperature is a key factor for most species, but rainfall, through its influence on soil moisture, is also important [78,79]. Finally, in subtropical desert ecosystems, phenological variations resulting from climate change may occur due to changes in the timing and quantity of rainfall and increases in temperature [80]. Therefore, phenology varies over geographic gradients according to climate zone and vegetation type.
In contrast to direct phenological observations, satellite data have the advantage of high spatial coverage. It allows global, regional, and local studies. NDVI and climatic time series provide key information in the study of vegetation. The analysis of vegetation requires a spatial and temporal approach to explain changes and dynamics in vegetation. Temporal, seasonal and phenological metrics obtained by vegetation index series from remote sensing scenes provide useful information to understand vegetation conditions at regional scales. Taking into account previous research, the climatic response of grasslands varies geographically. Previous studies have shown that a warming climate has a positive impact on greening in the northern hemisphere, relying on the continuously warming air and spatially uneven trends of precipitation [81,82]. However, increased temperatures could not be the primary control of grassland. Rainfall is a factor that plays a relevant role in grassland. In the Biobio region, temperature increase from 2001 to 2020 has extended the period where conditions are favorable to plant growth; however, the grasslands have not benefitted from these conditions. Instead, grasslands in Biobio have shown little change at SOS and EOS. However, photosynthetic capacity, measured by NDVI, has increased significantly in all the agro-climatic regions of Biobio. These trends show how grasslands grow faster and more vigorously as temperatures increase. On the other hand, the incapacity of grasslands to exploit thermal resources can be potentially explained by a mismatch between current vegetation and climate conditions, which could be explained by the novel conditions caused by climate change, such as the reduction in rainfall. These results presented here align with other research, such as those in the Tibetan Plateau [83].
Grasslands are very sensitive to rainfall, so more arid conditions would lead to lower productivity. However, a decrease in winter frosts through higher minimum temperatures, higher temperatures in general and higher solar radiation would compensate for the negative effect of lower rainfall and may even result in a slight increase in productivity [84]. This contradictory situation could also be explained by the beneficial effect of high CO 2 concentrations in the atmosphere, which stimulates net plant photosynthesis [85], improving their growth and productivity [86,87], and water use efficiency [88]. CO 2 concentration can be measured through the use of onboard sensors on satellite platforms such as Sentinel 5P, GOSAT or OCO-2. These platforms have only recently been launched, however, and there are not enough data to create robust time-series data to analyze; thus, future work should evaluate climate and atmospheric variables together. Therefore, recent warming trends have been associated with earlier vegetation activity in spring and an overall extension in the length of the growing season [89]. On the other hand, there is less of a consensus on how climate change is affecting end-of-season phenology, although many biological events are observed to be occurring later [90].
While no comparisons have been made in this study between ground observations and the remote sensing-based results, this type of activity is needed to more fully understand and validate remote sensing-based observations of large-scale phenology in terrestrial vegetation.
On the other hand, variations in climate conditions have an important impact on water balance and material circulation in watersheds. Thus, the influence of climate change on water quality has been previously documented [91,92]. It is assumed that climate change can influence water environments by affecting water flow, pollutant transformation and migration, as well as toxicity [93,94]. Moreover, transformation and migration patterns of toxic pollutants in the water environment are complex processes, and the influencing mechanism of climate change on heavy metals is still unknown. Recently, a few studies demonstrated the impact of meteorological factors on heavy metals using model simulation methods or short-term water quality data [95,96]. High heavy metal concentrations in soils were cited as one of the factors limiting vegetation establishment and growth in areas suffering high anthropogenic impact [97][98][99]. Thus, heavy metal toxicity first affects the growth of vegetation roots because they are more sensitive to heavy metal stress. In addition, heavy metal pollution limits root growth and results in decreased root weight [100,101]. Based on previous research in the study area [102][103][104][105], future projects should analyze the influence of heavy metal and contaminant concentrations in soil, together with climatic variables and vegetation indices, to study changes in the phenology of vegetation, including grasslands.
Based on article 4.8 of the United Nations Framework Convention on Climate Change [106], Chile meets 7 of the 9 vulnerability characteristics to be identified as a region in the world affected by the alteration of the current global climate pattern. The map of grassland phenology in the Biobio region can support decision-making in monitoring vulnerability to climate variations and therefore serves as a useful tool for assessing measures taken against climate change.

Conclusions
In a scenario where climatic conditions are changing, monitoring tools are required to understand impacts on vegetation better. Therefore, an understanding of the drivers involved in phenology is needed even when detailed field observations are lacking, and meteorological stations are scarce. This result makes remote sensing datasets for vegetation monitoring a particularly useful tool under these scenarios. In this study, we evaluated the timing of phenology markers and their response to climate conditions in six agro-climatic regions in Biobio. Our results suggest that climatic conditions in Biobio have become more conducive to grassland growth over the 2001-2020 period. During this period, rainfall has decreased; however, favorable thermal conditions have extended. In this context, the grasslands have shown more intensive photosynthetic activity but without extending the period of activity. The six areas analyzed show that between 2001 and 2010, the evolution of NDVI values was stable. However, between 2010 and 2020, there was an increase in NDVI values to a greater or lesser extent depending on the area analyzed. The Precordillera region shows the greatest increase, while Secano Costero is the most stable.
In general, the results from this study highlight the relationship between grassland phenology and their response to the climatic variables of accumulated rainfall and temperature, both maximum and minimum, which has implications for developing policy frameworks for grassland management and protection and its relationship with future climate change.