Cluster-Segregate-Perturb (CSP): A Model-agnostic Explainability Pipeline for Spatiotemporal Land Surface Forecasting Models

## Abstract

Abstract Satellite images have become increasingly valuable for modelling regional climate change effects. Earth surface forecasting represents one such task that integrates satellite images with meteorological data to capture the joint evolution of regional climate change effects. However, understanding the complex relationship between specific meteorological variables and land surface evolution poses a significant challenge. In light of this challenge, our paper introduces a pipeline that integrates principles from both perturbation-based explainability techniques like LIME and global marginal explainability techniques like PDP, besides addressing the constraints of using such techniques when applying them to high-dimensional spatiotemporal deep models. The proposed pipeline simplifies the undertaking of diverse investigative analyses, such as marginal sensitivity analysis, marginal correlation analysis, lag analysis, etc., on complex land surface forecasting models In this study we utilised Convolutional Long Short-Term Memory (ConvLSTM) as the surface forecasting model and did analyses on the Normalized Difference Vegetation Index (NDVI) of the surface forecasts, since meteorological variables like temperature, pressure, and precipitation significantly influence it. The study area encompasses various regions in Europe. Our analyses show that precipitation exhibits the highest sensitivity in the study area, followed by temperature and pressure. Pressure has little to no direct effect on NDVI. Additionally, interesting nonlinear correlations between meteorological variables and NDVI have been uncovered.

###### Index Terms:

## I Introduction

The Earth’s climate is undergoing a significant transformation, posing a substantial threat to human existence. Its detrimental effects on terrestrial surfaces, which sustain most life on our planet, are becoming increasingly evident. From the depletion of Arctic sea ice [] to the intensification of fire incidents [], the repercussions of climate change manifest across diverse and variable geographic regions.
Studying the effects of meteorological variables like temperature, precipitation, and pressure is central to climate change analysis.

Over the past decade, there has been a notable increase in satellite sensors, leading to the availability of Earth observation data on an unprecedented scale. Initiatives like the Copernicus program [] offer high-resolution data with enhanced temporal coverage, enabling the generation of dense predictions and analyses that were previously unattainable.
Land surface forecasting using spatiotemporal forecasting models plays a crucial role in predicting changes in surface conditions over time, such as vegetation growth [] [], soil moisture levels [] [], and land use patterns [] []. However, interpreting the outputs of these models and understanding the factors driving their predictions pose significant challenges, particularly in the context of complex spatiotemporal data and high-dimensional feature spaces.

Explainability has emerged as a critical requirement for ensuring the transparency, trustworthiness, and reliability of machine learning models, including those used in land surface forecasting. Traditional explainability techniques, such as Local Interpretable Model-agnostic Explanations (LIME) [] and partial dependence plots (PDP) [] etc, have limitations when applied to land surface forecasting models. These models often operate in high-dimensional spatiotemporal feature spaces, making it challenging to isolate the effects of individual variables on model predictions.

We introduce the Cluster-Segregate-Perturb (CSP) pipeline, a novel approach to explainability in land surface forecasting models. CSP pipelines clustering, segregation, and perturbation to enable comprehensive investigative analyses of model predictions. By clustering instances of meteorological variables, we identify unique patterns within each meteorological variable. After clustering, the data samples are segregated based on the earlier clusters. This segregation organizes the data into segments we termed, weather segments, representing smaller, more homogeneous subsets of the dataset based on distinct sets of meteorological conditions. Finally, by perturbing the meteorological variables of data samples within each segment, the pipeline creates artificial samples representing variations from the original data. These artificial samples make the analyses more robust.

In this paper, we present the development and evaluation of the CSP pipeline for investigative analyses on land surface forecasting models. We demonstrate its effectiveness through empirical evaluations in uncovering the relationships between meteorological variables and land surface evolution via NDVI.

## II Challenges in Explainability of Spatiotemporal Land Surface Forecasting Models

- •
The main challenge is dealing with the temporal nature of the meteorological variables, Unlike global explainability techniques like PDP, which rely on simple artificially generated samples, to give insights into the marginal effect of a feature on the predicted output, generating high-dimensional natural spatiotemporal data poses significant difficulties. Despite all efforts, the curse of dimensionality emerges. Even after successfully generating spatiotemporal samples that appear natural, the reliability of the pre-trained model’s output may be compromised if it hasn’t been trained on data with similar natural patterns. This is a common challenge when dealing with complex models that possess a high-dimensional latent space, the approximation of the distribution is often sparse and hard to locate in such a latent space.
To overcome these challenges, we utilized the original training data and applied perturbations to the meteorological variables. While introducing variations, these perturbations remained within the range of natural weather patterns. By doing so, we effectively generated artificial data points that retained the characteristics of the observed patterns, thereby enabling the model to make reasonable and robust predictions besides eliminating the extra processing required to create artificial samples from scratch.
- •
Another challenge is handling weather-based variations. As we know, within a season, different temperature, pressure, and precipitation patterns will affect the land surface evolution differently. Therefore, any analysis without properly segregating these natural weather patterns might yield incorrect aggregated values.
To address this limitation, we segregated the samples into weather segments, representing smaller, more homogeneous subsets of the dataset based on distinct sets of meteorological conditions. By organizing the data in this manner, we were able to aggregate the effects within each segment, providing a more comprehensive and accurate understanding of the model’s behaviour across different weather scenarios. This is similar to perturbation-based explainability techniques like LIME, which operate at the instance level. Instead, here we work at the segment level.

## III Related Work

### III-A Explainable AI

Explainable AI (XAI) techniques have been developed to address the lack of interpretability in deep learning models, including Convolutional Neural Networks (CNNs) and LSTM models. These techniques aim to provide insights into the decision-making process of these models. One approach is to use saliency maps, feature attribution, and local interpretable model-agnostic explanations (LIME) to enhance transparency and trustworthiness, another approach involves modifying the architecture of CNNs to improve interoperability, such as in [] Habib et al. uses sinc-convolution layers and explanation vectors to identify domain-specific insights and in [] authors proposed a modification of the LSTM architecture called HydroLSTM, which enhances the interpretability of the model by representing internal system processes in a manner analogous to a hydrological reservoir.

