Automated Input Variable Selection for Analog Methods Using Genetic Algorithms

Analog methods (AMs) have long been used for precipitation prediction and climate studies. However, they rely on manual selections of parameters, such as predictor variables and analogy criteria. Previous work showed the potential of genetic algorithms (GAs) to optimize most of the AM parameters. This research goes one step further and investigates the potential of GAs for automating the selection of the input variables and the analogy criteria (distance metric between two data fields) in AMs. Our study focuses on the prediction of daily precipitation in central Europe, specifically Switzerland, as a representative case. Comparative analysis against established methods demonstrates the superiority of GA‐optimized AMs in terms of predictive accuracy. The selected input variables exhibit strong associations with key meteorological processes that influence the generation of precipitation. Further, we identify a new analogy criterion inspired by the Teweles‐Wobus criterion, which consistently performs better than other Euclidean distances and could be used in classic AMs. In contrast to conventional stepwise selection approaches, GA‐optimized AMs display a preference for a flatter structure characterized by a single level of analogy and an increased number of variables. Overall, our study demonstrates the successful application of GAs in automating input variable selection for AMs, with potential implications for application in diverse locations and data exploration to predict alternative predictands. In a broader context, GAs could be used to perform input variable selection in other data‐driven methods, opening perspectives for a broad range of applications.

