diff --git a/extreme_data/meteo_france_data/scm_models_data/visualization/plot_utils.py b/extreme_data/meteo_france_data/scm_models_data/visualization/plot_utils.py
index c16465d2ac9ed1f09f0a736ba8893d3a3602bed2..b029c576a3653017a0e35e7c8c77d2b5c2f851cd 100644
--- a/extreme_data/meteo_france_data/scm_models_data/visualization/plot_utils.py
+++ b/extreme_data/meteo_france_data/scm_models_data/visualization/plot_utils.py
@@ -20,7 +20,7 @@ def plot_against_altitude(altitudes, ax, massif_id, massif_name, values):
     ax.plot(altitudes, values, color=color, linewidth=2, label=massif_name, linestyle=linestyle)
 
 
-def load_plot(cmap, graduation, label, massif_name_to_value, altitude, fit_method):
+def load_plot(cmap, graduation, label, massif_name_to_value, altitude, fit_method, add_x_label=True):
     max_abs_change = max([abs(e) for e in massif_name_to_value.values()])
     ticks, labels = ticks_values_and_labels_for_percentages(graduation=graduation, max_abs_change=max_abs_change)
     min_ratio = -max_abs_change
@@ -47,5 +47,6 @@ def load_plot(cmap, graduation, label, massif_name_to_value, altitude, fit_metho
                                   )
     ax.get_xaxis().set_visible(True)
     ax.set_xticks([])
-    ax.set_xlabel('Altitude = {}m'.format(altitude), fontsize=15)
-    ax.set_title('Fit method is {}'.format(fit_method))
+    if add_x_label:
+        ax.set_xlabel('Altitude = {}m'.format(altitude), fontsize=15)
+        ax.set_title('Fit method is {}'.format(fit_method))
diff --git a/extreme_data/meteo_france_data/scm_models_data/visualization/study_visualizer.py b/extreme_data/meteo_france_data/scm_models_data/visualization/study_visualizer.py
index ed322e2158d8b39515011f207d3562d9d0f6d55b..34174e08676015ea7f00839da2d87e32fba4b7b8 100644
--- a/extreme_data/meteo_france_data/scm_models_data/visualization/study_visualizer.py
+++ b/extreme_data/meteo_france_data/scm_models_data/visualization/study_visualizer.py
@@ -708,18 +708,19 @@ class StudyVisualizer(VisualizationParameters):
 
     # PLot functions that should be common
 
-    def plot_map(self, cmap, fit_method, graduation, label, massif_name_to_value, plot_name):
-        load_plot(cmap, graduation, label, massif_name_to_value, self.study.altitude, fitmethod_to_str(fit_method))
+    def plot_map(self, cmap, fit_method, graduation, label, massif_name_to_value, plot_name, add_x_label=True):
+        load_plot(cmap, graduation, label, massif_name_to_value, self.study.altitude, fitmethod_to_str(fit_method),
+                  add_x_label=add_x_label)
         self.plot_name = plot_name
         # self.show_or_save_to_file(add_classic_title=False, tight_layout=True, no_title=True, dpi=500)
         self.show_or_save_to_file(add_classic_title=False, no_title=True, dpi=500)
         plt.close()
 
-    def plot_abstract(self, massif_name_to_value, label, plot_name, fit_method='', graduation=10.0, cmap=plt.cm.bwr):
+    def plot_abstract(self, massif_name_to_value, label, plot_name, fit_method='', graduation=10.0, cmap=plt.cm.bwr, add_x_label=True):
         # Regroup the plot by altitudes
         plot_name1 = '{}/{}'.format(self.study.altitude, plot_name)
         # Regroup the plot by type of plot also
         plot_name2 = '{}/{}'.format(plot_name.split()[0], plot_name)
         for plot_name in [plot_name1, plot_name2]:
-            self.plot_map(cmap, fit_method, graduation, label, massif_name_to_value, plot_name)
+            self.plot_map(cmap, fit_method, graduation, label, massif_name_to_value, plot_name, add_x_label)
 
diff --git a/extreme_fit/model/margin_model/linear_margin_model/linear_margin_model.py b/extreme_fit/model/margin_model/linear_margin_model/linear_margin_model.py
index 538fa81612f4aedeb2b17054aad521aaa13350f5..a8f5af9974140904a9c0034dbc8a891d966464d0 100644
--- a/extreme_fit/model/margin_model/linear_margin_model/linear_margin_model.py
+++ b/extreme_fit/model/margin_model/linear_margin_model/linear_margin_model.py
@@ -17,10 +17,9 @@ class LinearMarginModel(ParametricMarginModel):
     def load_margin_function(self, param_name_to_dims=None):
         assert param_name_to_dims is not None, 'LinearMarginModel cannot be used for sampling/fitting \n' \
                                                'load_margin_functions needs to be implemented in child class'