Quantifying the global relationship between input features and model predictions in time-series image forecasting presents a significant challenge and remains an active area of research and the attention given to model explainability in time-series applications has not been as significant as in the fields of computer vision or natural language processing []. Huang et al. in [] suggest two techniques for explainability in the spatiotemporal predictive learning task (SPLT): first, the synthesis of multiple independent components to analyze how the features contribute to the prediction; and second, a state decomposition and expansion technique to separate intertwined signals in the spatiotemporal dynamical system helping in exploring the mechanisms underlying motion formation. they concluded that a collaboration mechanism, namely, extending the present and erasing the past (EPEP), explains the motion formation in SPLT. Duckham et al. in [] utilises the Simple Event Model and PROV-O ontologies to enable queries not just about reasoner inferences but also about explanations for specific conclusions reached by the system. This capability is embedded in the NEXUS system, which combines multiple reasoning components that can support a wide range of spatiotemporal queries. Moosburger in [] proposed an interpretable and modular framework for unsupervised and weakly-supervised probabilistic topic modelling of time-varying data. This framework merges generative statistical models with computational geometric techniques. Pham et al. in [] proposed Temporally Weighted Spatiotemporal Explainable Neural Network for Multivariate Time Series (TSEM) merges RNN and CNN capabilities by using RNN hidden units to weigh the temporal axis of CNN feature maps. It matches STAM’s accuracy and fulfils several interpretability standards, including causality, fidelity, and spatiotemporality.

### III-B Exploring Meteorological Variables through NDVI Analysis

Several studies have examined the sensitivity, correlation, lag, etc of meteorological variables such as precipitation and temperature concerning NDVI. Li and Guo, in [], analyzed the response characteristics and sensitivity of NDVI to climatic factors in Tianjin, China. They found that NDVI increases with rising temperatures but gradually decreases with increasing precipitation. Moreover, the impact of temperature on NDVI was more pronounced in spring and autumn, while precipitation had a dominant effect in spring, autumn, and early summer. Shun et al. [], concluded that NDVI exhibits a higher correlation with air temperature in high-altitude alpine and plateau areas, whereas it correlates more strongly with precipitation in grassland and desert grassland regions.
Yujie et al. in [], concluded that precipitation was identified as the most significant factor affecting vegetation evolution, followed by temperature, land cover change, population, elevation, and nightlight.
Hao et al. [] studied the link between climatic variables and NDVI in the upper stream of the Yellow River and a strong ‘correlation between NDVI and precipitation for grassland and forest. Their results suggest that higher precipitation levels lead to elevated NDVI values. Additionally, the monthly highest temperature and precipitation significantly affected NDVI.
Similarly, Wang et al. in [] explored the spatial distribution and year-to-year changes in NDVI across the central Great Plains and found a clear link between NDVI patterns and average annual precipitation. A strong correlation was observed between NDVI and precipitation deviations during dry years, like 1989, following another dry year. Conversely, a weak correlation was observed during wet years, such as 1993, one of the wettest on record. These findings indicate higher correlation coefficients during or after dry periods.
In [], Jianming et al. . studied the time accumulation effect of meteorological variables on NDVI. They observed a positive correlation between NDVI and accumulated temperature, accumulated precipitation, and effective accumulated precipitation.

## IV Assumptions

- •
Availability of the original dataset used while training the model.
Understanding what the model learns is crucial for adjusting its architecture or addressing data-related biases. Additionally, there may be subsets of the samples that the model struggles to approximate. By solely examining training scores, shortcomings in these cases may not be apparent. The CSP pipeline can offer insights into such scenarios by organizing samples into weather segments, representing smaller, more similar subsets of the dataset based on distinct meteorological conditions.
- •
Input features are independent of each other.
In marginal feature analysis, feature independence simplifies the process by isolating the effect of one variable on the outcome while holding all others constant. This ensures that changes in one feature do not influence the distribution or behaviour of other features, leading to a clearer understanding of their relationship with the outcome.

## V CSP Pipeline

- •
Cluster
Identify the inherent and distinct spatiotemporal patterns within each input feature by clustering.
- •
Segregate
Partition the samples into segments containing smaller, more homogeneous subsets of the dataset based on distinct sets of spatiotemporal patterns of the features, By doing this, any analysis performed on these groups reflects the behaviour of a specific local environment, allowing for more focused and insightful interpretations.
- •
Perturb
Perturbation creates artificial samples in a segment. This facilitates many exploratory studies besides sensitivity analysis and makes the aggregated values more robust.

To demonstrate the application of the CSP pipeline to a land surface forecasting model, we employed a ConvLSTM model trained on the EarthNet21 dataset. It’s worth noting that the pipeline is model-agnostic, meaning it can be seamlessly integrated with any other model architecture.

### V-A Dataset, Model and Metric Overview

#### V-A 1 EarthNet2021 - Dataset

EarthNet2021 [] is a large-scale dataset and challenge for Earth surface forecasting, which involves predicting satellite imagery conditioned on future weather.
The Dataset consists of 32,000 samples within the European region, each comprising a series of 30 Sentinel-2 images, each captured at intervals of 5 days. These images contain four bands (red, green, blue, and near-infrared) with a spatial resolution of 128x128px or 2.56 km^2 and a ground resolution of 20m. Moreover, accompanying weather-related meteorological data is included, such as precipitation, sea level pressure, and temperature (minimum, maximum, and mean), each comprising a series of 150 images for 150 days, at a coarser spatial resolution of 80x80px or 102.4 km^2, sourced from the observational dataset E-OBS [].
Before the model training, as a preprocessing step, the spatiotemporal resolution of meteorological variables is matched to that of the Sentinel-2 images.

This study used the IID (In-Domain) train set fraction of the EarthNet2021 dataset denoted by $D$ containing 23904 samples denoted by $N$ in our analyses.

$$ $\begin{split}D=&\{x_{1},\ x_{2},\ x_{3},\ \ldots\ x_{N}\}\\ x_{i}=&(T_{\text{avg}_{i}},\ T_{\text{min}_{i}},\ T_{\text{max}_{i}},\ P_{i}, \ R_{i})\end{split}$ (1) $$

$x_{i}$ is a data sample $T_{\text{avg}_{i}},\ T_{\text{min}_{i}},\ T_{\text{max}_{i}},\ P_{i},\ R_{i}$ are the 30 timesteps spatiotemporal channels representing average temperature, minimum temperature, maximum temperature, pressure and precipitation respectively each $\in\mathbb{R}^{(30*128*128)}$. Additionally, other channels representing $DEM$, $red$, $blue$, $green$, $near-infrared$, $cloud\ mask$, $scene\ classification\ label$, and $data\ quality\ mask$ are also utilized during model training and inference. However, since they are not the main focus of this study, they are not explicitly notated.

