diff --git a/extreme_fit/model/abstract_model.py b/extreme_fit/model/abstract_model.py
index 331d275365931a3604c8889057b4616b785df772..cd6e735537e487583060702452dc7bcc9427ed9d 100644
--- a/extreme_fit/model/abstract_model.py
+++ b/extreme_fit/model/abstract_model.py
@@ -1,8 +1,11 @@
 class AbstractModel(object):
 
-    def __init__(self, params_sample=None):
-        self.default_params = None
-        self.user_params_sample = params_sample
+    def __init__(self, params_user=None):
+        self.user_params_sample = params_user
+
+    @property
+    def default_params(self):
+        return None
 
     @property
     def params_sample(self) -> dict:
diff --git a/extreme_fit/model/margin_model/abstract_margin_model.py b/extreme_fit/model/margin_model/abstract_margin_model.py
index 8358b41edb49334559e7f252e61946e3c76d8639..d87146e67cfeeac778ffea7d8ebca00837437629 100644
--- a/extreme_fit/model/margin_model/abstract_margin_model.py
+++ b/extreme_fit/model/margin_model/abstract_margin_model.py
@@ -18,9 +18,9 @@ class AbstractMarginModel(AbstractModel, ABC):
         -margin_function
     """
 
-    def __init__(self, coordinates: AbstractCoordinates, params_sample=None,
+    def __init__(self, coordinates: AbstractCoordinates, params_user=None,
                  params_class=GevParams):
-        super().__init__(params_sample)
+        super().__init__(params_user)
         assert isinstance(coordinates, AbstractCoordinates), type(coordinates)
         self.coordinates = coordinates
         self.margin_function = None  # type: AbstractMarginFunction
diff --git a/extreme_fit/model/margin_model/linear_margin_model/abstract_temporal_linear_margin_model.py b/extreme_fit/model/margin_model/linear_margin_model/abstract_temporal_linear_margin_model.py
index 5e1d147b8fa6a899c2d25666fad6c7351c61d27c..18c2afd59aad785e39bec63d4298d753ffa50564 100644
--- a/extreme_fit/model/margin_model/linear_margin_model/abstract_temporal_linear_margin_model.py
+++ b/extreme_fit/model/margin_model/linear_margin_model/abstract_temporal_linear_margin_model.py
@@ -20,13 +20,13 @@ class AbstractTemporalLinearMarginModel(LinearMarginModel):
     """Linearity only with respect to the temporal coordinates"""
 
     def __init__(self, coordinates: AbstractCoordinates,
-                 params_sample=None, starting_point=None,
+                 params_user=None, starting_point=None,
                  fit_method=MarginFitMethod.is_mev_gev_fit,
                  nb_iterations_for_bayesian_fit=5000,
                  params_initial_fit_bayesian=None,
                  type_for_MLE="GEV",
                  params_class=GevParams):
-        super().__init__(coordinates, params_sample, starting_point,
+        super().__init__(coordinates, params_user, starting_point,
                          params_class)
         self.type_for_mle = type_for_MLE
         self.params_initial_fit_bayesian = params_initial_fit_bayesian
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 f6091c29514a851159dea019cf26b8700c7ab69b..e3c86ebfecf17c99e51094a762e35c93e3b993f6 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
@@ -12,15 +12,11 @@ class LinearMarginModel(ParametricMarginModel):
         for param_name, coef_list in param_name_to_coef_list.items():
             for idx, coef in enumerate(coef_list, -1):
                 params[(param_name, idx)] = coef
-        return cls(coordinates, params_sample=params, params_class=params_class, **kwargs)
+        return cls(coordinates, params_user=params, params_class=params_class, **kwargs)
 
     def load_margin_functions(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 default params (with a dictionary format to enable quick replacement)
-        # IMPORTANT: Using a dictionary format enable using the default/user params methodology
-        self.default_params = self.default_param_name_and_dim_to_coef
-
         # Load sample coef
         coef_sample = self.param_name_to_linear_coef(param_name_and_dim_to_coef=self.params_sample)
         self.margin_function = LinearMarginFunction(coordinates=self.coordinates,
@@ -30,7 +26,7 @@ class LinearMarginModel(ParametricMarginModel):
                                                     params_class=self.params_class)
 
     @property
-    def default_param_name_and_dim_to_coef(self) -> dict:
+    def default_params(self) -> dict:
         default_intercept = 1
         default_slope = 0.01
         param_name_and_dim_to_coef = {}
diff --git a/extreme_fit/model/margin_model/linear_margin_model/temporal_linear_margin_models.py b/extreme_fit/model/margin_model/linear_margin_model/temporal_linear_margin_models.py
index 1f8c1e85187be003fb95527694f9bff71fdc0afc..b84dcd7f83dc75c6c699ac56a2464e9f5f72c6c5 100644
--- a/extreme_fit/model/margin_model/linear_margin_model/temporal_linear_margin_models.py
+++ b/extreme_fit/model/margin_model/linear_margin_model/temporal_linear_margin_models.py
@@ -71,9 +71,9 @@ class NonStationaryLocationAndScaleTemporalModel(AbstractTemporalLinearMarginMod
 class GumbelTemporalModel(StationaryTemporalModel):
 
     def __init__(self, coordinates: AbstractCoordinates,
-                 params_sample=None, starting_point=None, fit_method=MarginFitMethod.is_mev_gev_fit,
+                 params_user=None, starting_point=None, fit_method=MarginFitMethod.is_mev_gev_fit,
                  nb_iterations_for_bayesian_fit=5000, params_initial_fit_bayesian=None):
-        super().__init__(coordinates, params_sample, starting_point, fit_method,
+        super().__init__(coordinates, params_user, starting_point, fit_method,
                          nb_iterations_for_bayesian_fit, params_initial_fit_bayesian, type_for_MLE="Gumbel")
 
 
diff --git a/extreme_fit/model/margin_model/parametric_margin_model.py b/extreme_fit/model/margin_model/parametric_margin_model.py
index 3e638248a11397b69e498b3e9b7944a530b58ca9..6aba0d054ba4807be604a151a70f4a0fc566df41 100644
--- a/extreme_fit/model/margin_model/parametric_margin_model.py
+++ b/extreme_fit/model/margin_model/parametric_margin_model.py
@@ -17,7 +17,7 @@ from spatio_temporal_dataset.coordinates.abstract_coordinates import AbstractCoo
 class ParametricMarginModel(AbstractMarginModel, ABC):
 
     def __init__(self, coordinates: AbstractCoordinates,
-                 params_sample=None, starting_point=None, params_class=GevParams,
+                 params_user=None, starting_point=None, params_class=GevParams,
                  fit_method=MarginFitMethod.spatial_extremes_mle):
         """
         :param starting_point: starting coordinate for the temporal trend