-        # Load sample coef
-        coef_sample = self.param_name_to_linear_coef(param_name_and_dim_to_coef=self.params_sample)
+        param_name_to_coef = self.param_name_to_linear_coef(param_name_and_dim_to_coef=self.params_sample)
         return LinearMarginFunction(coordinates=self.coordinates,
-                                    param_name_to_coef=coef_sample,
+                                    param_name_to_coef=param_name_to_coef,
                                     param_name_to_dims=param_name_to_dims,
                                     starting_point=self.starting_point,
                                     params_class=self.params_class)
diff --git a/projects/contrasting_trends_in_snow_loads/snowfall_versus_time_and_altitude/main_snowfall_article.py b/projects/contrasting_trends_in_snow_loads/snowfall_versus_time_and_altitude/main_snowfall_article.py
index cbc4c2af911af72196d53a379b1956e64df5bcd3..33c73e08ecd5391ef601eb1d3d51482280ac8ebc 100644
--- a/projects/contrasting_trends_in_snow_loads/snowfall_versus_time_and_altitude/main_snowfall_article.py
+++ b/projects/contrasting_trends_in_snow_loads/snowfall_versus_time_and_altitude/main_snowfall_article.py
@@ -71,7 +71,7 @@ def intermediate_result(altitudes, massif_names=None,
             _ = compute_minimized_aic(visualizer)
 
     # Plots
-    # validation_plot(altitude_to_visualizer)
+    validation_plot(altitude_to_visualizer)
     plot_snowfall_mean(altitude_to_visualizer)
 
 
diff --git a/projects/contrasting_trends_in_snow_loads/snowfall_versus_time_and_altitude/snowfall_plot.py b/projects/contrasting_trends_in_snow_loads/snowfall_versus_time_and_altitude/snowfall_plot.py
index 5d17a9583ea1f72bc49f81aeb65a1ad29bec47fa..dfbbf81ed3123f5a5b5581050027cd1bfc61a4ef 100644
--- a/projects/contrasting_trends_in_snow_loads/snowfall_versus_time_and_altitude/snowfall_plot.py
+++ b/projects/contrasting_trends_in_snow_loads/snowfall_versus_time_and_altitude/snowfall_plot.py
@@ -24,16 +24,23 @@ def plot_snowfall_mean(altitude_to_visualizer: Dict[int, StudyVisualizerForMeanV
     visualizer = list(altitude_to_visualizer.values())[0]
     # Plot the curve for the evolution of the mean
     massif_name_to_a, massif_name_to_b, massif_name_to_r2_score = plot_mean(altitude_to_visualizer, derivative=False)
-    # Plot map with the coefficient a
-    visualizer.plot_abstract_fast(massif_name_to_a, label='a')
-    visualizer.plot_abstract_fast(massif_name_to_b, label='b')
-    visualizer.plot_abstract_fast(massif_name_to_r2_score, label='r2')
+    # Augmentation every km
+    massif_name_to_augmentation_every_km = {m: a * 1000 for m, a in massif_name_to_a.items()}
+    visualizer.plot_abstract_fast(massif_name_to_augmentation_every_km,
+                                  label='Augmentation of mean annual maxima for every km of elevation ()',
+                                  add_x_label=False)
+    # Value at 2000 m
+    massif_name_to_mean_at_2000 = {m: a * 2000 + massif_name_to_b[m] for m, a in massif_name_to_a.items()}
+    visualizer.plot_abstract_fast(massif_name_to_mean_at_2000, label='Mean annual maxima of at 2000 m ()',
+                                  add_x_label=False)
+    # R2 score
+    visualizer.plot_abstract_fast(massif_name_to_r2_score, label='r2', graduation=0.1, add_x_label=False)
 
 
 def plot_mean(altitude_to_visualizer: Dict[int, StudyVisualizerForMeanValues], derivative=False):
+    ax = plt.gca()
     massif_name_to_linear_regression_result = {}
 
-    altitudes = list(altitude_to_visualizer.keys())
     visualizers = list(altitude_to_visualizer.values())
     visualizer = visualizers[0]
     study = visualizer.study
@@ -51,22 +58,22 @@ def plot_mean(altitude_to_visualizer: Dict[int, StudyVisualizerForMeanValues], d
             else:
                 moment = 'mean'
                 values = [t.mean_value(year=year) for t in trend_tests]
-            massif_name_to_linear_regression_result[massif_name] = fit_linear_regression(altitudes, values)
-            plot_values_against_altitudes(altitudes_massif, massif_id, massif_name, moment, study, values, visualizer)
+            massif_name_to_linear_regression_result[massif_name] = fit_linear_regression(altitudes_massif, values)
+            plot_values_against_altitudes(ax, altitudes_massif, massif_id, massif_name, moment, study, values, visualizer)
+    plt.close()
 
     return [{m: t[i][0] if i == 0 else t[i]
              for m, t in massif_name_to_linear_regression_result.items()} for i in range(3)]
 
 
-def plot_values_against_altitudes(altitudes, massif_id, massif_name, moment, study, values, visualizer):
-    ax = plt.gca()
+def plot_values_against_altitudes(ax, altitudes, massif_id, massif_name, moment, study, values, visualizer):
     plot_against_altitude(altitudes=altitudes, ax=ax, massif_id=massif_id, massif_name=massif_name, values=values)
     plot_name = '{} annual maxima of {}'.format(moment, SCM_STUDY_CLASS_TO_ABBREVIATION[type(study)])
     ax.set_ylabel('{} ({})'.format(plot_name, study.variable_unit), fontsize=15)
     ax.set_xlabel('altitudes', fontsize=15)
-    lim_down, lim_up = ax.get_ylim()
-    lim_up += (lim_up - lim_down) / 3
-    ax.set_ylim([lim_down, lim_up])
+    # lim_down, lim_up = ax.get_ylim()
+    # lim_up += (lim_up - lim_down) / 3
+    # ax.set_ylim([lim_down, lim_up])
     ax.tick_params(axis='both', which='major', labelsize=13)
     visualizer.plot_name = plot_name
     visualizer.show_or_save_to_file(dpi=500, add_classic_title=False)
diff --git a/projects/contrasting_trends_in_snow_loads/snowfall_versus_time_and_altitude/study_visualizer_for_mean_values.py b/projects/contrasting_trends_in_snow_loads/snowfall_versus_time_and_altitude/study_visualizer_for_mean_values.py
index 334ec2ca1d2c2d25316a870b124a048f30b7c8d7..70e5d768e149f0762c0c2fe1458291cccae86048 100644
--- a/projects/contrasting_trends_in_snow_loads/snowfall_versus_time_and_altitude/study_visualizer_for_mean_values.py
+++ b/projects/contrasting_trends_in_snow_loads/snowfall_versus_time_and_altitude/study_visualizer_for_mean_values.py
@@ -52,6 +52,6 @@ class StudyVisualizerForMeanValues(StudyVisualizerForNonStationaryTrends):
             massif_name_to_relative_difference[massif_name] = relative_diference
         return massif_name_to_relative_difference
 
-    def plot_abstract_fast(self, massif_name_to_value, label, graduation=10.0, cmap=plt.cm.coolwarm):
-        super().plot_abstract(massif_name_to_value, label, label, self.fit_method, graduation, cmap)
+    def plot_abstract_fast(self, massif_name_to_value, label, graduation=10.0, cmap=plt.cm.coolwarm, add_x_label=True):
+        super().plot_abstract(massif_name_to_value, label, label, self.fit_method, graduation, cmap, add_x_label)
 
diff --git a/projects/contrasting_trends_in_snow_loads/snowfall_versus_time_and_altitude/validation_plot.py b/projects/contrasting_trends_in_snow_loads/snowfall_versus_time_and_altitude/validation_plot.py
index 6b6302525b5b472601d79606411b2cfb912fef4e..88a5a5b5f178b4acc503226f8eedc0eff0c46c4f 100644
--- a/projects/contrasting_trends_in_snow_loads/snowfall_versus_time_and_altitude/validation_plot.py
+++ b/projects/contrasting_trends_in_snow_loads/snowfall_versus_time_and_altitude/validation_plot.py
@@ -27,7 +27,7 @@ def plot_shoe_relative_differences_distribution(altitude_to_relative_differences
     width = 150
     ax.boxplot([altitude_to_relative_differences[a] for a in altitudes], positions=altitudes, widths=width)
     ax.set_xlim([min(altitudes) - width, max(altitudes) + width])
-    ylabel = 'Relative difference between empirical mean and parametric mean (%)'
+    ylabel = 'Relative difference between empirical mean and parametric mean (\%)'
     ax.set_ylabel(ylabel)
     ax.set_xlabel('Altitude (m)')
     ax.legend()