diff --git a/data/fieldlist_CMIP.jsonc b/data/fieldlist_CMIP.jsonc index 19417846c..8f29b95cb 100644 --- a/data/fieldlist_CMIP.jsonc +++ b/data/fieldlist_CMIP.jsonc @@ -451,6 +451,11 @@ "realm": "ocean", "units": "m", "ndim": 3 + }, + "mrro": { + "standard_name": "runoff_flux", + "units": "kg m-2 s-1", + "ndim": 3 } }, "env_vars" : { diff --git a/diagnostics/runoff_sensitivities/doc/Figure1.png b/diagnostics/runoff_sensitivities/doc/Figure1.png new file mode 100644 index 000000000..02fce09ef Binary files /dev/null and b/diagnostics/runoff_sensitivities/doc/Figure1.png differ diff --git a/diagnostics/runoff_sensitivities/doc/Figure2.png b/diagnostics/runoff_sensitivities/doc/Figure2.png new file mode 100644 index 000000000..e2254c8a9 Binary files /dev/null and b/diagnostics/runoff_sensitivities/doc/Figure2.png differ diff --git a/diagnostics/runoff_sensitivities/doc/runoff_sensitivities.rst b/diagnostics/runoff_sensitivities/doc/runoff_sensitivities.rst new file mode 100644 index 000000000..d5d93e038 --- /dev/null +++ b/diagnostics/runoff_sensitivities/doc/runoff_sensitivities.rst @@ -0,0 +1,161 @@ +Runoff Sensitivities Diagnostic Documentation +============================================= + +Last update: 04/07/2025 + +Introduction to runoff sensitivity +---------------------------------- + +Runoff projections are essential for future water security. While Earth System Models (ESMs) are now being used to inform water resource assessments, the substantial model uncertainty has undermined the reliability of ESM projections. The primary source of the runoff projection uncertainty is the model's climate forcing of the land surface, namely inter-model differences in precipitation (P) and surface air temperature (T) trends. + +However, runoff projections are generally more uncertain than either P and T projections, implying that additional uncertainties arise at the terrestrial interface (:ref:`Lehner et al., 2019 `; :ref:`Wang et al., 2022 `). These additional uncertainties stem from differences in the sensitivity of runoff to climate change - how runoff responds to a given set of climate forcings. Previous studies have measured this using runoff sensitivity, a metric that quantifies runoff (Q) changes in response to changes in precipitation (P sensitivity; δQ/δP, %/%) and temperature (T sensitivity; δQ/δT, %/K) (:ref:`Tang and Lettenmaier, 2012 `; :ref:`Hoerling et al., 2019 `; :ref:`Lehner et al., 2019 `; :ref:`Milly and Dunne, 2020 `). + +Uncertainty in this runoff sensitivity, inherent to each ESM, can contribute to the projection uncertainty as much as the uncertainty in meteorological forcings (Kim et al. 2025, *in preparation*). In addition, the runoff sensitivity is often biased in ESMs, such that ESMs tend to underestimate future runoff declines in several global basins (:ref:`Zhang et al., 2023 `; :ref:`Douville et al., 2024 `; Kim et al. 2025, *in preparation*). Therefore, reducing the runoff sensitivity bias is important for producing more reliable projections of terrestrial climate change. + +Calculation of runoff sensitivity +--------------------------------- + +In this diagnostics package, we quantify the runoff sensitivity of a target model using multiple linear regression. The historical timeseries (1945-2014) are averaged for each water year (Oct.-Sep.) and 131 global river basins. +Then, the annual timeseries is smoothed with a 5-year moving window to minimize the storage effect of runoff from previous years. Eventually, the runoff sensitivity is calculated as: + +.. math:: + + \delta{Q} = {\alpha}\delta{P} + {\beta}\delta{T} + {\gamma}\delta{P}\delta{T} + +where α and β are estimates of P and T sensitivity respectively, δ represents 5-year averaged anomalies relative to long-term mean, δQ and δP are percent anomalies in runoff and precipitation, δT is temperature anomalies, δPδT is an interaction term. Note that the interaction term is usually negligible and does not affect the values of α and β in CMIP5/6 models. + +The calculated runoff sensitivity for a target model is then compared to a pre-calculated observational estimate and a set of CMIP5/6 models. For the observational estimate, we used the machine learning-based global runoff reanalysis (GRUN-ENSEMBLE; Ghiggi et al. 2021). We used 100 realizations of GRUN-ENSEMBLES to sample observational uncertainty. For CMIP models, we used historical simulation from 1945 to 2014. For CMIP5, we extended the timeseries to 2014 using the RCP4.5 scenario. The diagnostic page shows the ranges from observational uncertainty, inter-model spread, and uncertainty of regression coefficients for specific river basins. + +Version & Contact info +---------------------- + +- Version/revision information: version 2 (04/07/2025) +- PI: Flavio Lehner, Cornell University, flavio.lehner@cornell.edu +- Developer/point of contact: Hanjun Kim, Cornell University, hk764@cornell.edu +- Other contributors: Andy wood, David Lawrence, Katie Dagon, Sean Swenson, Nathan Collier + +Open source copyright agreement +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +The MDTF framework is distributed under the LGPLv3 license (see LICENSE.txt). + +Functionality +------------- + +The main driver code (runoff_sensitivities_diag.py) includes all functions, codes for calculations and figures. + +0) Functions are pre-defined for further analysis. +1) The code will load data and calculate averages for each water year and 131 global river basins (Basin Mask source: GRDC Major River Basins of the World). +2) Water budget closure (precipitation - evapotranspiration == runoff) in a long-term average sense is checked. +3) Using the pre-defined function "runoff_sens_reg", runoff sensitivities are calculated. +4) Calculated runoff sensitivities for target models are saved as .nc files. +5) The diagnostic figures will be plotted with the pre-calculated OBS and CMIP data. + + +Required programming language and libraries +------------------------------------------- + +Python 3 is used to calculate and draw the figures. + +- All libraries used in this diagnostic are available in MDTF conda environment "_MDTF_python3_base". +- Used libraries: "scipy", "numpy", "matplotlib", "netCDF4", "cartopy", "sklearn" +- To deal with the shapefile, "cartopy.io.shapereader" and "matplotlib.path" are used. +- For multi-linear regression, "sklearn.linear_model" is used. + +**Caution**: In Oct. 2023, the diagnostic did not work after an update in "pydantic" library. +Below commands for downgrading "pydantic" solved the problem for us. + +.. code-block:: restructuredtext + + conda activate _MDTF_base + conda install -c conda-forge pydantic==1.10.9 + + +Required model output variables +------------------------------- + +The monthly mean output from historical simulations, including period 1945-2014, are needed. +(Model outputs are assumed to have the same naming convention as CMOR-ized CMIP output.) + +Target variables: + - ``tas`` (surface air temperature, K), [time, lat, lon] + - ``pr`` (precipitation, kg m-2 s-1), [time, lat, lon] + - ``hfls`` (latent heat flux, W m-2), [time, lat, lon] + - ``mrro`` (runoff, kg m-2 s-1), [time, lat, lon] + +Lon-lat grids for all variables have to be the same. In CMIP, there are some models in which grids are slightly different between land and atmospheric variables. Checking and interpolation to a common grid is recommended. + + +More about this diagnostic +-------------------------- +Runoff sensitivity biases +^^^^^^^^^^^^^^^^^^^^^^^^^ +.. _my-figure-tag: + +.. figure:: Figure1.png + :align: center + :width: 100 % + + Fig. 1. CMIP6 multi-model median (MMM) (a) P sensitivity and (b) T sensitivity. Observational estimate of (c) P sensitivity and (d) T sensitivity. (e,f) The CMIP6 MMM sensitivity biases. Note that all sensitivities are calculated for the historical period from 1947 to 2017. + +The ESM's P sensitivity is positive across all basins, ranging from 1%/% to 3.5%/% with consistent values among ESMs (Fig. 1a). T sensitivity is negative for most basins and models in the northern mid-latitudes, while in the tropics and Southern Hemisphere its sign is varied and not consistent across models (Fig. 1b). A negative T sensitivity indicates more evapotranspiration per amount of warming and thus lower runoff, triggered by greater energy input at the surface. Additional factors, such as the phase shift from snow to rain and increased water use due to vegetation greening, can further contribute to negative T sensitivity. On the other hand, factors like enhanced stomatal closure, increased extreme precipitation events, and glacier melt can contribute to relatively more positive T sensitivity. These competing effects and their imperfect representations in models result in a wide range of T sensitivities across models. + +The observed P and T sensitivities generally range from 1%/% to 3.5%/% and -30%/K to slightly positive, respectively (Fig. 1c,d). CMIP6 models show significantly different P sensitivities in 106 out of 131 basins, overestimating the sensitivities in 102 of those 106 basins (Fig. 1e). On the other hand, CMIP6 models’ T sensitivities differ significantly from the observations in 90 out of 131 basins, underestimating the negative sensitivities in 75 of those 90 basins (Fig. 1f). **Overall, climate models tend to exhibit more positive P sensitivity and less negative T sensitivity compared to the observational estimates.** + + +Impact on the future runoff projections +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +.. _my-figure-tag: + +.. figure:: Figure2.png + :align: center + :width: 100 % + + Fig. 2. (a) Algorithm for testing the statistical significance of an observational constraint based on the models' runoff sensitivity biases. Observational constraining effect on MMM runoff projection in (b) CMIP6 SSP2-4.5 scenario and (c) CMIP5 RCP4.6 scenario (% change, 2030-2070 vs. 1947-2017). MMM runoff projection (d) before and (e) after correction by the observational constraint. + +How does such runoff sensitivity bias affect future projections? This motivates an observational constraint: the use of the observed runoff sensitivity to constrain the ESM runoff projections. We do this by convolving the observational sensitivity with the ESM's projection in T and P: + +.. math:: + + \Delta{Q}_{obs} = {\alpha}_{obs}\Delta{P}_{ESM} + {\beta}_{obs}\Delta{T}_{ESM} + +where {\alpha}_{obs} and {\beta}_{obs} are observational estimate of P and T sensitivity, \Delta{P}_{ESM} and \Delta{T}_{ESM} represent the ESM projection in P and T, and \Delta{Q}_{obs} is the observationally-constrained runoff projection. + +The difference between this observationally-constrained projection and the unconstrained ESM projection is the observational constraining effect. We developed a systematic approach to test the significance and robustness of the constraining effect (Fig.2a; see details in Kim et al., 2025, *in preparation*). The observationally-constrained projections indicate a drier future than the unconstrained projections (Fig. 2b,c). The downward correction is significant for 41 of 131 global river basins. Specifically, in CMIP6, the modest runoff reductions in regions like the Amazon, Eastern Europe, and Australia are intensified (Fig. 1d,e). Moreover, the originally projected runoff increases over North Africa and Southeast Asia become weaker or are even reversed. The downward correction mainly arises from the T sensitivity bias (Fig. 2b,c). **In other words, the ESMs underestimate the future runoff decline mainly because negative T sensitivity is too weak compared to observations.** This systematic underestimation of runoff decline is consistent with previous studies that have adjusted runoff projections downward using other observational datasets and statistical methods (:ref:`Zhang et al., 2023 `; :ref:`Douville et al., 2024 `). + +The causes of the models’ runoff sensitivity biases are not identified here and warrant further investigation. A sensitivity analysis of mean state climate variables for each basin shows that P sensitivity generally exhibits an inverse relationship with mean runoff and runoff ratio (Q/P). In contrast, temperature sensitivity generally displays no systematic inter-model correlation with any mean state variables (Kim et al. 2025, *in preparation*). Depending on the basins, a significant correlation to certain mean state variables exists, but it is difficult to identify a more global culprit. **This suggests that traditional modeling approaches, which focus on improving mean state biases, may improve predictions related to P sensitivity but are unlikely to resolve the more critical T sensitivity biases.** These challenges motivated the development of the metrics package, which can be used and expanded by the broader scientific community to diagnose model biases and aid model development. + + +References +---------- +.. _ref-Lehner: + +1. Lehner et al. (2019): The potential to reduce uncertainty in regional runoff projections from climate models. *Nature Climate Change*, **9** (12), 926-933, `doi:10.1038/s41558-019-0639-x `__. + +.. _ref-Wang: + +2. Wang et al. (2022): Future Changes in Global Runoff and Runoff Coefficient From CMIP6 Multi‐Model Simulation Under SSP1‐2.6 and SSP5‐8.5 Scenarios. *Earth’s Future*, **10** (12), e2022EF002910, `doi:10.1029/2022EF002910 `__. + +.. _ref-Tang: + +3. Tang, Q., & Lettenmaier, D. P. (2012): 21st century runoff sensitivities of major global river basins. *Geophysical Research Letters*, **39** (6), 2011GL050834, `doi:10.1029/2011GL050834 `__. + +.. _ref-Hoerling: + +4. Hoerling et al. (2019): Causes for the Century-Long Decline in Colorado River Flow. *Journal of Climate*, **32** (23), 8181–8203, `doi:10.1175/JCLI-D-19-0207.1 `__. + +.. _ref-Zhang: + +5. Zhang et al. (2023): Future global streamflow declines are probably more severe than previously estimated. *Nature Water*, **1** (3), 261–271, `doi:10.1038/s44221-023-00030-7 `__. + +.. _ref-Milly: + +6. Milly, P. C. D., & Dunne, K. A. (2020): Colorado River flow dwindles as warming-driven loss of reflective snow energizes evaporation. *Science*, **367** (6483), 1252–1255, `doi:10.1126/science.aay9187 `__. + +.. _ref-Douville: + +7. Douville, H. (2024): Observational Constraints on Basin-Scale Runoff: A Request for Both Improved ESMs and Streamflow Reconstructions. * Geophysical Research Letters*, **51** (13), e2024GL108824, `doi:10.1029/2024GL108824 `__. + +.. _ref-Ghiggi: + +8. Ghiggi et al. (2021): G‐RUN ENSEMBLE: A multi‐forcing observation‐based global runoff reanalysis. *Water Resources Research*, **57** (5), e2020WR028787, `doi:10.1029/2020WR028787 `__. diff --git a/diagnostics/runoff_sensitivities/preprocessing_code_at_local_cluster.py b/diagnostics/runoff_sensitivities/preprocessing_code_at_local_cluster.py new file mode 100644 index 000000000..4edffc757 --- /dev/null +++ b/diagnostics/runoff_sensitivities/preprocessing_code_at_local_cluster.py @@ -0,0 +1,708 @@ +#################################################################### +## import libraries +import os +import numpy as np +import cartopy.io.shapereader as shpreader +import netCDF4 as nc +import matplotlib.pyplot as plt +import matplotlib.patches as patches +import matplotlib.colors as colors +from sklearn.linear_model import LinearRegression +import scipy.stats as stats +import scipy.io as sio + +#################################################################### +## Define functions +## calculate moving average! (smoothing) +def moving_avg(x_1d,n_mw): + nt=len(x_1d) + weight = np.ones((n_mw))/n_mw + n_delete=int((n_mw)/2); + smoothed=np.convolve(x_1d,weight,mode='same') + if n_delete != 0: + smoothed[0:n_delete]=np.nan + smoothed[-n_delete:nt]=np.nan + return smoothed + +#################################################################### +## file path +dpath="C:/runoff_mrb_grdc3/" +figpath="C:/Users/reapr/내 드라이브/Research/runoff_project/MDTF/ver241202/runoff_sensitivities"; +os.makedirs(figpath, exist_ok=True) + +## load basin info +basin_path=dpath + 'basins_pp.mat' +matfile=sio.loadmat(basin_path) +basins=matfile['basins'] +basin_names_raw=matfile['basin_names'][0] +nbasin=len(basins)-4 +for b in range(nbasin): + point=[(x,y) for x,y in zip(list(basins[b][0][2][0]),list(basins[b][0][3][0]))] + if b==0: + basin_points=[point[0::10]] + basin_names=[basin_names_raw[0][0]] + else: + basin_points.extend([point[0::10]]) + basin_names.extend([basin_names_raw[b][0]]) + +## function for runoff sensitivity calculation +def runoff_sens_reg(r, p, t, mw=1, alpha=0.1): + nan_val=np.mean(np.isnan(r))+np.mean(np.isnan(p))+np.mean(np.isnan(t)) + if nan_val == 0: + # Ensure input vectors are column vectors + if r.ndim == 1: + r = r[:, np.newaxis] + if p.ndim == 1: + p = p[:, np.newaxis] + if t.ndim == 1: + t = t[:, np.newaxis] + # Create the regression matrix and do normalization + Xraw = np.column_stack((p, t, p * t)) + Xstd = np.nanstd(Xraw,axis=0) + rstd = np.nanstd(r,axis=0) + Xnorm = (Xraw - np.nanmean(Xraw,axis=0)) / Xstd + rnorm = (r - np.nanmean(r,axis=0)) / rstd + # Perform linear regression + model = LinearRegression() + model.fit(Xnorm, rnorm) + # Get regression coefficients + a1 = np.squeeze(model.coef_[0][0] * rstd / Xstd[0]) + b1 = np.squeeze(model.coef_[0][1] * rstd / Xstd[1]) + c1 = np.squeeze(model.coef_[0][2] * rstd / Xstd[2]) + # Calculate the standard errors for coefficients + y_pred = a1*Xraw[:,0] + b1*Xraw[:,1] + c1*Xraw[:,2] + residuals = r - y_pred[:,np.newaxis] + n, k = Xraw.shape + # autocorr=get_autocorr(r,1); + # dof2 = (n-k)*(1-autocorr)/(1+autocorr) + dof=np.ceil((n-k)/mw) + mse = np.sum(residuals ** 2) / (dof) + var_b = mse * np.linalg.pinv(np.dot(Xraw.T, Xraw)) + se_a1 = np.sqrt(var_b[0, 0]) + se_b1 = np.sqrt(var_b[1, 1]) + se_c1 = np.sqrt(var_b[2, 2]) + # Calculate the t-statistic for a given confidence level (e.g., alpha = 0.05) + t_critical = stats.t.ppf(1 - alpha / 2, dof) + # Calculate the confidence intervals + a2 = np.array([a1 - t_critical * se_a1, a1 + t_critical * se_a1]) + b2 = np.array([b1 - t_critical * se_b1, b1 + t_critical * se_b1]) + c2 = np.array([c1 - t_critical * se_c1, c1 + t_critical * se_c1]) + # R-squared value + corr, _ = stats.pearsonr(np.squeeze(r), np.squeeze(y_pred)) + r2 = corr**2 + return a1, a2, b1, b2, c1, c2, r2 + else: + a1=np.nan; a2=(np.nan,np.nan); b1=np.nan; b2=(np.nan,np.nan); c1=np.nan; c2=(np.nan,np.nan); r2=np.nan; + return a1, a2, b1, b2, c1, c2, r2 + + +###################################################################### +# start and end year for sensitivity calculation +syr_sens_cri=1945 +eyr_sens_cri=2014 + +###################################################################### +# target moving windows +n_mw=[1,3,5,9] +tmw=2 # mw=5 year is target moving window + +##################### +## CMIP6 ## +##################### +############################## +## HIST6 + SSP245 ## +############################## +## load variables +matfile = sio.loadmat(dpath + "tasb_wy_c2f6_cmip6.mat") +tasb_wy_c2f6 = matfile["tasb_wy_c2f6"][:,:,:] +matfile = sio.loadmat(dpath + "prb_wy_c2f6_cmip6.mat") +prb_wy_c2f6 = matfile["prb_wy_c2f6"][:,:,:] +matfile = sio.loadmat(dpath + "evspsblb_wy_c2f6_cmip6.mat") +evspsblb_wy_c2f6 = matfile["evspsblb_wy_c2f6"][:,:,:] +matfile = sio.loadmat(dpath + "mrrob_wy_c2f6_cmip6.mat") +mrrob_wy_c2f6 = matfile["mrrob_wy_c2f6"][:,:,:] +model_cmip6 = matfile["model_cmip6"][:][0] +syr_c2f6 = matfile["syr_c2f6"][0] +eyr_c2f6 = matfile["eyr_c2f6"][0] +nmodel_cmip6 = len(model_cmip6) +matfile = sio.loadmat(dpath + "prb_grun_1903_2012_intp.mat") +area_basins = np.squeeze(matfile["area_basins"]) + +################################################################ +# water budget closure for historical period (picontrol is similar) +# use the model with all variables we need +hval_available=np.isnan(np.mean(mrrob_wy_c2f6[:,:,3],axis=1))|np.isnan(np.mean(prb_wy_c2f6[:,:,3],axis=1))|np.isnan(np.mean(evspsblb_wy_c2f6[:,:,3],axis=1))|np.isnan(np.mean(tasb_wy_c2f6[:,:,3],axis=1)); +nanind_model_c2f6=np.nonzero(hval_available) +tasb_wy_c2f6[nanind_model_c2f6,:,:]=np.nan +prb_wy_c2f6[nanind_model_c2f6,:,:]=np.nan +evspsblb_wy_c2f6[nanind_model_c2f6,:,:]=np.nan +mrrob_wy_c2f6[nanind_model_c2f6,:,:]=np.nan + +# check the errors +pmeb_model = np.mean(prb_wy_c2f6[:,0:165,0:nbasin] - evspsblb_wy_c2f6[:,0:165,0:nbasin], axis=1) +error_val_c2f6 = (pmeb_model - np.mean(mrrob_wy_c2f6[:,0:165,0:nbasin], axis=1)) / np.mean(mrrob_wy_c2f6[:,0:165,0:nbasin], axis=1) +error_val_c2f6[nanind_model_c2f6,:]=np.nan +pval = 0.1 +closed_fraction_c2f6=np.nansum( np.abs(error_val_c2f6)0.001)[0] +nmodel_c2f6_closed=len(mind_c2f6) +nmodel_c2f6_available=len(mind_c2f6_available) +print(f"c2f6: {nmodel_c2f6_closed}/{nmodel_c2f6_available}") +print(f"negative values in c2f6 closed: {len(np.intersect1d(mind_negative_c2f6,mind_c2f6))}") + +# incorporate HIST/SSP model info +nanind_cmip6=nanind_c2f6 +mind_cmip6_available=mind_c2f6_available +mind_cmip6=mind_c2f6 +mind_cmip6=np.setdiff1d(mind_cmip6,69) # 69-TaiESM1's IA R2 is very low. P cannot explain R variability among 58/87 basins +mind_cmip6=np.setdiff1d(mind_cmip6,34) # 34-FGOALS-g3's IA R2 is very low. P cannot explain R variability among 29/87 basins, and R2 values are in generall too weak +mind_cmip6=np.setdiff1d(mind_cmip6,57) # 57-MIROC-ES2L has starnge value for northern sierras - due to low resolution, coastal regions have strange value +nmodel_cmip6_available=len(mind_cmip6_available) +nmodel_cmip6_closed=len(mind_cmip6) +print(f"HIST6/SSP245: {nmodel_cmip6_closed}/{nmodel_cmip6_available}") +closed_fraction_cmip6=closed_fraction_c2f6 +closed_area_fraction_cmip6=closed_area_fraction_c2f6 +closed_nmodel_cmip6=closed_nmodel_c2f6 +############################################ +temp2=negative_fraction_c2f6>0.1 +ind_c2f6_nm=np.nonzero(temp2)[0] +ind_c2f6_nb=np.nonzero(temp2)[1] + +for mind,bind in zip(ind_c2f6_nm,ind_c2f6_nb): + # display(f'c2f6 {model_cmip6[mind]}({mind}); {basin_names[bind]}({bind})') + # display(f'c2f6 mean Q: {np.nanmean(mrrob_wy_c2f6[mind,0:150,bind])}') + tasb_wy_c2f6[mind,:,bind]=np.nan + prb_wy_c2f6[mind,:,bind]=np.nan + evspsblb_wy_c2f6[mind,:,bind]=np.nan + mrrob_wy_c2f6[mind,:,bind]=np.nan + +# assign nan to unclosed models +tasb_wy_c2f6[nanind_cmip6,:,:]=np.nan +prb_wy_c2f6[nanind_cmip6,:,:]=np.nan +evspsblb_wy_c2f6[nanind_cmip6,:,:]=np.nan +mrrob_wy_c2f6[nanind_cmip6,:,:]=np.nan + + +############################################################################ +def duplicated_fraction(data): + ndiff_frac=(len(data)-len(set(data)))/len(data) + return ndiff_frac # No duplicates found + +s=0 +for m in range(nmodel_cmip6): + for b in range(nbasin): + val=duplicated_fraction(mrrob_wy_c2f6[m,:,b]) + if val>0.1: + tasb_wy_c2f6[m,:,b]=np.nan + prb_wy_c2f6[m,:,b]=np.nan + evspsblb_wy_c2f6[m,:,b]=np.nan + mrrob_wy_c2f6[m,:,b]=np.nan + s=s+1 + +############################################################################ +## calculate sensitivity +# moving average +nmw=len(n_mw) +tasb_wym_c2f6 = np.full(tasb_wy_c2f6.shape+(len(n_mw),),np.nan); +prb_wym_c2f6 = np.full(prb_wy_c2f6.shape+(len(n_mw),),np.nan); +evspsblb_wym_c2f6 = np.full(evspsblb_wy_c2f6.shape+(len(n_mw),),np.nan); +mrrob_wym_c2f6 = np.full(mrrob_wy_c2f6.shape+(len(n_mw),),np.nan); +for m in mind_cmip6: + for b in range(nbasin): + for n in range(len(n_mw)): + tasb_wym_c2f6[m,:,b,n] = moving_avg(tasb_wy_c2f6[m,:,b],n_mw[n]) + prb_wym_c2f6[m,:,b,n] = moving_avg(prb_wy_c2f6[m,:,b],n_mw[n]) + evspsblb_wym_c2f6[m,:,b,n] = moving_avg(evspsblb_wy_c2f6[m,:,b],n_mw[n]) + mrrob_wym_c2f6[m,:,b,n] = moving_avg(mrrob_wy_c2f6[m,:,b],n_mw[n]) +yrs=np.arange(syr_c2f6,eyr_c2f6) +inds1=int(np.nonzero(yrs==syr_sens_cri)[0]); inde1=int(np.nonzero(yrs==eyr_sens_cri)[0])+1; + +## get long-term average for % anomaly +prb_hist6_mw=np.mean(np.tile(prb_wy_c2f6[:,inds1:inde1,:,np.newaxis],(1,1,1,nmw)),axis=1); +tasb_hist6_mw=np.mean(np.tile(tasb_wy_c2f6[:,inds1:inde1,:,np.newaxis],(1,1,1,nmw)),axis=1); +mrrob_hist6_mw=np.mean(np.tile(mrrob_wy_c2f6[:,inds1:inde1,:,np.newaxis],(1,1,1,nmw)),axis=1); +evspsblb_hist6_mw=np.mean(np.tile(evspsblb_wy_c2f6[:,inds1:inde1,:,np.newaxis],(1,1,1,nmw)),axis=1); + +# yearly anomalies for regression analysis +nt=prb_wym_c2f6.shape[1] +aprb_wym_c2f6 = prb_wym_c2f6[:,:,:,:] - np.transpose(np.tile(prb_hist6_mw[:,:,:,np.newaxis],(1,1,1,nt)),(0,3,1,2)); +atasb_wym_c2f6 = tasb_wym_c2f6[:,:,:,:] - np.transpose(np.tile(tasb_hist6_mw[:,:,:,np.newaxis],(1,1,1,nt)),(0,3,1,2)); +amrrob_wym_c2f6 = mrrob_wym_c2f6[:,:,:,:] - np.transpose(np.tile(mrrob_hist6_mw[:,:,:,np.newaxis],(1,1,1,nt)),(0,3,1,2)); +aevspsblb_wym_c2f6 = evspsblb_wym_c2f6[:,:,:,:] - np.transpose(np.tile(evspsblb_hist6_mw[:,:,:,np.newaxis],(1,1,1,nt)),(0,3,1,2)); + +# percent anomaly +aprb_wym_c2f6_pct = aprb_wym_c2f6 / np.transpose(np.tile(prb_hist6_mw[:,:,:,np.newaxis],(1,1,1,nt)),(0,3,1,2)) * 100; +amrrob_wym_c2f6_pct = amrrob_wym_c2f6 / np.transpose(np.tile(mrrob_hist6_mw[:,:,:,np.newaxis],(1,1,1,nt)),(0,3,1,2)) * 100; +aevspsblb_wym_c2f6_pct = aevspsblb_wym_c2f6 / np.transpose(np.tile(evspsblb_hist6_mw[:,:,:,np.newaxis],(1,1,1,nt)),(0,3,1,2)) * 100; + +## get sensitivities +# initailizations +psens_hist6=np.full((len(mind_cmip6),nbasin,3),np.nan) +tsens_hist6=np.full((len(mind_cmip6),nbasin,3),np.nan) +intsens_hist6=np.full((len(mind_cmip6),nbasin,3),np.nan) +r2_ia_int_hist6=np.full((len(mind_cmip6),nbasin),np.nan) +r2_ia_hist6=np.full((len(mind_cmip6),nbasin),np.nan) +# calculate runoff sensitivity +yrs=np.arange(syr_c2f6,eyr_c2f6) +inds=int(np.nonzero(yrs==syr_sens_cri)[0]); inde=int(np.nonzero(yrs==eyr_sens_cri)[0])+1; +m2=0 +for m in mind_cmip6: + for b in range(nbasin): + a1, a2, b1, b2, c1, c2, r2 = runoff_sens_reg(amrrob_wym_c2f6_pct[m,inds:inde,b,tmw], aprb_wym_c2f6_pct[m,inds:inde,b,tmw], atasb_wym_c2f6[m,inds:inde,b,tmw], mw=n_mw[tmw]) + + psens_hist6[m2,b,0] = a1 + psens_hist6[m2,b,1:3] = a2 + tsens_hist6[m2,b,0] = b1 + tsens_hist6[m2,b,1:3] = b2 + intsens_hist6[m2,b,0] = c1 + intsens_hist6[m2,b,1:3] = c2 + r2_ia_int_hist6[m2,b] = r2 + + if not np.isnan(a1): + pred = a1*aprb_wym_c2f6_pct[m,inds:inde,b,tmw] + b1*atasb_wym_c2f6[m,inds:inde,b,tmw] + corr, _ = stats.pearsonr(pred,amrrob_wym_c2f6_pct[m,inds:inde,b,tmw]) + r2_ia_hist6[m2,b]= corr**2 + m2=m2+1 + + + +###################################################################### +# start and end year for sensitivity calculation +syr_sens_cri=1945 +eyr_sens_cri=2014 + +##################### +## CMIP5 ## +##################### +############################## +## HIST5 + RCP45 ## +############################## +## load variables +matfile = sio.loadmat(dpath + "tasb_wy_c2f5_cmip5.mat") +tasb_wy_c2f5 = matfile["tasb_wy_c2f5"][:,:,:] +matfile = sio.loadmat(dpath + "prb_wy_c2f5_cmip5.mat") +prb_wy_c2f5 = matfile["prb_wy_c2f5"][:,:,:] +matfile = sio.loadmat(dpath + "evspsblb_wy_c2f5_cmip5.mat") +evspsblb_wy_c2f5 = matfile["evspsblb_wy_c2f5"][:,:,:] +matfile = sio.loadmat(dpath + "mrrob_wy_c2f5_cmip5.mat") +mrrob_wy_c2f5 = matfile["mrrob_wy_c2f5"][:,:,:] +model_cmip5 = matfile["model_cmip5"][:][0] +syr_c2f5 = matfile["syr_c2f5"][0] +eyr_c2f5 = matfile["eyr_c2f5"][0] +nmodel_cmip5 = len(model_cmip5) + +################################################################ +# water budget closure for historical period (picontrol is similar) +# use the model with all variables we need +hval_available=np.isnan(np.mean(mrrob_wy_c2f5[:,:,3],axis=1))|np.isnan(np.mean(prb_wy_c2f5[:,:,3],axis=1))|np.isnan(np.mean(evspsblb_wy_c2f5[:,:,3],axis=1))|np.isnan(np.mean(tasb_wy_c2f5[:,:,3],axis=1)); +nanind_model_c2f5=np.nonzero(hval_available) +tasb_wy_c2f5[nanind_model_c2f5,:,:]=np.nan +prb_wy_c2f5[nanind_model_c2f5,:,:]=np.nan +evspsblb_wy_c2f5[nanind_model_c2f5,:,:]=np.nan +mrrob_wy_c2f5[nanind_model_c2f5,:,:]=np.nan + +# check the errors +pmeb_model = np.mean(prb_wy_c2f5[:,0:165,0:nbasin] - evspsblb_wy_c2f5[:,0:165,0:nbasin], axis=1) +error_val_c2f5 = (pmeb_model - np.mean(mrrob_wy_c2f5[:,0:165,0:nbasin], axis=1)) / np.mean(mrrob_wy_c2f5[:,0:165,0:nbasin], axis=1) +error_val_c2f5[nanind_model_c2f5,:]=np.nan +pval = 0.1 +closed_fraction_c2f5=np.nansum( np.abs(error_val_c2f5)0.001)[0] +nmodel_c2f5_closed=len(mind_c2f5) +nmodel_c2f5_available=len(mind_c2f5_available) +print(f"c2f5: {nmodel_c2f5_closed}/{nmodel_c2f5_available}") +print(f"negative values in c2f5 closed: {len(np.intersect1d(mind_negative_c2f5,mind_c2f5))}") + +# incorporate PI/HIST/SSP model info +nanind_cmip5=nanind_c2f5 +mind_cmip5_available=mind_c2f5_available +mind_cmip5=mind_c2f5 +mind_cmip5=np.setdiff1d(mind_cmip5,12) # 12-CSIRO-Mk3-6-0, continuous 0 values for many basins +nmodel_cmip5_available=len(mind_cmip5_available) +nmodel_cmip5_closed=len(mind_cmip5) +print(f"HIST6/SSP245: {nmodel_cmip5_closed}/{nmodel_cmip5_available}") +closed_fraction_cmip5=closed_fraction_c2f5 +closed_area_fraction_cmip5=closed_area_fraction_c2f5 +closed_nmodel_cmip5=closed_nmodel_c2f5 +############################################ +temp2=negative_fraction_c2f5>0.1 +ind_c2f5_nm=np.nonzero(temp2)[0] +ind_c2f5_nb=np.nonzero(temp2)[1] +for mind,bind in zip(ind_c2f5_nm,ind_c2f5_nb): + # display(f'c2f5 {model_cmip5[mind]}({mind}); {basin_names[bind]}({bind})') + # display(f'c2f5 mean Q: {np.nanmean(mrrob_wy_c2f5[mind,0:150,bind])}') + tasb_wy_c2f5[mind,:,bind]=np.nan + prb_wy_c2f5[mind,:,bind]=np.nan + evspsblb_wy_c2f5[mind,:,bind]=np.nan + mrrob_wy_c2f5[mind,:,bind]=np.nan + +# assign nan to unclosed models +tasb_wy_c2f5[nanind_cmip5,:,:]=np.nan +prb_wy_c2f5[nanind_cmip5,:,:]=np.nan +evspsblb_wy_c2f5[nanind_cmip5,:,:]=np.nan +mrrob_wy_c2f5[nanind_cmip5,:,:]=np.nan + +############################################################################ +def duplicated_fraction(data): + ndiff_frac=(len(data)-len(set(data)))/len(data) + return ndiff_frac # No duplicates found + +s=0 +for m in range(nmodel_cmip5): + for b in range(nbasin): + val=duplicated_fraction(mrrob_wy_c2f5[m,:,b]) + if val>0.1: + tasb_wy_c2f5[m,:,b]=np.nan + prb_wy_c2f5[m,:,b]=np.nan + evspsblb_wy_c2f5[m,:,b]=np.nan + mrrob_wy_c2f5[m,:,b]=np.nan + s=s+1 + +############################################################################ +## calculate sensitivity +# moving average +nmw=len(n_mw) +tasb_wym_c2f5 = np.full(tasb_wy_c2f5.shape+(len(n_mw),),np.nan); +prb_wym_c2f5 = np.full(prb_wy_c2f5.shape+(len(n_mw),),np.nan); +evspsblb_wym_c2f5 = np.full(evspsblb_wy_c2f5.shape+(len(n_mw),),np.nan); +mrrob_wym_c2f5 = np.full(mrrob_wy_c2f5.shape+(len(n_mw),),np.nan); +for m in mind_cmip5: + for b in range(nbasin): + for n in range(len(n_mw)): + tasb_wym_c2f5[m,:,b,n] = moving_avg(tasb_wy_c2f5[m,:,b],n_mw[n]) + prb_wym_c2f5[m,:,b,n] = moving_avg(prb_wy_c2f5[m,:,b],n_mw[n]) + evspsblb_wym_c2f5[m,:,b,n] = moving_avg(evspsblb_wy_c2f5[m,:,b],n_mw[n]) + mrrob_wym_c2f5[m,:,b,n] = moving_avg(mrrob_wy_c2f5[m,:,b],n_mw[n]) +yrs=np.arange(syr_c2f5,eyr_c2f5) +inds1=int(np.nonzero(yrs==syr_sens_cri)[0]); inde1=int(np.nonzero(yrs==eyr_sens_cri)[0])+1; + +## get long-term average for % anomaly +prb_hist5_mw=np.mean(np.tile(prb_wy_c2f5[:,inds1:inde1,:,np.newaxis],(1,1,1,nmw)),axis=1); +tasb_hist5_mw=np.mean(np.tile(tasb_wy_c2f5[:,inds1:inde1,:,np.newaxis],(1,1,1,nmw)),axis=1); +mrrob_hist5_mw=np.mean(np.tile(mrrob_wy_c2f5[:,inds1:inde1,:,np.newaxis],(1,1,1,nmw)),axis=1); +evspsblb_hist5_mw=np.mean(np.tile(evspsblb_wy_c2f5[:,inds1:inde1,:,np.newaxis],(1,1,1,nmw)),axis=1); + +# yearly anomalies for regression analysis +nt=prb_wym_c2f5.shape[1] +aprb_wym_c2f5 = prb_wym_c2f5[:,:,:,:] - np.transpose(np.tile(prb_hist5_mw[:,:,:,np.newaxis],(1,1,1,nt)),(0,3,1,2)); +atasb_wym_c2f5 = tasb_wym_c2f5[:,:,:,:] - np.transpose(np.tile(tasb_hist5_mw[:,:,:,np.newaxis],(1,1,1,nt)),(0,3,1,2)); +amrrob_wym_c2f5 = mrrob_wym_c2f5[:,:,:,:] - np.transpose(np.tile(mrrob_hist5_mw[:,:,:,np.newaxis],(1,1,1,nt)),(0,3,1,2)); +aevspsblb_wym_c2f5 = evspsblb_wym_c2f5[:,:,:,:] - np.transpose(np.tile(evspsblb_hist5_mw[:,:,:,np.newaxis],(1,1,1,nt)),(0,3,1,2)); + +# percent anomaly +aprb_wym_c2f5_pct = aprb_wym_c2f5 / np.transpose(np.tile(prb_hist5_mw[:,:,:,np.newaxis],(1,1,1,nt)),(0,3,1,2)) * 100; +amrrob_wym_c2f5_pct = amrrob_wym_c2f5 / np.transpose(np.tile(mrrob_hist5_mw[:,:,:,np.newaxis],(1,1,1,nt)),(0,3,1,2)) * 100; +aevspsblb_wym_c2f5_pct = aevspsblb_wym_c2f5 / np.transpose(np.tile(evspsblb_hist5_mw[:,:,:,np.newaxis],(1,1,1,nt)),(0,3,1,2)) * 100; + +## get sensitivities +# initailizations +psens_hist5=np.full((len(mind_cmip5),nbasin,3),np.nan) +tsens_hist5=np.full((len(mind_cmip5),nbasin,3),np.nan) +intsens_hist5=np.full((len(mind_cmip5),nbasin,3),np.nan) +r2_ia_int_hist5=np.full((len(mind_cmip5),nbasin),np.nan) +r2_ia_hist5=np.full((len(mind_cmip5),nbasin),np.nan) +# calculate runoff sensitivity +yrs=np.arange(syr_c2f5,eyr_c2f5) +inds=int(np.nonzero(yrs==syr_sens_cri)[0]); inde=int(np.nonzero(yrs==eyr_sens_cri)[0])+1; +m2=0 +for m in mind_cmip5: + for b in range(nbasin): + a1, a2, b1, b2, c1, c2, r2 = runoff_sens_reg(amrrob_wym_c2f5_pct[m,inds:inde,b,tmw], aprb_wym_c2f5_pct[m,inds:inde,b,tmw], atasb_wym_c2f5[m,inds:inde,b,tmw], mw=n_mw[tmw]) + + psens_hist5[m2,b,0] = a1 + psens_hist5[m2,b,1:3] = a2 + tsens_hist5[m2,b,0] = b1 + tsens_hist5[m2,b,1:3] = b2 + intsens_hist5[m2,b,0] = c1 + intsens_hist5[m2,b,1:3] = c2 + r2_ia_int_hist5[m2,b] = r2 + + if not np.isnan(a1): + pred = a1*aprb_wym_c2f5_pct[m,inds:inde,b,tmw] + b1*atasb_wym_c2f5[m,inds:inde,b,tmw] + corr, _ = stats.pearsonr(pred,amrrob_wym_c2f5_pct[m,inds:inde,b,tmw]) + r2_ia_hist5[m2,b]= corr**2 + m2=m2+1 + + +###################################################################### +# start and end year for sensitivity calculation +syr_sens_cri=1945 +eyr_sens_cri=2014 + +##################################################### +## OBS (GRUN + [HADCRUT, GPCC, ERA5-land]) ## +##################################################### +## load data +filename=f'{dpath}tasb_wy_obs2.mat' +matfile=sio.loadmat(filename) +tasb_wy_obs2=matfile['tasb_wy_obs2'][0:5,:,:] +syr_obs=matfile['syr_obs'][0] +eyr_obs=matfile['eyr_obs'][0] + +filename=f'{dpath}prb_wy_obs2.mat' +matfile=sio.loadmat(filename) +prb_wy_obs2=matfile['prb_wy_obs2'][0:5,:,:] + +filename=f'{dpath}mrrob_wy_obs2_ens.mat' +matfile=sio.loadmat(filename) +mrrob_wy_obs2=matfile['mrrob_wy_obs2_ens'][0:5,:,:,:] +nfoc=mrrob_wy_obs2.shape[0] +nens=mrrob_wy_obs2.shape[1] + +## do moving average +nmw=len(n_mw) +tasb_wym_obs2 = np.full(tasb_wy_obs2.shape+(len(n_mw),),np.nan); +prb_wym_obs2 = np.full(prb_wy_obs2.shape+(len(n_mw),),np.nan); +for f in range(nfoc): + for b in range(nbasin): + for w in range(len(n_mw)): + tasb_wym_obs2[f,:,b,w] = moving_avg(tasb_wy_obs2[f,:,b],n_mw[w]) + prb_wym_obs2[f,:,b,w] = moving_avg(prb_wy_obs2[f,:,b],n_mw[w]) + +mrrob_wym_obs2 = np.full(mrrob_wy_obs2.shape+(len(n_mw),),np.nan); +for f in range(nfoc): + for e in range(nens): + for b in range(nbasin): + for w in range(len(n_mw)): + mrrob_wym_obs2[f,e,:,b,w] = moving_avg(mrrob_wy_obs2[f,e,:,b],n_mw[w]) + +## get % anomalies +nyr=prb_wy_obs2.shape[1] +yrs=np.arange(syr_obs,eyr_obs) +inds=int(np.nonzero(yrs==syr_sens_cri)[0]); inde=int(np.nonzero(yrs==eyr_sens_cri)[0])+1; + +prb_wym_obs2_mean=np.tile(np.nanmean(prb_wy_obs2[:,inds:inde,:],axis=1)[:,:,np.newaxis],(1,1,nmw)); +tasb_wym_obs2_mean=np.tile(np.nanmean(tasb_wy_obs2[:,inds:inde,:],axis=1)[:,:,np.newaxis],(1,1,nmw)); +mrrob_wym_obs2_mean=np.tile(np.nanmean(mrrob_wy_obs2[:,:,inds:inde,:],axis=2)[:,:,:,np.newaxis],(1,1,nmw)); + +aprb_wym_gobs = prb_wym_obs2[:,:,:,:] - np.transpose(np.tile(prb_wym_obs2_mean[:,:,:,np.newaxis],(1,1,1,nyr)),(0,3,1,2)); +atasb_wym_gobs = tasb_wym_obs2[:,:,:,:] - np.transpose(np.tile(tasb_wym_obs2_mean[:,:,:,np.newaxis],(1,1,1,nyr)),(0,3,1,2)); +amrrob_wym_gobs = mrrob_wym_obs2[:,:,:,:,:] - np.transpose(np.tile(mrrob_wym_obs2_mean[:,:,:,:,np.newaxis],(1,1,1,1,nyr)),(0,1,4,2,3)); +aprb_wym_gobs_pct = aprb_wym_gobs / np.transpose(np.tile(prb_wym_obs2_mean[:,:,:,np.newaxis],(1,1,1,nyr)),(0,3,1,2)) * 100; +amrrob_wym_gobs_pct = amrrob_wym_gobs / np.transpose(np.tile(mrrob_wym_obs2_mean[:,:,:,:,np.newaxis],(1,1,1,1,nyr)),(0,1,4,2,3)) * 100; + + +psens_obs=np.full((nfoc-1,nens,nbasin,3),np.nan) +tsens_obs=np.full((nfoc-1,nens,nbasin,3),np.nan) +intsens_obs=np.full((nfoc-1,nens,nbasin,3),np.nan) +r2_ia_int_obs=np.full((nfoc-1,nens,nbasin),np.nan) +r2_ia_obs=np.full((nfoc-1,nens,nbasin),np.nan) +# calculate runoff sensitivity +f2=0 +for f in [0,1,3,4]: # skip GSWP3 + print(f'forcing: {f}') + for e in range(nens): + for b in range(nbasin): + a1, a2, b1, b2, c1, c2, r2 = runoff_sens_reg(amrrob_wym_gobs_pct[f,e,inds:inde,b,tmw], aprb_wym_gobs_pct[f,inds:inde,b,tmw], atasb_wym_gobs[f,inds:inde,b,tmw], mw=n_mw[tmw]) + + psens_obs[f2,e,b,0] = a1 + psens_obs[f2,e,b,1:3] = a2 + tsens_obs[f2,e,b,0] = b1 + tsens_obs[f2,e,b,1:3] = b2 + intsens_obs[f2,e,b,0] = c1 + intsens_obs[f2,e,b,1:3] = c2 + r2_ia_int_obs[f2,e,b] = r2 + + if not np.isnan(a1): + pred = a1*aprb_wym_gobs_pct[f,inds:inde,b,tmw] + b1*atasb_wym_gobs[f,inds:inde,b,tmw] + corr, _ = stats.pearsonr(pred,amrrob_wym_gobs_pct[f,e,inds:inde,b,tmw]) + r2_ia_obs[f2,e,b]= corr**2 + f2=f2+1 + + + +################################################################ +## save netcdf file ####### +################################################################ +import numpy as np +from netCDF4 import Dataset + +# define the dimensions +nb = nbasin +nfoc = 4 +nmodel6 = len(mind_cmip6) +nmodel5 = len(mind_cmip5) + +# change name for r2 values (training accuracy) +r2_obs_pred = r2_ia_obs +r2_hist5_pred = r2_ia_hist5 +r2_hist6_pred = r2_ia_hist6 + +sens_types = ['sensitivity value', 'lower bound (95% confidence interval of reg. coeff.)', 'upper bound (95% confidence interval of reg. coeff.)'] +foc_names=['CRUTSv4.04', 'PGFv2', 'GSWP3-EWEMBI', 'GSWP3-W5E5'] + +## save OBS data +# Create a NetCDF file +fpath='C:/runoff_mdtf2/' +ncfile = Dataset(fpath + 'runoff_sensitivity_obs.nc', 'w', format='NETCDF4') +# Define the dimensions +ncfile.createDimension('foc', nfoc) +ncfile.createDimension('ens', nens) +ncfile.createDimension('basin', nb) +ncfile.createDimension('sens_type', 3) +ncfile.createDimension('string_length', max([len(name) for name in sens_types])) # For basin_names +# Create variables +psens_obs_var = ncfile.createVariable('psens_obs', np.float32, ('foc', 'ens', 'basin', 'sens_type')) +tsens_obs_var = ncfile.createVariable('tsens_obs', np.float32, ('foc', 'ens', 'basin', 'sens_type')) +intsens_obs_var = ncfile.createVariable('intsens_obs', np.float32, ('foc', 'ens', 'basin', 'sens_type')) +r2_obs_pred_var = ncfile.createVariable('r2_obs_pred', np.float32, ('foc', 'ens', 'basin')) +basin_names_var = ncfile.createVariable('basin_names', 'S1', ('basin', 'string_length')) +sens_type_names_var = ncfile.createVariable('sens_type_names', 'S1', ('sens_type', 'string_length')) +foc_names_var = ncfile.createVariable('foc_names', 'S1', ('foc', 'string_length')) +# Write data to variables +psens_obs_var[:,:,:,:] = psens_obs +tsens_obs_var[:,:,:,:] = tsens_obs +intsens_obs_var[:,:,:,:] = intsens_obs +r2_obs_pred_var[:,:,:] = r2_obs_pred +for i, name in enumerate(basin_names): + basin_names_var[i, :len(name)] = list(name) +for i, type in enumerate(sens_types): + sens_type_names_var[i, :len(type)] = list(type) +for i, name in enumerate(foc_names): + foc_names_var[i, :len(name)] = list(name) +# additional informations +ncfile.description = "Observational estimation of runoff sensitivity from GRUN dataset (ref: https://doi.org/10.1029/2020WR028787)." +ncfile.ensemble = "Different ensembles are results of different sub-sampling during the trianing of machine-learning algorithm. Refer to above publication for details." +ncfile.methods = "Sensitivities are estimated as lienar regression coefficients of runoff onto the precipitation and temperature, using 5-year averaged water-year anomalies in 1945-2014 period." +ncfile.contact = "Hanjun Kim (hanjunkim0617@gmail.com)" +psens_obs_var.units = "%/%" +tsens_obs_var.units = "%/K" +intsens_obs_var.description = "Sensitivity of interaction term, which often has negligible role compared to other two sensitivities." +r2_obs_pred_var.description = "R2 values between predictions (using psens and tsens) and original runoff variation, representing the singificance of the regression coefficients." +basin_names_var.description = "Names of major river basins. Masks used for basin-average are from GRDC (https://mrb.grdc.bafg.de/)." +sens_type_names_var.description = "Sensitivity value is provided with the lower and upper bounds of regression coefficients using t-test with 95% confidence interval." +foc_names_var.description = "The atmospheric forcings used to generate GRUN runoff dataset (ref: https://doi.org/10.1029/2020WR028787)." +# Close the NetCDF file +ncfile.close() + +## save hist6 data +# Create a NetCDF file +fpath='C:/runoff_mdtf2/' +ncfile = Dataset(fpath + 'runoff_sensitivity_hist6.nc', 'w', format='NETCDF4') +# Define the dimensions +ncfile.createDimension('model', nmodel6) +ncfile.createDimension('basin', nb) +ncfile.createDimension('sens_type', 3) +# Create variables +psens_hist6_var = ncfile.createVariable('psens_hist6', np.float32, ('model', 'basin', 'sens_type')) +tsens_hist6_var = ncfile.createVariable('tsens_hist6', np.float32, ('model', 'basin', 'sens_type')) +intsens_hist6_var = ncfile.createVariable('intsens_hist6', np.float32, ('model', 'basin', 'sens_type')) +r2_hist6_pred_var = ncfile.createVariable('r2_hist6_pred', np.float32, ('model', 'basin')) +# Write data to variables +psens_hist6_var[:,:,:] = psens_hist6 +tsens_hist6_var[:,:,:] = tsens_hist6 +intsens_hist6_var[:,:,:] = intsens_hist6 +r2_hist6_pred_var[:,:] = r2_hist6_pred +# Close the NetCDF file +ncfile.close() + +## save hist5 data +# Create a NetCDF file +fpath='C:/runoff_mdtf2/' +ncfile = Dataset(fpath + 'runoff_sensitivity_hist5.nc', 'w', format='NETCDF4') +# Define the dimensions +ncfile.createDimension('model', nmodel5) +ncfile.createDimension('basin', nb) +ncfile.createDimension('sens_type', 3) +# Create variables +psens_hist5_var = ncfile.createVariable('psens_hist5', np.float32, ('model', 'basin', 'sens_type')) +tsens_hist5_var = ncfile.createVariable('tsens_hist5', np.float32, ('model', 'basin', 'sens_type')) +intsens_hist5_var = ncfile.createVariable('intsens_hist5', np.float32, ('model', 'basin', 'sens_type')) +r2_hist5_pred_var = ncfile.createVariable('r2_hist5_pred', np.float32, ('model', 'basin')) +# Write data to variables +psens_hist5_var[:,:,:] = psens_hist5 +tsens_hist5_var[:,:,:] = tsens_hist5 +intsens_hist5_var[:,:,:] = intsens_hist5 +r2_hist5_pred_var[:,:] = r2_hist5_pred +# Close the NetCDF file +ncfile.close() + + + + +################################################################ +## calculation check: draw maps with filled river basins ####### +################################################################ +def plot_and_save_basin_filled(values, basin_points, color_bins, color_bins2, plt_colormap, plt_unit, plt_title, plt_path, coast_path): + ## data needed for plotting + # assign color index to each values + values_ind=np.digitize(values, color_bins2)-1 + # make colormap + custom_cmap=colors.LinearSegmentedColormap.from_list('custom_colormap',plt_colormap, N=len(color_bins)-1) + # load coastline + lonmap = nc.Dataset(coast_path)['lonmap'][0:9850] + latmap = nc.Dataset(coast_path)['latmap'][0:9850] + ## draw figure + fig, ax = plt.subplots(figsize=(12, 4.8)) + plt.rcParams.update({'font.size': 12}) + # coastline + ax.plot(lonmap, latmap, color=[0.8, 0.8, 0.8, 0], linewidth=0.5) + ax.add_patch(patches.Polygon(np.column_stack([lonmap, latmap]), closed=True, facecolor=(0.8, 0.8, 0.8), edgecolor=(0.5, 0.5, 0.5), linestyle='none')) + # fill basins with colors corresponding to target values + nb=len(values) + for b in range(nb): + X = [item[0] for item in basin_points[b]] + Y = [item[1] for item in basin_points[b]] + X = X[0:-int(np.floor(len(X)/100))] + Y = Y[0:-int(np.floor(len(Y)/100))] + if not np.isnan(values_ind[b]): + ax.add_patch(patches.Polygon(np.column_stack([X, Y]), closed=True, facecolor=custom_cmap(values_ind[b]), edgecolor=(0.5, 0.5, 0.5), linewidth=0.5)) + else: + ax.plot(X, Y, color=(0.5, 0.5, 0.5), linewidth=0.5) + # Set colormap and colorbar + cb = plt.colorbar(plt.cm.ScalarMappable(norm=plt.Normalize(0, 1), cmap=custom_cmap)) + cb.set_ticks(np.linspace(0, 1, len(color_bins))) + clabels = [str(b) for b in color_bins] + cb.set_ticklabels(clabels) + cb.set_label(plt_unit, fontsize=12) + # Customize the ticks and labels + ax.set_xticks(range(-180, 181, 30)) + ax.set_xticklabels(['', '', '120°W', '', '60°W', '', '0°', '', '60°E', '', '120°E', '', '']) + ax.set_yticks(range(-90, 91, 15)) + ax.set_yticklabels(['', '', '60°S', '', '30°S', '', '0°', '', '30°N', '', '60°N', '', '']) + ax.set_xlim([-180, 180]) + ax.set_ylim([-60, 80]) + # Set title + ax.set_title(plt_title, fontsize=17) + # Save the figure as a eps + plt.savefig(plt_path, bbox_inches='tight') + # plt.close() + +## figures for T sensitivity ## +bins = [-30,-15,-12,-9,-6,-3,0,3,6,9,12,15,30] +bins2 = [-60,-15,-12,-9,-6,-3,0,3,6,9,12,15,60] +plot_colormap = [(0.4000, 0, 0, 1),(0.7706, 0, 0, 1),(0.9945, 0.0685, 0.0173, 1),(0.9799, 0.2483, 0.0627, 1),(0.9715, 0.4442, 0.0890, 1),(0.9845, 0.6961, 0.0487, 1),(0.9973, 0.9480, 0.0083, 1),(1.0000, 1.0000, 0.3676, 1),(1.0000, 1.0000, 1.0000, 1),(1.0000, 1.0000, 1.0000, 1),(0.6975, 0.8475, 0.9306, 1),(0.4759, 0.7358, 0.8797, 1),(0.2542, 0.6240, 0.8289, 1),(0.0436, 0.5130, 0.7774, 1),(0.0533, 0.4172, 0.7138, 1),(0.0630, 0.3215, 0.6503, 1),(0.0411, 0.1760, 0.5397, 1),(0, 0, 0.4000, 1)] +plot_unit='[%/K]' +coast_path = dpath + "/coastline.nc" +# OBS +values=np.nanmean(np.nanmean(tsens_obs[:,:,:,0],axis=1),axis=0) +# values=np.nanmean(tsens_hist5[:,:,0],axis=0) +# values=np.nanmean(tsens_hist6[:,:,0],axis=0) +plot_title='T sensitivity' +plot_path = figpath + 'mdtf_check.png' +plot_and_save_basin_filled(values, basin_points, bins, bins2, \ + plot_colormap, plot_unit, plot_title, plot_path, coast_path) diff --git a/diagnostics/runoff_sensitivities/runoff_sensitivities.html b/diagnostics/runoff_sensitivities/runoff_sensitivities.html new file mode 100644 index 000000000..a8e859860 --- /dev/null +++ b/diagnostics/runoff_sensitivities/runoff_sensitivities.html @@ -0,0 +1,365 @@ + +runoff sensitivity diagnostic + +

+Runoff sensitivity: The sensitivity of runoff to surface air temperature and precipitation +

+ +

+General discription
+The sensitivity of runoff (Q) to temperature (T sensitivity; %/K) and precipitation (P sensitivity; %/%) +are calculated using the multiple linear regression, using the 5-year averaged annual timeseries of the historical period during 1945-2014. +The sensitivity is calculated respectively for 131 global river basins. +Before calculating the sensitivity, the variables are averaged for water year (Oct. to Sep.) and each basin using the GRDC Major River Basin masks.

+ +The calculated runoff sensitivity is then compared with the pre-calculated sensitivities for observations (GRUN-ENSEMBLE), CMIP5/6 historical simulations. +GRUN-ENSEMBLE is machine-learning based global runoff reanalysis, providing 100 realizations to sample to observational uncertainty (Ghiggi et al. 2021). +For CMIP models, we use historical simulation from 1945 to 2014. For CMIP5, we extend the timeseries to 2014 using RCP4.5 scenario. +The ranges from observational uncertainty, inter-model spread, and uncertainty of regression coefficients are shown for individual river basins.

+ +In CMIP5 and CMIP6, the uncertainty of T sensitivity can contribute to the runoff projection uncertainty as much as the different P projections among models. +The runoff sensitivity in climate model is often biased; the negative T sensitivity is often too weak in climate models, indicating the underestimation of future runoff declines.

+ +While the P sensitivity is generally correlated with the mean state runoff ratio (Q/P), the T sensitivity exhibits no systematic relationship with mean state variables. +Hence, the traditional modeling approaches, which focus on improving mean state biases, may not resolve the biases in T sensitivity. +Therefore, the new diagnostic metric related to the runoff sensitivity is needed for future land model development. +

+ +

+References
+Calculation and interpretation of runoff sensitivities:
+Lehner, F., Wood, A. W., Vano, J. A., Lawrence, D. M., Clark, M. P., and Mankin, J. S. (2019). +The potential to reduce uncertainty in regional runoff projections from climate models. Nature Climate Change, 9(12), 926-933.
+Kim, H., Lehner, F., Dagon K., Lawrence, D. M., Swenson, S., Wood, A. W. (2025). +Constraining climate model projections with observations amplifies future runoff declines. in revision.
+Observational estimate (GRUN-ENSEMBLE):
+Ghiggi, G., Humphrey, V., Seneviratne, S. I., & Gudmundsson, L. (2021). +G‐RUN ENSEMBLE: A multi‐forcing observation‐based global runoff reanalysis. Water Resources Research, 57(5), e2020WR028787. +

+

+Caution
+For some models, the surface water budget (Precipitation - Evaporation = Runoff) is not closed in the long-term average.
+To evaluate the water budget closure (WBC) for the period of {{FIRSTYR}}-{{LASTYR}}, we specify river basins in which WBC is kept within 10% error. The fraction of closed river basins compared to all river basins is used as criteria.
+For example, WBC = 0.91 for CESM2, and it means that WBC is kept for 91% of the area of global river basins.
+if WBC < 0.6, we consider the model's runoff data to be inappropriate for runoff sensitivity calculations. We do not know why the water balance is not closeable in some models; it might have more to do with the data they provide in the CMIP archive rather than the model itself.
+WBC value is shown in the plots of MODEL runoff sensitivities (title). +

+

+Casename: {{CASENAME}} +

+ +


+ + + + + + + + +
Runoff sensitivities + OBS + MODEL + MODEL-OBS + CMIP6-OBS +
T sensitivity (%/K) + plot + plot + plot + plot +
P sensitivity (%/%) + plot + plot + plot + plot +
+ +


+ + + + + + + + +
Summary +
Prediction accuracy + T sensitivity + P sensitivity +
plot + plot + plot +
+ +


+ + + + + + + + + + + + + + + + +
Global map of river basins with significant biases +
In the above map, the sensitivity bias is significant if the model value falls outside the 5th–95th percentile range of the observational sensitivity. +
Basins outlined in black indicate regions where the T sensitivity bias is significant. +
Basins with black numbering indicate regions where either the T or P sensitivity bias is significant. +
Note that the criterion of being “outside the 5th–95th percentile range” is more stringent than conventional significance tests, such as the t-test. +
Basins with hatchings indicate regions where the regression model for runoff sensitivity does not have significant R2 relative to the original variability; this means that the sensitivity is not reliable for those basins in the target model. +
+ + +


+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
For individual river basins:  +Sensitivities with uncertainty ranges +
  +1:CONGO   +2:CUANZA   +3:DRAA   +4:KUNENE   +5:NIGER   +
  +6:NILE   +7:OGOOUE   +8:ORANGE   +9:SANAGA   +10:SENEGAL   +
  +11:VOLTA   +12:JUBBA   +13:LIMPOPO   +14:RUFIJI   +15:RUVUMA   +
  +16:SAVE   +17:TANA   +18:ZAMBEZI   +19:ABHE BAD   +20:CHAD   +
  +21:CUVELAI   +22:GUIR   +23:LAKE TURKANA   +24:MELRHIR   +25:OKAVANGO   +
  +26:RHARSA   +27:INDIGIRKA   +28:KHATANGA   +29:KOLYMA   +30:LENA   +
  +31:OB   +32:OLENEK   +33:PUR   +34:PYASINA   +35:TAZ   +
  +36:YANA   +37:YENISEY   +38:BRAHMAPUTRA   +39:GANGES   +40:GODAVARI   +
  +41:INDUS   +42:IRRAWADDY   +43:KRISHNA   +44:MAHANADI   +45:MEGHNA   +
  +46:SALWEEN   +47:SHATT AL ARAB   +48:AMUR   +49:ANADYR   +50:CHAO PHRAYA   +
  +51:HUAI HE   +52:LIAO HE   +53:MEKONG   +54:RED   +55:YANGTZE   +
  +56:YELLOW RIVER   +57:HAI HE   +58:ZHU   +59:ARAL SEA   +60:CHUY   +
  +61:HAMUN   +62:HAMUN-I-MASHKEL   +63:HAR US NUUR   +64:LAKE BALKHASH   +65:SARYSU   +
  +66:TARIM HE   +67:TORGHAY   +68:URAL   +69:AMAZON   +70:CHUBUT   +
  +71:COLORADO   +72:MAGDALENA   +73:NEGRO   +74:ORINOCO   +75:PARANA   +
  +76:PARNAIBA   +77:RIO DE LA PLATA   +78:SAO FRANCISCO   +79:TOCANTINS   +80:URUGUAY   +
  +81:MAR CHIQUITA   +82:ALBANY   +83:BACK   +84:CANIAPISCAU   +85:CHURCHILL   +
  +86:HAYES   +87:LA GRANDE   +88:MACKENZIE   +89:MOOSE   +90:NELSON   +
  +91:SEVERN   +92:THELON   +93:BRAZOS   +94:COLORADO   +95:MISSISSIPPI   +
  +96:MOBILE RIVER   +97:RIO GRANDE   +98:SAINT LAWRENCE   +99:USUMACINTA   +100:COLORADO   +
  +101:COLUMBIA   +102:FRASER   +103:KUSKOKWIM   +104:SANTIAGO   +105:YUKON   +
  +106:GREAT SALT LAKE   +107:FITZROY   +108:MURRAY   +109:BURDEKIN   +110:FITZROY   +
  +111:FLINDERS   +112:VICTORIA   +113:LAKE EYRE   +114:LAKE GAIRDNER   +115:PECHORA   +
  +116:SEVERNAYA DVINA   +117:DANUBE   +118:DNIEPER   +119:DON   +120:DOURO   +
  +121:ELBE   +122:LOIRE   +123:NEVA   +124:ODER   +125:RHINE   +
  +126:RHONE   +127:VISTULA   +128:KURA   +129:VOLGA   +130:SACRAMENTO   +
  +131:SAN JOAQUIN   +
+ +

+For each river basin, the uncertainty ranges for OBS, MODEL, CMIP5, and CMIP6 are defined as below:
+OBS (black errorbar) - 5%-95% confidence interval of regression coefficients using t-test, averaged across 100 observational members (black dots)
+OBS (gray errorbar) - 5%-95% percentiles using 100 members. Note that 5%-95% confidence interval using t-test becomes too narrow with sample size of 100
+ +MODEL - 5%-95% confidence interval of regression coefficients using t-test
+CMIP5 - 5%-95% confidence interval of inter-model spread using t-test (22 models, green dots)
+CMIP6 - 5%-95% confidence interval of inter-model spread using t-test (28 models, blue dots)
+

+ diff --git a/diagnostics/runoff_sensitivities/runoff_sensitivities_diag.py b/diagnostics/runoff_sensitivities/runoff_sensitivities_diag.py new file mode 100644 index 000000000..041b5a56a --- /dev/null +++ b/diagnostics/runoff_sensitivities/runoff_sensitivities_diag.py @@ -0,0 +1,1019 @@ +# ================================================================================ +# MDTF Diagnostic POD +# ================================================================================ +# +# Runoff sensitivities Diagnostic POD +# +# Last update: 12/02/2024 +# +# Synopsis +# +# Runoff projections are essential for future water security. +# The Earth System Models (ESMs) are increasingly being utilized for future water resource risk assessment. +# However, the runoff projections based on ESMs are highly uncertain, in part due to differences in the land process representation among ESMs. +# +# In this diagnostics, we try to measure the land process biases in ESMs. +# First, the land processes related to runoff projections in each ESM can be statistically emulated by +# quantifying the inter-decadal sensitivity of runoff to temperature (temperature sensitivity) and +# precipitation (precipitation sensitivity) using multiple linear regression (Lehner et al. 2019; Kim et al. 2024). +# To represent the land process biases, the runoff senstivities for each ESM are compared to +# observational estimations, which is pre-calculated using same regression method. +# For the observational estimation, we used the GRUN-ENSEMBLE data +# - global reanalysis of monthly runoff using the machine learning technique (Ghiggi et al. 2021). +# The runoff sensitivities from CMIP5/6 models are also prepared to facilitate the comparison for new model development. +# The uncertainty ranges from internal variability, observational uncertainty, and inter-moel spread +# are provided for specific river basins to assess the significance of ESM biases. +# +# Version & Contact info +# +# - Version/revision information: version 3 (12/02/2024) +# - PI: Flavio Lehner, Cornell University, flavio.lehner@cornell.edu +# - Developer/point of contact: Hanjun Kim, Cornell University, hk764@cornell.edu +# - Other contributors: +# +# Open source copyright agreement +# +# The MDTF framework is distributed under the LGPLv3 license (see LICENSE.txt). +# +# Functionality +# +# The main driver code (runoff_sensitivities_diag.py) include all functions, codes for calculations and figures. +# 0) Functions are pre-defined for further analysis. +# 1) The codes will load data and do area average for 131 global river basins for each water year (Mask source - GRDC Major River Basins of the World). +# 2) Water budget clousre (precipitaton - evapotranspiration == runoff) in long term average is checked. +# 3) Using the pre-defined function "runoff_sens_reg", runoff sensitivities are calculated. +# 4) Calculated runoff sensitivities for target models are saved as .nc file. +# 5) The diagnostic figures will be plotted with the pre-calculated OBS and CMIP data. +# +# Programm language summary +# +# Python 3 is used to calculate and draw the figures. +# All libraries used in this diagnostic are available in MDTF conda environment "_MDTF_python3_base". +# - Used libraries: "scipy", "numpy", "matplotlib", "netCDF4", "cartopy", "sklearn" +# - To deal with the shapefile, "cartopy.io.shapereader" and "matplotlib.path" is utilized. +# - For multi-linear regression, "sklearn.linear_model" is utilized. +# +# Required model output variables +# +# - The monthly historical simulations including period 1940-2014 are needed (Model outputs are assumed to be same with CMIP6 output). +# - Target variables +# - tas (surface air temperature, K), [time, lat, lon] +# - pr (precipitaiton, kg m-2 s-1), [time, lat, lon] +# - hfls (latent heat flux, W m-2), [time, lat, lon] +# - mrro (runoff, kg m-2 s-1), [time, lat, lon] +# - lon-lat grids for 4 variables have to be same +# (In CMIP, there are some cases where grids are slightly different between land and atm variables. Checking/interpolation is recommended) +# +# +# References +# +# Lehner, F., Wood, A. W., Vano, J. A., Lawrence, D. M., Clark, M. P., & Mankin, J. S. (2019). +# The potential to reduce uncertainty in regional runoff projections from climate models. Nature Climate Change, 9(12), 926-933. +# doi: 10.1038/s41558-019-0639-x +# +# Kim, H., Lehner, F., Dagon, K., Lawrence, D. M., Swenson, S., & Wood, A. W. (2024). +# Constraining climate model projections with observations amplifies future runoff declines. In preparation +# +# Ghiggi, G., Humphrey, V., Seneviratne, S. I., & Gudmundsson, L. (2021). +# G‐RUN ENSEMBLE: A multi‐forcing observation‐based global runoff reanalysis. Water Resources Research, 57(5), e2020WR028787. +# doi:10.1029/2020WR028787 +# ================================================================================ + +import os +import warnings +import time as time +import numpy as np # python library for dealing with data array +import cartopy.io.shapereader as shpreader # cartopy library for loading with shapefiles +import matplotlib +# matplotlib.use('Agg') # non-X windows backend +import matplotlib.pyplot as plt # python library we use to make plots +import matplotlib.patches as patches # python library for filling the polygons in plots +import matplotlib.colors as colors # python library for customizing colormaps +import matplotlib.path as mpltPath # matplotlib library to deal with polygons +from sklearn.linear_model import LinearRegression # python library we use to do multi-linear regression +import netCDF4 as nc # python library for handling nc files +import scipy.stats as stats # python library we use to do basic statistics + +############################################################################################## +### 0) pre-defined functions for calculations: ############################################### +############################################################################################## +## function to get area weight +def area_weight(lat, lon): + a = 6371 * 1e3 + nlat = len(lat) + nlon = len(lon) + lon2d, lat2d = np.meshgrid(lon, lat) + dx = a * np.outer(np.cos(np.deg2rad(lat)), np.gradient(np.deg2rad(lon))) + dy = np.tile(a * np.gradient(np.deg2rad(lat)), (nlon, 1)).T + aw = dx * dy + return aw + +################################################################ +## function for runoff sensitivity calculation (least square) ## +################################################################ +from sklearn.linear_model import LinearRegression +import scipy.stats as stats +import numpy as np +def runoff_sens_reg(r, p, t, mw=1, alpha=0.1): + nan_val=np.mean(np.isnan(r))+np.mean(np.isnan(p))+np.mean(np.isnan(t)) + if nan_val == 0: + # Ensure input vectors are column vectors + if r.ndim == 1: + r = r[:, np.newaxis] + if p.ndim == 1: + p = p[:, np.newaxis] + if t.ndim == 1: + t = t[:, np.newaxis] + + # Create the regression matrix and do normalization + Xraw = np.column_stack((p, t, p * t)) + Xstd = np.nanstd(Xraw,axis=0) + rstd = np.nanstd(r,axis=0) + Xnorm = (Xraw - np.nanmean(Xraw,axis=0)) / Xstd + rnorm = (r - np.nanmean(r,axis=0)) / rstd + + # Perform linear regression + model = LinearRegression() + model.fit(Xnorm, rnorm) + + # Get regression coefficients + a1 = np.squeeze(model.coef_[0][0] * rstd / Xstd[0]) + b1 = np.squeeze(model.coef_[0][1] * rstd / Xstd[1]) + c1 = np.squeeze(model.coef_[0][2] * rstd / Xstd[2]) + + # Calculate the standard errors for coefficients + y_pred = a1*Xraw[:,0] + b1*Xraw[:,1] + c1*Xraw[:,2] + residuals = r - y_pred[:,np.newaxis] + n, k = Xraw.shape + # autocorr=get_autocorr(r,1); + # dof2 = (n-k)*(1-autocorr)/(1+autocorr) + dof=np.ceil((n-k)/mw) + mse = np.sum(residuals ** 2) / (dof) + var_b = mse * np.linalg.pinv(np.dot(Xraw.T, Xraw)) + se_a1 = np.sqrt(var_b[0, 0]) + se_b1 = np.sqrt(var_b[1, 1]) + se_c1 = np.sqrt(var_b[2, 2]) + + # Calculate the t-statistic for a given confidence level (e.g., alpha = 0.05) + t_critical = stats.t.ppf(1 - alpha / 2, dof) + + # Calculate the confidence intervals + a2 = np.array([a1 - t_critical * se_a1, a1 + t_critical * se_a1]) + b2 = np.array([b1 - t_critical * se_b1, b1 + t_critical * se_b1]) + c2 = np.array([c1 - t_critical * se_c1, c1 + t_critical * se_c1]) + + # R-squared value + corr, _ = stats.pearsonr(np.squeeze(r), np.squeeze(y_pred)) + r2 = corr**2 + + return a1, a2, b1, b2, c1, c2, r2 + else: + a1=np.nan; a2=np.nan; b1=np.nan; b2=np.nan; c1=np.nan; c2=np.nan; r2=np.nan; + return a1, a2, b1, b2, c1, c2, r2 + +## function for calculate moving average! +def moving_avg(x_1d,n_mw): + nt=len(x_1d) + weight = np.ones((n_mw))/n_mw + n_delete=int((n_mw)/2); + smoothed=np.convolve(x_1d,weight,mode='same') + if n_delete != 0: + del_ind=np.concatenate((np.arange(n_delete),np.arange(nt-n_delete,nt))) + smoothed_final=np.delete(smoothed,del_ind) + return smoothed_final + +## function for ignoring NaN value for correlation +def corr_nan(temp1,temp2): + indnan=np.nanmean(temp1)+np.nanmean(temp2) + if ~np.isnan(indnan): + tempind=np.nonzero(np.isnan(temp1)|np.isnan(temp2)) + temp1=np.delete(temp1,tempind) + temp2=np.delete(temp2,tempind) + corr,_=stats.pearsonr(temp1,temp2) + return corr + else: + return np.nan + +## Calculate confidence interval of the data by t-test +def get_confidence_interval(data,alpha): + # alpha=0.9 --> 5%-95% confidence interval + # data is assumed to be not containing nan-value + mean = np.mean(data) + std_dev = np.std(data, ddof=1) # ddof=1 for sample standard deviation + # Number of samples + n = len(data) + # dof for t-distributions + dof = n - 1 + # Calculate the standard error of the mean + se = std_dev / np.sqrt(n) + # Calculate the t-score for the confidence level + t_score = stats.t.ppf((1 + alpha) / 2, dof) + # Calculate the margin of error + margin_of_error = t_score * se + return margin_of_error + +## find statistically significant r, given n +import numpy as np +import scipy.stats as stats +def significant_pearsonr(n,alpha): + # Degrees of freedom + df = n - 2 + # Critical t-value (two-tailed) + t_critical = stats.t.ppf(1 - alpha/2, df) + # Calculate critical r + r_critical = t_critical / np.sqrt(df + t_critical**2) + r_critical = np.abs(r_critical) # Consider absolute value for significance + return r_critical + +################################################################################ +### 1) Loading model data files and doing area average for each water-year ##### +################################################################################ +# +# The framework copies model data to a regular directory structure of the form +# //...nc +# Here and frequency are requested in the "varlist" part of +# settings.json. + +# The following command sets input_path to the value of the shell environment +# variable called TAS_FILE. This variable is set by the framework to let the +# script know where the locally downloaded copy of the data for this variable +# (which we called "tas") is. + +## load model deta +# data path from framework +pr_path = os.environ["PR_FILE"] +tas_path = os.environ["TAS_FILE"] +hfls_path = os.environ["HFLS_FILE"] +mrro_path = os.environ["MRRO_FILE"] + +# model grids (pr,tas,hfls,mrro must have same lon,lat grids) +lat_coor_name = os.environ["lat_coord"] +lon_coor_name = os.environ["lon_coord"] +time_coord_name = os.environ["time_coord"] +latm = nc.Dataset(pr_path).variables[lat_coor_name][:] +lonm = nc.Dataset(pr_path).variables[lon_coor_name][:] +time = nc.Dataset(pr_path).variables[time_coord_name][:] +nlat = len(latm) +nlon = len(lonm) +nt = len(time) +nyr = int(nt/12) +syr_model = int(os.environ["startdate"]) +eyr_model = int(os.environ["enddate"]) + +# variables +pr_var_name = os.environ["pr_var"] +tas_var_name = os.environ["tas_var"] +hfls_var_name = os.environ["hfls_var"] +mrro_var_name = os.environ["mrro_var"] +prm = nc.Dataset(pr_path).variables[pr_var_name][:,:,:].filled(fill_value=np.nan) +tasm = nc.Dataset(tas_path).variables[tas_var_name][:,:,:].filled(fill_value=np.nan) +hflsm = nc.Dataset(hfls_path).variables[hfls_var_name][:,:,:].filled(fill_value=np.nan) +mrrom = nc.Dataset(mrro_path).variables[mrro_var_name][:,:,:].filled(fill_value=np.nan) + +## change lon & variables to western degree, since the basin mask is not closed in eastern degree. +if np.mean(lonm) > 100: + # change lon to western degree + lonw = lonm - 180 + lonw[lonw < 0] = 1000 + ind180 = np.argmin(lonw) + lonwm = np.concatenate((lonm[ind180:] - 360, lonm[0:ind180])) + western_ind=np.concatenate((range(ind180, nlon), range(0, ind180))) + prm = prm[:, :, western_ind] + tasm = tasm[:, :, western_ind] + mrrom = mrrom[:, :, western_ind] + hflsm = hflsm[:, :, western_ind] + +## water year average (Oct-Sep water year for every basin) +## do time average and then do basin average (save computation time a lot) +syr_sens=1945 +eyr_sens=2014 +nyr = eyr_sens - syr_sens +inds=(syr_sens-syr_model)*12 + 10 - 1 +inde=inds + (nyr)*12 +prm_wy = np.nanmean(np.reshape(prm[inds:inde,:,:],(12,nyr,nlat,nlon),order="F"), axis=0) +tasm_wy = np.nanmean(np.reshape(tasm[inds:inde,:,:],(12,nyr,nlat,nlon),order="F"), axis=0) +mrrom_wy = np.nanmean(np.reshape(mrrom[inds:inde,:,:],(12,nyr,nlat,nlon),order="F"), axis=0) +hflsm_wy = np.nanmean(np.reshape(hflsm[inds:inde,:,:],(12,nyr,nlat,nlon),order="F"), axis=0) + +## area average for specific river basins +# load basin masks from shapefile provided by GRDC +shapefile_path = "{OBS_DATA}/mrb_basins.shp".format(**os.environ) +records=list(list(shpreader.Reader(shapefile_path).records())) + +# Select large basins based on the pre-defined indices (this reduces computation time) +bind = np.array([9, 12, 14, 19, 28, 29, 33, 34, 39, 42, 46, 54, 55, 62, 64, 65, 66, 71, 72, 73, 76, 79, 82, 83, 84, 85, 89, 90, 91, \ + 92, 94, 95, 99, 100, 101, 102, 103, 106, 110, 111, 113, 114, 116, 117, 119, 123, 124, 128, 129, 130, 138, 144, 148, 153, 160, 161, 163, 164, \ + 166, 169, 170, 172, 173, 174, 176, 178, 179, 181, 187, 194, 196, 211, 216, 218, 222, 223, 227, 229, 231, 232, 270, 276, 280, 283, 284, 293, 297, \ + 298, 299, 300, 306, 307, 312, 318, 328, 329, 340, 347, 358, 363, 364, 368, 372, 383, 393, 394, 399, 403, 409, 415, 416, 439, 443, 445, 452, 454, \ + 461, 463, 465, 467, 470, 485, 494, 495, 500, 501, 517, 519, 520, 379, 381]) + +# get basin masks/names for the slected ones +bind = bind - 1; +basin_points=[] +basin_names=[] +for i in bind: + geom = records[i].geometry + if geom.geom_type == 'Polygon': + coords = list(geom.exterior.coords) + basin_points.append(coords[0::10]) + elif geom.geom_type == 'MultiPolygon': + npoly=len(geom.geoms) + coords = list(geom.geoms[npoly-1].exterior.coords) + basin_points.append(coords[0::10]) + basin_names.append(records[i].attributes['RIVER_BASI']) + +nb=len(basin_points) + +# basin mask for each basin +basin_maskm = np.zeros((nlat, nlon, nb)) +lon2d, lat2d = np.meshgrid(lonwm, latm) +lon_points=np.reshape(lon2d,(nlat*nlon,1),order='F') +lat_points=np.reshape(lat2d,(nlat*nlon,1),order='F') +for b in range(nb): + path = mpltPath.Path(basin_points[b]) + inside=path.contains_points(np.column_stack([lon_points,lat_points])) + inside2d=np.reshape(inside,(nlat,nlon),order='F') + basin_maskm[:, :, b] = inside2d + +# area weight for each basin +aw = area_weight(latm, lonwm) +aw_basins = np.tile(aw[:,:,np.newaxis], (1, 1, nb)) * basin_maskm +# print(np.nansum(np.nansum(aw_basins,axis=1),axis=0)) +## note that basin 114 (PECHORA) is not captured by CESM2 grid -> area is 0 + +# basin average +awm = np.transpose(np.tile(aw_basins[:,:,:,np.newaxis],(1,1,1,nyr)),(3,0,1,2)) + +prb_wy = np.nansum(np.nansum( \ + np.tile(prm_wy[:, :, :, np.newaxis], (1, 1, 1, nb)) * awm \ + , axis=1), axis=1) / np.nansum(np.nansum(awm, axis=1), axis=1) + +tasb_wy = np.nansum(np.nansum( \ + np.tile(tasm_wy[:, :, :, np.newaxis], (1, 1, 1, nb)) * awm \ + , axis=1), axis=1) / np.nansum(np.nansum(awm, axis=1), axis=1) + +mrrob_wy = np.nansum(np.nansum( \ + np.tile(mrrom_wy[:, :, :, np.newaxis], (1, 1, 1, nb)) * awm \ + , axis=1), axis=1) / np.nansum(np.nansum(awm, axis=1), axis=1) + +hflsb_wy = np.nansum(np.nansum( \ + np.tile(hflsm_wy[:, :, :, np.newaxis], (1, 1, 1, nb)) * awm \ + , axis=1), axis=1) / np.nansum(np.nansum(awm, axis=1), axis=1) + + + +################################################################################ +### 2) check water budget closure, negative values, consistent values: ############################################# +################################################################################ +## sanity check to ensure that runoff output is not strange + +## Note that CESM is almost free from this problem when we analyzed CEMS1 and CESM2 + +## some CMIP model's mrro (runoff) shows some erronous behaviour for some basins +# 1. water budget is not closed; P(precip)-E(evap) is way different from Q(runoff) +# 2. runoff is negative in some basins, leading to negative or very small climatological runoff +# 3. same and small runoff value is repeated for long times in some dry basins + +## 1. check the water budget closure +pmeb_model = np.nanmean(prb_wy - hflsb_wy/(2.5*1e6), axis=0) +error_val = (pmeb_model - np.nanmean(mrrob_wy, axis=0)) / np.nanmean(mrrob_wy, axis=0) +pval = 0.05 +closed_fraction=np.nansum( np.abs(error_val)<0.05 ,axis=0) / nb +if closed_fraction<0.8: + raise ValueError('water budget is not closed more than 80% of ' + str(nb) + 'basins') + + + +############################################################################## +## 3) calculate runoff sensitivities ########## +############################################################################## +# initialization +psens_model = np.full((nb, 3),np.nan) +tsens_model = np.full((nb, 3),np.nan) +intsens_model = np.full((nb, 3),np.nan) +r2_mrro_wy_int = np.full((nb),np.nan) +r2_model_int = np.full((nb),np.nan) +r2_model_pred = np.full((nb),np.nan) +r2_model_pred_int = np.full((nb),np.nan) + +# water year average (Oct-Sep water year for every basin) +# after 5-year average, the definition of differnet water year for each basin does not matter much +mw=5 # 5-year moving window +for b in range(nb): + # specify start and end indices + + # smoothing (5-year moving average) + mrrobym = moving_avg(mrrob_wy[:,b],mw) + prbym = moving_avg(prb_wy[:,b],mw) + tasbym = moving_avg(tasb_wy[:,b],mw) + + # get anomalies + amrrobym = mrrobym - np.nanmean(mrrobym) + aprbym = prbym - np.nanmean(prbym) + atasbym = tasbym - np.nanmean(tasbym) + + # get % anomalies + amrrobym_pct = amrrobym / np.nanmean(mrrobym) * 100 + aprbym_pct = aprbym / np.nanmean(prbym) * 100 + + # calculate runoff sensitivities using def runoff_sens_reg(r, p, t): + a1, a2, b1, b2, c1, c2, r2 = runoff_sens_reg(amrrobym_pct, aprbym_pct, atasbym, mw=5, alpha=0.1) + psens_model[b,0] = a1 + psens_model[b,1:3] = a2 + tsens_model[b,0] = b1 + tsens_model[b,1:3] = b2 + intsens_model[b,0] = c1 + intsens_model[b,1:3] = c2 + r2_mrro_wy_int[b] = r2 + + # prediction and resulting R2 (training accuracy) + amrrobym_pct_pred_model = a1 * aprbym_pct + b1 * atasbym + amrrobym_pct_pred_int_model = a1 * aprbym_pct + b1 * atasbym + c1 * atasbym * aprbym_pct + corr = corr_nan(amrrobym_pct_pred_model, amrrobym_pct_pred_int_model) + r2_model_int[b] = corr**2 + corr = corr_nan(amrrobym_pct_pred_model, amrrobym_pct) + r2_model_pred[b] = corr**2 + corr = corr_nan(amrrobym_pct_pred_int_model, amrrobym_pct) + r2_model_pred_int[b] = corr**2 + +# display warning if the interaction term is affecting the regression models +# role of interaction term (P*T) is usually negligible, consistent to Lehner et al. 2019 +condition=np.asarray(r2_model_int<0.8) +# condition=np.asarray(np.abs(r2_model_pred_int-r2_model_pred)>0.2) +for b in range(b): + if condition[b]: + warning_message = f'Warning: Interaction term may matter for model + basin #{b+1} (r={r2_model_int[b]:.2f})' + warnings.warn(warning_message) + + + +##################################################### +### 4) Saving output data: ########################## +##################################################### +## save runoff sensitivity data in netCDF4 file +out_path = "{WORK_DIR}/model/netCDF/runoff_sensitivities_{CASENAME}.nc".format(**os.environ) +# Create a NetCDF dataset +dataset = nc.Dataset(out_path, 'w') +# Create dimensions +basin_dim = tsens_model.shape[0] +sens_dim = tsens_model.shape[1] +dataset.createDimension("basin", basin_dim) +dataset.createDimension("sens_type", sens_dim) +string_length = 45 +dataset.createDimension("string_len", string_length) +# Create variables +basin_var = dataset.createVariable('basin', 'S1', ('basin', 'string_len')) +tsens_model_var = dataset.createVariable('tsens_model', np.float64, ('basin', 'sens_type'), fill_value=-9999) +psens_model_var = dataset.createVariable('psens_model', np.float64, ('basin', 'sens_type'), fill_value=-9999) +r2_model_pred_var = dataset.createVariable('r2_model_pred', np.float64, ('basin'), fill_value=-9999) +# assign attributes +basin_var.long_name = 'target basins' +basin_var.reference = 'ref: https://grdc.bafg.de/GRDC/EN/02_srvcs/22_gslrs/221_MRB/riverbasins_node.html' +tsens_model_var.long_name = f'temperature sensitivity ({syr_sens}-{eyr_sens})' +tsens_model_var.units = '%/K' +tsens_model_var.sens_type_index = '1: sensitivity value / 2,3: 95% confidence interval' +psens_model_var.long_name = f'precipitation sensitivity ({syr_sens}-{eyr_sens})' +psens_model_var.units = '%/%' +psens_model_var.sens_type_index = '1: sensitivity value / 2,3: 95% confidence interval' +r2_model_pred_var.long_name = f'regression model accuracy (R^2) ({syr_sens}-{eyr_sens})' +r2_model_pred_var.units = 'no units' +# assign variables +for i, s in enumerate(basin_names): + basin_var[i, :len(s)] = np.array(list(s), dtype='S1') +tsens_model_var[:] = tsens_model +psens_model_var[:] = psens_model +r2_model_pred_var[:] = r2_model_pred +# close and save +dataset.close() + + +################################################################################ +### 5) Loading digested data & plotting obs figures: ########################## +################################################################################ +## load pre-calculated obs and cmip data +# load obs +obs_path = "{OBS_DATA}/runoff_sensitivity_obs.nc".format(**os.environ) +tsens_obs = nc.Dataset(obs_path)['tsens_obs'][:,:,0:nb,:] +psens_obs = nc.Dataset(obs_path)['psens_obs'][:,:,0:nb,:] +r2_obs_pred = nc.Dataset(obs_path)['r2_obs_pred'][:,:,0:nb] +# load cmip5 +hist5_path = "{OBS_DATA}/runoff_sensitivity_hist5.nc".format(**os.environ) +tsens_hist5 = nc.Dataset(hist5_path)['tsens_hist5'][:,0:nb,:] +psens_hist5 = nc.Dataset(hist5_path)['psens_hist5'][:,0:nb,:] +r2_hist5_pred = nc.Dataset(hist5_path)['r2_hist5_pred'][:,0:nb] +# load cmip6 +hist6_path = "{OBS_DATA}/runoff_sensitivity_hist6.nc".format(**os.environ) +tsens_hist6 = nc.Dataset(hist6_path)['tsens_hist6'][:,0:nb,:] +psens_hist6 = nc.Dataset(hist6_path)['psens_hist6'][:,0:nb,:] +r2_hist6_pred = nc.Dataset(hist6_path)['r2_hist6_pred'][:,0:nb] + + +############################################# +## draw maps with filled river basins ####### +############################################# +def plot_and_save_basin_filled(values, basin_points, color_bins, color_bins2, plt_colormap, plt_unit, plt_title, plt_path, coast_path): + ## data needed for plotting + # assign color index to each values + values_ind=np.digitize(values, color_bins2)-1 + # make colormap + custom_cmap=colors.LinearSegmentedColormap.from_list('custom_colormap',plt_colormap, N=len(color_bins)-1) + # load coastline + lonmap = nc.Dataset(coast_path)['lonmap'][0:9850] + latmap = nc.Dataset(coast_path)['latmap'][0:9850] + ## draw figure + fig, ax = plt.subplots(figsize=(12, 4.8)) + plt.rcParams.update({'font.size': 12}) + # coastline + ax.plot(lonmap, latmap, color=[0.8, 0.8, 0.8, 0], linewidth=0.5) + ax.add_patch(patches.Polygon(np.column_stack([lonmap, latmap]), closed=True, facecolor=(0.8, 0.8, 0.8), edgecolor=(0.5, 0.5, 0.5), linestyle='none')) + # fill basins with colors corresponding to target values + for b in range(nb): + X = [item[0] for item in basin_points[b]] + Y = [item[1] for item in basin_points[b]] + X = X[0:-int(np.floor(len(X)/100))] + Y = Y[0:-int(np.floor(len(Y)/100))] + if not np.isnan(values_ind[b]): + ax.add_patch(patches.Polygon(np.column_stack([X, Y]), closed=True, facecolor=custom_cmap(values_ind[b]), edgecolor=(0.5, 0.5, 0.5), linewidth=0.5)) + else: + ax.plot(X, Y, color=(0.5, 0.5, 0.5), linewidth=0.5) + # Set colormap and colorbar + cb = plt.colorbar(plt.cm.ScalarMappable(norm=plt.Normalize(0, 1), cmap=custom_cmap)) + cb.set_ticks(np.linspace(0, 1, len(color_bins))) + clabels = [str(b) for b in color_bins] + cb.set_ticklabels(clabels) + cb.set_label(plt_unit, fontsize=12) + # Customize the ticks and labels + ax.set_xticks(range(-180, 181, 30)) + ax.set_xticklabels(['', '', '120°W', '', '60°W', '', '0°', '', '60°E', '', '120°E', '', '']) + ax.set_yticks(range(-90, 91, 15)) + ax.set_yticklabels(['', '', '60°S', '', '30°S', '', '0°', '', '30°N', '', '60°N', '', '']) + ax.set_xlim([-180, 180]) + ax.set_ylim([-60, 80]) + # Set title + ax.set_title(plt_title, fontsize=17) + # Save the figure as a eps + plt.savefig(plt_path, bbox_inches='tight') + plt.close() + + + +## figures for T sensitivity ## +bins = [-30,-15,-12,-9,-6,-3,0,3,6,9,12,15,30] +bins2 = [-60,-15,-12,-9,-6,-3,0,3,6,9,12,15,60] +plot_colormap = [(0.4000, 0, 0, 1),(0.7706, 0, 0, 1),(0.9945, 0.0685, 0.0173, 1),(0.9799, 0.2483, 0.0627, 1),(0.9715, 0.4442, 0.0890, 1),(0.9845, 0.6961, 0.0487, 1),(0.9973, 0.9480, 0.0083, 1),(1.0000, 1.0000, 0.3676, 1),(1.0000, 1.0000, 1.0000, 1),(1.0000, 1.0000, 1.0000, 1),(0.6975, 0.8475, 0.9306, 1),(0.4759, 0.7358, 0.8797, 1),(0.2542, 0.6240, 0.8289, 1),(0.0436, 0.5130, 0.7774, 1),(0.0533, 0.4172, 0.7138, 1),(0.0630, 0.3215, 0.6503, 1),(0.0411, 0.1760, 0.5397, 1),(0, 0, 0.4000, 1)] +plot_unit='[%/K]' +coast_path = "{OBS_DATA}/coastline.nc".format(**os.environ) + +# OBS +obs_val=np.nanmean(np.nanmean(tsens_obs[:,:,:,0],axis=0),axis=0) +plot_title='T sensitivity: OBS' +plot_path = "{WK_DIR}/obs/PS/tsens_obs.eps".format(**os.environ) +plot_and_save_basin_filled(obs_val, basin_points, bins, bins2, \ + plot_colormap, plot_unit, plot_title, plot_path, coast_path) + +# model +values=tsens_model[:,0] +plot_title=f'T sensitivity: MODEL (WBC={closed_fraction:0.2f})' +plot_path= "{WK_DIR}/model/PS/tsens_model.eps".format(**os.environ) +plot_and_save_basin_filled(values, basin_points, bins, bins2, \ + plot_colormap, plot_unit, plot_title, plot_path, coast_path) + +# model bias +values=tsens_model[:,0]-obs_val +plot_title='T sensitivity bias: MODEL - OBS' +plot_path= "{WK_DIR}/model/PS/tsens_model_bias.eps".format(**os.environ) +plot_and_save_basin_filled(values, basin_points, bins, bins2, \ + plot_colormap, plot_unit, plot_title, plot_path, coast_path) + +# hist6 bias +values=np.nanmean(tsens_hist6[:,:,0],axis=0)-obs_val +plot_title='T sensitivity bias: CMIP6 - OBS' +plot_path= "{WK_DIR}/obs/PS/tsens_hist6_bias.eps".format(**os.environ) +plot_and_save_basin_filled(values, basin_points, bins, bins2, \ + plot_colormap, plot_unit, plot_title, plot_path, coast_path) + + +## figures for P sensitivity ## +bins=[-3.5, -3, -2.5, -2, -1.6, -1.2, -0.8, -0.4, 0, 0.4, 0.8, 1.2, 1.6, 2, 2.5, 3, 3.5]; +bins2=[-100, -3, -2.5, -2, -1.6, -1.2, -0.8, -0.4, 0, 0.4, 0.8, 1.2, 1.6, 2, 2.5, 3, 100]; +plot_colormap = [(0.4000, 0, 0, 1), (0.7316, 0, 0, 1), (0.9975, 0.0306, 0.0078, 1), (0.9845, 0.1915, 0.0484, 1), (0.9715, 0.3525, 0.0890, 1), (0.9776, 0.5636, 0.0700, 1), (0.9891, 0.7889, 0.0338, 1), (1.0000, 1.0000, 0.0263, 1), (1.0000, 1.0000, 0.4408, 1), (1.0000, 1.0000, 1.0000, 1), (1.0000, 1.0000, 1.0000, 1), (0.7325, 0.8651, 0.9386, 1), (0.5342, 0.7652, 0.8931, 1), (0.3359, 0.6652, 0.8476, 1), (0.1375, 0.5652, 0.8020, 1), (0.0477, 0.4727, 0.7506, 1), (0.0563, 0.3870, 0.6938, 1), (0.0651, 0.3013, 0.6370, 1), (0.0368, 0.1575, 0.5250, 1), (0, 0, 0.4000, 1)] +plot_colormap = [(0.4000, 0, 0, 1), (0.7316, 0, 0, 1), (0.9975, 0.0306, 0.0078, 1), (0.9845, 0.1915, 0.0484, 1), (0.9715, 0.3525, 0.0890, 1), (0.9776, 0.5636, 0.0700, 1), (0.9891, 0.7889, 0.0338, 1), (1.0000, 1.0000, 0.0263, 1), (1.0000, 1.0000, 0.4408, 1), (1.0000, 1.0000, 1.0000, 1), (1.0000, 1.0000, 1.0000, 1), (0.7325, 0.8651, 0.9386, 1), (0.5342, 0.7652, 0.8931, 1), (0.3359, 0.6652, 0.8476, 1), (0.1375, 0.5652, 0.8020, 1), (0.0477, 0.4727, 0.7506, 1), (0.0563, 0.3870, 0.6938, 1), (0.0651, 0.3013, 0.6370, 1), (0.0368, 0.1575, 0.5250, 1), (0, 0, 0.4000, 1)] +plot_unit='[%/%]' + +# OBS +obs_val=np.nanmean(np.nanmean(psens_obs[:,:,:,0],axis=0),axis=0) +plot_title='P sensitivity: OBS' +plot_path= "{WK_DIR}/obs/PS/psens_obs.eps".format(**os.environ) +plot_and_save_basin_filled(obs_val, basin_points, bins, bins2, \ + plot_colormap, plot_unit, plot_title, plot_path, coast_path) + +# model +values=psens_model[:,0] +plot_title=f'P sensitivity: MODEL (WBC={closed_fraction:0.2f})' +plot_path= "{WK_DIR}/model/PS/psens_model.eps".format(**os.environ) +plot_and_save_basin_filled(values, basin_points, bins, bins2, \ + plot_colormap, plot_unit, plot_title, plot_path, coast_path) + +# model bias +values=psens_model[:,0]-obs_val +plot_title='P sensitivity bias: MODEL - OBS' +plot_path= "{WK_DIR}/model/PS/psens_model_bias.eps".format(**os.environ) +plot_and_save_basin_filled(values, basin_points, bins, bins2, \ + plot_colormap, plot_unit, plot_title, plot_path, coast_path) + +# hist6 bias +values=np.nanmean(psens_hist6[:,:,0],axis=0)-obs_val +plot_title='P sensitivity bias: CMIP6 - OBS' +plot_path= "{WK_DIR}/obs/PS/psens_hist6_bias.eps".format(**os.environ) +plot_and_save_basin_filled(values, basin_points, bins, bins2, \ + plot_colormap, plot_unit, plot_title, plot_path, coast_path) + + +########################################### +## summary figures for target variables ## +########################################### +def plot_and_save_basin_filled_summary(values_obs, values_model, values_hist5, values_hist6, basin_points, color_bins, color_bins2, plt_colormap, plt_unit, plt_title, plt_path, coast_path, closed_fraction): + # load coastline + lonmap = nc.Dataset(coast_path)['lonmap'][0:9850] + latmap = nc.Dataset(coast_path)['latmap'][0:9850] + # colormap + custom_cmap=colors.LinearSegmentedColormap.from_list('custom_colormap',plt_colormap, N=len(bins)-1) + + ## draw figure + fig = plt.figure(figsize=(17, 10)) + plt.rcParams.update({'font.size': 12}) # Change the font size to 12 + # OBS + values_ind=np.digitize(values_obs,color_bins2)-1 + h1 = fig.add_axes([0.06, 0.58, 0.43, 0.32]) + h1.plot(lonmap, latmap, color=[0.8, 0.8, 0.8, 0], linewidth=0.5) + h1.add_patch(patches.Polygon(np.column_stack([lonmap, latmap]), closed=True, facecolor=(0.8, 0.8, 0.8), edgecolor=(0.5, 0.5, 0.5), linestyle='none')) + for b in range(nb): + X = [item[0] for item in basin_points[b]] + Y = [item[1] for item in basin_points[b]] + X = X[0:-int(np.floor(len(X)/100))] + Y = Y[0:-int(np.floor(len(Y)/100))] + if not np.isnan(values_ind[b]): + h1.add_patch(patches.Polygon(np.column_stack([X, Y]), closed=True, facecolor=custom_cmap(values_ind[b]), edgecolor=(0.5, 0.5, 0.5), linewidth=0.5)) + else: + h1.plot(X, Y, color=(0.5, 0.5, 0.5), linewidth=0.5) + h1.set_xticks(range(-180, 181, 30)) + h1.set_xticklabels(['', '', '120°W', '', '60°W', '', '0°', '', '60°E', '', '120°E', '', '']) + h1.set_yticks(range(-90, 91, 15)) + h1.set_yticklabels(['', '', '60°S', '', '30°S', '', '0°', '', '30°N', '', '60°N', '', '']) + h1.set_xlim([-180, 180]) + h1.set_ylim([-60, 80]) + h1.set_title('OBS', fontsize=17) + plt.text(195,105,plt_title, fontsize=20, ha='center', fontweight='bold') + # model + values_ind=np.digitize(values_model,color_bins2)-1 + h2 = fig.add_axes([0.54, 0.58, 0.43, 0.32]) + h2.plot(lonmap, latmap, color=[0.8, 0.8, 0.8, 0], linewidth=0.5) + h2.add_patch(patches.Polygon(np.column_stack([lonmap, latmap]), closed=True, facecolor=(0.8, 0.8, 0.8), edgecolor=(0.5, 0.5, 0.5), linestyle='none')) + for b in range(nb): + X = [item[0] for item in basin_points[b]] + Y = [item[1] for item in basin_points[b]] + X = X[0:-int(np.floor(len(X)/100))] + Y = Y[0:-int(np.floor(len(Y)/100))] + if not np.isnan(values_ind[b]): + h2.add_patch(patches.Polygon(np.column_stack([X, Y]), closed=True, facecolor=custom_cmap(values_ind[b]), edgecolor=(0.5, 0.5, 0.5), linewidth=0.5)) + else: + h2.plot(X, Y, color=(0.5, 0.5, 0.5), linewidth=0.5) + h2.set_xticks(range(-180, 181, 30)) + h2.set_xticklabels(['', '', '120°W', '', '60°W', '', '0°', '', '60°E', '', '120°E', '', '']) + h2.set_yticks(range(-90, 91, 15)) + h2.set_yticklabels(['', '', '60°S', '', '30°S', '', '0°', '', '30°N', '', '60°N', '', '']) + h2.set_xlim([-180, 180]) + h2.set_ylim([-60, 80]) + h2.set_title(f'MODEL (WBC={closed_fraction:0.2f})', fontsize=17) + # hist5 + values_ind=np.digitize(values_hist5,color_bins2)-1 + h3 = fig.add_axes([0.06, 0.17, 0.43, 0.32]) + h3.plot(lonmap, latmap, color=[0.8, 0.8, 0.8, 0], linewidth=0.5) + h3.add_patch(patches.Polygon(np.column_stack([lonmap, latmap]), closed=True, facecolor=(0.8, 0.8, 0.8), edgecolor=(0.5, 0.5, 0.5), linestyle='none')) + for b in range(nb): + X = [item[0] for item in basin_points[b]] + Y = [item[1] for item in basin_points[b]] + X = X[0:-int(np.floor(len(X)/100))] + Y = Y[0:-int(np.floor(len(Y)/100))] + if not np.isnan(values_ind[b]): + h3.add_patch(patches.Polygon(np.column_stack([X, Y]), closed=True, facecolor=custom_cmap(values_ind[b]), edgecolor=(0.5, 0.5, 0.5), linewidth=0.5)) + else: + h3.plot(X, Y, color=(0.5, 0.5, 0.5), linewidth=0.5) + h3.set_xticks(range(-180, 181, 30)) + h3.set_xticklabels(['', '', '120°W', '', '60°W', '', '0°', '', '60°E', '', '120°E', '', '']) + h3.set_yticks(range(-90, 91, 15)) + h3.set_yticklabels(['', '', '60°S', '', '30°S', '', '0°', '', '30°N', '', '60°N', '', '']) + h3.set_xlim([-180, 180]) + h3.set_ylim([-60, 80]) + h3.set_title('CMIP5 MMM (21 models)', fontsize=17) + # hist6 + values_ind=np.digitize(values_hist6,color_bins2)-1 + h4 = fig.add_axes([0.54, 0.17, 0.43, 0.32]) + h4.plot(lonmap, latmap, color=[0.8, 0.8, 0.8, 0], linewidth=0.5) + h4.add_patch(patches.Polygon(np.column_stack([lonmap, latmap]), closed=True, facecolor=(0.8, 0.8, 0.8), edgecolor=(0.5, 0.5, 0.5), linestyle='none')) + for b in range(nb): + X = [item[0] for item in basin_points[b]] + Y = [item[1] for item in basin_points[b]] + X = X[0:-int(np.floor(len(X)/100))] + Y = Y[0:-int(np.floor(len(Y)/100))] + if not np.isnan(values_ind[b]): + h4.add_patch(patches.Polygon(np.column_stack([X, Y]), closed=True, facecolor=custom_cmap(values_ind[b]), edgecolor=(0.5, 0.5, 0.5), linewidth=0.5)) + else: + h4.plot(X, Y, color=(0.5, 0.5, 0.5), linewidth=0.5) + h4.set_xticks(range(-180, 181, 30)) + h4.set_xticklabels(['', '', '120°W', '', '60°W', '', '0°', '', '60°E', '', '120°E', '', '']) + h4.set_yticks(range(-90, 91, 15)) + h4.set_yticklabels(['', '', '60°S', '', '30°S', '', '0°', '', '30°N', '', '60°N', '', '']) + h4.set_xlim([-180, 180]) + h4.set_ylim([-60, 80]) + h4.set_title('CMIP6 MMM (26 models)', fontsize=17) + # Colorbar + cax = fig.add_axes([0.06, 0.085, 0.91, 0.025]) + sm=plt.cm.ScalarMappable(norm=plt.Normalize(0, 1), cmap=custom_cmap) + cb = plt.colorbar(sm,cax=cax, cmap=custom_cmap, orientation='horizontal') + cb.set_ticks(np.linspace(0, 1, len(color_bins))) + cb.set_ticklabels([str(b) for b in color_bins]) + cb.ax.set_xlabel(plt_unit, fontsize=12) + # Save the figure + plt.savefig(plt_path) + plt.close() + +# summary for training accuracy +bins=[0, 0.16, 0.36, 0.49, 0.64, 0.81, 1] +bins2=[0, 0.16, 0.36, 0.49, 0.64, 0.81, 1] +plt_colormap=[(1.0000, 1.0000, 1.0000), (0.9961, 0.9020, 0), (1.0000, 0.6980, 0), (1.0000, 0.2980, 0.0039), (0.8471, 0, 0.0039), (0.5882, 0.0196, 0.0196)] +values_obs=np.nanmean(np.nanmean(r2_obs_pred[:,:,:], axis=0), axis=0) +values_model=r2_model_pred[:] +values_hist5=np.nanmean(r2_hist5_pred[:,:], axis=0) +values_hist6=np.nanmean(r2_hist6_pred[:,:], axis=0) +plt_title='Summary: training accuracy' +plt_unit='[R$^2$]' +plt_path="{WK_DIR}/model/PS/summary_prediction_accuracy.eps".format(**os.environ) +plot_and_save_basin_filled_summary(values_obs, values_model, values_hist5, values_hist6, \ + basin_points, bins, bins2, plt_colormap, plt_unit, plt_title, plt_path,\ + coast_path, closed_fraction) + +# summary for T sensitivity +bins = [-30,-15,-12,-9,-6,-3,0,3,6,9,12,15,30] +bins2 = [-60,-15,-12,-9,-6,-3,0,3,6,9,12,15,60] +plt_colormap=[(0.4000, 0, 0, 1),(0.7706, 0, 0, 1),(0.9945, 0.0685, 0.0173, 1),(0.9799, 0.2483, 0.0627, 1),(0.9715, 0.4442, 0.0890, 1),(0.9845, 0.6961, 0.0487, 1),(0.9973, 0.9480, 0.0083, 1),(1.0000, 1.0000, 0.3676, 1),(1.0000, 1.0000, 1.0000, 1),(1.0000, 1.0000, 1.0000, 1),(0.6975, 0.8475, 0.9306, 1),(0.4759, 0.7358, 0.8797, 1),(0.2542, 0.6240, 0.8289, 1),(0.0436, 0.5130, 0.7774, 1),(0.0533, 0.4172, 0.7138, 1),(0.0630, 0.3215, 0.6503, 1),(0.0411, 0.1760, 0.5397, 1),(0, 0, 0.4000, 1)] +values_obs=np.nanmean(np.nanmean(tsens_obs[:, :, :, 0], axis=0), axis=0) +values_model=tsens_model[:, 0] +values_hist5=np.nanmean(tsens_hist5[:, :, 0], axis=0) +values_hist6=np.nanmean(tsens_hist6[:, :, 0], axis=0) +plt_title='Summary: T sensitivity' +plt_unit='[%/K]' +plt_path="{WK_DIR}/model/PS/summary_tsens.eps".format(**os.environ) +plot_and_save_basin_filled_summary(values_obs, values_model, values_hist5, values_hist6, \ + basin_points, bins, bins2, plt_colormap, plt_unit, plt_title, plt_path,\ + coast_path, closed_fraction) + +# summary for P sensitivity +bins=[-3.5, -3, -2.5, -2, -1.6, -1.2, -0.8, -0.4, 0, 0.4, 0.8, 1.2, 1.6, 2, 2.5, 3, 3.5]; +bins2=[-100, -3, -2.5, -2, -1.6, -1.2, -0.8, -0.4, 0, 0.4, 0.8, 1.2, 1.6, 2, 2.5, 3, 100]; +plt_colormap=[(0.4000, 0, 0, 1), (0.7316, 0, 0, 1), (0.9975, 0.0306, 0.0078, 1), (0.9845, 0.1915, 0.0484, 1), (0.9715, 0.3525, 0.0890, 1), (0.9776, 0.5636, 0.0700, 1), (0.9891, 0.7889, 0.0338, 1), (1.0000, 1.0000, 0.0263, 1), (1.0000, 1.0000, 0.4408, 1), (1.0000, 1.0000, 1.0000, 1), (1.0000, 1.0000, 1.0000, 1), (0.7325, 0.8651, 0.9386, 1), (0.5342, 0.7652, 0.8931, 1), (0.3359, 0.6652, 0.8476, 1), (0.1375, 0.5652, 0.8020, 1), (0.0477, 0.4727, 0.7506, 1), (0.0563, 0.3870, 0.6938, 1), (0.0651, 0.3013, 0.6370, 1), (0.0368, 0.1575, 0.5250, 1), (0, 0, 0.4000, 1)] +values_obs=np.nanmean(np.nanmean(psens_obs[:, :, :, 0], axis=0), axis=0) +values_model=psens_model[:, 0] +values_hist5=np.nanmean(psens_hist5[:, :, 0], axis=0) +values_hist6=np.nanmean(psens_hist6[:, :, 0], axis=0) +plt_title='Summary: P sensitivity' +plt_unit='[%/%]' +plt_path="{WK_DIR}/model/PS/summary_psens.eps".format(**os.environ) +plot_and_save_basin_filled_summary(values_obs, values_model, values_hist5, values_hist6, \ + basin_points, bins, bins2, plt_colormap, plt_unit, plt_title, plt_path,\ + coast_path, closed_fraction) + + + + +########################################### +## significant test ## +########################################### +## Are the runoff sensitivities reliable? (Prediction skill of the regression model itself) +# observations, CMIP6 models, target model +nfoc=r2_obs_pred.shape[0]; nens=r2_obs_pred.shape[1]; +hval_r2_ia_gobs=np.full(nb,False); +hval_r2_ia_hist6=np.full(nb,False); +hval_r2_ia_model=np.full(nb,False); +r2_ia_sig=significant_pearsonr(np.floor(70/5),0.05)**2 +hval_r2_ia_gobs = np.sum(np.reshape(r2_obs_pred,(nfoc*nens,nb),order='F')>r2_ia_sig,axis=0)==nfoc*nens +nmodel6=r2_hist6_pred.shape[0] +for b in range(nb): + hval_r2_ia_hist6[b] = np.sum(r2_hist6_pred[:,b]>r2_ia_sig,axis=0) == np.nansum(~np.isnan(r2_hist6_pred[:,b])) +hval_r2_ia_model = r2_model_pred > r2_ia_sig +hval_r2_ia = hval_r2_ia_gobs & hval_r2_ia_hist6 & hval_r2_ia_model +ind_r2_sig = np.nonzero(hval_r2_ia)[0] + +## Are the runoff sensitivity bias large enough to overcome the uncertainty? +# confidence interval of observational uncertainty +pct90_tsens_obs=np.full((nb,2),np.nan) +pct90_psens_obs=np.full((nb,2),np.nan) +nobs=nfoc*nens +for b in range(nb): + tsens=np.reshape(tsens_obs[:,:,b,0],(nobs),order='F') + pct90_tsens_obs[b,:] = (np.percentile(tsens,5),np.percentile(tsens,95)) + psens=np.reshape(psens_obs[:,:,b,0],(nobs),order='F') + pct90_psens_obs[b,:] = (np.percentile(psens,5),np.percentile(psens,95)) + +# compare the observational uncertainty to the pre-calculated uncertainty in regression coefficient +def range_difference(a1,a2,b1,b2): + return (a2