Understanding Intra-Annual Dynamics of Ecosystem Services Using Satellite Image Time Series

Landscape processes fluctuate over time, influencing the intra-annual dynamics of ecosystem services. However, current ecosystem service assessments generally do not account for such changes. This study argues that information on the dynamics of ecosystem services is essential for understanding and monitoring the impact of land management. We studied two regulating ecosystem services (i. erosion prevention, ii. regulation of water flows) and two provisioning services (iii. provision of forage, iv. biomass for essential oil production) in thicket vegetation and agricultural fields in the Baviaanskloof, South Africa. Using models based on Sentinel-2 data, calibrated with field measurements, we estimated the monthly supply of ecosystem services and assessed their intra-annual variability within vegetation cover types. We illustrated how the dynamic supply of ecosystem services related to temporal variations in their demand. We also found large spatial variability of the ecosystem service supply within a single vegetation cover type. In contrast to thicket vegetation, agricultural land showed larger temporal and spatial variability in the ecosystem service supply due to the effect of more intensive management. Knowledge of intra-annual dynamics is essential to jointly assess the temporal variation of supply and demand throughout the year to evaluate if the provision of ecosystem services occurs when most needed.


Introduction
Agricultural land use, including arable lands, permanent meadows, and pastures, represents 37% of the ice-free global land surface [1]. These landscapes provide and simultaneously depend on ecosystem services that are essential for human wellbeing [2,3]. However, land degradation affects 40% of the agricultural land on earth, reducing the provision of ecosystem services and resulting in adverse environmental, social, and economic consequences [4][5][6]. Given the increased pressure on ecosystem services, numerous projects and research groups are working to improve the understanding and management of these services [7]. The ecosystem services concept is included in local, national, and international policies and their supply is a prerequisite to meet the Sustainable Development Goals [8]. During the last decades, and especially since the publication of the Millennium Ecosystem Assessment [9], there has been a noticeable increase in the use of the ecosystem service approach as an explicit decision and policy-making tool [10][11][12][13].
The ecosystem service supply is not static but depends on the dynamic structures and functions of ecosystems [14]. Weather fluctuations throughout the year affect the biophysical conditions that determine the intra-annual supply of ecosystem services, especially those that depend directly on green biomass production, such as the provision of forage [15] and crop productivity [16,17]. Soil occasionally fall below 0 °C [58].
In the study area, we focused on the subtropical arid thicket areas, where spekboom (Portulacaria afra) is a dominant and heavily grazed species [59]. We also included pastures and fields planted with rosemary for essential oil production. In this area, thicket has been heavily degraded due to unsustainable pastoralism [60,61]. Since spekboom is a succulent species that propagates vegetatively [62], spontaneous recovery does not occur in heavily degraded sites [63,64]. Land degradation has resulted in severe and widespread soil erosion [58]. The reduction of the natural vegetation, which is the common source of food for the extensively farmed goat and sheep in the area, has also contributed to a dramatic decline in agricultural returns in recent years. Also, the degradation of succulent thicket affects water infiltration by decreasing the proportion of the soil surface covered with plant litter [63]. For more than a decade, the planting of spekboom cuttings has been implemented as a practical method of restoration in the area [61,[65][66][67]. Several farmers in the study area are transitioning from extensive goat and sheep farming to more sustainable farming activities, such as essential oil production and tourism. This transition is made in partnership with three local and international non-governmental organizations that are active in the study area, i.e., (1) Living Lands (https://livinglands.co.za), (2) Grounded (https://www.grounded.co.za), and (3) Commonland (https://www.commonland.com), which are local and international non-governmental organizations. These organizations support large-scale and long-term restoration and sustainable land use initiatives. Essential oil production is considered a more sustainable farming practice in the area as it requires limited water and fertilizer inputs and needs less land to be profitable, compared to goat farming. In the study area, we focused on the subtropical arid thicket areas, where spekboom (Portulacaria afra) is a dominant and heavily grazed species [59]. We also included pastures and fields planted with rosemary for essential oil production. In this area, thicket has been heavily degraded due to unsustainable pastoralism [60,61]. Since spekboom is a succulent species that propagates vegetatively [62], spontaneous recovery does not occur in heavily degraded sites [63,64]. Land degradation has resulted in severe and widespread soil erosion [58]. The reduction of the natural vegetation, which is the common source of food for the extensively farmed goat and sheep in the area, has also contributed to a dramatic decline in agricultural returns in recent years. Also, the degradation of succulent thicket affects water infiltration by decreasing the proportion of the soil surface covered with plant litter [63]. For more than a decade, the planting of spekboom cuttings has been implemented as a practical method of restoration in the area [61,[65][66][67]. Several farmers in the study area are transitioning from extensive goat and sheep farming to more sustainable farming activities, such as essential oil production and tourism. This transition is made in partnership with three local and international non-governmental organizations that are active in the study area, i.e., (1) Living

Workflow
We selected four ecosystem services (Table 1) based on recent management objectives that aim to overcome local environmental challenges related to land degradation. We considered two regulating services, i.e., erosion prevention and regulation of water flows, and two provisioning services, i.e., the provision of forage and biomass for essential oil production. Thicket vegetation was evaluated for the supply of erosion prevention, regulation of water flows, and provision of forage. Even with the current transition from extensive livestock farming towards essential oil farming, forage is still a valuable ecosystem service in the area, not only for small livestock but also for wildlife species. We excluded irrigated agricultural land for the estimation of the regulation of water flows in order to prevent disruptions in the calculated infiltration levels due to heterogeneous initial soil moisture.
The general workflow of this study is presented in Figure 2. The first stage consisted of (1) collecting field, RS, and soil and terrain data; (2) calibrating Sentinel-2 models; and (3) estimating the supply of ecosystem services. This first stage is summarized in Section 2.2. and described in a previous study aiming to compare restoration interventions [68]. The second stage involved a temporal variability assessment in relation to the demand side of the studied ecosystem services.    To estimate ecosystem services indicators and build the models, we measured canopy dimensions, canopy cover, and infiltration in 30 plots distributed within the study area ( Table 2). Plots were sampled in pairs, ensuring a similar slope angle, orientation, and geographical vicinity to avoid wide variations in soil and weather conditions. To validate the models, we used five-fold cross-validation, which we repeated 100 times [69]. Given the low number of plots, separating the 30 samples into training and validation samples would not be a valid approach, thus justifying the use of k-fold cross-validation. We combined the field-measured fractional vegetation cover of various forms of vegetation to calculate the stratified vegetation cover and quantify erosion prevention [41]. We measured the soil infiltration rates under different vegetation cover to assess the potential regulation of water flows, as infiltration levels in successfully revegetated landscapes are expected to increase over time due to higher soil macropores [70]. We used a previously developed allometric equation to estimate the biomass of rosemary based on the measured canopy dimensions [68], and to estimate the green biomass based on the measured vegetation cover for grasses and shrubs [71]. We used the calculated fresh plant biomass as an indicator for the production of essential oil derivatives and the green biomass to quantify the production of forage. Using the field-based ecosystem service estimations, we calibrated and tested the ability of the models to estimate ecosystem services based on 13 indices derived from Sentinel-2 level 2A images and a digital elevation model (see Section 2.3). The second stage consisted of calculating the supply of ecosystem services for the whole study area using the models in Table 2 and estimating the variation of their supply for the years 2017 and 2018. To quantify how much the ecosystem service supply varies over time at the pixel level, we calculated per pixel the temporal standard deviation of each ecosystem service, separately for the years 2017 and 2018. We used a vegetation map as input to assess the supply of ecosystem services for all thicket Remote Sens. 2020, 12, 710 6 of 19 vegetation in the area [72]. In some cases, built-up, agricultural areas, or vegetation along the river were misclassified as thicket on the vegetation map. To improve the vegetation map, the thicket cover class was manually corrected using vector files of rivers, roads, and agricultural lands provided by Living Lands (Table 3). We used Google Earth Pro to correct misclassified thicket vegetation in built-up areas.
We applied the ecosystem service models under the assumption that the RS models calibrated with the field data of 30 sampled plots can be spatially and temporally extrapolated to areas with the same vegetation type, similar land management, landscape, and weather characteristics. Using the collected data of rosemary fields from 2017 and 2018, we tested how well the model calibrated for 2018 was able to estimate the ecosystem service supply in 2017. The mean fresh aboveground biomass (fAGB) measured in the field in 2017 was 196.1 g m −2 . Our RS-model predictions underestimated the provision of fAGB by 29 g m −2 on average (15%), which is within the reasonable limits for temporal extrapolation of the equations.
Our RS-based method also allows identification of the differences in ecosystem service provision between and within land cover classes. To explore the intra-annual variation within a single land cover type, we grouped each vegetation cover into five clusters. The clustering was calculated using as input the time series of estimated ecosystem service indicators derived from the 60 cloud-free Sentinel-2 images acquired in 2017 and 2018 for the study area. The clustering was performed separately for each ecosystem service, and based on how similar pixels were in their temporal behavior of the RS-derived ecosystem service values. This was achieved using ISODATA unsupervised classification [73]. We used a maximum of 50 iterations and a convergence threshold of 1. For each ecosystem service, we derived five clusters. This number was selected somewhat arbitrarily but was generally guided by our knowledge of the area to distinguish the main location clusters with differential temporal behavior for that ecosystem service. The resulting clusters were then used to spatially illustrate (1) the annual variation measured as the annual standard deviation of the ecosystem service indicators at the pixel level, and (2) describe the monthly changes of these clusters and relate them with the demand side of the assessed ecosystem services. The latter is described in the following section.
To obtain one mid-month ecosystem service indicator value, for each cluster and indicator, we interpolated the mean supply using the Akima's univariate interpolation and smoothing method [74], which uses a non-linear algorithm consisting of a set of third-degree polynomials applied to successive intervals of given points. The interpolated values were calculated using the 'akimaInterp' function of the R Package 'pracma' for regular or irregular gridded input data [75]. We plotted the monthly supply of ecosystem services to illustrate the intra-annual behavior of each ecosystem service between the selected trend classes during the years 2017 and 2018. To compare the variation of different ecosystem services and to understand the intra-annual relative variability of the resulting clusters, we calculated the coefficient of determination for the years 2017 and 2018.

Intra-Annual Dynamics of the Supply and Demand of Ecosystem Services
In this study, the demand of ecosystem services is defined as the moment at which ES are benefitted from in a particular area over a given time period [76]. By estimating the monthly demand for ecosystem services throughout the year, we analyzed the relationships between the timing of ecosystem services' peak availability and the period with the highest demand. The data source for the estimation of the ecosystem service demand included biophysical data, a literature review, and expert consultation (Table 3).
To illustrate when water erosion prevention is mostly needed, we averaged the daily rain records of six stations located across the study area and extracted the monthly maximum and the cumulative rainfall (mm) [56]. Pastures and rosemary fields are located in the lower parts of the valley and are prone to heavy wind erosion and soil deterioration when the finer particles are blown away [77]. We considered the monthly maximum wind speed (m s −1 ) as a source quantitative proxy of the main agent of wind erosion [78]. For the demand of the regulation of water flows, we used the cumulative and maximum monthly rainfall since higher infiltration rates are needed in months of higher rainfall. Regarding the demand for the provision of forage, we considered the average dry matter intake (DMI) in kg per day for an individual kudu (Tragelaphus strepsiceros) and Angora goat (Capra hircus aegagrus), representing one of the key wildlife species and small livestock species farmed in the study area, respectively. We did not consider the quality and palatability of the forage. In the case of the kudu, we used the seasonal DMI in proportion to body mass [79] and assumed a body weight of 250 kg for males and 200 kg for females. Across species of ungulates, the last gestating trimester and lactation are the biological stages when daily energy costs are the highest for females [80][81][82][83]. This energy demanding period occurs from December to February in the study area [84]. For our estimations, we assumed that calving occurred in January for both years. We considered increases of DMI for a gestating female of between 50% and 70% compared to a non-pregnant female between the last trimester of gestation and two months post-partum [82,85]. For goats, we assumed a bodyweight of 40 kg for a male and 35 kg for a female. We considered a constant DMI of 1.1 for a male goat and 1 kg/day for a female, assuming a diet containing 9 MJ/kg DM of metabolizable energy and 12% crude protein, gaining 25 g/day of non-fiber tissue, and producing 15 g/day of clean mohair fiber [86]. For a gestating goat, our estimations considered a peak DMI of 2.5 kg/day at two months after parturition and an increase DMI of 85% of the maximum value during the first month of lactation [87].
To illustrate the demand of materials for essential oil production accurately for the local conditions and management, the production manager from the Baviaanskloof Development Company for essential oil production (Devco) provided expert knowledge to understand the production goals of their rosemary fields and how they planned to manage their harvest timing throughout the year. Table 3 provides an overview of the data used in the two stages described in Section 2.2. For the first stage, we calculated spectral indices using a Sentinel-2 image from 24/06/2017 acquired over tile 34HGH, corresponding to the middle of the fieldwork period (May to July 2017). For the estimation of the biomass for essential oil production of rosemary, we collected additional field data in September-October 2018 to (1) build allometric equations to calculate the fresh biomass, and (2) to calibrate the models due to the lack of repetitions in 2017. To build the model for fresh rosemary biomass, we followed the same procedure described above using a Sentinel-2 image from 7/10/2018. In addition to the RS data, we extracted the slope (degrees), altitude (m), and aspect (north, east, south, west) from a 12.5-m resolution ALOS PALSAR-derived DEM [88].

Data Description
As input for the intra-annual variability assessment, we selected 60 Sentinel-2 images that constituted all available cloud-free images over the entire study area during the years 2017 and 2018. We used 24 images for the year 2017 and 36 for 2018. The ESA Sen2Cor processor, available in the Sentinel Application Platform (SNAP) version 6.0.2, was used for the atmospheric and topographic correction of the Sentinel-2 top-of-atmosphere level 1C images [89], i.e., to generate level 2A products. All the bands from Sentinel-2 were resampled to 10 m before calculating the spectral indices. Finally, for the illustration of the demand, we collected biophysical information from weather stations, literature reviews, and expert consultation.  Figure 3 shows the annual standard deviation per pixel of (a) erosion prevention, (b) regulation of water flows, and c) provision of forage calculated by the RS-based models for 36 moments in 2018. Higher standard deviations (red colors in the maps) indicate larger variability in the ecosystem service supply within 2018. The insets in the plots show the five clusters for thicket vegetation. Clusters' numbers were ordered according to their level of provision of ecosystem services over time, with cluster 1 representing locations with low supply, and cluster 5 high supply. Heavily degraded areas with low or absent vegetation cover (cluster 1) presented a smaller variation in the ecosystem service supply in thicket vegetation throughout the year. On the other hand, areas having relatively higher erosion prevention (Figure 3a, cluster 5) also showed larger annual fluctuations in the ecosystem service supply. The results for 2017 are presented in the Supplementary Materials ( Figures S3 and S4).

RS to Capture Ecosystem Service Supply over Space and Time
The distribution of the variation of the rosemary fields during 2018 (Figure 4a, b, c) as well as their respective clusters (Figure 4d-f). The clusters denote a large heterogeneity of biomass supply areas with similar temporal behavior within the agricultural fields. Field-specific causes exist for the within-field heterogeneity, such as concentrated runoff (Figure 4a) and the presence of different combinations of weeds and cover crops (Figure 4b).      Figure 5 shows the monthly supply and demand for ecosystem services. The ecosystem service supply is shown per cluster in each graph. The right-side axes indicate the proxies used to describe the variability in demand. These graphs do not indicate whether the total supply and demand for ecosystem services in the area match, because the units between the two y-axes are different. Nonetheless, they do indicate if a moment of high demand coincides with a relatively higher level of provision. The coefficient of variation of each ecosystem service and vegetation cover type is described in Table 4 to compare the degree of annual variability between ecosystem service supplies, years, and demands.

Intra-Annual Dynamics of Supply and the Demand of Ecosystem Services
Comparing thicket vegetation (Figure 5a,d,e) to agricultural fields (Figure 5b,c,f-h), intra-annual patterns of the ecosystem services supply in thicket vegetation showed a larger similarity between clusters. Agricultural fields also showed drastic changes in ecosystem service supply throughout the year, whereas gradual variations were present in thicket vegetation. Particularly low intra-annual variability of supply and similar clusters patterns were observed in the regulation of water flows. Clusters with overall larger provision of ecosystem services (clusters 4 and 5) generally presented more intra-annual variability of the ecosystem service supply than clusters with low ecosystem service supply. The timing in both the distribution and timing of peaks of the ecosystem service supply appeared to differ between the years 2017 and 2018. Clusters with the lowest provision of ecosystem services represented the larger proportion of the studied area. For example, thicket vegetation clusters 1 and 2 together accounted for 47% of the total area for erosion prevention, 73% for the regulation of water flows, and 66% of the provision of forage.
For erosion prevention in thicket, a mismatch in the peak ecosystem service supply and demand moments was observed during months of heavy rainfall (high demand) following dry periods of low supply of erosion prevention (e.g., January and September 2018, Figure 5a), especially for clusters with low vegetation cover (cluster 1,2). The prevention of wind erosion in agricultural fields has a relatively constant demand, so the temporal mismatches are characterized only by the supply of ecosystem services. Based on our assumptions and estimations, the periods of the high demand for forage by gestating kudu and Angora goat availability (spring and summer) did not match the peaks of supply in the thicket vegetation. With regard to the demand for fresh biomass for essential oil production, the local producers aim to harvest four metric tons per hectare (400 g m −2 ) per year. They currently harvest once a year, cutting the branches higher than 20-30 cm from the ground, which, depending on the plant size, represents between 50% and 70% of the total fresh biomass. The producers are evaluating if harvesting two or three times a year could improve productivity. The variation between the ecosystem service supply clusters indicates that the best harvest moment varies across areas.

RS to Capture Ecosystem Service Supply over Space and Time
The used Sentinel-2 images time series showed spatially diverse intra-annual changes in vegetation that affect the temporal supply of ecosystem services. This RS approach provided insights into where ecosystem services are supplied and how they vary throughout the year. We identified clusters with specific intra-annual ecosystem service behavior within a single vegetation type. These clusters not only reflected different levels of supply of ecosystem services in one moment but also denoted how (un)stable the supply of ecosystem services are in particular areas. The presented intraannual assessment provides a better understanding of the spatial and temporal distribution of ecosystem services, helping to improve monitoring of agroecosystems for more precise and timely land management according to specific contexts. The commonly used static and aggregated approach

RS to Capture Ecosystem Service Supply over Space and Time
The used Sentinel-2 images time series showed spatially diverse intra-annual changes in vegetation that affect the temporal supply of ecosystem services. This RS approach provided insights into where ecosystem services are supplied and how they vary throughout the year. We identified clusters with specific intra-annual ecosystem service behavior within a single vegetation type. These clusters not only reflected different levels of supply of ecosystem services in one moment but also denoted how (un)stable the supply of ecosystem services are in particular areas. The presented intra-annual assessment provides a better understanding of the spatial and temporal distribution of ecosystem services, helping to improve monitoring of agroecosystems for more precise and timely land management according to specific contexts. The commonly used static and aggregated approach that only considers the averages of ecosystem services per land cover per year misses this valuable information. For example, forage provision in the thicket in our study area ranged between 14.9 and 45.4 kg m −2 compared to an aggregated yearly average of 27.8 kg m −2 .
In contrast to the resilience goals where little ecosystem service supply variability is regarded as a characteristic of more resilient ecosystems [91], in this study, we found that a large proportion of the assessed thicket was severely degraded, showing low supplies of ecosystem services and a low absolute intra-annual variation. To better understand this variation in ecosystem services, it is essential to consider that parts of the most heavily degraded areas have lost most of their soil. Even under optimal weather conditions, little vegetation can grow in areas where the soil is lost, and the ground is constituted mainly by parental rock and stones. Most of the observed intra-annual changes are produced by sporadic herbaceous growth and, secondly, by the regrowth of shrubs leaves. Therefore, the degree of variation of a particular cluster of erosion prevention and forage availability in thicket vegetation is an indicator of the presence of herbaceous vegetation and/or shrub regrowth. On the other hand, areas with a higher provision of ecosystem services presented a highly unstable provision of ecosystem services. These variations are mainly caused due to the growth of plants with annual phenological cycles, affecting the general vegetation composition throughout the year, and consequently, the retrieved RS index signal [92].
Unlike a previous study carried out in Spain, where higher infiltration rates were recorded in summer than in autumn due to the higher initial soil moisture content [19], based on our model we did not find evident seasonal changes in the supply of water flow regulation services. This could indicate that our simple model may not accurately represent variability in infiltration, given that our RS-based models were built using only field measurements in a dry period. Regardless of the estimated low intra-annual variation in the supply of water flow regulation services within thicket classes, the changes were even smaller in drier (2017) than in wetter (2018) years. In the model used to estimate the infiltration rate, we considered vegetation cover as the main factor affecting infiltration rates. However, the possible occurrence of soil crusting can drastically decrease infiltration rates [93]. We indirectly considered the effect of this soil surface crust by including the slope as an additional variable, where flat areas showed a higher tendency to develop an impermeable layer. However, the model for predicting the regulation of water flows can only explain 60% of the variability in infiltration rates. We also observed soils under greener vegetation do not necessarily have higher infiltration rates, limiting our current estimation of infiltration through Sentinel-2 indices.
The assessed agricultural lands showed substantial differences between their clusters due to diverse management practices and local conditions. The creation of clusters in agricultural fields, supported with field knowledge, allowed us to identify the proportion of crops, cover crop, and weeds during the assessed period. Knowing the composition of each cluster per field is crucial for the correct interpretation of results and could help to overcome the challenge of accurately estimating fresh biomass in heterogeneous and intercropped fields with Sentinel-2 images. Therefore, the creation of clusters could improve the estimation of ecosystem services that relate to particular vegetation types, such as the provision of forage by cover crops and biomass for essential oil production by rosemary plants.

Describing Intra-Annual Dynamics of Supply and the Demand of Ecosystem Services
This study illustrates the importance of considering the time of the year and the intra-annual variability of the supply of ecosystem services in relation to their dynamic demand. When the demands for ecosystem services vary, special attention is needed to identify possible mismatches of demand and supply. Land management and agricultural practices can help to minimize this temporal mismatch. For example, fodder storage can work as a buffer for periods of low supply to: (i) Directly complement the livestock requirements of forage during the gestating period; (ii) allow for the recovery of vegetation cover; and (iii) improve soil erosion control and regulation of water flows.
Several studies have shown how sustainable land management and ecological restoration can be oriented towards the promotion of ecosystem service supply in degraded ecosystems [55,[94][95][96][97][98][99][100]. Considering the increasingly erratic global [101][102][103] and local weather behavior [104], additional efforts are needed to ascertain a sufficient supply of ecosystem services during periods of high demand. Grouping areas based on different temporal behavior in ecosystem service supply into clusters can assist in locating areas for specific management actions. Even though RS can be used to detect the location, size, and status of beneficiaries of ecosystem services [105][106][107], we focused the demand analysis on when ecosystem services are needed throughout the year. Regarding the spatial dimension, we assumed that ecosystem services are needed everywhere within the analyzed vegetation type. For example, vegetation cover is required within thicket vegetation to prevent soil erosion and regulate water flows. Still, weather conditions, such as rainfall and wind, would determine when vegetation cover is more critical. In the case of provisioning ecosystem services (forage and biomass for essential oil), a higher spatial resolution would be needed to capture animal movement or harvest activities in intercropped fields. Further attention is required on the feedback between demand and supply, including the effect of seasonal factors on the supply and demand of ecosystem services, and between different types of demand [34]. However, since RS is a physical-based approach for recording object and feature characteristics, it is generally more suitable for estimating the ecosystem service supply than their demand [24].
In this study, we assumed that our RS-models that were calibrated for 30 plots in 2017 could be spatially and temporally scaled-up to areas with similar vegetation characteristics. Potential error is introduced by these assumptions, in addition to the ecosystem service variation not captured by our RS models. In addition, we emphasize that this study did not intend to determine if particular levels of the supply of ecosystem services are sufficient for their demand at a given moment. We can only identify the temporal mismatch between the two. Moreover, because different units are used to assess the supply and demand of ecosystem services, absolute values of supply and demand cannot be directly compared. Additional fieldwork and data would be required to quantify the minimum supply for erosion prevention, regulation of water flows, fodder availability, or biomass for production that is required to satisfy specific demands. Finally, it is important to keep in mind that demand for the provision of forage in the study was estimated using several assumptions (e.g., the constant protein content of vegetation for the provision of forage, arbitrary animal species, and their body weight, calving months, diet composition). Therefore, even when our estimations provide an idea of forage demand variability, these values most likely differ from reality. For more accurate estimates, it would be necessary to include management details at the field or sub-field level and animal density information at the farm level.

Conclusions
This study aimed to illustrate the use and relevance of satellite time series to capture intra-annual variation in ecosystem services. Sentinel-2 satellite time series data can capture vegetation dynamics and, as such, can be used to assess the spatial and temporal variability of the ecosystem service supply. The consideration of the intra-annual dynamics of supply and demand provides a more realistic overview of the state of ecosystem services. It allows us to identify across locations the periods when the ecosystem service supply shows a mismatch with the demand. In addition, clustering locations based on temporal trajectories can help to capture heterogeneity within one land cover class. Understanding and accounting for this spatial and temporal variability can improve ecosystem service estimates compared to static land cover-based assessments. This study is a first step in RS-based approaches to assess the intra-annual dynamics of ecosystem service supply and demand. There are still several challenges to solve in the future related to improving the estimates of ecosystem service supply, especially when the ecosystem service is related to a specific species within a complex vegetation composition. Agricultural lands have both large intra-annual variability and large spatial differences in ecosystem service supply, highlighting the need for information on intra-annual and spatial variability in ecosystem service assessment. Quantitative assessments of the intra-annual dynamics of ecosystem service supplies related to their demand can support more effective monitoring and timely ecosystem-based management of agricultural landscapes. For their wellbeing, people need sufficient provision of ecosystem services and temporally reliable levels of ecosystem service supply that match their demand. Funding: None to be declared.