Stepwise AMs for predicting precipitation commonly have a first level of analogy based on the atmospheric circulation.The variable of interest is the geopotential height (Z) at various pressure levels and specific times throughout the day (Table 2; Obled et al., 2002;Horton et al., 2018).Bontron (2004) introduced a second level of analogy based on a moisture index that is the product of the relative humidity at 850 hPa and the total precipitable water (RM3 method in Table 2).Other consecutive studies selected different pressure levels (Horton et al., 2018, method RM4 in Table 2) or added a wind component to the moisture index (Marty, 2010).Ben Daoud et al. (2016) inserted an additional level of analogy between the circulation and the moisture analogy based on the vertical velocity at 850 hPa (methods RM6 in Table 2) and named it "SANDHY" for Stepwise Analog Downscaling method for Hydrology (Ben Daoud et al., 2016;Caillouet et al., 2016).
To calibrate the method, a semi-automatic sequential procedure (Ben Daoud et al., 2016;Bontron, 2004;Radanovics et al., 2013) has often been used to optimize the size of the domain and the number of analogs.However, predictor variables, vertical levels, temporal windows (time of day), and analogy criteria were manually selected.This manual selection requires the comparison of numerous combinations and a comprehensive assessment of some parameter ranges.Moreover, the sequential calibration procedure successively calibrates the different levels of analogy, and thus it does not handle parameters inter-dependencies.Considering these limitations, Horton et al. (2017) introduced a global optimization of the AM using genetic algorithms (GAs).Using this approach, an automatic and objective selection of the temporal windows, the vertical levels, the domains, and the number of analogs became possible, improving the prediction skills of the method (Horton et al., 2018).A weighting of the predictor variables has also been introduced.The only parameters left for manual selection were the meteorological variables and the analogy criteria.
Selecting predictors for precipitation prediction with AMs in Europe has been the focus of multiple studies aiming to improve prediction skills (Ben Daoud et al., 2016;Bontron, 2004;Gibergans-Báguena & Llasat, 2007;Obled et al., 2002;Radanovics et al., 2013).Thus, the relevant predictors are likely to be known nowadays and supported by expert knowledge.However, transferring AMs to a region with different climatic conditions or to another predictand would involve reconsidering the selected meteorological variables.This work aims to test a fully automatic optimization of all AM parameters, including the selection of the meteorological variables and even the analogy criteria, using GAs.GAs have already been used for input variable selection (IVS) in other contexts (D'heygere et al., 2003;Huang et al., 2007;Cateni et al., 2010;Gobeyn et al., 2017).
GAs have also been used in the context of AMs for other tasks, such as the selection of optimal vertices in an unstructured grid approach to reduce computational resources when working with high-resolution data (Hu & Cervone, 2019).An alternative approach to IVS, proposed by Hu et al. (2023), is to compress multiple predictors into latent features using a deep learning network and then select the analogs in this latent space.This approach eliminates the need for the prior selection of predictors; however, it sacrifices the advantage of interpretability provided by an analogy computed on the original variables.
Here, we seek to assess the potential of GAs for input variable selection in the context of the analog method.Moreover, we want to test the GAs' ability to jointly select the distance metric in addition, that is, the analogy criterion.To compare with well-established AMs, daily precipitation in central Europe, specifically in Switzerland, has been chosen as predictand.Also, as is often the case, the AMs were optimized in the perfect prognosis framework, using predictors from reanalyses.This work focuses mainly on the proof of concept of automatic IVS for AMs rather than the details of the case study.
The paper is organized as follows.Section 2 describes the data sets, the fundamentals of AMs, the characteristics of the GAs implementation, the software used, and the details of the experiment setup.Section 3 presents the results of different analyses, such as the selection of the best predictor variable, the relevance of various AM structures, and the accuracy of the optimized methods.Section 4 discusses some findings of the work.Finally, Section 5 summarizes the main contributions of the work and open perspectives for applications of the developed approach.

Data
The target variable (predictand) is daily precipitation derived from the RhiresD gridded data set from MeteoSwiss (2021).It is a daily aggregation (from 06 UTC of day D to 06 UTC of day D + 1) at a 1 km resolution with data from 1961 onward.It is produced using an interpolation scheme between gauging stations (Frei & Schär, 1998).The gridded data were spatially aggregated across 25 catchments of about 200 km 2 (Table 1).These catchments were chosen to cover the different climatic regions of Switzerland (Schüepp & Gensler, 1980), as illustrated in Figure 1.
As is often done in the context of the perfect prognosis framework, we used variables provided by global reanalyses.Although most reanalysis provides good quality data in Europe, differences still exist, and the choice of the reanalysis data set can impact the accuracy of the AM even more substantially than the choice of the predictor variables (Horton & Brönnimann, 2019).Thus, it was considered advisable to test some of the following analyses with another reanalysis to assess the robustness of the selected variables.The main reanalysis used in this work is ERA-Interim (ERA-I, Dee et al., 2011), which was produced by the European Center for Medium-Range Weather Forecasts (ECMWF) and covers the period from 1979 to 2019.The forecast model uses a hybrid sigma-pressure vertical coordinate on 60 layers and has a T255 horizontal resolution Z850@12hr MI700@24hr Z700@24hr MI600@12hr Z400@12hr RM5 T925@36hr Z1000@12hr MI925@12 + 24hr Ben Daoud et al. ( 2016) T600@12hr Z500@24hr MI700@12 + 24hr RM6 T925@36hr Z1000@12hr W850@06-24hr MI925@12 + 24hr Ben Daoud et al. ( 2016) T600@12hr Z500@24hr MI700@12 + 24hr Note.The analogy criterion is S 1 for Z and RMSD for the other variables.The predictors are described by meteorological variables (e.g., Z), pressure levels (e.g., 1,000), and time windows (e.g., 12hr UTC).Multiple time steps can be considered.For example, 'MI850@12 + 24hr' indicates that the moisture index is being considered at 12 and 24hr UTC.In the case of 'W850@06-24hr', all 6-hourly time steps between 6 and 24hr UTC are being used.Z, geopotential height; T, air temperature; W, vertical velocity; MI, moisture index.

Water Resources Research
10.1029/2023WR035715 (about 79 km) and a 30 min time step.The output variables have a grid resolution of 0.75°.This work started before the release of ERA5, the successor of ERA-I.
The Climate Forecast System Reanalysis (CFSR, Saha et al., 2010a), provided by NCEP, was used for the first experiment to compare the results obtained with ERA-I.The model used to produce CFSR has a horizontal resolution of T382 (about 38 km) and 64 levels on sigma-pressure hybrid vertical coordinates.The period covered is 1979 to August 2019, and the output variables have a spatial resolution of 0.5°.
Finally, ERA5 (Hersbach et al., 2019) was used for the last analysis.ERA5 provides more variables and a higher spatial (0.25°, but used here at 0.5°) and temporal resolution (hourly, but used here at a 3-hourly time step).ERA5 assimilates considerably more data than ERA-I and provides, among others, more consistent sea surface temperature and sea ice, improved representation of tropical cyclones, better balance of evaporation and precipitation, and improved soil moisture.ERA5 also relies on more appropriate radiative forcing and boundary conditions (e.g., changes in greenhouse gases, aerosols, SST, and sea ice) (Hersbach et al., 2019).

Analog Methods
AMs are based on the rationale that two similar synoptic situations may produce similar local weather (Lorenz, 1956(Lorenz, , 1969)).It thus consists of extracting past atmospheric situations similar to a target date.Selected predictor fields define this similarity.The analogy is defined by.
2. The vertical levels at which the predictors are selected.
3. The spatial windows (domains) over which the predictors are compared.4. The hours of the day at which the predictors are considered.5.The analogy criteria (distance metric to rank candidate situations).6. Possible weights between the predictors.7. The number of analog situations N i to select for the level of analogy i.
AMs usually start with a seasonal preselection to cope with seasonal effects (Lorenz, 1969).The seasonal preselection is often implemented as a moving window of 120 days centered around the target date (Ben Daoud et al., 2016;Bontron, 2004;Horton et al., 2012;Marty et al., 2012).Alternatively, the candidate dates can be preselected based on similar air temperatures at the nearest grid point (Ben Daoud et al., 2016, methods RM5 and RM6 in Table 2).In this work, we used the temporal moving window to reduce the number of potential candidate dates and, thus, the computing time.
The first level of analogy in AMs for precipitation is often based on the atmospheric circulation using the geopotential height (Z) at different pressure levels and hours of the day (Table 2).The distance (analogy criterion) between two Z fields is calculated on the vector components of the gradient, that is, using the difference between adjacent grid cells, rather than comparing absolute values.The Teweles-Wobus criterion (S 1 , Equation 1, Teweles & Wobus, 1954;Drosdowsky & Zhang, 2003) was identified as the most suitable by different studies (Bontron, 2004;Guilbaud & Obled, 1998;Wilson & Yacowar, 1980;Woodcock, 1980).It is defined as follows: where Δẑ i is the gradient component between the ith pair of adjacent points from the geopotential field of the target situation, and Δz i is the corresponding gradient component in the candidate situation.The gradient components are computed in both the latitude and longitude directions.S 1 ranges from 0 to 200.The smaller the S 1 values, the more similar the shape of the pressure fields and therefore the atmospheric circulation.S 1 was developed to verify prognostic charts (Teweles & Wobus, 1954).It was computed using pressure differences between stations arranged in north-south and east-west lines.The "difficulty coefficient" (the denominator) reduces the influence of the seasons and the strength of the weather systems on the score.
After selecting a certain number of analog dates based on Z, subsequent steps can be added to subsample a lower number of analog situations based on other predictors.The method developed by Ben Daoud et al. (2016) has, for example, three levels of analogy (and a preliminary level for the preselection) where each level subselects a smaller number of analog situations from the candidates provided by the previous level (RM6 method in Table 2).For other predictors than the geopotential height (e.g., for moisture variables), classic criteria representing Euclidean distances between grid point values are used.
The output of the AM is a probabilistic prediction for the target day.It is provided by the empirical conditional distribution of the N i predictand values corresponding to the N i dates selected at the last level of analogy.

Genetic Algorithms
Genetic Algorithm (GA) is a global optimization technique inspired by genetics and natural selection (Holland, 1992).It belongs to the family of evolutionary algorithms and comprises different operators such as natural selection, couples selection, chromosome crossover, mutation, and elitism.These operators act on parameter sets of the problem to optimize by mixing, combinations, and random modifications.GA aims to combine, over time, the strength of different parameter sets and to explore the parameter space while converging toward the global optimum.The optimization starts here with 2,000 random parameter sets (as defined in Section 2.2) and is stopped when the best parameter set cannot be improved after 30 iterations.
A variant of GA has been tailored to optimize AMs by Horton et al. (2017).All the parameters of the method, except the meteorological predictor variables and the analogy criteria, have already been successfully optimized using GAs (Horton et al., 2018).All parameters were optimized jointly on the different levels of analogy.The use of GAs provided for the first time an objective and global optimization of AMs, resulting in gains in prediction accuracy.To bring the optimization further, the selection of the predictor variables and the analogy criteria were performed here by GAs.
The reason why the predictor variables and analogy criteria were left out in the previous GA-AM setup by Horton et al. ( 2017) is the different nature of these variables.The parameters optimized so far by Horton et al. (2017) were quantitative variables, that is, numerical values (e.g., location and size of the spatial windows or the number of analogs), which have a notion of continuity.However, the parameters characterizing the selected predictors or analogy criteria are entries in a list of possible variables/criteria that can be considered as categorical variables as there is no relationship among entries.They are treated as arrays of independent values by the algorithm.Therefore, the mutation operator relying on a search radius in the parameter space (Horton et al., 2017) cannot be applied.Instead, a simple random sampling was used for these parameters when selected for mutation.In addition to the increased difficulty due to the higher number of parameters to optimize, this aspect will likely slow down the optimization.
In GAs, the mutation operator changes a parameter value (gene) if this parameter was selected to mutate (all parameters have a certain mutation probability).The new value assigned depends on the rules of the mutation operator applied.This operator enables the optimization to explore new areas of the parameter space and was shown to have the greatest impact on the success of AM optimizations (Horton et al., 2017).Thus, as suggested in Horton et al. (2017), five variants of this operator were used in parallel optimizations (see details in Appendix B): three variants of the non-uniform mutation (Michalewicz, 1996), the multiscale mutation, and the chromosome of adaptive search radius.The non-uniform mutation aims to reduce the magnitude of the search in the parameter space with the evolution of the population to transition from the exploration of the whole parameter space to the exploitation of local solutions.This operator has three controlling variables, which makes it difficult to adjust, and thus is used with three different configurations.The multiscale mutation considers both exploration and exploitation in parallel.It has no controlling parameters and no evolution during the optimization.The chromosome of adaptive search radius was introduced by Horton et al. (2017) and is inspired by the non-uniform mutation.It takes an auto-adaptive approach by adding two chromosomes, one for the mutation rate and one for controlling the search magnitude (see details in Horton et al., 2017).Therefore, it has no controlling parameters, is thus easier to use, and automatically transitions from the exploration phase to exploitation.

Software
The optimization of AMs with GAs is implemented in the open-source AtmoSwing software (https://atmoswing.org/) (Horton, 2019a) that has been used for this work.AtmoSwing is written in object-oriented C++ and has been optimized for computational performance.It scales well on HPC infrastructures, as the different members of the GAs populations, that is, the various parameter sets, can be assessed in parallel using multiple independent Water Resources Research 10.1029/2023WR035715 threads.However, due to the increasingly large number of assessments needed by GAs with the increasing complexity of the problem, a further reduction in computing time became necessary.Indeed, while applying AMs to perform a prediction for a single target date is a very fast and light process, GAs require a substantial amount of parameter assessments over long calibration periods.
Despite being simple methods, AMs require many comparisons of gridded fields during the calibration phase.For example, this work used a 24-year calibration period.For each target day, a gridded predictor needs to be compared to about 2,820 candidate situations (24*120-60, using a 120-day temporal window minus 60 days in the target year that are excluded).Over the entire calibration period, this amounts to about 2.47⋅10 7 (24*365*2,820) field comparisons per predictor of the first level of analogy.Here, one optimization required, on average, about 200 generations made of 2,000 individuals, which brings the average number of grid comparisons to about 1 ⋅ 10 13 per predictor of the first level of analogy.The comparison of the gridded predictors-that is, the calculation of the analogy criteria-was identified by profilers as the most time-consuming task, despite using the efficient linear algebra library Eigen 3 (Guennebaud et al., 2010).
A first attempt to reduce the computing time was based on storing the whole history of the optimization in memory and looking up for similar already-assessed parameters to a newly generated parameter set.However, this approach turned out to be even more time-consuming after several generations and led to memory issues for long optimizations.Finally, computation using graphics processing units (GPUs) was implemented for this study in a new version of AtmoSwing, v.2.1.2(Horton, 2019b).The calculation of the analogy criteria has been written using NVIDIA's CUDA.The details of the implementation and the results of a benchmark experiment can be found in Appendix A. When optimizing the methods using ERA5 at a 3-hourly time step and 0.5°resolution, the difference is substantial.One generation (2,000 evaluations) took 8 to more than 10 hr using 20 CPU threads, while 50-80 min were needed using 3 CPU threads and 3 GPU devices (NVIDIA GeForce703 RTX 2080).

Experimental Setup
The experiments were carried out over a 30-year period, from 1981 to 2010, divided into a calibration period (CP) and an independent validation period (VP).An additional test period (TP), covering the years 2011-2017, was allocated to evaluate the accuracy of the optimized methods on unseen data (Section 3.4).To reduce the impact of potential inhomogeneities in the time series, the selection of the validation period (VP) was evenly distributed over the entire series (as in Ben Daoud, 2010).A total of 6 years was used for the VP by selecting one year out of every five (explicitly: 1985, 1990, 1995, 2000, 2005, 2010).The archive period (AP), where the analog dates are being retrieved, is the same as the CP.The VP is also excluded from the AP (days from the VP were never used as candidate situations for the selection of analogs), as well as a period of ±30 days around the target date to exclude potential dependent meteorological situations.Unless stated otherwise, all results are presented for the VP.
The GAs optimized all parameters of the method.Only the AM structure (number of analogy levels and predictors) was not optimized.Different structures were tested in Section 3.2.For each level of analogy and each predictor, the following parameters were optimized within the corresponding ranges.[ 10,20].The spatial windows differ between predictors, even within the same level of analogy. 5. Analogy criterion: see Section 2.5.2. 6. Weight [0, 1] with a precision of 0.01 (0.05 for Experiment 2).The optimizer can turn off a variable by setting its weight to zero.7. Number of analogs: varies according to the structure, but with an overall range of [5, 300] and a step of 5.The optimizer can turn off a level of analogy by setting its number of analogs to the same value as the previous level of analogy.
The CRPS (Continuous Ranked Probability Score; Brown, 1974;Matheson & Winkler, 1976;Hersbach, 2000) was used to assess the accuracy of the predictions and is the objective function used by the GAs.It evaluates the

Water Resources Research
10.1029/2023WR035715 predicted cumulative distribution functions F(y), here of the precipitation values y associated with the analog situations, compared to the single observed value y 0 for a day i: where H y y 0 i ) is the Heaviside function that is null when y y 0 i < 0, and one otherwise; the better the prediction, the lower the score.The CRPS was here computed using the rectangle approximation based on the discrete precipitation values provided by the analog dates.

Meteorological Variables
The meteorological variables were considered for different types of vertical levels: surface or entire atmosphere (to capture, for example, the moisture content of an entire air column), pressure levels (1,000, 950, 900, 850, 800, 700, 600, 500, 400, 300, 200 hPa, to capture the vertical structure), potential temperature levels (290,300,310,320,330,350,400 K, necessary to include potential vorticity), and potential vorticity levels.The selected variables are listed in Table 3.The optimization can pick any variable on any level type and value, as long as it is available.Precipitation variables from reanalyses were not considered potential predictors.Precipitation is often considered as a predictor in AMs used in a post-processing context (where the same precipitation product is used for training and then predicting).However, AMs targeting downscaling tasks or alternative forecasts to NWP models do not rely on precipitation, as a method developed in the perfect prognosis context (using reanalyses data sets that can potentially assimilate precipitation data) would then be difficult to use in other contexts (using other model outputs) due to the high uncertainties and the biases associated with precipitation predicted by an NWP or a climate model.
The variables were standardized (using the overall climatology) on-the-fly by AtmoSwing when loaded from files.Standardization has no impact on the selection of analog situations for a single predictor, but it makes the combination of predictors within one level of analogy more balanced, as they might have very different orders of magnitude and units.It allows for a more effective optimization of the weights between predictors.

Analogy Criteria
During the optimization, the GAs can select different analogy criteria for the different predictors, with the aim of identifying the best criteria per predictor.These are then combined using weights to provide a single analogy distance value.The most common analogy criteria in AMs are the Root Mean Square Deviation (RMSD) and the Teweles-Wobus criterion (S 1 , see Section 2.2).Other criteria were made available to the GAs in order to explore potential new characterizations of the analogy metrics.Two of these criteria are new and are derived from S1.The potential criteria made available to the GAs are the following.
1. RMSD: the Root Mean Square Deviation.2. MD: the Mean Absolute Difference, or Mean Absolute Error.It differs from the RMSD in that the differences are not squared.3. S 1 : the Teweles-Wobus index as defined in Equation 1 from Section 2.2.It consists of a comparison of the gradients, mainly used for the geopotential height.4. S 2 : inspired by the Teweles-Wobus index, we introduced a new criterion based on the second derivative of the fields instead of the gradients: where ∇ 2 xi is the second derivative between the ith triplet of adjacent points from the predictor field of the target situation, and ∇ 2 x i is the corresponding second derivative in the candidate situation.Please note that it differs from the S 2 index from Teweles and Wobus (1954).
where xi is the ith point from the predictor field of the target situation, and x i is the corresponding point in the candidate situation.The reason for adding such a criterion was accidental, as it was an erroneous implementation of S 2 .However, it turned out to be relevant (see Sections 3 and 4).
6. DSD: difference in standard deviation over the spatial window.It is a non-spatial criterion, as the location of the features does not matter.7. DMV: absolute difference in mean value.It is also non-spatial, as the means are computed over the spatial window before comparison.

Design of Experiments
The input variable selection with GAs has been assessed in sequential steps.First, GAs were used to identify the single best predictor variables and their associated analogy criteria for each catchment (Section 3.1).The objective was to assess the consistency of the selected variables in the most straightforward configuration.Then, since AMs can be made of different levels of analogy with multiple predictors, the second experiment assessed the accuracy associated with different structures and the ability of GAs to deal with these, using a limited number of catchments (Section 3.2).Based on these results, the third experiment performed the input variable selection for each catchment (Section 3.3).For all experiments, the GAs used the CRPS of the precipitation prediction as the objective function.The selection of the meteorological variables and the analogy criteria, along with the other parameters, is thus done to improve the accuracy of the AM.

Best Single Variables
The first experiment evaluates the use of GAs to select a single predictor variable and analogy criterion for each catchment.The selection has been performed using ERA-I (Figure 2) but also CFSR for comparison (Figure 3), with six optimizations per catchment and data set.The six optimizations were based on different mutation operators (the five variants but twice the chromosome of adaptive search radius).The purpose of using two reanalyses is to assess the consistency and possible differences in the variable selection between two data sets.
One of the first elements that can be seen for both data sets is the dominance of the S 0 criterion, selected 60% of the time for ERA-I and more than 55% of the time for CFSR, along with the other criteria based on Teweles-Wobus (Figure 4).The other analogy criteria were rarely selected, if at all.The same applies to the RMSD, commonly used in analog methods.The GAs could better predict using S 0 as a metric for the Euclidian distance between the predictor fields.This result is further discussed in Section 4.
The variable selection results show some variability per catchment but similar accuracy.Although GAs can, in theory, identify the global optimum, this search is highly time consuming for such complex problems, and we have to stop the optimizations at a good-enough solution.These factors explain the variability that can be observed in the results.Nevertheless, this variability provides information about alternative variables with almost the same predictive accuracy.
Figures 2 and 3 demonstrate that the optimal variables can vary across different regions.Figure 5 illustrates this information spatially for the ERA-I variables.In terms of similarities, the vertical velocity (W) at 700 and 800 hPa is the most frequently selected for both data sets and is quantified using the S 0 criteria.Upward vertical winds at these levels are typically associated with precipitation generation.Within the southern Alpine climatic region (catchments 19, 20, 21), Z (based on the S 1 criterion) emerges as the best single predictor for ERA-I, which is not so clear with CFSR.Heavy precipitation events in this region are primarily the result of orographic effects related to sustained southerly advection of moisture-laden air masses (Massacand et al., 1998).Other regional clusters can be observed using ERA-I, such as the meridional wind V (with S 1 ) in the eastern part of Switzerland, also likely related to the southerly advection, STR(D) (surface net thermal radiation and surface thermal radiation downwards) in northern Switzerland (see the discussion in Section 4), and the second derivative of Z (with S 2 ) for several catchments at similar latitudes.The second derivative of Z is also frequently selected for CFSR.While cloud water (CWAT) is often chosen from CFSR, it is not available directly in ERA-I.3 for the variables abbreviations) from ERA-I for the 25 catchments (abscissa).The colors represent the analogy criteria, and the size of the dots is proportional to the accuracy of the resulting method (the larger the dots, the better), within a range of 5% of the best result (those with lower accuracy are hidden).

