Code for "Propagation of precipitation trends into the terrestrial water cycle"
Creators: name: Stephan, Ruth
Publisher: name: Zenodo
Year: 2025

description: The here provided code was used for the study "Propagation of precipitation trends into the terrestrial water cycle". We used Python 3.11.5 and the environment py_rues that is described in detail with PythonEnvironment_py_rues.txt. Please be aware that you need to adapt the paths.

01: Preprocessing with 

01a_Preprocessing_MaskOceans.py imports the mask for terrestrial regions and all the nc-files from the different models and available variables for r1i1p1f1 and ssp585. This script saves the data in data/CMIP6/model_individuals/

01b_Preprocessing_CombineVariables.py imports the masked variables files from all models and combines them in one nc-file. Furthermore the Residual (P-ET-R) is calculated saved as data_var called PPETR. The _combined.nc file is stored per model in data/CMIP6/model_individuals/.

01c_Preprocessing_EnsembleMean.py imports the data from all models that provide tas, pr, evspsbl, mrro, mrso and calculates the mean to establish one nc-file as CMIP6 Mean. Thus this script calculated P mean, which is used by other scripts to identify whether the CMIP6 Mean shows P increase or decrease in the for future. Furthermore, this script computes the masks YearlyPrecip_Mask.nc (only gridcells with P >= 200mm/yr based on the longterm mean between 2000-2100), YearlyPrecipChange_Mask.nc (only gridcells where P change between 2000-2025 and 2075-2100 is greater than 5 %), and SeasonsPrecip_Mask.nc (only gridcells where seasonal P change between 2000-2025 and 2075-2100 is greater than 5 %), all of which are used later to mask out gridcells that do not match those criterias.

02: Analyse Ratios with

02a_AnalyseRatios_ModelIndividuals.py calculates for each model the change of the model variables between the time period 2075-2100 and 2000-2025 and stores those changes in a nc-file called delta_ds. The change is first calculated as a difference and then as a ratio. Later in the analysis, we only use the ratios, since they provide more information. For the changes the function calc_vardiff (for the annual mean) and calc_vardiff_seasons (for the seasons) are used. The ratios are stored in the delta_ds as ratio_ETP (∆ET/∆P), ratio_RP (∆Q/∆P), and ratio_PETRP (∆RES/∆P). The changes are first masked by different criterias (see 01c) and then mapped. In addition this script identifies the ratio closest to 1 per gridcell (highest similarity) among ∆ET/∆P, ∆Q/∆P and ∆RES/∆P. This highest ratios are stored as data_var called index_arr with 1 meaning ∆ET/∆P is closest to 1, 2 meaning ∆Q/∆P is closest to 1 and 3 meaning ∆RES/∆P is closest to 1. In a next step, the identified ratios are masked with Pmean in a way that the index shows whether the gridcell shows P increase or P decrease in the future. This is done for the annual mean (delta.nc) and for all seasons (sdelta.nc). The final files for each model is stored with the corresponding model name and as a list of nc-files called all_models_delta.pkl (annual mean) and all_models_sdelta.pkl (seasons).

02b_AnalyseRatios_AcrossModels.py imports the list of delta.nc and sdelta.nc files created by 02a and evaluates the model agreement. Different agreements are plotted. Agreement on NaN values, P surplus, P deficit, highest similarities. Those agreements are then the result of the CMIP6 Ensemble. For many maps, gridcells are blacked out when less than 50% of the models agree. In addition, the script identifies the model that shows the highest similarity with the CMIP6 Ensemble. This scripts imports many user-functions from Functiony.py to fasten the computation.

02c_Ratios_ModelMean.py does the same analysis for the CMIP6 Mean as described in 02a.

03: Analyse Correlations with 

03a1_GetTrends_ModelIndividuals.py calculates the longterm trends of all variables per model and therefore uses the corresponding _combined.nc files. The trends are calculated with LOWESS-filtering using a moving window that captures 25% of the data. The trends are calculated per gridcell and for the annual mean and for all seasons. The trend.nc (annual mean) and strend.nc (seasons) is stored per model and as a list of all models (all_models_trend.pkl and all_models_strend.pkl). Despite user-functions imported from Functions.py, the calculation of the trends is computational and took some time.

03a2_AnalyseCorr_ModelIndividuals.py analysis the trends of each model and therefore imports the lists of nc-files created by 03a1. Per gridcell, the trends of P and ET, Q , RES are correlated with each other to identify how similar they behave. These correlations are stored as corr_PET (Cor(P,ET)), corr_PR (Cor(P,Q)) and corr_PPETR (Cor(P,RES)). The highest correlation out of those three is computed and stores as data_var index_ arr with 1 meaning Cor(P, ET) is highest, 2 meaning Cor(P, Q) is highest and 3 meaning Cor(P, RES) is highest. In a next step, the identified ratios are masked with Pmean in a way that the index shows whether the gridcell shows P increase or P decrease in the future. This is done for the annual mean (corr.nc) and for all seasons (scorr.nc). The final files for each model is stored with the corresponding model name and as a list of nc-files called all_models_corr.pkl (annual mean) and all_models_scorr.pkl (seasons).

