diff --git a/src/xtgeo/surface/_regsurf_gridding.py b/src/xtgeo/surface/_regsurf_gridding.py
index 0f760ac4f..3bb6d5445 100644
--- a/src/xtgeo/surface/_regsurf_gridding.py
+++ b/src/xtgeo/surface/_regsurf_gridding.py
@@ -1,8 +1,10 @@
 # -*- coding: utf-8 -*-
 """Do gridding from 3D parameters"""
-
+from __future__ import annotations
 
 import warnings
+from collections.abc import Callable
+from typing import TYPE_CHECKING
 
 import numpy as np
 import numpy.ma as ma
@@ -12,6 +14,9 @@
 import xtgeo
 from xtgeo.common import null_logger
 
+if TYPE_CHECKING:
+    from xtgeo.surface import RegularSurface
+
 logger = null_logger(__name__)
 
 # Note: 'self' is an instance of RegularSurface
@@ -322,21 +327,25 @@ def surf_fill(self, fill_value=None):
         logger.info("Do fill... DONE")
 
 
-def smooth_median(self, iterations=1, width=1):
-    """Smooth a surface using a median filter.
+def _smooth(
+    self: RegularSurface,
+    window_function: Callable[[np.ndarray], np.ndarray],
+    iterations: int = 1,
+) -> None:
+    """
+    Smooth a RegularSurface using a window function.
 
-    .. versionadded:: 2.1
+    Original mask (undefined values) is stored before applying
+    smoothing on a filled array. Subsequently the original mask
+    is used to restore the undefined values in the output.
     """
 
     mask = ma.getmaskarray(self.values)
-    tmpv = ma.filled(self.values, fill_value=np.nan)
 
-    for _itr in range(iterations):
-        tmpv = scipy.ndimage.median_filter(tmpv, width)
+    self.fill()
 
-    tmpv = ma.masked_invalid(tmpv)
+    smoothed_values = self.values
+    for _ in range(iterations):
+        smoothed_values = window_function(smoothed_values, mode="nearest")
 
-    # seems that false areas of invalids (masked) may be made; combat that:
-    self.values = tmpv
-    self.fill()
-    self.values = ma.array(self.values, mask=mask)
+    self.values = ma.array(smoothed_values, mask=mask)
diff --git a/src/xtgeo/surface/regular_surface.py b/src/xtgeo/surface/regular_surface.py
index 6a2e5142c..b86a11653 100644
--- a/src/xtgeo/surface/regular_surface.py
+++ b/src/xtgeo/surface/regular_surface.py
@@ -49,6 +49,7 @@
 import numpy as np
 import numpy.ma as ma
 import pandas as pd
+import scipy.ndimage
 
 import xtgeo
 import xtgeo.common.sys as xtgeosys
@@ -2065,18 +2066,41 @@ def fill(self, fill_value=None):
         """
         _regsurf_gridding.surf_fill(self, fill_value=fill_value)
 
-    def smooth(self, method="median", iterations=1, width=1):
+    def smooth(
+        self,
+        method: Literal["median", "gaussian"] = "median",
+        iterations: int = 1,
+        width: float = 1,
+    ) -> None:
         """Various smoothing methods for surfaces.
 
         Args:
-            method: Smoothing method (median)
+            method: Smoothing method (median or gaussian)
             iterations: Number of iterations
-            width: Range of influence (in nodes)
+            width:
+                - If method is 'median' range of influence is in nodes.
+                - If method is 'gaussian' range of influence is standard
+                  deviation of the Gaussian kernel.
 
         .. versionadded:: 2.1
         """
+
         if method == "median":
-            _regsurf_gridding.smooth_median(self, iterations=iterations, width=width)
+            _regsurf_gridding._smooth(
+                self,
+                window_function=functools.partial(
+                    scipy.ndimage.median_filter, size=int(width)
+                ),
+                iterations=iterations,
+            )
+        elif method == "gaussian":
+            _regsurf_gridding._smooth(
+                self,
+                window_function=functools.partial(
+                    scipy.ndimage.gaussian_filter, sigma=width
+                ),
+                iterations=iterations,
+            )
         else:
             raise ValueError("Unsupported method for smoothing")
 
diff --git a/tests/test_surface/test_regular_surface.py b/tests/test_surface/test_regular_surface.py
index 09f887552..b8fa70397 100644
--- a/tests/test_surface/test_regular_surface.py
+++ b/tests/test_surface/test_regular_surface.py
@@ -1148,20 +1148,35 @@ def test_fill():
 
 
 def test_smoothing():
-    """Smooth the the surface"""
+    """Smooth the surface using median/gaussian filter"""
 
     srf = xtgeo.surface_from_file(TESTSET1, fformat="irap_binary")
 
     mean1 = srf.values.mean()
     assert mean1 == pytest.approx(1698.65, abs=0.1)
 
-    srf.smooth(iterations=1, width=5)
-
-    mean2 = srf.values.mean()
+    # test median smoothing methoid
+    srf2 = srf.copy()
+    srf2.smooth(method="median", iterations=1, width=5)
+    mean2 = srf2.values.mean()
     assert mean2 == pytest.approx(1698.65, abs=0.3)  # smoothed ~same mean
-
     assert mean1 != mean2  # but not exacly same
 
+    # test gaussian smoothing methoid
+    srf3 = srf.copy()
+    srf3.smooth(method="gaussian", iterations=1, width=5)
+    mean3 = srf3.values.mean()
+    assert mean3 == pytest.approx(1698.65, abs=0.3)  # smoothed ~same mean
+    assert mean1 != mean3  # but not exacly same
+
+    # check that the three surfaces have different min values
+    min1 = srf.values.min()
+    min2 = srf2.values.min()
+    min3 = srf3.values.min()
+    assert min1 == pytest.approx(1547.20, abs=0.1)
+    assert min2 == pytest.approx(1553.62, abs=0.1)
+    assert min3 == pytest.approx(1566.91, abs=0.1)
+
 
 def test_loadvalues_before_remove_deprecated(default_surface):
     """Test that load_values() has the claimed effect before deprecating __init__"""