Assessment of AM Structures
The analysis of different AM structures (Section 2.5.3)aims to identify the best performing structures, that is, the optimal number of analogy levels and predictors.We first considered one to four levels of analogy, with one to four predictors per level.Five optimizations were performed for each of these 16 structures with the different mutation operators.As this assessment requires 80 optimizations, it was performed on only four catchments (L'Allaine (1), Sitter (15), Doveria (19), Flaz ( 25)).These were selected to maximize the diversity of climatic conditions represented.A complementary analysis was performed on two catchments (L'Allaine (1) and Doveria ( 19)) to explore the use of up to eight predictors on one and two levels of analogy.These experiments also allowed comparing the performance of the mutation operators for different problem complexities.
Even though the structure is provided to the GAs, it can still evolve to a simpler version by assigning a zero weight to some predictors or by setting the same number of analogs for two successive levels of analogy.This simplification often happened, such as that no solution ended up with the structure 4 × 4 (four levels of analogy with four predictors each).The best performing methods in the validation period were always made of one or two levels of analogy (Figures 6 and 7).Although some AMs have up to four levels of analogy (Section 2.2), the use of normalized variables and weights might here favor their combination in the same level of analogy.Methods with fewer levels of analogy present less of a hierarchy among the predictors.However, not having a systematic constraint by the atmospheric circulation, as in most AMs, results in more influence from other variables.Although atmospheric circulation is often of primary importance for heavy precipitation events, there may be situations where it is preferable to relax these constraints.Nevertheless, we cannot conclude that two levels of analogy are the maximum to be considered, as the optimizer might have failed to optimize complex structures satisfactorily.
The results also show notable performance differences between the mutation operators (Section 2.3).The chromosome of adaptive search radius (option #1) provides the best performing parameter sets 76.3% of the time for the calibration period and 62.5% of the time for the validation period (Figure B1).The second best is the nonuniform mutation with a mutation probability (p mut ) of 0.1 (option #4), which is the best option for 11.3% of the optimizations for the calibration period and 21.3% for the validation period.However, the same operator with a mutation probability (p mut ) of 0.2 (option #5; G m,r = 100) is the worst-performing option, with a success rate of  1.3% for the calibration period and 2.5% for the validation period.It quite well illustrates the difficulty of tuning such operators and the risk of a badly configured mutation operator, and thus the benefit of an auto-adaptive option such as the chromosome of adaptive search radius with no controlling parameters.Moreover, the latter usually performed better for more complex AM structures.

Full Optimization
The third experiment used different AM structures to perform the full input variable selection for each catchment.
Only the chromosome of adaptive search radius has been used because of its higher performance.

Using Variables From ERA-I
Based on the previous results, three AM structures were selected: 1 level of analogy with 8 (1 × 8) or 12 predictors (1 × 12), and 2 levels with 6 predictors (2 × 6) (Section 2.5.3).Two optimizations were performed by structure and catchment.The structure with two levels of analogy (2 × 6) turned out to be simplified by the GAs to a single level of analogy (1 × 6) for several catchments.Consequently, this structure resulted in lower accuracy as fewer predictors were used.Therefore, only structures with a single level of analogy (1 × 8 and 1 × 12) are further analyzed here.Figure 8 shows the different variables selected for each catchment along with the analogy criteria (color) and the weights (size).Figure 9 synthesizes the 30 most often selected variables and the associated analogy criteria, temporal windows, and spatial windows across catchments.These results again show a strong dominance of the S 0 , S 1 , and S 2 analogy criteria, the others being only rarely selected, including RMSD.S 0 is most often selected.The properties of S 0 are further discussed in Section 4.
Vertical velocity (W) at 700 hPa (and sometimes at 600 or 800 hPa) is the most frequently selected variable, also for the catchments that previously selected another best single variable (Section 3.1).Those with higher elevations and located in the southern part of the country additionally selected W at 500 hPa or even higher.
The surface solar radiation downward (SSRD) is the second most selected variable and is mainly relevant when compared in terms of gradients (S 1 ) rather absolute values.Other radiation variables occupy the fourth and  fifth ranks, such as surface thermal radiation downwards (STRD) and surface net thermal radiation (STR).These are mainly relevant when compared in terms of absolute values (S 0 ), although there is a non-negligible representation of the S 1 criteria (see discussion on radiation variables in Section 4).
CAPE is the third most selected variable, and the total column water (TCW) is the sixth variable.In the ninth position comes the meridional wind at 10 m, but using S 1 or even S 2 .The derivative of the wind can be informative on the location of frontal systems and convergence or divergence zones.Then comes the meridional wind on the PV level.The 2 m temperature has the twelfth position and is compared in terms of gradients (S 1 ), which can reflect the position of fronts.Follows the geopotential height (Z) at 700 and 600 hPa compared primarily using the second derivatives of the fields (S 2 ).The curvature of the geopotential height helps to identify and characterize synoptic-scale features, such as ridges and troughs in the atmosphere.A bit further down on the list, SLP is also compared in terms of its second derivative.Other variables such as RH, PV, D, and U also populate the 30 best variables.
The optimal spatial windows (Figure 9) cover Switzerland most of the time, with different extents depending on the variables.For example, while the medians of the optimal domains for W and CAPE are slightly larger than Switzerland, PV is here considered over a larger domain.The 2 m temperature (T2m) is characterized by unusual longitudinally extended domains, with the main body in southern Switzerland extending to the northern Mediterranean.Thus, it likely represents information at a synoptic scale, such as the location of fronts, rather than local conditions.Note that SST was also in the pool of potential variables, but has never been selected as relevant.
The optimal temporal windows (time of day) show substantial variability between the predictor variables.At the lower end of the range is TCW, which is considered better at the beginning of the precipitation accumulation period (06 UTC).The top of the range (06 UTC the next day, corresponding to the end of the accumulation period) was favored by the divergence (D at 285°K) and some low-level W (W900 and W950) or Z (Z900).It should be noted here that the radiation variables used were cumulative variables that were not decomposed prior to the analysis.Therefore, most of the selected temporal windows correspond to the beginning of the accumulation period, that is, 15 UTC.

Using Variables From ERA5
A similar experiment has been carried out using ERA5 and a single method structure (1 × 12).ERA5 has been used at a 3-hourly time step, which may be more relevant than 6-hourly when considering radiation variables, and at a 0.5°spatial resolution.The potential analogy criteria were limited to S 0 , S 1 and S 2 and the spatial domains were slightly reduced (latitudes = [39, 55], longitudes = [ 4, 20]).If previously the weights could be null for a predictor, a minimum of 0.01 was enforced here to force the GAs to select a relevant predictor.Finally, some predictors, often selected in the previous experiment, were forced: W700 (with S 0 criterion), CAPE (with S 0 criterion), TCW (with S 0 or S 1 criteria); leaving nine predictors unconstrained.
In addition, only the variables found relevant when using ERA-I were selected as potential predictors, thus decreasing the pool of variables.Also, potential temperature levels and PV levels were not considered further.However, cloud cover variables were added to the potential predictors to assess whether the radiation variables served as a proxy for cloud cover.Therefore, this experiment should not be considered a complete exploration of ERA5 as it builds on the results obtained for ERA-I.
The selected variables from ERA5 are shown in Figures 10 and 11.When compared with the ERA-I results, TCW gained importance, as it was the most selected variable here.Similarly, the relative humidity at 1,000 and 1850 hPa increased in importance as if its relevance improved in ERA5.There were also changes in the radiation variables, with the added top (top-of-atmosphere) net thermal radiation (TTR) taking the fourth slot and being completed by other ones in the top 30 variables: top net solar radiation (TSR), surface latent heat flux (SLHF), surface net thermal radiation (STR), surface solar radiation downwards (SSRD), and surface net solar radiation (SSR).These variables are likely highly correlated, and the selection could be reduced.It can also be noted that these variables are still often considered in terms of gradient (using S 1 ), even though cloud cover variables were made available (see discussion on radiation variables in Section 4).As for cloud cover variables, different ones were selected in the top 30: low cloud cover (LCC) and cloud cover (CC) at 600, 1,000, and 1500 hPa.Although LCC was most often considered in terms of gradients, the absolute values of the other cloud cover variables were mostly selected.The importance of low-level PV also increased compared to ERA-I.Conversely, the geopotential

Water Resources Research
10.1029/2023WR035715 height was only selected at 500 hPa in the top 30 predictors, SLP is no longer among the best, and the presence of the divergence variables also decreased.
The optimal spatial domains are comparable with those selected for ERA-I, including the 2-m temperature extension to the south.Regarding the temporal windows, TCW is again mainly selected between 6 and 12 UTC, and RH at different times of the day.PV is often selected at the end of the day, along with W at 1,000 hPa, the surface latent heat flux (SLHF), and the 2-m temperature (T2m).The other variables are mainly selected during the daytime.

Accuracy of the Optimized Methods
To assess the relevance of the methods optimized in this work, their accuracy has been compared to the benchmark methods (Section 2.2). Figure 12 shows the CRPS skill score and the Brier skill scores for different thresholds (1 mm, 10 mm, 50 mm) using the simplest RM1 method (for the CP) as reference.The best optimization result per catchment was selected based on the VP score.The scores for the test period (TP) were then calculated from unseen data for these selected parameter sets.
The skill scores are shown for the first single variable selection from ERA-I (ERA-I GAS 1 × 1), and the full optimizations using ERA-I (ERA-I GAS 1 × 8, 1 × 12) or ERA5 (ERA5 GAS 1 × 12).One can see in Figure 12 that the selection of a single best variable (GAS 1 × 1) shows similar accuracy to the RM1 method.Obviously, the  3 for the variables abbreviations) from ERA5 for the 1 × 12 structure for the different catchments.The variables that were forced into the AM are marked with a rectangle.The colors represent the analogy criteria, and the size of the dots is proportional to the weight given to the predictor within the range [0.02, 0.2].Variables that were never selected with a weight equal to or larger than 0.05 are not represented.