#### V-A 2 ConvLSTM

ConvLSTM (Convolutional Long Short-Term Memory) networks have been proposed as effective methods for remote sensing time series analysis. These networks leverage the temporal and spatial contextual information present in time series images to improve classification accuracy. It was first used for precipitation nowcasting in [] since then it has been used for tasks such as land cover classification, change detection, and time series reconstruction.

In this study, we utilized the ConvLSTM model, as described by Diaconu et al. []. We trained the model on the IID (In-Domain) split of the EarthNet2021 dataset for 60 epochs. During this training phase, we attained an EarthNetScore of 0.3257, closely aligning with the score reported in the original paper, 0.3266. Here the EarthNetScore (ENS) is a composite evaluation metric used for assessing the performance of Earth surface prediction models it is the harmonic mean of the four components (MAD, OLS, EMD, SSIM), scaled between 0 (worst) and 1 (best) as described in [].

Let $M_{\theta}$ denote the ConvLSTM model having $\theta$ as pre-trained weights. then the prediction $y_{1}$ is equated as:

$$ $y_{i}=M_{\theta}\left(x_{i}\right)$ (2) $$

here $y_{i}=(r_{i},\ g_{i},\ b_{i},\ nir_{i})$, where $r_{i}$, $g_{i}$, $b_{i}$ and $nir_{i}$ represent the red, blue, green and near-infrared channels of the output respectively, and each channel $\in\mathbb{R}^{(20*128*128)}$.

#### V-A 3 Soft-DTW (Soft Dynamic Time Warping)

Soft-DTW [] is an extension of Dynamic Time Warping, which is a method for measuring similarity between two sequences that may vary in time or speed. DTW is particularly useful when comparing time series data where there might be temporal distortions or differences in the pacing of events.

In our investigation, soft-DTW served as the distance metric within k-means for clustering meteorological time series variables. For this purpose, we employed tslearn, a Python library specialized in time series data analysis, featuring pre-implemented soft-DTW functionalities [].
The formal definition of soft-DTW :

$$ $soft-DTW^{\gamma}(\tau,\tau^{\prime})={\min}^{\gamma}_{\pi\in\phi(\tau,\tau^{ \prime})}\sum_{(i,j)\in\pi}f(\tau_{i},\tau^{\prime}_{j})^{2}$ (3) $$

$$ $\min{}^{\gamma}(q_{1},\dots,q_{n})=-\gamma\log\sum_{i}e^{-q_{i}/\gamma}$ (4) $$

Here $\min{}^{\gamma}$ is the soft-min operator parametrized by a smoothing factor $\gamma$ this makes it differentiable everywhere.
$\phi(\tau,\tau^{\prime})$ represents the set of all possible alignments between the two input sequences $\tau$ and $\tau^{\prime}$. This set contains all the possible ways the elements of the two sequences can be aligned to each other. Each alignment, denoted by $\pi$, is a set of pairs $(i,j)$ where $i$ is an index from sequence $\tau$ and $j$ is an index from sequence $\tau^{\prime}$. These pairs represent the correspondences between elements of the two sequences in a particular alignment.
$f(\tau_{i},\tau^{\prime}_{j})$ represents the distance or dissimilarity between the $i$-th element of sequence $\tau$ and the $j$-th element of sequence $\tau^{\prime}$, in eq 4 it is denoted as $q_{i}$.
In the case of tslearn [], this distance of dissimilarity is Euclidean distance.

### V-B Applying the Stages of the CSP Pipeline

#### V-B 1 Clustering the Meteorological Variables

In the EarthNet21 dataset, each sample covers a small geographical area of only 2.56 km^2, resulting in minimal variability of meteorological variables like temperature, precipitation, and pressure across the spatial resolution at any given time step. Leveraging this characteristic, we simplified the clustering process from spatiotemporal to temporal clustering. This involved downsampling the meteorological variables to a single value for each time step, followed by time-series k-means clustering using soft-DTW on the downsampled meteorological variables.

To downsample the meteorological variables of each data sample we computed the average of the pixel values of each timestep image. The spatial resolution of each image is (128,128) px representing 2.56 km^2 ground resolution, we average it down to a single value. This approach condenses the meteorological variables to 30 values each representing one timestep. This downsampling is justified as the pixel values of these variables exhibit minimal variations within the channel images, and this even aligns with the amount of variance of a meteorological variable in a small geographical area.
Following the downsampling we conducted k-means clustering on each meteorological variable across different cluster sizes and obtained the optimal cluster size $K$ by selecting the one with the highest cluster $GoodnessScore$.

Figure: Figure 1: We’ve identified $9$ unique clusters (representing base temporal patterns) for $precipitation$ and $pressure$, and $7$ unique clusters for $temperature$. The $frequency$ displayed on each plot indicates the number of data samples assigned to each cluster. In the pressure column the $2^{nd}$, $4^{th}$ and $8^{th}$ plots exhibited noise patterns in the pressure data. Therefore, we can safely disregard these clusters.
Refer to caption: extracted/5786296/images/uniqueClusters.png

$$ $GoodnessScore_{K}=\frac{interCentroidScore_{K}}{intraClusterScore_{K}}$ (5) $$

We performed clustering for cluster sizes $\in\{2,3....15\}$. The goal is to identify a set of unique cluster centroids for each meteorological variable exhibiting a low $intraClusterScore$ and a high $interCentroidScore$.

$$ ${intraClusterScore_{K}}=\sum_{i=1}^{K}\frac{\sum_{j\in d_{i}}{e^{DTW(centroid_ {i}^{K},j)}{N_{i}}}}{N}$ (6) $$

After clustering for a specific cluster size $K$, $d_{i}$ is the subset of samples assigned to the cluster with centroid $centroid_{i}^{K}$ where $i\in\{1,2....K\}$.
In Equation 6, we calculate the fractional average exponential DTW value between the centroid and the samples assigned to that centroid where $\frac{N_{i}}{N}$ is the weight for the fractional average. $N$ refers to the total number of data samples, and $N_{i}$ refers to the count of samples belonging to cluster $i$. The function $DTW(x,y)$ is unbounded from above and its lower bound is zero. When the two temporal signals are more similar, the value of the function approaches zero so ultimately a lower value of $intraClusterScore$ is desirable.

