diff --git a/Notebooks/eval_precipitation.ipynb b/Notebooks/eval_precipitation.ipynb
index b1ea5e217ce0f05df12783555122bab52e44ddc2..4bdac8d9b968b1fcbeee519b34fb67684df33631 100644
--- a/Notebooks/eval_precipitation.ipynb
+++ b/Notebooks/eval_precipitation.ipynb
@@ -59,7 +59,10 @@
     "SPREDICTORS = 'p'\n",
     "DEM = 'dem'\n",
     "DEM_FEATURES = ''\n",
-    "DOY = 'doy'"
+    "DOY = 'doy'\n",
+    "# WET_DAY_THRESHOLD = 1\n",
+    "# LOSS = 'MSELoss'\n",
+    "# LOSS = 'BernoulliGammaLoss'"
    ]
   },
   {
@@ -86,12 +89,13 @@
     "import xarray as xr\n",
     "import numpy as np\n",
     "import matplotlib.pyplot as plt\n",
+    "from mpl_toolkits.axes_grid1.inset_locator import inset_axes\n",
     "import scipy.stats as stats\n",
     "from IPython.display import Image\n",
     "from sklearn.metrics import r2_score, roc_curve, auc, classification_report\n",
     "\n",
     "# locals\n",
-    "from climax.main.io import ERA5_PATH, OBS_PATH, TARGET_PATH\n",
+    "from climax.main.io import ERA5_PATH, OBS_PATH, TARGET_PATH, DEM_PATH\n",
     "from pysegcnn.core.utils import search_files\n",
     "from pysegcnn.core.graphics import plot_classification_report"
    ]
@@ -193,7 +197,21 @@
     "PATTERN = '_'.join([PATTERN, SPREDICTORS]) if SPREDICTORS else PATTERN\n",
     "PATTERN = '_'.join([PATTERN, DEM]) if DEM else PATTERN\n",
     "PATTERN = '_'.join([PATTERN, DEM_FEATURES]) if DEM_FEATURES else PATTERN\n",