Water Resources Research
10.1029/2023WR035715 skill of a single variable remains lower than that of more complex AMs.The other optimized methods (GAs 1 × 8 or 1 × 12) show a higher CRPSS than the benchmark methods.Thus, despite having a single level of analogy, they outperform complex stepwise AMs in terms of CRPS.Brier skill scores for the prediction of the precipitation occurrence (threshold of 1 mm) of the optimized methods show values similar to those of RM6 when ERA-I is used and some further improvements when ERA5 is used.Brier skill scores of the optimized methods show similar skill to the best benchmark methods for a threshold of 10 mm, but lower values for a threshold of 50 mm.This can result from either an underestimation or an overestimation of the prediction.The GAs optimized the methods by minimizing the CRPS only, and a combined objective function that also accounts for the Brier score, for example, could be used instead to improve other properties of the resulting AMs.The gain obtained by using ERA5 instead of ERA-I may be due to higher spatial and temporal resolutions or better variables (Horton, 2021).Some differences can be observed between the three splits (CP, VP, TP), also for the benchmark methods.However, there is no clear trend, and the distributions remain relatively close.These differences can have multiple origins: the presence of stronger precipitation events in some splits, inhomogeneities in the quality of the predictor variables, or just natural variability.The three splits of the method optimized with ERA5 provide very similar results, which can be due to its higher quality, the variables selected, or just luck.Anyway, the selection of the predictor variables and the analogy criteria by GAs, along with all other parameters, provides AMs that prove relevant and consistent among different periods.No overfitting from the GAs can be observed.
Rank histograms have been computed for some catchments.Figure 13 shows such a plot for the Ergolz catchment.
Other catchments show similar results, that is, that the prediction by AMs tends to be over-dispersive without presenting a clear bias.This observation is not specific to these results, but is a common behavior of AMs.
An additional experiment has been attempted by forcing the predictor variables (along with their vertical level and their time) and the analogy criteria and letting the GAs optimize the weights between these variables, along with the spatial domains.To this end, 26 of the most commonly selected ERA5 variables were provided to the optimizer, organized in a single level of analogy (1 × 26).The results are shown in Appendix C.This approach did not provide the best accuracy (not shown), which can be due to non-optimal choices made to homogenize the vertical levels or times of the day, for example, In addition, this approach is not computationally efficient, as it requires loading variables that barely play a role in the selection of analog situations.Therefore, we do not recommend using this strategy.
The predictions provided by the optimized AMs for the strongest precipitation event of the test period are illustrated for some catchments in Figure 14.While most of these heavy precipitation events were captured satisfactorily by the optimized AMs, few were underestimated in some catchments as shown in Figure 14 (bottom), where the two worst predictions are shown.