$$ $\begin{aligned} {interCentroidScore_{K}}&=\frac{CentroidPairs(K)- SimilarCentroids(\lambda)^{2}}{CentroidPairs(K)}\\ CentroidPairs(K)&=\frac{K(K-1)}{2}\\ SimilarCentroids(\lambda)&=\sum_{i=1}^{K}\sum_{i=j+1}^{K}\mathbbm{1}\left(DTW( centroid_{i},centroid_{j})<\lambda\right)\end{aligned}$ (7) $$

For $interCentroidScore$ in the set of equations in 7 we penalise the cluster for having more similar centroids because we aim to find unique temporal patterns. Hence, a higher value of $interCentroidScore$ is desirable. Again we used the $DTW(x,y)$ function for calculating the similarity score, as discussed earlier this function is only bounded from below so we put a threshold $\lambda$ for classifying the centroids into two classes i.e., similar or sufficiently different to be deemed unique. through trial and error, we set its value to $0.4$.

Figure 1 illustrates the unique clusters identified for each meteorological variable.

After the clustering step, the prediction of the cluster-index $I_{\alpha}$ of a time series signal $Q$ is defined as

$$ $I_{\alpha}=m^{K}_{\alpha}(Q)$ (8) $$

where $m^{K}_{\alpha}$ denote the k-means model with the highest $GoodnessScore$ having $K$ unique clusters for meteorological variable $\alpha$ where $\alpha$ is a $\in$ $\{t_{avg},t_{min},t_{max},p,r\}$ shortened notation for meteorological variables.

#### V-B 2 Segregating Samples into Weather Segments

We partitioned the dataset into segments, containing smaller, more homogeneous subsets of distinct temporal meteorological patterns. These segments are recognized by a tuple $c_{i}$ having five cluster-index $I_{\alpha}$, one for each meteorological variable $\alpha$. Since a segment is a set of samples with similar weather conditions regardless of their geographical location, we refer to it as weather segment and notate this set as $S_{i}$.

$$ $c_{i}=(I_{r},I_{p},I_{t_{\text{avg}}},I_{t_{\text{min}}},I_{t_{\text{max}}})$ (9) $$

We identify the weather segment $S_{i}$ of a sample by determining the tuple $c_{i}$, by applying the k-means models in 8 to each meteorological variable of the sample.

$$ $\begin{split}S_{i}&=\{x_{j}:(m^{9}_{r}(R_{j}),\ m^{9}_{p}(P_{j}),\ m^{7}_{t_{ \text{avg}}}(T_{\text{avg}_{j}}),\\ &\ m^{7}_{t_{\text{min}}}(T_{\text{min}_{j}}),\ m^{7}_{t_{\text{max}}}(T_{ \text{max}_{j}}))=c_{i}\}\end{split}$ (10) $$

Thus the entire dataset is segregated into distinct segments, $D=S_{1}\cup S_{2}\cup S_{3}\ldots\cup S_{N_{c}}$ and through this process we uncovered a total of $N_{c}=784$ weather segments.

#### V-B 3 Perturbation

Figure: Figure 2: Illustrate average NDVI signals for different perturbations of the meteorological variables for a weather segment whose $c_{i}$ is (6,5,5,4,2), and frequency is 465 samples.
Refer to caption: extracted/5786296/images/perturbation.png

We created artificial samples by making small adjustments to one of the meteorological variables at a time while treating the others as constant across all the samples within a segment.

Since the samples originate from Central and Western Europe, primarily within the temperate zone, we have taken precautions to ensure that the variables remain within natural bounds. The average temperature of the dataset hovers around $20$ $\mathrm{\SIUnitSymbolCelsius}$. Therefore, the additive adjustment for temperature $\delta_{1}\in[-10,-5,0,+5,+10,+15,+20]$ $\mathrm{\SIUnitSymbolCelsius}$. Similarly, the average precipitation revolves around 1.5 $\mathrm{mm}$/30 timesteps and the average pressure revolves around 1020 $\mathrm{hPa}$ so the additive adjustment in precipitation and pressure are $\delta_{3}\in[-2,0,+2,+4,+6,+8,+10]$ mm/30 timesteps and $\delta_{2}\in[-30,-20,-10,0,+10,+20,+30]$ $\mathrm{hPa}$ respectively. The curves in fig 2 show a weather segment’s average perturbed NDVI signals.

Let $V^{(\delta_{1},\delta_{2},\delta_{3})}$ be the small additive adjustment done to the meteorological variables, defined as:

$$ $\begin{split}&V^{(\delta_{1},\delta_{2},\delta_{3})}=(v^{\delta_{1}}_{t_{\text {avg}}},v^{\delta_{1}}_{t_{\text{min}}},v^{\delta_{1}}_{t_{\text{max}}},v^{ \delta_{2}}_{p},v^{\delta_{3}}_{r})\\ v^{t}_{\alpha}&=A_{ijs}\mid\forall i\forall j\forall s\{A_{ijs}=t\}\text{ and} \in\mathbb{R}^{(30*128*128)}\end{split}$ (11) $$

eq: 11 simply suggests that adjustment $V^{(\delta_{1},\delta_{2},\delta_{3})}$ is a tuple similar in shape and order to a data sample defined in eq: 1. $\delta_{1},\delta_{2},\delta_{3}$ is the amount of additive adjustment for temperature, pressure and precipitation channels respectively.

Also, the meteorological channels in the EarthNet21 dataset have been normalized using the eq: 12 so we also normalized these adjustments $\delta_{1},\delta_{2},\delta_{3}$ before adding them to their respective meteorological channels.

$$ $R_{$\mathrm{mm}$}=50R,\quad P_{$\mathrm{hPa}$}=200P+900,\quad T_{$\mathrm{ \SIUnitSymbolCelsius}$}=50(2T-1)$ (12) $$

## VI Experiments

We conducted two investigative analyses: marginal sensitivity analysis and marginal correlation analysis between the meteorological variables and NDVI of the output of the ConvLSTM prediction.

We choose NDVI (Normalized Difference Vegetation Index) for this study since changes in meteorological variables like temperature, precipitation, and pressure greatly affect vegetation health and density.

$$ $NDVI\left(y_{i}\right)=\frac{nir_{i}-r_{i}}{nir_{i}+r_{i}}$ (13) $$

### VI-A Marginal Sensitivity Analysis

We conducted the marginal sensitivity analysis of each meteorological variable on the NDVI of the predicted output within each weather segment. This analysis is localized because each weather segment contains only a subset of the samples thus following a distinct pattern of these meteorological variables.

