Identifying Spatio-Temporal Drivers of Extreme Events

Necessary background knowledge and related work should be supplemented as thoroughly as possible.

Thank you for mentioning this issue. In the revised version, we will extend the introduction and related works sections to provide as much as possible background knowledge given the the page limit.

The anomaly detection part in Section 2 is redundant

In the paragraph "Anomaly detection algorithms" in Section 2, we wanted to give a broader overview why current methods for anomaly detection cannot be applied to the problem that is addressed in this work. We will revise this part in Section 2 and only focus on the most relevant anomaly detection approaches and emphasize the differences to our approach. In fact, none of these approaches is similar to our proposed approach.

Font size in figures

We will increase the font size in the figures.

We appreciate and thank the reviewer for recognizing and appreciating the efforts behind this work. Please see our responses to the questions below.

Code and datasets

Please note that we will release the code and datasets along with the documentations upon publication. We describe the datasets and benchmarks in the Appendix in Sections A, H, and I (pages 18-21, and 32-35). Please let us know if more information is required. We are happy to clarify any open questions and provide more details.

The contribution of the benchmark in comparison to the method

We appreciate that the addressed problem and proposed benchmarks are considered as an important contribution. Since the proposed problem has not been addressed before, we also present a novel approach that is not a simple adaptation of existing works to the problem. As you have acknowledged, the proposed approach outperforms adaptation of one-class unsupervised [25; 107; 37], reconstruction-based [44; 49], and multiple instance learning [92; 93; 94] approaches to the problem. We also included additional comparisons in the global response. We thus argue that both the benchmarks and the proposed method are a valuable and important contributions. However, we will change the order of the contributions in Section 1 as suggested.

The backbone and the global level analysis

We already compare 6 different backbones which differ in the number of parameters in the Appendix (Table 11, page 29). We agree that using Video Swin Transformer as backbone for the global level would be very expensive. However, the backbone can be replaced by more efficient backbones. We conducted an additional experiment where we replaced the attention block with a linear selective state space model (Gu et al. "Mamba: Linear-time sequence modeling with selective state spaces", 2023). The results are reported below for F1 scores on the validation set for anomalies/extremes detection:

Please note that we did not yet fully train the network with MAMBA as backbone, but the results already indicate that using Mamba instead of Swin Transformer achieves similar or better results when less parameters are used. In contrast to Swin Transformer, Mamba can be scaled to the global scale. Please note we already cover large areas at the continental scale (Table 17, page 33) and we think that identifying spatio-temporal relations at the continental scale is at the moment already challenging enough. We also mention this limitation already lines 345-348.

I have also listed several confusions and suggestions in the paper. If the author can address my concerns, I would be willing to support its publication here and recommend an increase in its score, even though it is more suitable for the DB track.

We appreciate your careful reading of our work and your suggestions. Please find below our responses to your questions that have not been already answered.

Line 113 and a simple illustrative example figure

We will revise the sentence. Please see the author rebuttal and the general response above. We have also included a figure in the PDF (Fig. 2).

The spatial resolution

We provide the spatial resolution in Section 4.2, page 5 and in Appendix Section I.I., page 34. Please note that depending on the coordinate system, the area on the Earth surface changes with respect to the location (i.e., it decreases toward the pole). The spatial resolution for ERA5-Land is $0.1^\circ \times 0.1^\circ$ on the regular latitude longitude grid. CERRA Reanalysis has a spatial resolution of $5.5$ km $\times~ 5.5$ km on its Lambert conformal conical grid. While the remote sensing data has a high resolution of $0.05^\circ \times 0.05^\circ$. We will make it clearer in the revision. We conducted an additional experiment regarding the impact of the spatial resolution:

This shows that the spatial resolution matters as expected.

Loss functions in Table 2(a)

Without the $\mathcal{L}\_{(anomaly)}$ loss, the detection of anomalies is not reliable since pixels at regions and intervals where no extreme event occurred can be assigned to $z_{q=1}$ (anomaly) as well. The other loss functions actually do not prevent that this is happening.