Discussion
The primary objective of this study was to assess the relevance of GAs in selecting input variables for AMs.The results demonstrated that GAs could identify pertinent predictors and analogy criteria.However, caution is due when extrapolating the use of these selected predictors to different contexts, as their applicability may not be universally optimal.In fact, the compilation of potential variables must be tailored to the specific requirements of the AM application.For example, in forecasting applications, only meteorological variables deemed reliably predicted should be included.In the context of climate impact studies, the selection is constrained by the limited  availability of meteorological variables compared to the extensive output provided by reanalysis and NWP models.Furthermore, it is crucial to exercise discretion in selecting variables that exhibit a causal relationship with the predictand of interest and avoid undesirable covariability.Essentially, adapting the pool of potential variables to the application at hand is fundamental for a robust use of the optimized AM.
Radiation variables were often selected as relevant predictors.When using ERA-I, SSRD is the second most selected variable, and STRD and STR are the fourth and fifth.When using ERA5, TTR is the fourth most important variable.As these variables were selected so often, they did provide useful information for precipitation prediction, but their role is not easily interpretable.Hereafter, we propose some hypotheses about the information that can potentially be retrieved from these variables.First, STR values for days with high precipitation values show positive anomalies, meaning that the long-wave radiation from the atmosphere towards the surface is anomalously high (vertical fluxes are positive downward).Thermal radiation emitted by clouds and the atmosphere contributes to the downward STR.It is possible that the thermal radiation flux towards the surface is increased due to a high concentration of water vapor in the lower atmosphere and/or the presence of low clouds (with a higher cloud base temperature).Therefore, it can be used as a proxy for the presence of low clouds.Low clouds can interact with the topography, and this interaction might not be reflected in the vertical motion in ERA-I due to the relatively coarse spatial resolution of the orography in the reanalysis.The information from STR would then compensate for missing local processes at some locations, which potentially have a better representation in ERA5.
Then, SSRD was selected as a relevant predictor, but with the analogy criteria comparing the gradients rather than the absolute values, meaning that the pattern of SSRD matters more than its values.Gradients in SSRD could be an indication of the presence of fronts or thunderstorm clouds.Finally, in ERA-5, TTR anomalies are selected.They can be a proxy for high cloud tops with lower temperatures and, therefore, might provide information on the cloud thickness.Further research is needed to explore these hypotheses.
The triplet S 0 , S 1 and S 2 dominate the selection of analogy criteria.The S 1 score originally developed by Teweles and Wobus (1954) to verify prognostic charts was then used because it penalizes forecasters who tend to be overly conservative by forecasting weak systems too often.The rationale behind this lies in the denominator, which is determined by the sum of the maximum gradients of either the forecast or the observation.Consequently, forecasting a weaker system incurs a greater penalty than forecasting a stronger one.However, it should be noted that this approach may lead to the opposite effect, as forecasters may find it safer to predict stronger systems with larger gradients, thereby inflating the denominator (Thompson & Carter, 1972).This can be transposed to the AM, where stronger gradients in Z from analog situations are preferred over weaker ones.
The S 0 and S 2 criteria share a key characteristic with S 1 by imposing heavier penalties on weaker fields.Consequently, the analog selection based on S 0 , S 1 , and S 2 exhibits asymmetry, favoring the selection of analog fields close to the target but tending towards greater rather than weaker values (see Appendix C).The inherent asymmetry of S 0 , S 1 , and S 2 proves advantageous for prediction.Optimal analog situations are skewed toward being slightly stronger than weaker.Considering that the CRPS is strongly influenced by heavy precipitation events, this suggests a hypothesis: given the potential underrepresentation of large precipitation events in the archive, AMs benefit from selecting stronger predictor fields, often associated with higher precipitation.This selection bias may function as a compensatory mechanism for the underrepresentation of intense precipitation events.These assumptions would need to be further investigated.