The average NDVI signal of a segment labelled with tuple $c_{i}$ having the sample set $S_{i}$ corresponding to the perturbations $(\delta_{1},\delta_{2},\delta_{3})$ in a meteorological variable is defined as:

$$ $\overline{NDVI}_{(\delta_{1},\delta_{2},\delta_{3})}^{c_{i}}=\frac{\sum_{x\in S _{i}}NDVI\left(M_{\theta}(x+V^{(\delta_{1},\delta_{2},\delta_{3})})\right)}{ \left|S_{i}\right|}$ (14) $$

The curves from the eq: 14 are visualized in fig:2.

Since we are doing marginal analysis, we set $\delta_{2}=\delta_{3}=0$ when computing the above metric 14 for temperature. Similarly for pressure $\delta_{1}=\delta_{3}=0$ and precipitation $\delta_{1}=\delta_{2}=0$. Hence local marginal sensitivity of temperature($t_{avg}$), pressure($p$) and precipitation($r$) for the weather segment labeled with tuple $c_{i}$ is given as:

$$ $\begin{split}Senstivity_{t_{avg}}^{c_{i}}=\sum_{a}\sum_{b\mid a\neq b}\frac{ \left|\overline{NDVI}_{(a,0,0)}^{c_{i}}-\overline{NDVI}_{(b,0,0)}^{c_{i}} \right|}{\left|a-b\right|}\\ Senstivity_{p}^{c_{i}}=\sum_{a}\sum_{b\mid a\neq b}\frac{\left|\overline{NDVI} _{(0,a,0)}^{c_{i}}-\overline{NDVI}_{(0,b,0)}^{c_{i}}\right|}{\left|a-b\right|} \\ Senstivity_{r}^{c_{i}}=\sum_{a}\sum_{b\mid a\neq b}\frac{\left|\overline{NDVI} _{(0,0,a)}^{c_{i}}-\overline{NDVI}_{(0,0,b)}^{c_{i}}\right|}{\left|a-b\right|} \end{split}$ (15) $$

After determining the marginal sensitivity of meteorological variables $\alpha\in\{t_{avg},p,r\}$ for individual weather segments, we found out that irrespective of the weather segment the sensitivity of the variables remained almost the same so we approximated the global sensitivity by computing the weighted average where weights are the cardinality of the sample sets of each weather segment, denoted by $\left|S_{i}\right|$.

$$ $Senstivity_{\alpha}=\sum_{i=1}^{N_{c}}(Senstivity_{\alpha}^{c_{i}}*\left|S_{i} \right|)$ (16) $$

#### Result

**TABLE I: Marginal NDVI sensitivity table**
| Variable | Sensitivity | SD | Unit |
| --- | --- | --- | --- |
| Precipitation | 0.0183 | 0.0043 | per $\mathrm{mm}$ |
| Temperature | 0.0034 | 0.0014 | per $\mathrm{\SIUnitSymbolCelsius}$ |
| Pressure | 0.0015 | 0.0007 | per $\mathrm{hPa}$ |

Table I demonstrates that within the region of study i.e., Europe, a unit change in precipitation has the most significant effect on the NDVI value, followed by temperature and pressure in sequence.
To quantify this impact, the weight of precipitation is approximately $12$ times greater than that of pressure and roughly $5$ times greater than that of temperature.

### VI-B Marginal Correlation Analysis

Figure: Figure 3: Correlation patterns between meteorological variables and NDVI of the predictions. Temperature is split into two curves to enhance visualization: lower NDVI scenes exhibit greater correlation curve curvature, which decreases as NDVI increases.
Refer to caption: extracted/5786296/images/correlation.png

We conducted a correlation analysis between meteorological variables and NDVI of the predictions. We aimed to identify a best-fitting correlation curve for each meteorological variable $\alpha\in\{t_{avg},p,r\}$ within each weather segment denoted by the tuple $c_{i}$. However, for curve fitting, we needed a set of points representing the aggregated curve for each meteorological variable $\alpha$ within a weather segment $c_{i}$:

$$ $Points_{\alpha}^{c_{i}}=\{(x_{b_{1}},y_{b_{1}}),(x_{b_{2}},y_{b_{2}}),\dots\}$ (17) $$

Here $y_{b_{i}}$ is the median value obtained by downsampling the average NDVI signal of the weather segment $c_{i}$ for the perturbation $b_{i}$ i.e., downsampling the curve obtained through the eq: 14 and $x_{b_{i}}$ is the mean value of the meteorological variable $\alpha$ of the weather segment $c_{i}$ added to the perturbation $b_{i}$. Value of the perturbation $b_{i}\in\delta_{1}$ in case of temperature, $\in\delta_{2}$ in case of pressure and $\in\delta_{3}$ in case of precipitation as given in 11 and $\forall i\forall j\{b_{i}\neq b_{j}\}$.

With the points set in hand, we conducted curve fitting using a range of linear and non-linear models, including polynomial, exponential, logarithmic, sinusoidal, and Gaussian curves.

$$ $coeff_{\alpha}^{c_{i}}=fit(eq,Points_{\alpha}^{c_{i}})$ (18) $$

Here the $fit$ function takes a curve equation and a set of points as parameters and returns the coefficients of the best-fitting curve.
It was discovered that second-degree polynomials gave the best fit for most weather segments across all meteorological variables. Figure 4 shows one instance of the fitted curves.

We plotted the curves from all the weather segments to visualize the underlying pattern in 3. We standardized the range of $x_{a_{i}}$ for different meteorological variables i.e., for temperature, the range was standardized to $[0,35]{$\mathrm{\SIUnitSymbolCelsius}$}$ similarly for pressure and precipitation it was $[990,1040]{$\mathrm{hPa}$}$ and $[-2,12]{$\mathrm{mm}$}$ respectively. Finally, we used the best-fitted curve to calculate the approximated value of $y_{a_{i}}$ for the ranges of $x_{a_{i}}$ as mentioned above.

#### Result

**TABLE II: Correlation curves**
| Curve | $a^{\prime}$ | $b^{\prime}$ | $c^{\prime}$ |
| --- | --- | --- | --- |
| Precipitation | -0.0034 | 0.0252 | 0.5554 |
| Temp-${NDVI}_{high}$ | 0.0001 | 0.0044 | 0.5805 |
| Temp-${NDVI}_{low}$ | -0.0002 | 0.0096 | 0.3750 |
| Pressure | 2.18e-05 | -0.0441 | 22.8677 |