In case of $\mathcal{L}\_{(extreme)}$ multi-heads, we observe that anomalies are identified in a small subset of variables because the network omits some variables if there is a correlation with other variables. Please see Fig. 4 in the rebuttal PDF. In case of a single head such flips occur less often, but they can occur. If the $\mathcal{L}\_{(anomaly)}$ loss is used, such flips cannot occur and the multi-head improves both F1-scores by a large margin. When comparing row 2 and 4, there is a slight decrease in extreme prediction but a large improvement in anomaly detection. Note that there is always a trade-off between extreme and anomaly detection. The anomalies generate a bottleneck of information. When more information goes through the bottleneck the better the extreme prediction gets. Without any anomaly detection, the extreme prediction is best as shown in Table 12 (page 29), but the increase in F1 score is only moderate. This is also visible in row 2 and 3 in Table 2(b). Cross-attention improves extreme prediction, but it hurts the detection of anomalies since the information is propagated between the variables.

Thank you for the detailed review and the thoughtful feedback. We are glad that you found our work interesting and important for future research.

Terminology

Thank you for this suggestion. We will follow your suggestion and define the terminology in the introduction section of the paper in the revised version. Please see the discussion in the global response above.

Albedo & Soil temperature

Please note that ERA5-Land does not use data assimilation directly. The evolution and the simulated land fields are controlled by the ERA5 atmospheric forcing. We conducted 4 more experiments on both CERRA and ERA5-Land where we trained models that take only one variable al/fal or stl as input and predict the extreme events directly without the anomaly detection step. In all of these experiments, the F1-score was very low. In the next experiment, we increased the threshold for VHI and trained new models to predict extremes directly. The results for the validation set are shown below:

The first potential reason to consider is that some land surface variables might deviate from the reality. Another reason might be that when training only on extremes (VHI < 26), there are not enough samples to learn the relations. Please note that VHI is a combination of both TCI and VCI. Most extremes (VHI < 26) might result from a deficiency in both stl/t2m and vsw. This might also explain why stl and albedo cannot be that informative to predict very extreme events. We will discuss this issue.

Soil moisture in ERA5-Land

It is true that the volumetric soil water variable in ERA5-land has some biases. One solution is to use satellite observations for the top layer. However, our experiments showed that the model relates vsw anomalies with the extremes in VHI and provides reasonable predictions. Although we do not consider this as a major issue, we will mention it in the revised version.

Anomalies in precipitation

Thank you for this notice. At the moment, we treat each input variable in the same way and do not apply any pre-processing that is specific to a single variable. The proposed approach to pre-process precipitation is indeed an interesting direction. We will mention this in the discussion in the revised version.

Connecting anomalies in atmospheric variables to the anomalies in land surface states (VHI)

We agree that this is the mid-term goal, but this is beyond the scope of the paper. The purpose of the paper is to present a novel approach that addresses this very important problem and a benchmark that allows to systematically evaluate approaches for this new task. The impact will be two-fold. First, methods can be further developed that improve the detection of drivers on the synthetic datasets. Second, the anomalies that are detected by the method in atmospheric variables can be further investigated by statistically approaches.

The performance on the real-world data

We respectfully disagree regarding this point. The performance depends on the type and ratio of extremes, spatio-temporal resolution, the quality and consistency between the remote sensing and the reanalysis data. The performance is consistent with other recent works for predicting extremes on real-world data (Nearing et al. "Global prediction of extreme floods in ungauged watersheds", Nature, 2024). Note that the F1 scores substantially increase when the threshold on VHI is increased (see previous answer regarding stl1). We include a figure in the rebuttal PDF (Fig. 1), which shows the predictions. Given the quantitative and qualitative results, we think that the model provides reasonable predictions. Note that it is not required to predict all extremes in order to learn some relations from the predicted events.

Baseline as an interpretable forecasting

Please see the general author rebuttal and response above for the requested experiment.

The spatial aspect