Conclusions
The objective of the work was to assess the ability of GAs to select the input variables of the analog method along with the analogy criteria.The experiment was successful, as the selected predictors provided better accuracy (in terms of CRPS that was used as the objective function for the optimizations) than the benchmark methods, without overfitting.In addition, most of the selected variables can be related to the meteorological processes involved in precipitation generation.For example, among the most selected variables are: the vertical velocity (W) at 700 hPa (along with other levels), the total column water (TCW), the convective available potential energy (CAPE), radiation variables, the potential vorticity (PV), the relative humidity (RH), cloud cover variables, wind components, the geopotential height, air temperature, and the divergence.

Water Resources Research
10.1029/2023WR035715 The selection of analogy criteria also proved fruitful, as there were clear trends toward a dominant criterion for a given variable.The unexpected result was the success of the criterion S 0 , inspired by the Teweles-Wobus criterion.This new S 0 turned out to be the most often selected analogy criterion for the characterization of Euclidean distances.Three analogy criteria were most often selected and all are derived from the Teweles-Wobus criterion; one is based on the raw point values, another on the gradients, and the third on the second derivative of the fields.All of them are normalized by the sum of the largest point (pair)wise values from the target or the candidate fields.This normalization makes the criteria asymmetrical, so that higher values are preferred to lower ones.These new criteria should be further investigated and could be used in classic AMs.
Another unexpected result is the preferred structure for analog methods.While most benchmark methods build on a stepwise selection of predictors with successive levels of analogy subsampling from the previous one by using different predictors, here, the GAs preferred a flatter structure, mainly with a single level of analogy, but more variables.The benchmark methods most often start with selecting candidate analogs using the geopotential height and then narrowing down the selection using vertical velocity or moisture variables.A primary difference with the benchmark methods is that the variables are standardized here and weights are used (and optimized) to combine them in a given level of analogy.These two elements make the combination of variables with different value ranges easier.However, it cannot be excluded that deeper structures can provide better results, but that GAs did not find these solutions.
Such optimization is computationally intensive.The new GPU-based computations brought notable time improvement, particularly for high-resolution data.Other approaches could be considered to decrease the computation time, such as a faster exploration of the data set using a smaller period for data pre-screening, or the division of the whole period into smaller batches.An alternative could be to reduce the number of days with small precipitation amounts, as they have a small impact on the CRPS, while weighting their contributions by using a weighted CRPS approach.
for the analogy criterion S 1 , with gradients pre-processed using CPUs only (counted in the total time).The other analogy criteria showed similar results.The task consisted of extracting analogs for 32 years using the other 31 years as archives for candidate situations within a 120-day temporal window.It makes a total of 43.5 ⋅ 10 6 field comparisons per predictor of the first level of analogy.
The experiment was carried out on the UBELIX cluster of the University of Bern, using the same node for the whole benchmark and processing on a single NVIDIA GeForce RTX 2080 graphics card.The CPU processingusing the linear algebra library Eigen 3 (Guennebaud et al., 2010)-was done on a single thread.Although AtmoSwing can parallelize the calculation of the analogy criteria on multiple CPU threads, it uses a single thread for this task when optimizing with GAs because it parallelizes the evaluation of the different individuals on multiple threads.With GPUs, it still assesses the individuals on multiple CPU threads, each of them being able to use a different GPU device to calculate the analogy criteria.It is thus parallelizing both on CPUs and GPUs.
The benchmark (Figure A1) shows that GPU computations are systematically faster than those on the CPU, and this difference increases with the number of grid points.The GPU computations were 13 times faster on average and up to 38 times faster (5.2 s instead of 3.3 min) when using 2,048 points.NWP model outputs and reanalyses show an increase in spatial resolution; therefore, the impact on the computation time will become increasingly important.When using CPU only, adding a predictor in the first level of analogy has a much higher impact on time than adding a second level of analogy.It is explained by the fact that it needs to process the analogy criteria for the whole archive for each predictor of the first level of analogy, while the second level has only a few candidate situations to assess.
Where p mut is the mutation probability, G m,r is the maximum number of generations (G) during which the magnitude of the research varies, and w is a threshold chosen to maintain a minimum search magnitude when G > G m,r .
Figure B1 shows the performance of these five mutation operators for different AM structures and the different catchments considered in Section 3.2.Overall, the chromosome of adaptive search radius has a success rate of 76.25% in calibration and 62.5% in validation, the multiscale mutation 7.5%, and 8.75% respectively, and the nonuniform mutation with its different options: (3) 11.25% and 10%, (4) 11.25% and 21.25%, and (5) 1.25% and 2.5% respectively.
Thus, it is quite clear that the chromosome of adaptive search radius obtains the best results, all the more so with more complex structures, that is, more predictor variables.Although its success rate decreases slightly in validation, it remains much larger than the other options.The non-uniform mutation shows notable variability of performance depending on its options.