@@ -26,7 +26,7 @@ class ParametricMarginModel(AbstractMarginModel, ABC):
         self.starting_point = starting_point
         self.margin_function = None  # type: ParametricMarginFunction
         self.drop_duplicates = True
-        super().__init__(coordinates, params_sample, params_class)
+        super().__init__(coordinates, params_user, params_class)
 
     def fitmargin_from_maxima_gev(self, data: np.ndarray, df_coordinates_spat: pd.DataFrame,
                                   df_coordinates_temp: pd.DataFrame) -> ResultFromSpatialExtreme:
diff --git a/extreme_fit/model/margin_model/spline_margin_model.py b/extreme_fit/model/margin_model/spline_margin_model.py
index f9fec6a4209463bb63624e38a10d3f10f08ac7b3..83addd3c7a9a1f1bba4365638459e2a7372ab6f3 100644
--- a/extreme_fit/model/margin_model/spline_margin_model.py
+++ b/extreme_fit/model/margin_model/spline_margin_model.py
@@ -12,8 +12,8 @@ from spatio_temporal_dataset.coordinates.abstract_coordinates import AbstractCoo
 class SplineMarginModel(ParametricMarginModel):
 
     def __init__(self, coordinates: AbstractCoordinates,
-                 params_sample=None):
-        super().__init__(coordinates, params_sample)
+                 params_user=None):
+        super().__init__(coordinates, params_user)
 
     def load_margin_functions(self, param_name_to_dims: Dict[str, List[int]] = None,
                               param_name_to_coef: Dict[str, AbstractCoef] = None,
diff --git a/extreme_fit/model/max_stable_model/abstract_max_stable_model.py b/extreme_fit/model/max_stable_model/abstract_max_stable_model.py
index 8e7cb652f55d7ce08b80f59189527813b1e5216d..09513822dda16ee94d5f5020261b5804c758cdac 100644
--- a/extreme_fit/model/max_stable_model/abstract_max_stable_model.py
+++ b/extreme_fit/model/max_stable_model/abstract_max_stable_model.py
@@ -14,8 +14,8 @@ from spatio_temporal_dataset.coordinates.abstract_coordinates import AbstractCoo
 
 class AbstractMaxStableModel(AbstractModel):
 
-    def __init__(self, params_sample=None):
-        super().__init__(params_sample)
+    def __init__(self, params_user=None):
+        super().__init__(params_user)
         self.cov_mod = None
 
     @property
@@ -105,12 +105,15 @@ class CovarianceFunction(Enum):
 
 class AbstractMaxStableModelWithCovarianceFunction(AbstractMaxStableModel):
 
-    def __init__(self, params_sample=None,
+    def __init__(self, params_user=None,
                  covariance_function: CovarianceFunction = None):
-        super().__init__(params_sample)
+        super().__init__(params_user)
         assert covariance_function is not None
         self.covariance_function = covariance_function
-        self.default_params = {
+
+    @property
+    def default_params(self):
+        return {
             'range': 3,
             'smooth': 0.5,
             'nugget': 0.5
diff --git a/extreme_fit/model/max_stable_model/max_stable_models.py b/extreme_fit/model/max_stable_model/max_stable_models.py
index bcc3ad8332b1bdb4de97a9f0003e50bcad8e094b..7ba235478f393c9d4c00db77121b7119caa28a0e 100644
--- a/extreme_fit/model/max_stable_model/max_stable_models.py
+++ b/extreme_fit/model/max_stable_model/max_stable_models.py
@@ -9,7 +9,10 @@ class Smith(AbstractMaxStableModel):
     def __init__(self, *args, **kwargs):
         super().__init__(*args, **kwargs)
         self.cov_mod = 'gauss'
-        self.default_params = {
+
+    @property
+    def default_params(self):
+        return {
             'var': 1,
             'cov11': 1,
             'cov12': 0,
@@ -29,7 +32,10 @@ class BrownResnick(AbstractMaxStableModel):
     def __init__(self, *args, **kwargs):
         super().__init__(*args, **kwargs)
         self.cov_mod = 'brown'
-        self.default_params = {
+
+    @property
+    def default_params(self):
+        return {
             'range': 3,
             'smooth': 0.5,
         }
@@ -40,7 +46,6 @@ class Schlather(AbstractMaxStableModelWithCovarianceFunction):
     def __init__(self, *args, **kwargs):
         super().__init__(*args, **kwargs)
         self.cov_mod = self.covariance_function.name
-        self.default_params.update({})
 
 
 class Geometric(AbstractMaxStableModelWithCovarianceFunction):
@@ -48,7 +53,12 @@ class Geometric(AbstractMaxStableModelWithCovarianceFunction):
     def __init__(self, *args, **kwargs):
         super().__init__(*args, **kwargs)
         self.cov_mod = 'g' + self.covariance_function.name
-        self.default_params.update({'sigma2': 0.5})
+
+    @property
+    def default_params(self):
+        d = super().default_params
+        d.update({'sigma2': 0.5})
+        return d
 
 
 class ExtremalT(AbstractMaxStableModelWithCovarianceFunction):
@@ -56,7 +66,12 @@ class ExtremalT(AbstractMaxStableModelWithCovarianceFunction):
     def __init__(self, *args, **kwargs):
         super().__init__(*args, **kwargs)
         self.cov_mod = 't' + self.covariance_function.name
-        self.default_params.update({'DoF': 2})
+
+    @property
+    def default_params(self):
+        d = super().default_params
+        d.update({'DoF': 2})
+        return d
 
 
 class ISchlather(AbstractMaxStableModelWithCovarianceFunction):
@@ -64,4 +79,9 @@ class ISchlather(AbstractMaxStableModelWithCovarianceFunction):
     def __init__(self, *args, **kwargs):
         super().__init__(*args, **kwargs)
         self.cov_mod = 'i' + self.covariance_function.name
-        self.default_params.update({'alpha': 0.5})
+
+    @property
+    def default_params(self):
+        d = super().default_params
+        d.update({'alpha': 0.5})
+        return d
diff --git a/test/test_extreme_fit/test_model/test_margin_model.py b/test/test_extreme_fit/test_model/test_margin_model.py
index 626a624b3ff4ca9fd9d71fc612d60c5086c3accf..b26149b5e73f06b29213eefca1174e24bac467a3 100644
--- a/test/test_extreme_fit/test_model/test_margin_model.py
+++ b/test/test_extreme_fit/test_model/test_margin_model.py
@@ -19,7 +19,7 @@ class TestVisualizationLinearMarginModel(unittest.TestCase):
 
     def test_example_visualization_1D(self):
         coordinates = LinSpaceSpatialCoordinates.from_nb_points(nb_points=self.nb_points)
-        self.margin_model = self.margin_model_class(coordinates=coordinates, params_sample={(GevParams.SHAPE, 0): 0.02})
+        self.margin_model = self.margin_model_class(coordinates=coordinates, params_user={(GevParams.SHAPE, 0): 0.02})
 
     def test_example_visualization_2D_spatial(self):
         spatial_coordinates = LinSpaceSpatial2DCoordinates.from_nb_points(nb_points=self.nb_points)
@@ -42,7 +42,7 @@ class TestVisualizationSplineMarginModel(unittest.TestCase):
 
     def test_example_visualization_1D_spline(self):
         coordinates = LinSpaceSpatialCoordinates.from_nb_points(nb_points=self.nb_points, start=0.0)
-        self.margin_model = self.margin_model_class(coordinates=coordinates, params_sample={(GevParams.SHAPE, 1): 0.02})
+        self.margin_model = self.margin_model_class(coordinates=coordinates, params_user={(GevParams.SHAPE, 1): 0.02})
 
     def test_example_visualization_2D_spatial_spline(self):
         spatial_coordinates = LinSpaceSpatial2DCoordinates.from_nb_points(nb_points=self.nb_points)
diff --git a/test/test_unitary_r_packages/test_spatial_extreme/test_rmaxstab/test_rmaxstab_without_margin.py b/test/test_unitary_r_packages/test_spatial_extreme/test_rmaxstab/test_rmaxstab_without_margin.py
index 5e32c74009af2be983b6d841110818e7fdd6eb63..92e9cee552f938a905fc5e686f52eb0fed99d945 100644
--- a/test/test_unitary_r_packages/test_spatial_extreme/test_rmaxstab/test_rmaxstab_without_margin.py
+++ b/test/test_unitary_r_packages/test_spatial_extreme/test_rmaxstab/test_rmaxstab_without_margin.py
@@ -26,7 +26,7 @@ class TestRMaxStab(TestUnitaryAbstract):
         coordinates = AbstractSpatialCoordinates.from_df(df)
         # Load max stable model
         params_sample = {'range': 3, 'smooth': 0.5, 'nugget': 0}
-        max_stable_model = Schlather(covariance_function=CovarianceFunction.whitmat, params_sample=params_sample)
+        max_stable_model = Schlather(covariance_function=CovarianceFunction.whitmat, params_user=params_sample)
         return coordinates, max_stable_model
 
     @property