diff --git a/step4_plot_regional_sealevel.py b/step4_plot_regional_sealevel.py
index 55b34e9..2507a41 100644
--- a/step4_plot_regional_sealevel.py
+++ b/step4_plot_regional_sealevel.py
@@ -12,7 +12,7 @@
 
 from config import settings
 from tide_gauge_locations import extract_site_info, find_nearest_station_id
-from slr_pkg import abbreviate_location_name  # found in __init.py__
+from slr_pkg import abbreviate_location_name, choose_montecarlo_dir  # found in __init.py__
 from surge import tide_gauge_library as tgl
 from directories import read_dir, makefolder
 from plotting_libraries import location_string, scenario_string, \
@@ -29,12 +29,13 @@ def compute_uncertainties(df_r_list, scenarios, tg_years, tg_amsl):
     :return: years, scenario uncertainty, model uncertainty, internal
     variability
     """
-    allpmid = np.zeros([3, 94])
-    allunc = np.zeros([3, 94])
+    nyrs = np.arange(2007, settings["projection_end_year"] + 1).size
+    allpmid = np.zeros([3, nyrs])
+    allunc = np.zeros([3, nyrs])
 
     # Estimate internal variability from de-trended gauge data
     tg_years_arr = np.array(tg_years, dtype='int')
-    vals = np.ma.masked_values(tg_amsl, -99999.)
+    vals = np.ma.masked_invalid(tg_amsl)
     IntV = compute_variability(tg_years_arr / 1000., vals / 1.)
 
     # Get the values of the multi-index: years and percentile values
@@ -53,7 +54,7 @@ def compute_uncertainties(df_r_list, scenarios, tg_years, tg_amsl):
     # Model Uncertainty (already expressed as 90% confidence interval)
     UncM = np.mean(allunc, axis=0)
 
-    return years, UncS, UncM, IntV
+    return years, np.squeeze(UncS), np.squeeze(UncM), IntV
 
 
 def compute_variability(x, y, factor=1.645):
@@ -64,10 +65,8 @@ def compute_variability(x, y, factor=1.645):
     :param factor: multiplication factor
     :return: internal variability component of uncertainty
     """
-    mask = np.ma.getmask(y)
-    index = np.where(mask is True)
-    new_x = np.delete(x, index)
-    new_y = np.delete(y, index)
+    new_x = x[~y.mask]
+    new_y = y[~y.mask]
     fit = np.polyfit(new_x, new_y, 1)
     fit_data = new_x * fit[0] + fit[1]
     stdev = np.std(new_y - fit_data)
@@ -711,6 +710,37 @@ def plot_figure_seven(g_df_list, r_df_list, site_name, scenarios, fig_dir):
     plt.close()
 
 
+def plot_figure_eight(df_r_list, tg_years, tg_amsl, site_name, scenarios, fig_dir):
+    years, UncS, UncM, IntV = compute_uncertainties(df_r_list, scenarios, tg_years, tg_amsl)
+    
+    IntV = np.full(UncS.shape, IntV)
+    
+    UncT = np.sum([UncS, UncM, IntV], axis=0)
+    
+    fig = plt.figure(figsize=(4.5, 4.5))
+    
+    ax = fig.add_subplot(111)
+    matplotlib.rcParams['font.size'] = 9
+    
+    ax.fill_between(years, 0, UncM/UncT, color='blue', label='Model')
+    ax.fill_between(years, UncM/UncT, (UncM+UncS)/UncT, color='green', label='Scenario')
+    ax.fill_between(years, (UncM + UncS)/UncT, 1, color='orange', label='Variability')
+    
+    ax.legend(loc='lower right', fancybox=True, framealpha=1)
+    ax.set_title(f'{location_string(site_name)}')
+    ax.set_xlabel('Year')
+    ax.set_ylabel('Fraction of total variance')
+    
+    ax.set_xlim(2000, years[-1])
+    ax.set_ylim(0, 1)
+    
+    fig.tight_layout()
+    outfile = f'{fig_dir}08_{site_name}_uncertainty_budget.png'
+    plt.savefig(outfile, dpi=200, format='png')
+    plt.show()
+    plt.close()
+
+
 def plot_tg_data(ax, nflag, flag, tg_years, non_missing, tg_amsl, tg_name):
     """
     Plot the annual mean sea levels from the tide gauge data.
@@ -787,7 +817,9 @@ def read_IPCC_AR5_Levermann_proj(scenarios, refname='sum'):
     print('running function read_IPCC_AR5_Levermann_proj')
 
     # Directory of Monte Carlo time series for new projections
-    mcdir = settings["montecarlodir"]
+    mcdir = choose_montecarlo_dir()
+    
+    nyrs = np.arange(2007, settings["projection_end_year"] + 1).size
 
     ar5_low = []
     ar5_mid = []
@@ -795,11 +827,11 @@ def read_IPCC_AR5_Levermann_proj(scenarios, refname='sum'):
 
     for sce in scenarios:
         reflow = iris.load_cube(
-            os.path.join(mcdir, sce + '_' + refname + 'lower.nc')).data
+            os.path.join(mcdir, sce + '_' + refname + 'lower.nc'))[:nyrs].data
         refmid = iris.load_cube(
-            os.path.join(mcdir, sce + '_' + refname + 'mid.nc')).data
+            os.path.join(mcdir, sce + '_' + refname + 'mid.nc'))[:nyrs].data
         refupp = iris.load_cube(
-            os.path.join(mcdir, sce + '_' + refname + 'upper.nc')).data
+            os.path.join(mcdir, sce + '_' + refname + 'upper.nc'))[:nyrs].data
 
         ar5_low.append(reflow)
         ar5_mid.append(refmid)
@@ -964,6 +996,12 @@ def main():
         # Same style as Figure 4 Howard and Palmer (2019)
         plot_figure_seven(g_df_list, r_df_list, df_loc, rcp_scenarios,
                           sealev_fdir)
+        
+        # Figure 8
+        # Uncertainty budget for local sea level projections, 
+        # based on Hawkins & Sutton (2009)
+        plot_figure_eight(r_df_list, tg_years, tg_amsl, df_loc, 
+                          rcp_scenarios, sealev_fdir)
 
 
 if __name__ == '__main__':