Figure 1 .
Figure 1.Location of the 25 selected catchments in Switzerland along with the climatic regions (dashed lines) and the river network (source: SwissTopo, HADES).

Figure 2 .
Figure 2. Best single variable selected (ordinate; see Table3for the variables abbreviations) from ERA-I for the 25 catchments (abscissa).The colors represent the analogy criteria, and the size of the dots is proportional to the accuracy of the resulting method (the larger the dots, the better), within a range of 5% of the best result (those with lower accuracy are hidden).

Figure 3 .
Figure 3. Same as Figure 2 but for CFSR.

Figure 4 .
Figure 4. Frequency of the criteria selection for both reanalysis data sets.

Figure 5 .Figure 6 .
Figure 5. Map of the best variables for ERA-I for each catchment.

Figure 7 .
Figure 7. Same as Figure 6 for different AM structures with up to two levels of analogy and eight variables per level for two catchments in Switzerland.

Figure 8 .
Figure 8. Selected variables (see Table3the variables abbreviations) from ERA-I for the 1 × 8 and 1 × 12 structures for the different catchments.The colors represent the analogy criteria, and the size of the dots is proportional to the weight given to the predictor within the range [0.02, 0.2].Variables that were never selected with a weight equal to or larger than 0.05 are not represented.

Figure 9 .
Figure9.Statistics of the 30 most selected variables from ERA-I for the 1 × 8 and 1 × 12 structures for the different catchments (100 optimizations) along with the analogy criteria, the temporal window (30 = next day at 06 UTC; some radiation variables were considered at 15 UTC), and the spatial windows (longitudes and latitudes).The extent of Switzerland is shown in gray on the plots of the spatial windows.