We present a general approach that is not limited to droughts and specific input variables. We used droughts only as a real-world example. For instance, variables over sea regions can impact variables over land regions.

VHI is a general vegetation condition index

Thank you for this suggestion. We mention this issue in lines 987-989. We followed the general association of VCI and VHI with agricultural droughts, see i.e., (Hao et al. "Seasonal drought prediction: Advances, challenges, and future prospects.", Reviews of Geophysics, 2018). We will discuss this in the revision.

Predicting the exact values of VHI

Since we will include the new baselines based on forecasting and integrated gradients (see previous answer), we will briefly discuss the mentioned methods on vegetation forecasting. The mentioned work by Shams Eddin et al. ("Focal-TSMP: deep learning for vegetation health prediction and agricultural drought assessment from a regional climate simulation", GMD 2024) has shown that it is hard to predict VHI directly. Instead, they predict NDVI and BT and then normalize the predicted values to estimate VHI. Another reason is that some extremes cannot be derived from satellite products but are stored as binary or discrete variables in databases. Having such scenarios in mind as well, we decided to present extremes as binary variables. Extending the approach to continuous variables is a potential future direction.

Typos

Thank you for the careful reading. We will fix this typo and check the paper for any other typos.

How do you ensure your models are not overfitting?

We follow the common practice in climate science where we define different time periods for the training/validation/test sets (Table 17 in the Appendix, page 33). As also shown in the Appendix (Table 11, page 29), increasing the model parameters still does not show a sign of overfitting.

We appreciate the reviewer’s feedback for this work and that the reviewer recognized the novelty of this work. We respond to the questions below.

The method only shows results on droughts. Have you tried the method on other extreme events other than droughts?

We tested the algorithm on 9 types of extreme events; real-world agricultural droughts (Table 14 in Appendix, page 30), synthetic CERRA extreme events (Tables 1 and 3, pages 6 and 20), 5 variations of the synthetic CERRA extreme events (please see Fig. 3 where we changed the coupling between the climate variables and consequently changed the type of the events), synthetic NOAA extreme events (Tables 4 and 6, pages 20 and 27), and synthetic artificial extreme events (Tables 5 and 7, pages 21 and 27). The synthetic datasets show that the approach is not limited to droughts and can be applied to other extremes. The main difficulty is that suitable datasets for extreme events are rare, i.e., the dataset should have a high resolution and a long-term and large-scale coverage. Building a dataset from existing sources is thus a major effort (see Section H and I in the Appendix, pages 32-35). In the future, we aim to test the model on more types of events like flood, but this requires to prepare the data first.

The model is dependent on temporal resolution, which might not be well documented in all parts of the world. What if there is discontinuity in terms of temporal data?

We assume that you are referring to the temporal gaps in the reanalysis and remote sensing data. In fact, this is an issue for real-world data. To tackle this issue for the remote sensing data, a temporal decomposition was conducted to remove some discontinuity and aggregate the data into a weekly product. However, some pixels will still be empty. To tackle this issue, we first check if the pixel was covered by another satellite. If it was not the case, we flag the pixel as invalid and we discard it from the training and evaluation. Regarding the input reanalysis data, we first normalize the data using the pre-computed statistics and then replace the invalid pixels with zero values. We will add more details to the Appendix.

Binary masks tend to oversimplify real-world events.

Please note that we only use the binary masks as flags where anomalies or extremes are detected. Using the binary approach allows to use the approach in the future for other extremes that cannot be derived from satellite products but that are stored in a binary format in databases. We agree that an analysis regarding the physical implication requires the continuous values where binary masks only indicate the existence of anomalies or extremes, see paragraph "Physical consistency" in Section 5.2, page 9.

The method seems to bluntly connect anomalies with extremes without specific theoretical reasoning.