Table II illustrates the correlation curves. The analysis suggests that the correlations exhibit nonlinear behaviour and can be described with the equation of the parabolic curve:

$$ $y=a^{\prime}x^{2}+b^{\prime}x+c^{\prime}$ $$

.
Figure 3 illustrates that increasing precipitation and temperature lead to higher NDVI values. However, with precipitation, this rise reaches a threshold somewhere around 4-5 $\mathrm{mm}$, beyond which NDVI declines. Notably, when NDVI is low, the correlation curves for both temperature and precipitation exhibit more pronounced curvature. This curvature diminishes as NDVI increases, suggesting a diminishing effect of additional precipitation on NDVI in scenes with already high vegetation indices. Additionally, the relationship between pressure and NDVI appears nearly linear, indicating the minimal impact of pressure changes on NDVI.

Figure: Figure 4: The figure highlights the curve fitting of the correlation graphs of a weather segment for different meteorological variables
Refer to caption: extracted/5786296/images/curveFitting.png

## VII Challenges

While the CSP pipeline is simple, model-agnostic and inherently generic, the approach has some challenges.

- •
Spatiotemporal clustering is still an active area of research. selecting the appropriate clustering method according to the type of data is crucial. In [] Ansari et al. categorised spatiotemporal clustering into six broad categories however the review focuses more on individual spatiotemporal points, such as events, geo-referenced data items, time series, moving objects, and trajectories unlike which most land surface prediction dataset contains sequences of satellite images and meteorological data.
For clustering entire spatiotemporal samples, the technique would need to be adjusted to emphasize the clustering of these samples as a whole, rather than focusing solely on the characteristics of individual points. Time complexity would be another crucial aspect to consider.
Another way to look at this challenge is video clustering. This approach aims to uncover inherent structures or themes within video datasets, including common topics, visual motifs, or recurrent motion patterns. While it can be effectively applied to clustering satellite images based on visual motifs, it is not well-suited for clustering meteorological features, which rely more on temporal patterns for detail.
- •
Maintaining temporal consistency and spatial coherence can be challenging when creating artificial data by perturbing high-dimensional spatiotemporal samples. The process differs vastly between perturbing satellite images, where spatial details stay pretty much the same or change at an overall constant rate over time. It is very different from dynamic media like music videos, where things change significantly. So, it’s important to pick the right perturbation method based on the dataset being worked with.

## VIII Conclusion and Future Work

This paper presented a pipeline to improve the explainability of complex land surface prediction models like ConvLSTM. The proposed methodology enables various investigative analyses, enhancing our understanding of the relationship between meteorological variables and model predictions. Our analysis revealed that NDVI exhibits the highest marginal sensitivity against precipitation, followed by temperature and pressure, with approximate ratios of 12:2:1. Moreover, we observed a nonlinear correlation between NDVI and meteorological variables, resembling a parabolic curve. Furthermore, as the average NDVI of the scene increases, the influence of precipitation and temperature on the curvature of the correlation curve diminishes. Additionally, it is concluded that pressure has little to no direct effect on NDVI.

In the future, we aim to pursue further advanced studies. This involves exploring the transitional impacts of meteorological variables by using techniques like Perlin noise [] which can ensure smoother interpolations, alongside conducting lag analysis. Nonetheless, the challenges in the clustering process need to be researched more, It is imperative to cluster samples based on factors other than meteorological variables, such as crop type, elevation, building density, etc and downsampling might not be the option every time, to address this obstacle, we intend to develop a spatiotemporal deep clustering method, thus enhancing the methodology’s adaptability to handle more diverse datasets.

## Peer Reviews

confidence: 4
contribution: 1
initial_rating: 3
presentation: 2
questions: 1.	Could you elaborate on your method’s contribution and explain how it outperforms existing methods? Have you compared CSP to other approaches?
2.	Though you claim CSP can be integrated with any other model architecture in the paper, have you tried any other models except for ConvLSTM? If so, do you have the corresponding results? 
3.	Have you evaluated CSP in scenarios where input features are correlated and do not satisfy the assumption?
soundness: 3
strengths: 1. The idea of designing a pipeline for the complex spatiotemporal data and high-dimensional feature spaces is intriguing.
2. The description of the challenges in explainability of spatiotemporal land surface forecasting is technically sound.
summary: This paper introduces the Cluster-Segregate-Perturb (CSP) pipeline as an approach for achieving explainability in land surface forecasting. The CSP pipeline includes clustering, segregation, and perturbation. This pipeline facilitates analyses such as marginal sensitivity, correlation, and lag analysis, etc for complex land forecasting models. CSP has been tested on EarthNet2021.
weaknesses: 1. The paper lacks novelty in its method design. The author may overstate the ability of CSP, for example, downsampling the images by averaging down to a single value is kind of rough. While it greatly reduces the complexity, it also loses the spatial information.
2. The paper’s presentation needs to be improved; the structure lacks coherence and unclear.
3.EARTHNET2021 is the only dataset being used. It is necessary to evaluate the method across various datasets in land surface forecasting fields. Also, consider using synthetic data to prove the functionality of CSP.
4. The experiments section mainly focuses on sensitivity and correlation analysis, lacking deeper and more detailed analysis. Additionally, there are no baselines for comparison.
5. The assumption of this method is not universally valid. Ignoring the relationships between input features may oversimplify the analysis since weather features are often correlated.

---

