diff --git a/ENDA.py b/ENDA.py
index 17368f29a77fc53ab50ab1d623ee76f5088bd6c7..4c573296f39f8167cf98000a3fc016642696a97a 100644
--- a/ENDA.py
+++ b/ENDA.py
@@ -20,6 +20,9 @@ from models import CBL
 
 verbose = True
 
+# Initiate a random number generator with preset seed, to ensure reproducibility
+RNG = np.random.default_rng(seed=240718)
+
 def transform_none(x,kind):
     if kind=='dir':
         return x
@@ -349,7 +352,8 @@ class cycle:
             # perturb the state here.
             if isinstance(nr,CBL):
                 if cbl_settings["simulate_error_growth"]:
-                    x0[:-cbl_settings["parameter_number"],:] += np.random.normal(scale=cbl_settings["error_growth_perturbations_amplitude"],size=nens)[None,:]
+                    x0[:-cbl_settings["parameter_number"],:] += RNG.normal(scale=cbl_settings["error_growth_perturbations_amplitude"],size=nens)[None,:]
+                    #x0[:-cbl_settings["parameter_number"],:] += np.random.normal(scale=cbl_settings["error_growth_perturbations_amplitude"],size=nens)[None,:]
 
             # Save initial conditions for next cycle
             da.update(x0)
@@ -442,7 +446,8 @@ class experiment:
                         variable = nr.history['v'][:,time_index]
                     truths[i,j] = observation_operator(variable,state_coordinates,obs_coordinate)
                     observations[i,j] = truths[i,j]+\
-                        np.random.normal(0,self.obs_error_sdev_generate[j])
+                        RNG.normal(0,self.obs_error_sdev_generate[j])
+                        #np.random.normal(0,self.obs_error_sdev_generate[j])
 
             # Store truths and observations
             self.truths = truths
diff --git a/models.py b/models.py
index 9261f94763cba46fe4aad34b7e61d97ea8806cfc..e9136eceed8117aadce7689f9cc6b708139f2ae3 100644
--- a/models.py
+++ b/models.py
@@ -5,6 +5,9 @@ from scipy.interpolate import CubicSpline
 
 verbose = True
 
+# Initiate a random number generator with preset seed, to ensure reproducibility
+RNG = np.random.default_rng(seed=240718)
+
 class CBL:
     def __init__(self,settings):
         for k, v in settings.items():
@@ -98,10 +101,13 @@ class CBL:
                 randomsize = 1
             
             if self.perturbations_type == "random" or self.perturbations_type == "uniform":
-                ppt = np.random.randn(randomsize,self.nens)
+                ppt = RNG.standard_normal((randomsize,self.nens))
+                #ppt = np.random.randn(randomsize,self.nens)
                 if self.is_bwind:
-                    ppu = np.random.randn(self.nens,randomsize)
-                    ppv = np.random.randn(self.nens,randomsize)
+                    ppu = RNG.standard_normal((self.nens,randomsize))
+                    ppv = RNG.standard_normal((self.nens,randomsize))
+                    #ppu = np.random.randn(self.nens,randomsize)
+                    #ppv = np.random.randn(self.nens,randomsize)
 
             # Smooth perturbations are slightly more complicated
             if self.perturbations_type == "smooth":
@@ -114,15 +120,18 @@ class CBL:
                 # Draw random perturbations, then interpolate
                 randomsize = self.perturbations_smooth_number
                 ppt = np.zeros((self.nz,self.nens))+np.nan
-                pert_t = np.random.randn(randomsize,self.nens)
+                pert_t = RNG.standard_normal((randomsize,self.nens))
+                #pert_t = np.random.randn(randomsize,self.nens)
                 for n in range(self.nens):
                     f = CubicSpline(ipert,pert_t[:,n])
                     ppt[:,n] = f(np.arange(self.nz))
                 if self.is_bwind:
                     ppu = np.zeros((self.nz,self.nens))+np.nan
                     ppv = np.zeros((self.nz,self.nens))+np.nan
-                    pert_u = np.random.randn(randomsize,self.nens)
-                    pert_v = np.random.randn(randomsize,self.nens)
+                    pert_u = RNG.standard_normal((randomsize,self.nens))
+                    pert_v = RNG.standard_normal((randomsize,self.nens))
+                    #pert_u = np.random.randn(randomsize,self.nens)
+                    #pert_v = np.random.randn(randomsize,self.nens)
                     for n in range(self.nens):
                         f = CubicSpline(ipert,pert_u[:,n])
                         ppu[:,n] = f(np.arange(self.nz))
@@ -212,7 +221,8 @@ class CBL:
         pp = np.zeros((self.parameter_number,self.nens))
 
         for k in range(-self.parameter_number,0):
-            dum = np.random.uniform(self.parameter_ensemble_min[k], self.parameter_ensemble_max[k], size=self.nens)
+            dum = RNG.uniform(self.parameter_ensemble_min[k], self.parameter_ensemble_max[k], size=self.nens)
+            #dum = np.random.uniform(self.parameter_ensemble_min[k], self.parameter_ensemble_max[k], size=self.nens)
             pp[k,:] = self.parameter_transform[k](dum , kind='dir')
 
         return pp
@@ -307,7 +317,8 @@ class CBL:
                 H0 = Hmax
             
             # Add random perturbations to the initial value
-            H0 += np.random.normal(scale=H0_perturbation_ampl_init)
+            H0 += RNG.normal(scale=H0_perturbation_ampl_init)
+            #H0 += np.random.normal(scale=H0_perturbation_ampl_init)
 
             # Set the surface momentum flux (ustar)
             ustar = 0
@@ -357,7 +368,8 @@ class CBL:
 
             # Add time-dependent surface perturbations
             # Then compute sensible heat flux and integrate T equation
-            H[0,j]     = H0 + np.random.normal(scale=H0_perturbation_ampl_time)
+            H[0,j]     = H0 + RNG.normal(scale=H0_perturbation_ampl_time)
+            #H[0,j]     = H0 + np.random.normal(scale=H0_perturbation_ampl_time)
             H[1:-1,j]  = -K[1:-1,j]*( (thetap[1:]-thetap[:-1])*rdz - gammac)
             H[nz,j]    = 2*H[nz-1,j]-H[nz-2,j]
             theta[:,j] = thetap[:]-dt*rdz*(H[1:,j]-H[:-1,j])