diff --git a/experiment/paper_past_snow_loads/result_trends_and_return_levels/main_result_trends_and_return_levels.py b/experiment/paper_past_snow_loads/result_trends_and_return_levels/main_result_trends_and_return_levels.py
index f0b21ab00c66bea03ae806ae6f954fa83948f011..796edeae09b7a27400473169eaa830dea82ee1f4 100644
--- a/experiment/paper_past_snow_loads/result_trends_and_return_levels/main_result_trends_and_return_levels.py
+++ b/experiment/paper_past_snow_loads/result_trends_and_return_levels/main_result_trends_and_return_levels.py
@@ -47,6 +47,9 @@ def intermediate_result(altitudes, massif_names=None,
     # Load altitude to visualizer
     altitude_to_visualizer = load_altitude_to_visualizer(altitudes, massif_names, non_stationary_uncertainty,
                                                          study_class, uncertainty_methods)
+    # Load variable object efficiently
+    for v in altitude_to_visualizer.values():
+        _ = v.study.year_to_variable_object
     # Plot trends
     visualizers = list(altitude_to_visualizer.values())
     if multiprocessing:
@@ -70,7 +73,7 @@ def major_result():
     uncertainty_methods = [ConfidenceIntervalMethodFromExtremes.my_bayes,
                            ConfidenceIntervalMethodFromExtremes.ci_mle][:]
     massif_names = None
-    for study_class in paper_study_classes[:2]:
+    for study_class in paper_study_classes[1:2]:
         if study_class == CrocusSnowLoadEurocode:
             non_stationary_uncertainty = [False]
         else:
@@ -79,12 +82,12 @@ def major_result():
 
 
 if __name__ == '__main__':
-    # major_result()
-    intermediate_result(altitudes=paper_altitudes, massif_names=['Maurienne'],
-                        uncertainty_methods=[ConfidenceIntervalMethodFromExtremes.my_bayes,
-                           ConfidenceIntervalMethodFromExtremes.ci_mle][:],
-                        non_stationary_uncertainty=[False, True][:],
-                        multiprocessing=False)
+    major_result()
+    # intermediate_result(altitudes=paper_altitudes, massif_names=['Maurienne'],
+    #                     uncertainty_methods=[ConfidenceIntervalMethodFromExtremes.my_bayes,
+    #                        ConfidenceIntervalMethodFromExtremes.ci_mle][1:],
+    #                     non_stationary_uncertainty=[False, True][1:],
+    #                     multiprocessing=True)
     # intermediate_result(altitudes=[900, 1200], massif_names=None)
     # intermediate_result(ALL_ALTITUDES_WITHOUT_NAN)
     # intermediate_result(paper_altitudes)
diff --git a/experiment/paper_past_snow_loads/result_trends_and_return_levels/plot_uncertainty_curves.py b/experiment/paper_past_snow_loads/result_trends_and_return_levels/plot_uncertainty_curves.py
index 9ba39fd55d086c64011321406d991eff5e3a4229..5ad1822272e170a9121befea9436ea83375ac2b2 100644
--- a/experiment/paper_past_snow_loads/result_trends_and_return_levels/plot_uncertainty_curves.py
+++ b/experiment/paper_past_snow_loads/result_trends_and_return_levels/plot_uncertainty_curves.py
@@ -84,8 +84,8 @@ def plot_single_uncertainty_massif_and_non_stationary_context(ax, massif_name, n
     altitudes = list(altitude_to_visualizer.keys())
     visualizer = list(altitude_to_visualizer.values())[0]
     alpha = 0.2
-    legend_size = 20
-    fontsize_label = 20
+    legend_size = 25
+    fontsize_label = 35
     # Display the EUROCODE return level
     eurocode_region = massif_name_to_eurocode_region[massif_name]()
 
@@ -103,13 +103,13 @@ def plot_single_uncertainty_massif_and_non_stationary_context(ax, massif_name, n
 
         # Plot bars of TDRL only in the non stationary case
         if j == 0 and non_stationary_context:
-            plot_tdrl_bars(altitude_to_visualizer, ax, massif_name, valid_altitudes, legend_size, fontsize_label)
+            plot_tdrl_bars(altitude_to_visualizer, ax, massif_name, valid_altitudes, legend_size, legend_size)
 
     ax.legend(loc=2, prop={'size': legend_size})
     # ax.set_ylim([-1, 16])
     ax.set_xlim([200, 1900])
     if massif_name == 'Maurienne':
-        ax.set_ylim([-1, 13])
+        ax.set_ylim([-1.5, 13])
     # add_title(ax, eurocode_region, massif_name, non_stationary_context)
     ax.set_xticks(altitudes)
     ax.tick_params(labelsize=fontsize_label)
@@ -155,7 +155,7 @@ def plot_tdrl_bars(altitude_to_visualizer, ax, massif_name, valid_altitudes, leg
     # ax2.legend(handles=legend_elements, bbox_to_anchor=(0.93, 0.7), loc='upper right')
     # ax2.annotate("Filled symbol = significant trend ", xy=(0.85, 0.5), xycoords='axes fraction', fontsize=7)
     ax2.legend(handles=legend_elements, loc='upper right', prop={'size': legend_size})
-    ax2.annotate("Filled symbol = significant trend ", xy=(0.5, 0.93), xycoords='axes fraction', fontsize=fontsize)
+    ax2.annotate("Filled symbol =\n significant trend ", xy=(0.6, 0.85), xycoords='axes fraction', fontsize=fontsize)
     ax2.set_yticks([])
 
 
diff --git a/experiment/paper_past_snow_loads/result_trends_and_return_levels/plot_uncertainty_histogram.py b/experiment/paper_past_snow_loads/result_trends_and_return_levels/plot_uncertainty_histogram.py
index c1d5f99b74da005a8a2302a080ba09111e2902c9..9fb61fd0b53d8d7fb5ca4614ac783d7cf1ab11b1 100644
--- a/experiment/paper_past_snow_loads/result_trends_and_return_levels/plot_uncertainty_histogram.py
+++ b/experiment/paper_past_snow_loads/result_trends_and_return_levels/plot_uncertainty_histogram.py
@@ -38,10 +38,13 @@ def plot_histogram(altitude_to_visualizer, non_stationary_context):
             bincenters = altitudes + offset
         width = 100
         plot_histogram_ci_method(visualizers, non_stationary_context, ci_method, ax, bincenters, width=width)
+    fontsize_label = 15
+    legend_size = 15
     ax.set_xticks(altitudes)
-    ax.legend(loc='upper left')
-    ax.set_ylabel('Massifs exceeding French standards (\%)')
-    ax.set_xlabel('Altitude (m)')
+    ax.tick_params(labelsize=fontsize_label)
+    ax.legend(loc='upper left', prop={'size': legend_size})
+    ax.set_ylabel('Massifs exceeding French standards (\%)', fontsize=fontsize_label)
+    ax.set_xlabel('Altitude (m)', fontsize=fontsize_label)
     ax.set_ylim([0, 100])
     ax.set_yticks([10 * i for i in range(11)])
     visualizer.plot_name = 'Percentages of exceedance with non_stationary={}'.format(non_stationary_context)