-    "PATTERN = '_'.join([PATTERN, DOY]) if DOY else PATTERN"
+    "PATTERN = '_'.join([PATTERN, DOY]) if DOY else PATTERN\n",
+    "# PATTERN = '_'.join([PATTERN, '{:0d}mm_{}'.format(WET_DAY_THRESHOLD, LOSS)])\n",
+    "PATTERN"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "ecaba394-f802-481f-8274-c44e2f5fdf1a",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# digital elevation model\n",
+    "dem = xr.open_dataset(search_files(DEM_PATH, 'eu_dem_v11_stt.nc').pop())\n",
+    "dem = dem.Band1.to_dataset().rename({'Band1': 'elevation'})"
    ]
   },
   {
@@ -207,7 +225,18 @@
    "source": [
     "# model predictions and observations NetCDF \n",
     "y_pred = xr.open_dataset(search_files(TARGET_PATH.joinpath(PREDICTAND), '.'.join([PATTERN, 'nc$'])).pop())\n",
-    "y_true = xr.open_dataset(search_files(OBS_PATH.joinpath(PREDICTAND), '.nc$').pop())"
+    "y_true = xr.open_dataset(search_files(OBS_PATH.joinpath(PREDICTAND), 'OBS_pr(.*).nc$').pop())"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "325f3086-85f0-4c28-b37d-7370d2d92405",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# reference dataset: ERA-5 precipitation\n",
+    "y_refe = xr.open_dataset(search_files(ERA5_PATH.joinpath('ERA5', 'total_precipitation'), '.nc$').pop())"
    ]
   },
   {
@@ -218,7 +247,9 @@
    "outputs": [],
    "source": [
     "# subset to time period covered by predictions\n",
-    "y_true = y_true.sel(time=y_pred.time)  "
+    "y_true = y_true.sel(time=y_pred.time)\n",
+    "y_refe = y_refe.sel(time=y_pred.time).drop_vars('lambert_azimuthal_equal_area')\n",
+    "y_refe = y_refe.rename({'tp': 'precipitation'})"
    ]
   },
   {
@@ -229,9 +260,22 @@
    "outputs": [],
    "source": [
     "# align datasets and mask missing values in model predictions\n",
-    "y_true, y_pred_pr, y_pred_prob = xr.align(y_true, y_pred.precipitation.to_dataset(), y_pred.prob.to_dataset(), join='override')\n",
+    "y_true, y_refe, y_pred_pr, y_pred_prob = xr.align(y_true, y_refe, y_pred.precipitation.to_dataset(), y_pred.prob.to_dataset(), join='override')\n",
     "y_pred_pr = y_pred_pr.where(~np.isnan(y_true.precipitation), other=np.nan)\n",
-    "y_pred_prob = y_pred_prob.where(~np.isnan(y_true.precipitation), other=np.nan)"
+    "y_pred_prob = y_pred_prob.where(~np.isnan(y_true.precipitation), other=np.nan)\n",
+    "y_refe = y_refe.where(~np.isnan(y_true.precipitation), other=np.nan)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "4e71997f-7808-463a-8b45-dcac639ebe88",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# align digital elevation model\n",
+    "_, dem = xr.align(y_true.precipitation.isel(time=0), dem, join='override')\n",
+    "dem = dem.where(~np.isnan(y_true.precipitation.isel(time=0)), other=np.nan)"
    ]
   },
   {
@@ -257,9 +301,9 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "# get predicted and observed values over entire time series and grid points\n",
-    "y_pred_values = y_pred_pr[NAMES[PREDICTAND]].groupby('time.month').mean(dim='time').values.flatten()\n",
-    "y_true_values = y_true[NAMES[PREDICTAND]].groupby('time.month').mean(dim='time').values.flatten()"
+    "# calculate monthly mean precipitation (mm / month)\n",
+    "y_pred_values = y_pred_pr[NAMES[PREDICTAND]].resample(time='1M').sum(skipna=False).groupby('time.month').mean(dim='time').values\n",
+    "y_true_values = y_true[NAMES[PREDICTAND]].resample(time='1M').sum(skipna=False).groupby('time.month').mean(dim='time').values"
    ]
   },
   {
@@ -283,7 +327,20 @@
    "outputs": [],
    "source": [
     "# calculate coefficient of determination\n",
-    "r2 = r2_score(y_true_values, y_pred_values)"
+    "r2 = r2_score(y_true_values, y_pred_values)\n",
+    "r2"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "a1b33431-9f2b-4e42-bbb0-f4fa258edd98",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# group timeseries by month and calculate mean over time and space\n",
+    "y_pred_ac = y_pred_pr.resample(time='1M').sum(skipna=False).groupby('time.month').mean(dim=('y', 'x', 'time'), skipna=True)\n",
+    "y_true_ac = y_true.resample(time='1M').sum(skipna=False).groupby('time.month').mean(dim=('y', 'x', 'time'), skipna=True)"
    ]
   },
   {
@@ -304,22 +361,33 @@
     "ax.plot(y_true_values, y_pred_values, 'o', alpha=.5, markeredgecolor='grey', markerfacecolor='none', markersize=3);\n",
     "\n",
     "# plot 1:1 mapping line\n",
-    "interval = np.arange(0, 11)\n",
+    "interval = np.arange(0, 300, 50)\n",
     "ax.plot(interval, interval, color='k', lw=2, ls='--')\n",
     "\n",
     "# add coefficient of determination: calculated on entire dataset!\n",
-    "ax.text(interval[-1], interval[0], s='Coefficient of determination R$^2$ = {:.2f}'.format(r2), ha='right', fontsize=14)\n",
+    "ax.text(interval[-1] - 2, interval[0] + 2, s='Coefficient of determination R$^2$ = {:.2f}'.format(r2), ha='right', fontsize=18)\n",
     "\n",
     "# format axes\n",
     "ax.set_ylim(interval[0], interval[-1])\n",
     "ax.set_xlim(interval[0], interval[-1])\n",
     "ax.set_xticks(interval)\n",
-    "ax.set_xticklabels(interval, fontsize=14)\n",
+    "ax.set_xticklabels(interval, fontsize=16)\n",
     "ax.set_yticks(interval)\n",
-    "ax.set_yticklabels(interval, fontsize=14)\n",
-    "ax.set_xlabel('Observed', fontsize=14)\n",
-    "ax.set_ylabel('Predicted', fontsize=14)\n",
-    "ax.set_title('Monthly mean {} (mm / day): 1991 - 2010'.format(NAMES[PREDICTAND]), fontsize=16, pad=10);\n",
+    "ax.set_yticklabels(interval, fontsize=16)\n",
+    "ax.set_xlabel('Observed', fontsize=18)\n",
+    "ax.set_ylabel('Predicted', fontsize=18)\n",
+    "ax.set_title('Monthly mean {} (mm / month)'.format(NAMES[PREDICTAND]), fontsize=20, pad=10);\n",
+    "\n",
+    "# add axis for annual cycle\n",
+    "axins = inset_axes(ax, width=\"30%\", height=\"40%\", loc=2, borderpad=0.25)\n",
+    "axins.plot(y_pred_ac[NAMES[PREDICTAND]].values, ls='--', color='k', label='Predicted')\n",
+    "axins.plot(y_true_ac[NAMES[PREDICTAND]].values, ls='-', color='k', label='Observed')\n",
+    "axins.legend(frameon=False, fontsize=12, loc='lower center');\n",
+    "axins.set_yticks(np.arange(0, 200, 50))\n",
+    "axins.set_yticklabels(np.arange(0, 200, 50), fontsize=12)\n",
+    "axins.yaxis.tick_right()\n",
+    "axins.set_xticks(np.arange(0, 12))\n",
+    "axins.set_xticklabels([calendar.month_name[i + 1] for i in np.arange(0, 12)], rotation=90, fontsize=12)\n",
     "\n",
     "# save figure\n",
     "fig.savefig('../Notebooks/Figures/{}_r2.png'.format(PREDICTAND), dpi=300, bbox_inches='tight')"
@@ -350,10 +418,13 @@
    "source": [
     "# yearly average bias over reference period\n",
     "y_pred_yearly_avg = y_pred_pr.groupby('time.year').mean(dim='time')\n",
+    "y_refe_yearly_avg = y_refe.groupby('time.year').mean(dim='time')\n",
     "y_true_yearly_avg = y_true.groupby('time.year').mean(dim='time')\n",
     "bias_yearly_avg = ((y_pred_yearly_avg - y_true_yearly_avg) / y_true_yearly_avg) * 100\n",
+    "bias_yearly_avg_ref = ((y_refe_yearly_avg - y_true_yearly_avg) / y_true_yearly_avg) * 100\n",
     "for var in bias_yearly_avg:\n",
-    "    print('Yearly average relative bias of {}: {:.2f}%'.format(var, bias_yearly_avg[var].mean().item()))"
+    "    print('(Model) Yearly average relative bias of {}: {:.2f}%'.format(var, bias_yearly_avg[var].mean().item()))\n",
+    "    print('(ERA-5) Yearly average relative bias of {}: {:.2f}%'.format(var, bias_yearly_avg_ref.mean().to_array().values.item()))"
    ]
   },
   {
@@ -364,9 +435,11 @@
    "outputs": [],
    "source": [
     "# mean absolute error over reference period\n",
-    "mae_avg = np.abs(y_pred_yearly_avg - y_true_yearly_avg).mean()\n",
+    "mae_avg = np.abs(y_pred_yearly_avg - y_true_yearly_avg)\n",
+    "mae_avg_ref = np.abs(y_refe_yearly_avg - y_true_yearly_avg)\n",
     "for var in mae_avg:\n",
-    "    print('Yearly average MAE of {}: {:.2f} '.format(var, mae_avg[var].item()) + 'mm / day')"
+    "    print('(Model) Yearly average MAE of {}: {:.2f} mm'.format(var, mae_avg[var].mean().item()))\n",
+    "    print('(ERA-5) Yearly average MAE of {}: {:.2f} mm'.format(var, mae_avg_ref.mean().to_array().values.item()))"
    ]
   },
   {
@@ -378,8 +451,10 @@
    "source": [
     "# root mean squared error over reference period\n",
     "rmse_avg = ((y_pred_yearly_avg - y_true_yearly_avg) ** 2).mean()\n",
+    "rmse_avg_ref = ((y_refe_yearly_avg - y_true_yearly_avg) **2).mean()\n",
     "for var in rmse_avg:\n",
-    "    print('Yearly average RMSE of {}: {:.2f} '.format(var, rmse_avg[var].item()) + 'mm / day')"
+    "    print('(Model) Yearly average RMSE of {}: {:.2f} mm / day'.format(var, rmse_avg[var].item()))\n",
+    "    print('(ERA-5) Yearly average RMSE of {}: {:.2f} mm / day'.format(var, rmse_avg_ref.mean().to_array().values.item()))"
    ]
   },
   {
@@ -400,6 +475,74 @@
     "print('Yearly average Pearson correlation coefficient for {}: {:.2f}'.format(var, np.asarray(r).mean()))"
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "ca683581-6b3b-4abd-ad46-157ebded19c6",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# plot yearly average MAE of reference vs. prediction\n",
+    "vmin, vmax = 0, 5\n",
+    "fig, axes = plt.subplots(1, 3, figsize=(24, 8), sharex=True, sharey=True)\n",
+    "\n",
+    "# plot bias of ERA-5 reference\n",
+    "reference = bias_yearly_avg_ref.mean(dim='year').to_array().squeeze()\n",
+    "im1 = axes[0].imshow(reference.values, origin='lower', cmap='RdBu_r', vmin=-40, vmax=40)\n",
+    "axes[0].text(x=reference.shape[0] - 2, y=2, s='Average: {:.1f}%'.format(reference.mean().item()), fontsize=14, ha='right')\n",
+    "\n",
+    "# plot MAE of model\n",
+    "prediction = bias_yearly_avg.mean(dim='year').to_array().squeeze()\n",
+    "im2 = axes[1].imshow(prediction.values, origin='lower', cmap='RdBu_r', vmin=-40, vmax=40)\n",
+    "axes[1].text(x=reference.shape[0] - 2, y=2, s='Average: {:.1f}%'.format(prediction.mean().item()), fontsize=14, ha='right')\n",
+    "\n",
+    "# plot topography\n",
+    "im_dem = axes[2].imshow(dem['elevation'].values, origin='lower', cmap='terrain', vmin=0, vmax=4000)\n",
+    "\n",
+    "# set titles\n",
+    "axes[0].set_title('ERA-5', fontsize=14, pad=10);\n",
+    "axes[1].set_title('DCEDN', fontsize=14, pad=10);\n",
+    "axes[2].set_title('Copernicus EU-DEM v1.1', fontsize=14, pad=10)\n",
+    "\n",
+    "# adjust axes\n",
+    "for ax in axes.flat:\n",
+    "    ax.axes.get_xaxis().set_ticklabels([])\n",
+    "    ax.axes.get_xaxis().set_ticks([])\n",
+    "    ax.axes.get_yaxis().set_ticklabels([])\n",
+    "    ax.axes.get_yaxis().set_ticks([])\n",
+    "    ax.axes.axis('tight')\n",
+    "    ax.set_xlabel('')\n",
+    "    ax.set_ylabel('')\n",
+    "    ax.set_axis_off()\n",
+    "\n",
+    "# adjust figure\n",
+    "# fig.suptitle('Average yearly mean absolute error: 1991 - 2010', fontsize=20);\n",
+    "fig.subplots_adjust(hspace=0, wspace=0, top=0.85)\n",
+    "\n",
+    "# add colorbar for bias\n",
+    "axes = axes.flatten()\n",
+    "cbar_ax_bias = fig.add_axes([axes[-1].get_position().x1 + 0.01, axes[-1].get_position().y0,\n",
+    "                             0.01, axes[-1].get_position().y1 - axes[-1].get_position().y0])\n",
+    "cbar_bias = fig.colorbar(im_dem, cax=cbar_ax_bias)\n",
+    "cbar_bias.set_label(label='Elevation (m)', fontsize=14)\n",
+    "cbar_bias.ax.tick_params(labelsize=14, pad=10)\n",
+    "\n",
+    "# add colorbar for predictand\n",
+    "cbar_ax_predictand = fig.add_axes([axes[0].get_position().x0, axes[0].get_position().y0 - 0.1,\n",
+    "                                   axes[-1].get_position().x0 - axes[0].get_position().x0,\n",
+    "                                   0.03])\n",
+    "cbar_predictand = fig.colorbar(im1, cax=cbar_ax_predictand, orientation='horizontal')\n",
+    "cbar_predictand.set_label(label='Relative mean error (%)', fontsize=14)\n",
+    "cbar_predictand.ax.tick_params(labelsize=14, pad=10)\n",
+    "\n",
+    "# add metrics: MAE and RMSE\n",
+    "#axes[1].text(x=ds.shape[0] - 2, y=2, s='MAE = {:.1f}'.format(mae_avg[NAMES[PREDICTAND]].item()) + 'mm day$^{-1}$', fontsize=14, ha='right')\n",
+    "#axes[1].text(x=ds.shape[0] - 2, y=12, s='RMSE = {:.1f}'.format(rmse_avg[NAMES[PREDICTAND]].item()) + 'mm$^2$ day$^{-2}$', fontsize=14, ha='right')\n",
+    "\n",
+    "# save figure\n",
+    "fig.savefig('../Notebooks/Figures/{}_rbias_ERA_vs_model.png'.format(PREDICTAND), dpi=300, bbox_inches='tight')"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,
@@ -409,7 +552,7 @@
    "source": [
     "# plot average of observation, prediction, and bias\n",
     "vmin, vmax = 0, 5\n",
-    "fig, axes = plt.subplots(len(y_pred_yearly_avg.data_vars), 3, figsize=(24, len(y_pred_yearly_avg.data_vars) * 6),\n",
+    "fig, axes = plt.subplots(len(y_pred_yearly_avg.data_vars), 3, figsize=(24, len(y_pred_yearly_avg.data_vars) * 8),\n",
     "                         sharex=True, sharey=True)\n",
     "axes = axes.reshape(len(y_pred_yearly_avg.data_vars), -1)\n",
     "for i, var in enumerate(y_pred_yearly_avg):\n",
@@ -457,8 +600,8 @@
     "cbar_predictand.ax.tick_params(labelsize=14)\n",
     "\n",
     "# add metrics: MAE and RMSE\n",
-    "axes[1].text(x=ds.shape[0] - 2, y=2, s='MAE = {:.1f}'.format(mae_avg[NAMES[PREDICTAND]].item()) + 'mm day$^{-1}$', fontsize=14, ha='right')\n",
-    "axes[1].text(x=ds.shape[0] - 2, y=12, s='RMSE = {:.1f}'.format(rmse_avg[NAMES[PREDICTAND]].item()) + 'mm$^2$ day$^{-2}$', fontsize=14, ha='right')\n",
+    "axes[1].text(x=ds.shape[0] - 2, y=2, s='MAE = {:.1f}'.format(mae_avg[NAMES[PREDICTAND]].mean().item()) + 'mm day$^{-1}$', fontsize=14, ha='right')\n",
+    "axes[1].text(x=ds.shape[0] - 2, y=12, s='RMSE = {:.1f}'.format(rmse_avg[NAMES[PREDICTAND]].mean().item()) + 'mm day$^{-1}$', fontsize=14, ha='right')\n",
     "\n",
     "# save figure\n",
     "fig.savefig('../Notebooks/Figures/{}_average_bias.png'.format(PREDICTAND), dpi=300, bbox_inches='tight')"
@@ -490,7 +633,22 @@
     "# group data by season: (DJF, MAM, JJA, SON)\n",
     "y_true_snl = y_true.groupby('time.season').mean(dim='time')\n",
     "y_pred_snl = y_pred_pr.groupby('time.season').mean(dim='time')\n",
-    "bias_snl = ((y_pred_snl - y_true_snl) / y_true_snl) * 100"
+    "y_refe_snl = y_refe.groupby('time.season').mean(dim='time')\n",
+    "bias_snl = ((y_pred_snl - y_true_snl) / y_true_snl) * 100\n",
+    "bias_snl_ref = ((y_refe_snl - y_true_snl) / y_true_snl) * 100"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "2a70e6ac-54fb-47b8-ae76-a1b51d9d1f17",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# print average bias per season: ERA-5\n",
+    "for var in bias_snl_ref.data_vars:\n",
+    "    for season in bias_snl_ref[NAMES[PREDICTAND]].season:\n",
+    "        print('(ERA-5) Average bias of mean {} for season {}: {:.1f}%'.format(var, season.values.item(), bias_snl_ref[var].sel(season=season).mean().item()))"
    ]
   },
   {
@@ -500,10 +658,10 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "# print average bias per season\n",
+    "# print average bias per season: model\n",
     "for var in bias_snl.data_vars:\n",
     "    for season in bias_snl[NAMES[PREDICTAND]].season:\n",
-    "        print('Average bias of mean {} for season {}: {:.1f}%'.format(var, season.values.item(), bias_snl[var].sel(season=season).mean().item()))"
+    "        print('(Model) Average bias of mean {} for season {}: {:.1f}%'.format(var, season.values.item(), bias_snl[var].sel(season=season).mean().item()))"
    ]
   },
   {
@@ -597,7 +755,8 @@
     "with warnings.catch_warnings():\n",
     "    warnings.simplefilter('ignore', category=RuntimeWarning)\n",
     "    y_pred_ex = y_pred_pr.groupby('time.year').quantile(quantile, dim='time')\n",
-    "    y_true_ex = y_true.groupby('time.year').quantile(quantile, dim='time')"
+    "    y_true_ex = y_true.groupby('time.year').quantile(quantile, dim='time')\n",
+    "    y_refe_ex = y_refe.groupby('time.year').quantile(quantile, dim='time')"
    ]
   },
   {
@@ -609,8 +768,31 @@
    "source": [
     "# calculate bias in extreme quantile for each year\n",
     "bias_ex = ((y_pred_ex - y_true_ex) / y_true_ex) * 100\n",
+    "bias_ex_ref = ((y_refe_ex - y_true_ex) / y_true_ex) * 100"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "9a610a70-bb94-47e6-b470-e01bf4a295c0",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# bias of extreme quantile: ERA-5\n",
+    "for var in bias_ex_ref:\n",
+    "    print('(ERA-5) Yearly average bias for P{:.0f} of {}: {:.1f}%'.format(quantile * 100, var, bias_ex_ref[var].mean().item()))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "423962d7-453d-4c86-b93e-a4571bfec1c4",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# bias of extreme quantile: Model\n",
     "for var in bias_ex:\n",
-    "    print('Yearly average bias for P{:.0f} of {}: {:.1f}%'.format(quantile * 100, var, bias_ex[var].mean().item()))"
+    "    print('(Model) Yearly average bias for P{:.0f} of {}: {:.1f}%'.format(quantile * 100, var, bias_ex[var].mean().item()))"
    ]
   },
   {
@@ -622,8 +804,31 @@
    "source": [
     "# mean absolute error in extreme quantile\n",
     "mae_ex = np.abs(y_pred_ex - y_true_ex).mean()\n",
+    "mae_ex_ref = np.abs(y_refe_ex - y_true_ex).mean()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "7331acdd-60bb-414e-a0e7-e81956c3a1bf",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# mae of extreme quantile: ERA-5\n",
+    "for var in mae_ex_ref:\n",
+    "    print('(ERA-5) Yearly average MAE for P{:.0f} of {}: {:.1f} mm / day'.format(quantile * 100, var, mae_ex_ref[var].item()))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "d7bf1a6c-e111-4c72-8774-64a3939cbe50",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# mae of extreme quantile: Model\n",
     "for var in mae_ex:\n",
-    "    print('Yearly average MAE for P{:.0f} of {}: {:.1f} mm / day'.format(quantile * 100, var, mae_ex[var].item()))"
+    "    print('(Model) Yearly average MAE for P{:.0f} of {}: {:.1f} mm / day'.format(quantile * 100, var, mae_ex[var].item()))"
    ]
   },
   {
@@ -633,10 +838,33 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "# root mean squared error over reference period\n",
+    "# root mean squared error in extreme quantile\n",
     "rmse_ex = ((y_pred_ex - y_true_ex) ** 2).mean()\n",
+    "rmse_ex_ref = ((y_refe_ex - y_true_ex) ** 2).mean()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "712a2109-14b9-45a3-9194-e2b917c5ba3f",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# rmse of extreme quantile: ERA-5\n",
+    "for var in rmse_ex_ref:\n",
+    "    print('(ERA-5) Yearly average RMSE for P{:.0f} of {}: {:.1f} mm / day'.format(quantile * 100, var, rmse_ex_ref[var].item()))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "1dfad687-c16f-4767-8ddd-96c3e56bd96a",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# rmse of extreme quantile: Model\n",
     "for var in rmse_ex:\n",
-    "    print('Yearly average RMSE for P{:.0f} of {}: {:.1f} mm / day'.format(quantile * 100, var, rmse_ex[var].item()))"
+    "    print('(Model) Yearly average RMSE for P{:.0f} of {}: {:.1f} mm / day'.format(quantile * 100, var, rmse_ex[var].item()))"
    ]
   },
   {
@@ -722,7 +950,8 @@
     "with warnings.catch_warnings():\n",
     "    warnings.simplefilter('ignore', category=RuntimeWarning)\n",
     "    y_true_ex_snl = y_true.groupby('time.season').quantile(quantile, dim='time')\n",
-    "    y_pred_ex_snl = y_pred_pr.groupby('time.season').quantile(quantile, dim='time')"
+    "    y_pred_ex_snl = y_pred_pr.groupby('time.season').quantile(quantile, dim='time')\n",
+    "    y_refe_ex_snl = y_refe.groupby('time.season').quantile(quantile, dim='time')"
    ]
   },
   {
@@ -733,7 +962,8 @@
    "outputs": [],
    "source": [
     "# compute relative bias in seasonal extremes\n",
-    "bias_ex_snl = ((y_pred_ex_snl - y_true_ex_snl) / y_true_ex_snl) * 100"
+    "bias_ex_snl = ((y_pred_ex_snl - y_true_ex_snl) / y_true_ex_snl) * 100\n",
+    "bias_ex_snl_ref = ((y_refe_ex_snl - y_true_ex_snl) / y_true_ex_snl) * 100"
    ]
   },
   {
@@ -743,10 +973,23 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "# print average bias in extreme per season\n",
+    "# print average bias in extreme per season: ERA-5\n",
+    "for var in bias_ex_snl_ref.data_vars:\n",
+    "    for season in bias_ex_snl_ref[NAMES[PREDICTAND]].season:\n",
+    "        print('(ERA-5) Average bias of P{:.0f} {} for season {}: {:.1f}%'.format(quantile * 100, var, season.values.item(), bias_ex_snl_ref[var].sel(season=season).mean().item()))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "0671b4f3-52c6-4167-b35b-03902dbe11a3",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# print average bias in extreme per season: Model\n",
     "for var in bias_ex_snl.data_vars:\n",
     "    for season in bias_ex_snl[NAMES[PREDICTAND]].season:\n",
-    "        print('Average bias of P{:.0f} {} for season {}: {:.1f}%'.format(quantile * 100, var, season.values.item(), bias_ex_snl[var].sel(season=season).mean().item()))"
+    "        print('(Model) Average bias of P{:.0f} {} for season {}: {:.1f}%'.format(quantile * 100, var, season.values.item(), bias_ex_snl[var].sel(season=season).mean().item()))"
    ]
   },
   {