03b_AnalyseCorr_AcrossModels.py imports the list of corr.nc and scorr.nc files created by 03a and evaluates the model agreement. Different agreements are plotted. Agreement on NaN values, P surplus, P deficit, highest similarities. Those agreements are then the result of the CMIP6 Ensemble. For many maps, gridcells are blacked out when less than 50% of the models agree or if the highest correlation is lower than 0.5. In addition, the script identifies the model that shows the highest similarity with the CMIP6 Ensemble. This scripts imports many user-functions from Functiony.py to fasten the computation.

03c_Ratios_ModelMean.py does the same analysis for the CMIP6 Mean as described in 03a1 and 03a2.

03a_GetAndAnalyseTrends_ModelIndividuals_withmrso.py, 03b_AnalyseSimilarTrends_ModelIndividuals_withmrso.py, 03a_AnalyseSimilarTrends_ModelMean_withmrso.py, Functions_withmrso: These scripts do the same analysis as described in 03: Analyse Correlations, but with soil moisture (mrso) instead of P-ET-Q (RES).

04 Attribution Analysis Correlations with 

4aa_PrepareData_AttributionAnalysis_Corr.py imports the masks by the CMIP6 Mean (P mean and P changes) for the annual mean and the seasons, as well as elevation data from NOAA. Then, this script computes for 6 ESMs the response and predictor data. For the responses Cor(P, ET), Cor(P, Q) and Cor(P, RES) the trend.nc files for each model are used, while for the predictors the combined.nc files for each are used. We calculated all predictors for each of the 6 ESMs shown in Figure 1, Step 3. The predictor “aridity” is based on radiation/pr. For few gridcells this resulted in extreme ratios that are excluded from the predictor data by only considering the ratios between the 5th and 95th percentile. Further the data is masked by the imported masked for P. In the next step, the script combined all predictor data in one dataframe and then excludes the columns for the three responses, the model name, lat, and lon to apply dimensional reduction for the random forest model. To Reduce collinearity by predictor selection, the script runs a cluster analysis to identify high correlated predictors and select one representative that shows the highest correlation sum to each predictor within the same cluster. The script stores the selected predictor, the corresponding cluster and the time scale, since this is done for the annual mean, across seasons and per season stored as Corr_Predictor_Selection.pkl. The predictor data for the annual mean is stored as Corr_PredictorData_yearly.pkl and for the seasons as Corr_PredictorData_seasons.pkl.

4a_AttributionAnalysis_Corr.py imports the response and predictor data created by 04aa. The script runs for 6 EMS separately a random forest with the same set of predictors and stores the OOB-Scores the HHI index and the SHAP-values for each predictor as *RF_output*. This is done for the annual mean, for each season in separate model, and across all seasons in one model. The calculation of the SHAP-values is computational and takes much time. Therefore, the following scripts are the same as 04a, but implemented for only 2 out of the 6 ESMs: 04a1, 04a2, 04a3.

4b_AttributionAnalysis_Evaluation_Corr.py imports the RF_output created by 04a and analysis the random forest models fit separately for 6 ESMs in summary plots. Therefore, the script is split in three parts The first focusing on the yearly time scale, the 2nd focusing on each season separately and the 3rd focusing on all seasons together. Each part loops over the different responses and created barplots and boxplot for the SHAP-values and OOB-Scores and HHI as inset plots.

04 Attribution Analysis Ratios with 

4cc_PrepareData_AttributionAnalysis_Ratio.py imports the masks by the CMIP6 Mean (P mean and P changes) for the annual mean and the seasons, as well as elevation data from NOAA. Then, this script computes for 6 ESMs the response and predictor data. For the responses ∆ET/∆P, ∆Q/∆P and ∆RES/∆P and the predictors the combined.nc files for each model are used. We calculated all predictors for each of the 6 ESMs shown in Figure 1, Step 3. The predictor “aridity” is based on radiation/pr. For few gridcells this resulted in extreme ratios that are excluded from the predictor data by only considering the ratios between the 5th and 95th percentile. Furthermore, the data is masked by the imported masked for P. In the next step, the script combined all predictor data in one dataframe and then excludes the columns for the three responses, the model name, lat, and lon to apply dimensional reduction for the random forest model. To Reduce collinearity by predictor selection, the script runs a cluster analysis to identify high correlated predictors and select one representative that shows the highest correlation sum to each predictor within the same cluster. The script stores the selected predictor, the corresponding cluster and the time scale, since this is done for the annual mean, across seasons and per season stored as Ratio_Predictor_Selection.pkl. The predictor data for the annual mean is stored as Ratio_PredictorData_yearly.pkl and for the seasons as Ratio_PredictorData_seasons.pkl.

4c_AttributionAnalysis_ Ratio.py imports the response and predictor data created by 04cc. The script runs for 6 EMS separately a random forest with the same set of predictors and stores the OOB-Scores the HHI index and the SHAP-values for each predictor as *RF_output*. This is done for the annual mean, for each season in separate model, and across all seasons in one model. The calculation of the SHAP-values is computational and takes much time. Therefore the following scripts are the same as 04a, but implemented for only 2 out of the 6 ESMs: 04c1, 04c2, 04c3.

4d_AttributionAnalysis_Evaluation_Ratio.py imports the RF_output created by 04c and analysis the random forest models fit separately for 6 ESMs in summary plots. Therefore, the script is split in three parts The first focusing on the yearly time scale, the 2nd focusing on each season separately and the 3rd focusing on all seasons together. Each part loops over the different responses and created barplots and boxplot for the SHAP-values and OOB-Scores (and HHI as inset plots).