diff --git a/models.py b/models.py
index 1eba7505483182b1bdbd8d5cc3aeee595beff195..bf98baaa2f34723472eecb46b0e71bfc4dcd2de2 100644
--- a/models.py
+++ b/models.py
@@ -240,6 +240,8 @@ class CBL:
         z0 = self.z0
         theta_0 = self.theta_0
         Hmax = self.Hmax
+        H0_perturbation_ampl_init = self.H0_perturbation_ampl_init
+        H0_perturbation_ampl_time = self.H0_perturbation_ampl_time
         exct = self.exct
         ug = self.ug
         vg = self.vg
@@ -311,6 +313,9 @@ class CBL:
             else:
                 H0 = Hmax
             
+            # Add random perturbations to the initial value
+            H0 += np.random.normal(scale=H0_perturbation_ampl_init)
+
             # Set the surface momentum flux (ustar)
             ustar = 0
 
@@ -357,8 +362,9 @@ class CBL:
             # Inverse grid spacing
             rdz = 1./self.dz
 
-            # Compute sensible heat flux and integrate T equation
-            H[0,j]     = H0
+            # 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[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])