diff --git a/notebooks/tutorials/massbalance_calibration.ipynb b/notebooks/tutorials/massbalance_calibration.ipynb
index aa0501dc..d9e61587 100644
--- a/notebooks/tutorials/massbalance_calibration.ipynb
+++ b/notebooks/tutorials/massbalance_calibration.ipynb
@@ -4,7 +4,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "# A look into the new mass balance calibration in OGGM v1.6"
+    "# A look into the mass balance calibration in OGGM v1.6"
    ]
   },
   {
@@ -43,7 +43,8 @@
     "\n",
     "from oggm import cfg, utils, workflow, tasks, graphics\n",
     "from oggm.core import massbalance\n",
-    "from oggm.core.massbalance import mb_calibration_from_scalar_mb, mb_calibration_from_geodetic_mb, mb_calibration_from_wgms_mb"
+    "from oggm.core.massbalance import mb_calibration_from_scalar_mb, mb_calibration_from_geodetic_mb, mb_calibration_from_wgms_mb, mb_calibration_to_rmsd\n",
+    "from oggm.utils import rmsd"
    ]
   },
   {
@@ -330,7 +331,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "# if you use the default calibration data from Hugonnet et al., 2021, \n",
+    "# if you use the default calibration data from Hugonnet et al., 2021,\n",
     "# we can get the geodetic data from any glacier from here:\n",
     "ref_mb_df = utils.get_geodetic_mb_dataframe().loc[gdir_hef.rgi_id]\n",
     "ref_mb_df = ref_mb_df.loc[ref_mb_df['period'] == cfg.PARAMS['geodetic_mb_period']].iloc[0]\n",
@@ -349,7 +350,7 @@
     "# overwrite_gdir has to be set to True,\n",
     "# because we want to overwrite the old calibration\n",
     "mb_calibration_from_scalar_mb(gdir_hef,\n",
-    "                              ref_mb = ref_mb, \n",
+    "                              ref_mb = ref_mb,\n",
     "                              ref_period=cfg.PARAMS['geodetic_mb_period'],\n",
     "                              overwrite_gdir=True);\n",
     "\n",
@@ -404,7 +405,7 @@
     "# overwrite_gdir has to be set to True,\n",
     "# because we want to overwrite the old calibration\n",
     "mb_calibration_from_scalar_mb(gdir_hef,\n",
-    "                              ref_mb = ref_mb, \n",
+    "                              ref_mb = ref_mb,\n",
     "                              ref_period=cfg.PARAMS['geodetic_mb_period'],\n",
     "                              calibrate_param1='temp_bias',\n",
     "                              overwrite_gdir=True)\n",
@@ -454,7 +455,7 @@
     "# overwrite_gdir has to be set to True,\n",
     "# because we want to overwrite the old calibration\n",
     "mb_calibration_from_scalar_mb(gdir_hef,\n",
-    "                              ref_mb = ref_mb, \n",
+    "                              ref_mb = ref_mb,\n",
     "                              ref_period=cfg.PARAMS['geodetic_mb_period'],\n",
     "                              calibrate_param1='prcp_fac',\n",
     "                              overwrite_gdir=True)\n",
@@ -565,7 +566,7 @@
    "source": [
     "mbdf_in_situ = gdir_hef.get_ref_mb_data()\n",
     "mbdf_in_situ['ref_mb'] = mbdf_in_situ['ANNUAL_BALANCE']\n",
-    "# check that every year between the beginning and the end has MB observations \n",
+    "# check that every year between the beginning and the end has MB observations\n",
     "assert len(mbdf_in_situ.index) == (mbdf_in_situ.index[-1] - mbdf_in_situ.index[0] + 1)\n",
     "ref_mb = mbdf_in_situ.ref_mb.mean()"
    ]
@@ -589,6 +590,15 @@
     "plt.xlabel('Year');"
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "mbdf_in_situ[['ref_mb','mod_mb']].mean()"
+   ]
+  },
   {
    "cell_type": "markdown",
    "metadata": {},
@@ -621,6 +631,231 @@
     "mbdf_in_situ[['ref_mb','mod_mb']].std()"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Calibrating on direct glaciological in-situ MB observations by minimising the RMSD\n",
+    "\n",
+    "In the previous section, we calibrated the model using in-situ observations via `mb_calibration_from_wgms_mb`, which internally relies on `mb_calibration_from_scalar_mb`. This approach minimises the *mean difference* between model and observations. While effective, it reduces the dataset to its average value and therefore ignores valuable information about year-to-year variability.\n",
+    "\n",
+    "However, WGMS mass balance observations provide *annual* values over the observation period. This richer dataset allows us to use calibration methods that account for interannual variability.\n",
+    "\n",
+    "OGGM now provides such a method: `mb_calibration_to_rmsd`. Instead of focusing only on the mean, it minimises the Root Mean Square Deviation (RMSD) between observed and modelled annual mass balance values.\n",
+    "\n",
+    "Another key difference is that this method performs a **multi-parameter optimisation**, whereas previous approaches typically optimise one parameter at a time.\n",
+    "\n",
+    "In the following, we compare this RMSD-based calibration to the previous approach, focusing on differences in mean, variability (standard deviation), and overall fit."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# The mass balance time series generated by calibrating melt_f\n",
+    "# to minimise the mean mass balance using mb_calibration_from_wgms_mb\n",
+    "print(f\"RMSD: {rmsd(mbdf_in_situ['ref_mb'], mbdf_in_situ['mod_mb']):.2f} kg m-2\")\n",
+    "\n",
+    "sd_diff = abs(mbdf_in_situ['ref_mb'].std() - mbdf_in_situ['mod_mb'].std())\n",
+    "print(f\"Std dev difference: {sd_diff:.2f} kg m-2\")\n",
+    "\n",
+    "mean_diff = abs(mbdf_in_situ['ref_mb'].mean() - mbdf_in_situ['mod_mb'].mean())\n",
+    "print(f\"Mean difference: {mean_diff:.2f} kg m-2\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The difference in mean mass balance is negligible. However, the differences in standard deviation and RMSD are much larger, indicating that important aspects of the variability are not well captured.\n",
+    "\n",
+    "Let us now consider an alternative approach, using the RMSD as the cost function.\n",
+    "\n",
+    "We begin by calibrating only the melt factor, while keeping the temperature bias and precipitation factor fixed at their default values. If no calibration parameters are specified, the function uses this configuration by default."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "mb_calibration_to_rmsd(gdir_hef, ref_df=mbdf_in_situ['ref_mb'], overwrite_gdir=True)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "mbmod_rmsd = massbalance.MonthlyTIModel(gdir_hef)\n",
+    "mbdf_in_situ['mod_mb_rmsd'] = mbmod_rmsd.get_specific_mb(\n",
+    "    h, w, year=mbdf_in_situ.index\n",
+    ")\n",
+    "\n",
+    "ref = mbdf_in_situ['ref_mb']\n",
+    "mod_mean = mbdf_in_situ['mod_mb']\n",
+    "mod_rmsd = mbdf_in_situ['mod_mb_rmsd']\n",
+    "\n",
+    "plt.plot(\n",
+    "    ref,\n",
+    "    color='grey',\n",
+    "    lw=3,\n",
+    "    label=(\n",
+    "        \"in-situ observations\\n\"\n",
+    "        f\"average: {ref.mean():.1f} \" + r\"kg m$^{-2}$\"\n",
+    "    ),\n",
+    ")\n",
+    "\n",
+    "plt.plot(\n",
+    "    mod_mean,\n",
+    "    label=(\n",
+    "        \"modelled mass balance\\n\"\n",
+    "        \"calibrated on mean\\n\"\n",
+    "        f\"average: {mod_mean.mean():.1f} \" + r\"kg m$^{-2}$\"\n",
+    "    ),\n",
+    ")\n",
+    "\n",
+    "plt.plot(\n",
+    "    mod_rmsd,\n",
+    "    label=(\n",
+    "        \"modelled mass balance\\n\"\n",
+    "        \"calibrated on RMSD\\n\"\n",
+    "        f\"average: {mod_rmsd.mean():.1f} \" + r\"kg m$^{-2}$\"\n",
+    "    ),\n",
+    ")\n",
+    "\n",
+    "plt.legend()\n",
+    "plt.ylabel(r\"specific mass balance (kg m$^{-2}$)\")\n",
+    "plt.xlabel(\"Year\")\n",
+    "plt.title(\"Calibrating melt_f using RMSD vs mean-based calibration\")\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "print(f\"RMSD: {rmsd(mbdf_in_situ['ref_mb'], mbdf_in_situ['mod_mb_rmsd']):.2f} kg m-2\")\n",
+    "\n",
+    "sd_diff = abs(mbdf_in_situ['ref_mb'].std() - mbdf_in_situ['mod_mb_rmsd'].std())\n",
+    "print(f\"Std dev difference: {sd_diff:.2f} kg m-2\")\n",
+    "\n",
+    "mean_diff = abs(mbdf_in_situ['ref_mb'].mean() - mbdf_in_situ['mod_mb_rmsd'].mean())\n",
+    "print(f\"Mean difference: {mean_diff:.2f} kg m-2\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "We can see a slight reduction in the RMSD and the differences in the standard deviation, however, now our differences in the means is far larger. Overall, all of these errors are still quite large. Let's see if we can reduce these errors further..."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Now let's try by calibrating all of the mass-balance parameters."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "mb_calibration_to_rmsd(gdir_hef, ref_df=mbdf_in_situ['ref_mb'], overwrite_gdir=True,\n",
+    "                       calibrate_params=('melt_f', 'prcp_fac', 'temp_bias'))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "mbmod_rmsd3 = massbalance.MonthlyTIModel(gdir_hef)\n",
+    "mbdf_in_situ['mod_mb_rmsd3'] = mbmod_rmsd3.get_specific_mb(h, w, year=mbdf_in_situ.index)\n",
+    "\n",
+    "ref = mbdf_in_situ['ref_mb']\n",
+    "mod_mean = mbdf_in_situ['mod_mb']\n",
+    "mod_rmsd = mbdf_in_situ['mod_mb_rmsd']\n",
+    "mod_rmsd3 = mbdf_in_situ['mod_mb_rmsd3']\n",
+    "\n",
+    "plt.plot(\n",
+    "    ref,\n",
+    "    color='grey',\n",
+    "    lw=3,\n",
+    "    label=(\n",
+    "        \"in-situ observations\\n\"\n",
+    "        f\"average: {ref.mean():.1f} \" + r\"kg m$^{-2}$\"\n",
+    "    ),\n",
+    ")\n",
+    "\n",
+    "plt.plot(\n",
+    "    mod_mean,\n",
+    "    label=(\n",
+    "        \"modelled mass balance\\n\"\n",
+    "        \"calibrated on mean\\n\"\n",
+    "        f\"average: {mod_mean.mean():.1f} \" + r\"kg m$^{-2}$\"\n",
+    "    ),\n",
+    ")\n",
+    "\n",
+    "plt.plot(\n",
+    "    mod_rmsd,\n",
+    "    label=(\n",
+    "        \"modelled mass balance\\n\"\n",
+    "        \"calibrated on RMSD\\n\"\n",
+    "        f\"average: {mod_rmsd.mean():.1f} \" + r\"kg m$^{-2}$\"\n",
+    "    ),\n",
+    ")\n",
+    "\n",
+    "plt.plot(\n",
+    "    mod_rmsd3,\n",
+    "    label=(\n",
+    "        \"modelled mass balance\\n\"\n",
+    "        \"calibrated on RMSD (3 params)\\n\"\n",
+    "        f\"average: {mod_rmsd3.mean():.1f} \" + r\"kg m$^{-2}$\"\n",
+    "    ),\n",
+    ")\n",
+    "\n",
+    "plt.legend()\n",
+    "plt.ylabel(r\"specific mass balance (kg m$^{-2}$)\")\n",
+    "plt.xlabel(\"Year\")\n",
+    "plt.title(\"Calibrating melt_f using RMSD vs mean-based calibration\")\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "print(f\"RMSD: {rmsd(mbdf_in_situ['ref_mb'], mbdf_in_situ['mod_mb_rmsd3']):.2f} kg m-2\")\n",
+    "\n",
+    "sd_diff = abs(mbdf_in_situ['ref_mb'].std() - mbdf_in_situ['mod_mb_rmsd3'].std())\n",
+    "print(f\"Std dev difference: {sd_diff:.2f} kg m-2\")\n",
+    "\n",
+    "mean_diff = abs(mbdf_in_situ['ref_mb'].mean() - mbdf_in_situ['mod_mb_rmsd3'].mean())\n",
+    "print(f\"Mean difference: {mean_diff:.2f} kg m-2\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "When calibrating all three mass-balance parameters simultaneously, we can identify more optimal parameter combinations that reduce the RMSD, while also improving the agreement in mean and standard deviation.\n",
+    "\n",
+    "This approach provides an alternative calibration strategy that captures multiple characteristics of the dataset, rather than focusing solely on the mean mass balance."
+   ]
+  },
   {
    "cell_type": "markdown",
    "metadata": {
@@ -726,10 +961,10 @@
    "outputs": [],
    "source": [
     "# Allowing another parameter to change is done by defining calibrate_param2\n",
-    "mb_calibration_from_scalar_mb(gdir_gol, ref_mb=ref_mb, \n",
+    "mb_calibration_from_scalar_mb(gdir_gol, ref_mb=ref_mb,\n",
     "                              write_to_gdir=False,\n",
-    "                              ref_period =cfg.PARAMS['geodetic_mb_period'], \n",
-    "                              calibrate_param1='melt_f', \n",
+    "                              ref_period =cfg.PARAMS['geodetic_mb_period'],\n",
+    "                              calibrate_param1='melt_f',\n",
     "                              calibrate_param2='temp_bias')"
    ]
   },
@@ -763,8 +998,8 @@
     "ref_period = '2000-01-01_2010-01-01'\n",
     "ref_mb = 2000 # Let's use an unrealistically positive mass-balance\n",
     "mb_calibration_from_scalar_mb(gdir_hef, ref_mb=ref_mb,\n",
-    "                              ref_period=ref_period, \n",
-    "                              overwrite_gdir=True, \n",
+    "                              ref_period=ref_period,\n",
+    "                              overwrite_gdir=True,\n",
     "                              calibrate_param2='prcp_fac');"
    ]
   },
@@ -899,20 +1134,20 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "# calibrate the melt_f and annual MB \n",
+    "# calibrate the melt_f and annual MB\n",
     "# Create a dataframe with precip factor going from 0.1 to 5.0 in 0.3 steps\n",
     "pd_prcp_fac_sens = pd.DataFrame(index=np.arange(0.1,5.0,0.3))\n",
     "# Calibrate the melt factor for each of these\n",
     "spec_mb_prcp_fac_sens_dict = {}\n",
     "for prcp_fac in pd_prcp_fac_sens.index:\n",
-    "    calib_param = mb_calibration_from_scalar_mb(gdir_hef, \n",
-    "                                                ref_mb=ref_mb, \n",
+    "    calib_param = mb_calibration_from_scalar_mb(gdir_hef,\n",
+    "                                                ref_mb=ref_mb,\n",
     "                                                ref_period=cfg.PARAMS['geodetic_mb_period'],\n",
-    "                                                prcp_fac=prcp_fac, \n",
+    "                                                prcp_fac=prcp_fac,\n",
     "                                                overwrite_gdir=True)\n",
     "    # Fill the dataframe with the calibrated parameters\n",
     "    pd_prcp_fac_sens.loc[prcp_fac, 'melt_f'] = calib_param['melt_f']\n",
-    "     \n",
+    "\n",
     "    # Check the modelled massbalance\n",
     "    mb_prcp_fac_sens = massbalance.MonthlyTIModel(gdir_hef)\n",
     "    # ok, we actually matched the new ref_mb\n",
@@ -958,11 +1193,11 @@
     "plt.figure()\n",
     "for j,prcp_fac in enumerate(pd_prcp_fac_sens.index):\n",
     "    plt.plot(np.arange(2000,2020,1),\n",
-    "             spec_mb_prcp_fac_sens_dict[prcp_fac], '-', \n",
+    "             spec_mb_prcp_fac_sens_dict[prcp_fac], '-',\n",
     "             color=colors_prcp_fac[j], label=f'{prcp_fac:.1f}')\n",
     "plt.plot(mbdf_in_situ.loc[2000:2019].index,\n",
     "         mbdf_in_situ.loc[2000:2019]['ANNUAL_BALANCE'],\n",
-    "         color='grey', lw=3, \n",
+    "         color='grey', lw=3,\n",
     "        label='observed in-situ')\n",
     "plt.ylabel(r'Annual mass-balance (kg m$^{-2}$)')\n",
     "plt.xlabel('Year')\n",
@@ -1012,9 +1247,9 @@
     "spec_mb_temp_b_sens_dict = {}\n",
     "\n",
     "for temp_bias in np.arange(-5,5.0,0.5):\n",
-    "    # for too negative temp_bias, no melt_f is found that matches the observations. We would need to \n",
+    "    # for too negative temp_bias, no melt_f is found that matches the observations. We would need to\n",
     "    # change the prcp_fac , but here we will just look\n",
-    "    # at those combinations where calibration works with a fixed prcp_fac. \n",
+    "    # at those combinations where calibration works with a fixed prcp_fac.\n",
     "    try:\n",
     "        calib_param = mb_calibration_from_scalar_mb(gdir_hef, ref_mb=ref_mb, ref_period=cfg.PARAMS['geodetic_mb_period'],\n",
     "                                                    temp_bias=temp_bias, overwrite_gdir=True)\n",
@@ -1039,17 +1274,17 @@
     "fig, axs = plt.subplots(1,2,figsize=(14,6))\n",
     "for j,temp_bias in enumerate(melt_f_dict_tb.keys()):\n",
     "    axs[0].plot(temp_bias, melt_f_dict_tb[temp_bias], 'o',\n",
-    "             color=colors_temp_bias(norm(temp_bias))) \n",
+    "             color=colors_temp_bias(norm(temp_bias)))\n",
     "axs[0].set_ylabel(r'melt_f (kg m$^{-2}$ day$^{-1}$ K$^{-1}$)')\n",
     "axs[0].set_xlabel('temp_bias (°C)');\n",
     "\n",
     "\n",
     "for temp_bias in melt_f_dict_tb.keys():\n",
     "    axs[1].plot(np.arange(2000,2020,1),\n",
-    "             spec_mb_temp_b_sens_dict[temp_bias], '-', \n",
+    "             spec_mb_temp_b_sens_dict[temp_bias], '-',\n",
     "             color=colors_temp_bias(norm(temp_bias)),\n",
     "             label=temp_bias)\n",
-    "axs[1].plot(mbdf_in_situ.loc[2000:2019].index, mbdf_in_situ.loc[2000:2019]['ANNUAL_BALANCE'], color='grey', lw=3, \n",
+    "axs[1].plot(mbdf_in_situ.loc[2000:2019].index, mbdf_in_situ.loc[2000:2019]['ANNUAL_BALANCE'], color='grey', lw=3,\n",
     "        label='observed in-situ')\n",
     "axs[1].set_ylabel(r'Annual mass-balance (kg m$^{-2}$)')\n",
     "axs[1].set_xlabel('Year')\n",
@@ -1083,9 +1318,9 @@
     "spec_mb_pf_temp_b_sens_dict = {}\n",
     "\n",
     "for temp_bias in np.arange(-5,5.0,0.5):\n",
-    "    # for too negative temp_bias, no prcp_fac is found that matches the observations. We would need to \n",
+    "    # for too negative temp_bias, no prcp_fac is found that matches the observations. We would need to\n",
     "    # change the melt_f , but here we will just look\n",
-    "    # at those combinations where calibration works with a fixed melt_f. \n",
+    "    # at those combinations where calibration works with a fixed melt_f.\n",
     "    try:\n",
     "        calib_param = mb_calibration_from_scalar_mb(gdir_hef, calibrate_param1='prcp_fac',\n",
     "                                                    ref_mb=ref_mb, ref_period=cfg.PARAMS['geodetic_mb_period'],\n",
@@ -1111,17 +1346,17 @@
     "fig, axs = plt.subplots(1,2,figsize=(14,6))\n",
     "for j,temp_bias in enumerate(pf_temp_b_dict_tb.keys()):\n",
     "    axs[0].plot(temp_bias, pf_temp_b_dict_tb[temp_bias], 'o',\n",
-    "             color=colors_temp_bias(norm(temp_bias))) \n",
+    "             color=colors_temp_bias(norm(temp_bias)))\n",
     "axs[0].set_ylabel(r'prcp_fac')\n",
     "axs[0].set_xlabel('temp_bias (°C)');\n",
     "\n",
     "\n",
     "for temp_bias in pf_temp_b_dict_tb.keys():\n",
     "    axs[1].plot(np.arange(2000,2020,1),\n",
-    "             spec_mb_pf_temp_b_sens_dict[temp_bias], '-', \n",
+    "             spec_mb_pf_temp_b_sens_dict[temp_bias], '-',\n",
     "             color=colors_temp_bias(norm(temp_bias)),\n",
     "             label=temp_bias)\n",
-    "axs[1].plot(mbdf_in_situ.loc[2000:2019].index, mbdf_in_situ.loc[2000:2019]['ANNUAL_BALANCE'], color='grey', lw=3, \n",
+    "axs[1].plot(mbdf_in_situ.loc[2000:2019].index, mbdf_in_situ.loc[2000:2019]['ANNUAL_BALANCE'], color='grey', lw=3,\n",
     "        label='observed in-situ')\n",
     "axs[1].set_ylabel(r'Annual mass-balance (kg m$^{-2}$)')\n",
     "axs[1].set_xlabel('Year')\n",
@@ -1176,8 +1411,8 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "workflow.execute_entity_task(mb_calibration_from_geodetic_mb, \n",
-    "                             gdir_hef, \n",
+    "workflow.execute_entity_task(mb_calibration_from_geodetic_mb,\n",
+    "                             gdir_hef,\n",
     "                             informed_threestep=True,   # \"informed_threestep\" is what OGGM used for 1.6\n",
     "                             overwrite_gdir=True)"
    ]
@@ -1373,7 +1608,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.11.3"
+   "version": "3.14.0"
   },
   "latex_envs": {
    "LaTeX_envs_menu_present": true,