Our aim is to investigate the relations between extreme events and their drivers from a data-driven perspective. The synthetic examples demonstrate that our proposed approach is able to achieve this. In contrast to statistical methods, our method does not require a prior hypothesis about drivers for extremes; instead, it generates hypotheses that can be verified by statistical methods in a second step. We believe that this is an important direction since climate reanalysis provides huge amount of data and it is infeasible to test all combinations. Data-driven approaches are therefore needed to generate potential candidates. Please also see a related discussion in the first paragraph "Anomalies and extreme events detection in climate data" in Section 2.

If possible, please add the reasoning behind feature representations and extremes.

We are not sure if the previous answers already answered the question.

Thank you for your time and for reviewing our work. We are glad that you found the task and the problem we address in this work important. In the following, we answer your questions.

Difference between anomaly and extreme in this paper is not clear. What is the difference between anomaly and extreme in this paper? Could you provide some examples to illustrate this?

We will define the terms more precisely in the revision. Please see the author rebuttal and the global response above for clarification.

The title is somehow misleading, actually what the paper does is more about extreme event prediction instead of learning spatial-temporal relations. The motivation of model design is not clear. For example, why do you need to detect the anomaly and then classify the extreme events instead of detecting extreme events directly? Such a designed pipeline will lead to more accumulated errors.

There seems to be some misunderstanding about the objective of this work, probably caused by the terms "anomaly" and "extreme". Our work is not about extreme event prediction in the first place but about identifying the anomalous drivers of these extreme events like droughts. Note that droughts can be observed, but it is unclear what are the anomalies in the atmospheric or hydrological state variables that are spatio-temporally connected with a drought. Note that the anomalies in the atmospheric or hydrological state variables can occur earlier than the drought and at a different spatial location.

We aim to identify these anomalies that are spatio-temporally connected with an observed extreme event. Since we only observe droughts and are only interested in atmospheric or hydrological state variables that are spatio-temporally connected with droughts, we design a network (Figure 1) that spatio-temporally identifies anomalies in the atmospheric or hydrological state variables and predicts from these the droughts in an end-to-end fashion. In other words, we force the network to reduce the input variables to spatio-temporal anomalies (quantization) that are sufficient to predict the drought. Due to the quantization, the prediction of a drought is lower compared to predicting droughts without identifying anomalies in the input as reported in lines 314-319 and in Table 12 (Appendix C.5, page 29). The drop of the F1-score on drought detection, however, is relatively small. Nevertheless, we are interested in detecting the driving anomalies that are spatio-temporally connected to a drought and not the drought itself. We will clarify this and revise the title if necessary.

Anomaly detection is a classification problem with severe class imbalance problems, how did the authors tackle this problem?

To solve the class imbalance issue, we utilized a weighted binary cross entropy for $\mathcal{L}_{(extreme)}$. Please see Appendix C.3 (page 28) and Table 10 (page 29). We will make this clearer in the method section in the revised version.

What is the ratio of extreme events to ordinary events?

The ratios of extremes are reported in Tables 3-5 (second last column, pages 20-21) in the Appendix. For the synthetic CERRA reanalysis (Table 3, page 20), the ratio is $1.16\%$ while the ratio of the correlated anomalies with these extremes is $1.69\%$ and the ratio of the random uncorrelated anomalies with extremes is $1.32\%$. For convenience, we summarize the ratios below:

Please note that there is no ground truth for anomalies in the real-world dataset. We only report the ratio of extreme events, which can be detected using remote sensing data:

When I read spatio-temporal relations, I thought this paper would build spatio-temporal graphs to describe the relations. Have you considered using graphs to tackle this problem?

Note that our real-world data like the CERRA dataset has a spatial resolution of $1069\times1069$ (Table 17 in Appendix, page 33) and we consider 8 time steps. A spatio-temporal graph would consists of $1069\times1069\times8=9,142,088$ nodes, which would be computationally very expensive. Even for visualizing the results shown in Figures 16, 18, 20, 22, 24, 26 and 28 (pages 37-43), a spatio-temporal graph would not be suitable.

The writing can be improved. For example, there are typos, such as 'MIL Is a weakly ...'.

Thank you for the careful reading of our paper. We will fix this typo and check the paper for any other typos.