Figure 10 .
Figure 10.Selected variables (see Table3for the variables abbreviations) from ERA5 for the 1 × 12 structure for the different catchments.The variables that were forced into the AM are marked with a rectangle.The colors represent the analogy criteria, and the size of the dots is proportional to the weight given to the predictor within the range [0.02, 0.2].Variables that were never selected with a weight equal to or larger than 0.05 are not represented.

Figure 11 .
Figure11.Statistics of the 30 most selected variables from ERA5 for the 1 × 12 structure for the different catchments (50 optimizations) along with the analogy criteria, the temporal window (30 = next day at 06 UTC), and the spatial windows (longitudes and latitudes).The extent of Switzerland is shown in gray on the plots of the spatial windows.

Figure 12 .
Figure12.Skill scores (CRPSS and Brier skill score) of the different benchmark and optimized methods on the calibration, validation, and test periods for the 25 catchments.The skill scores use the RM1 method on the CP as a reference.An LxP code represents the structures, with L being the number of levels of analogy and P being the number of predictors per level.

Figure 13 .
Figure13.Rank histogram for the prediction by the optimized AM for the Ergolz catchment on both the calibration and test periods.The frequencies were averaged over a boostrapping of 1,000 realizations to smooth out the effect of the random rank attribution of the zero precipitation cases.

Figure 14 .
Figure14.Illustration of predictions for the strongest precipitation event of the test period for several catchments.The predictions provided by the AMs are illustrated by their whole range as well as some quantiles often considered in operational forecasting.

Figure B1 .
Figure B1.Performance of the five mutation operators (Section 2.3) for different AM structures and the different catchments considered in Section 3.2.The values represent the number of optimizations for one mutation operator that resulted in the best performing AM.Results are shown for both calibration and validation.When multiple operators obtain the same accuracy, they all get a point.

Table 1
Characteristics of the 25 Selected Catchments in Switzerland

Table 2
Some Analog Methods Listed by Increasing Complexity

Table 3
Selected Variables for ERA-I, CFSR, and ERA5 for Different Types of Vertical Levels Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023WR035715by Schweizerische Akademie Der, Wiley Online Library on [13/04/2024].See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions)on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 5. S 0 : as with S 2 , this new criterion is derived from S 1 and is processed on the raw grid values.It differs from the MD mainly in that it is normalized by the sum of the maximum values instead of the number of points: Note.PL = pressure levels, PT = pot.temp.levels, PV = pot.vorticity levels, SC = single level, surface, or total column.a At 10 m. b At 2 m. c At mean sea level.d Moisture and temperature variable.Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023WR035715by Schweizerische Akademie Der, Wiley Online Library on [13/04/2024].See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions)on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023WR035715by Schweizerische Akademie Der, Wiley Online Library on [13/04/2024].See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions)on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License