confidence: 4
contribution: 1
initial_rating: 3
presentation: 2
questions: 1) I wonder if you could use climate analogues instead of the Perturbations? --> What if you switch the weather of one group of minicubes with that of another to compute the sensitivities? In that way you'd be more sure it is actually "physical" climate you are considering.
2) How do you match the resolution of the meteorological variables to Sentinel 2? Do you take a central cut-out that actually reflects the Sentinel 2 scale? Or do you interpolate, as in the spatial references do not match anymore?
3) Related to this, for the spatial averaging to do clustering of the meterological variables, are those time series representing 2.5km^2 or 100km^2?
4) Could you actually analyse if the relationships are different for different land cover types? E.g. Grasses and Trees should have very different relationships to weather, for instance under drought grasses can be brown within a few days, but evergreen trees take months before you visibly see changes in the NDVI.
5) The ConvLSTM uses memory, which may reflect some of the relationship between temperature and NDVI on a monthly time scale, in other words, the model might incorporate a trend for its meteorological variables or the response thereof in its memory - could that not potentially distort your analysis?
6) How do you account for lagged effects? The land surface can be much slower than the atmosphere, responding on very different time scales (slowly changing...)
7) Have you seen this recent pre-print? https://arxiv.org/abs/2410.01770 
8) More of a suggestion, I personally feel like the contribution in this paper may be much more on the "results" side of things, and less methodological. Hence, it might be much more suitable for a journal and not for ICLR.
soundness: 2
strengths: 1) The work is original in the sense that it tried working around the challenges the particular dataset (being high-dimensional and not global, and samples auto-correlated) brings, leading a highly tailored analysis for the EarthNet2021 dataset.
2) Reproducibility is given, as the work is based on open data and open source code of a model.
3) The paper highlights the challenging nature of the relationship between NDVI and weather, often being non-linear, and highly variable in space.
summary: This work tries to uncover the relationship of meteorological drivers to satellite-derived vegetation greenness through studying the partial sensitivities of a deep neural network to input perturbations. More specifically, they work on the EarthNet2021 dataset, which contains samples ("minicubes") of high resolution NDVI at many places in Europe co-located with meteorological variables. First, they cluster each individual meteorological variable. Then they identify similar minicubes, by grouping according to the cluster ids of five individual meteorological variables, resulting in ~800 meteorologically similar minicubes. For each group, they compute partial sensitivities of the ConvLSTM model to additive perturbations of the input meteorology, obtaining marginal response curves of NDVI to weather.
weaknesses: 1) My major concern is that the presented analysis does not achieve its goal, that is explaining the complex relationship between NDVI and weather at fine resolution. Here, i think a more careful presentation and interpretation of the results needs to be done, after reading this paper, the reader should have an increased understanding of the drivers of NDVI dynamics. 
2) The used dataset EarthNet2021 and the Diaconu ConvLSTM models both suffer from a few drawbacks. It would have been much better to use the more recent GreenEarthNet dataset, which is an improved version of EarthNet2021 and comes with many improved baseline models. https://openaccess.thecvf.com/content/CVPR2024/html/Benson_Multi-modal_Learning_for_Geospatial_Vegetation_Forecasting_CVPR_2024_paper.html 
3) The writing lacks clarity, I had to read the paper multiple times to understand the approach. Perhaps a figure that visually explains the approach could help here, but more importantly, many parts of the paper need to be rewritten to follow a clearer structure.
4) I am unsure if the clustering approach is meaningful for Precipitation data, which is generally exponentially distributed and has few spikes and else zeros. 
5) Introducing around ~800 groups of minicubes with similar weather does not necessarily increase interpretability, as one would still need to look at 800 different plots with 5 panels each (Fig. 2...) to fully understand the relationship between NDVI and weather.
6) The perturbations are done univariate, however the considered weather variables are not independent of each other.
7) The perturbations for precipitation seem physically not meaningful.
8) The related works section is poorly written. It should reflect what has been done and how that related to the presented study. Also, it lacks many related works. You could start by checking the references of the EarthNet2021 paper, which already contains many more especially regarding NDVI.
9) In Fig. 1, I am not sure the absolute x-coordinate makes so much sense. The EarthNet2021 minicubes start at random start dates, i.e. they are during different seasons, so patterns in meteorology should rather be translation-equivariant, i.e. considering relative coordinates, I am thinking in the direction of wavelets and other filters here.

---

confidence: 4
contribution: 2
initial_rating: 5
presentation: 2
questions: 1. What is the current state of spatiotemporal interpretability methods, and what challenges do they face?

2. What is the motivation behind the proposed CSP technique, and what implications does it have for future research?

3. What are the technical advantages of the proposed method compared to existing spatiotemporal interpretability methods, and is it more suitable for surface process spatiotemporal prediction problems?

4. How reasonable is the assumption of independent distribution among meteorological variables? In reality, meteorological variables typically exhibit significant correlations.

5. Can the authors provide validation results on synthetic datasets and quantitatively compare the accuracy of their method with other time series interpretability methods?

6. NDVI may experience saturation effects under strong solar radiation. Does the interpretability model account for the instability in independent-dependent variable relationships caused by this saturation phenomenon?

7. Can the authors demonstrate their method's performance on truly high-dimensional data to validate its advantages over existing methods in handling high-dimensional data?

8. Is there an issue with the y-axis labeling in Figure 1?
soundness: 3
strengths: This paper demonstrates good originality. The authors' analysis of existing PDP and LIME methods is reasonable, and they accurately identify the scientific challenges in these two approaches. To address these issues, the authors propose a sound scientific hypothesis: "explaining grouped spatiotemporally homogeneous data." Based on this hypothesis, they construct a three-step pipeline using existing methods to aggregate spatiotemporally homogeneous data and add perturbations that preserve temporal continuity. Using this approach, the authors discuss the influence of three main meteorological factors on NDVI and draw some interesting conclusions.
summary: This paper discusses how to develop an interpretability pipeline for spatiotemporal prediction problems in earth science. The authors point out that existing interpretability methods, such as PDP or LIME, may have two limitations when dealing with spatiotemporal prediction problems: 

1. These methods lack consideration of spatiotemporal variations in independent variables, and their randomly generated samples for explanatory analysis may lack temporal continuity.

2. These methods ignore that the relationships between independent and dependent variables may vary with seasons, lacking the ability to explain temporal dynamic mappings.

To address these issues, the authors propose an interpretability pipeline called CLUSTER-SEGREGATE-PERTURB (CSP). This method first uses Soft DWT+Kmeans to cluster the temporal variation curves of meteorological data into several patterns. Based on the patterns of meteorological attributes obtained from clustering, the authors divide the dataset into several subsets according to meteorological pattern similarity. Finally, the authors create perturbation sets for various variables based on climate patterns and data averages for subsequent interpretability analysis.

In the analysis phase, the authors selected three influence variables: temperature (min, max, avg), pressure, and precipitation. They trained a ConvLSTM model on the EarthNet21 dataset and used marginal sensitivity and marginal correlation analysis to interpret ConvLSTM on the segregation-perturbation sets generated by CSP. The authors analyzed the results of marginal sensitivity and marginal correlation. The results show that NDVI is most sensitive to precipitation, followed by temperature and pressure.
weaknesses: However, the paper has significant issues in its research foundation and experimental validation. Additionally, there are shortcomings in the literature review and paper organization.

Research Foundation:
The paper primarily builds upon improving PDP and LIME methods. While these methods indeed show limitations in explaining time series predictions, recent studies have already extended these interpretability methods to time series prediction:

[1] Shi, H., Yang, N., Yang, X., & Tang, H. (2023). Clarifying relationship between pm2. 5 concentrations and spatiotemporal predictors using multi-way partial dependence plots. Remote Sensing, 15(2), 358.

[2] Liu, J., & Zhang, X. (2022). ReX: A Framework for Incorporating Temporal Information in Model-Agnostic Local Explanation Techniques. arXiv preprint arXiv:2209.03798.

The authors' failure to consider and analyze these methods undermines the reliability of their research foundation.

Method Validation:
Given that this is not the only post-hoc interpretation method for time series data, it is essential to validate whether this method is more suitable for surface process spatiotemporal prediction than existing methods. Specifically, the paper lacks discussion of:

1. Quantitative evaluation of interpretability method accuracy (typically achieved using surface process simulation synthetic datasets)
2. Comparative analysis with other time series interpretation methods
3. Response to research hypotheses regarding existing methods' reliability with temporal and high-dimensional data

Additional Issues:
1. The related work section merely lists existing methods without analyzing their limitations or providing technical motivation

2. The paper's structure is overly segmented, mixing original methodology (5.2), previous research methods (5.1), and experimental results (5.2)

3. While claiming existing methods struggle with high-dimensional spatiotemporal data (lines 21 and 52), the experiments only use data with 3 attributes and 5 dimensions

4. There's an inconsistency between the criticism of PDP's inability to isolate variable effects (line 53) and the paper's assumption of variable independence in Section 4

---

confidence: 4
contribution: 2
initial_rating: 3
presentation: 2
questions: c.f. weaknesses above.
soundness: 1
strengths: The paper addresses a timely and important subject: can we understand how a certain variable affects what is observed in deep learning models, without explicitly modeling the physical phenomena. This topic is of utmost importance for deep learning, as they are often critiqued as black boxes. This paper thus proposes a method to disentangle this black box, and shed a light on the importance of each variable on the observed process.
summary: The paper proposes a pipeline to quantify the importance of different input variables for a Convolutional LSTM. The authors construct a model that aims at predicting the NDVI in a chip, based on different meteorological input variables. Once this model is trained, the authors propose a pipeline called Clustering, Segregation, Perturbation (CSP), that aims at identifying the importance of the different meteorological variables used in the ConvLSTM.
The pipeline starts by clustering in time: using k-means, the method identifies clusters of observations that can be grouped together over time. The second step is to segregate the clusters by k-means. The last step is to perturb the inputs, and observe the impact of the perturbation on the model output.
weaknesses: It is unclear to me what the clustering process adds here. The proposed perturbation resembles a simple sensitivity analysis. In table 1, the results are presented as averages over the clusters, and in figure 3 the different clusters are shown individually, but it is hard to understand weather these results would be different if the perturbation had been performed without the clustering.

Towards the end, the paper aims at demonstrating the impact of the different variables, and states an approximatively 12 times greater importance of the precipitation compared to the pressure, and 5 times greater than temperature. This seems significatif, but hard to judge given that we are not presented with a baseline. Comparing this method to other well known methods of sensitivity analysis, like simple perturbation or noise injection could inform about the significance of the result.

Overall, the paper was hard to read. The paper is not very well structured, and the methodology hard to follow. Multiple specific mistakes (c.f. specific comments below) made it hard to make sense of the methodology.

# Specific comments

Line 28: what interesting nonlinear correlations

Line 74: what is the curse of dimensionality?

Line 151: assumption: is this reasonable? Environmental variables are known to be not completely independent of each other

Line 177: so a total of 30 * 32000 images?

line 179: The spatial resolution is 20 meters, the dimension of the image is 128x128 pixels. I'm guessing the value 2.56 is the side of the image, so km, not an actual surface (20 m * 128 px / 1000 m/km = 2.56 km, the surface would be 6.5536 km2).

Line 181: "150 images for 150 days", so 1 image per day?

Line 181: unclear what the spatial resolution is. Assuming 102.4 is the surface (unit is km2, but this might be the same error as above), then the spatial resolution is 126.5 meters. But if it's the side, as above, then the spatial resolution is 1.28 km, which is probably more reasonable for a weather variable. Although E-OBS has a 10km resolution, so that would still be different.

Line 183: What is the resampling procedure? Probably nearest neighbor for spatial, but what about temporal?
How many pixels cover the same area as the Sentinel-2 chip? Assuming a 126.5 meters resolution, 20 pixels cover the chip, but if we take the value 1.28 km resolution, then 2 pixels per side (4 total) cover the chip, which means all the pixels in the resampled product basically have the same value. Even less if we consider the E-OBS resolution of 10km.

Line 185-196: not clear

Line 206: I'm not familiar with the EarthNetScore, but if it's scaled between 0 and 1, a value of 0.3257 seems very low. Wouldn't that very low modeling quality influence the interpretation of the results?

Line 216: Soft-DTW: is this a pre-processing step of the data, or is this already one of the clustering steps of the CSP? Unclear

Line 243-245: as we established above, 1 to 4 pixels basically cover the entirety of the Setinel-2 chip, so the lack of variability stems from the resampling of the dataset, not from the variables themselves.

Line 251: again wrong units? (km2?)

Line 251: given the spatial resolution of the meteorological dataset, what is the goal of first upscaling, then downscaling? Is the downscaled data used somewhere else?
Otherwise this step seems overkill to me, just take the average over the sentinel-2 extent.

Line 256-307: rephrase this section, the definitions are all over the place.

Figure 1: I'm not sure I understand the axis, the y-axis is all NDVI?

Line 326-334: introduce Si after equation 9, this is confusing.

Line 412: "irrespective of the weather segment the sensitivity of the variables remained almost the same". Does that mean that the clustering wasn't needed?

Line 414: "weighted by the cardinality of the sample sets". Not sure what that means.

# Minor comments

Figure 1: For each line, use same limits for y
Figure 2: use same limits for NDVI

# Grammar comments

Citations miss-match between automatic citations and manual:
117
131/132
145
... and many more

Author name twice:
107/108
113/114
119
128
134
... and many more

Line 24: makes it sound like EarthNet2021 is your previous study

Line 100: missing point between citation and reference

Line 101: "insights. In..." new sentence here
line 102: remove authors, "propose"

Line 128: replace "concerning" with "on"

Line 143: "during dry years, like 1989, following another dry year." unclear what that means