diff --git a/examples/optical_flow_methods_convergence.py b/examples/optical_flow_methods_convergence.py
index c5948ac4..9208af59 100644
--- a/examples/optical_flow_methods_convergence.py
+++ b/examples/optical_flow_methods_convergence.py
@@ -374,4 +374,22 @@ def plot_optflow_method_convergence(input_precip, optflow_method_name, motion_ty
 # ~~~~~~~~~~~~~~~~~
 plot_optflow_method_convergence(reference_field, "DARTS", "rotor")
 
+################################################################################
+# Farneback
+# ---------
+#
+# Constant motion x-direction
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~
+plot_optflow_method_convergence(reference_field, "farneback", "linear_x")
+
+################################################################################
+# Constant motion y-direction
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~
+plot_optflow_method_convergence(reference_field, "farneback", "linear_y")
+
+################################################################################
+# Rotational motion
+# ~~~~~~~~~~~~~~~~~
+plot_optflow_method_convergence(reference_field, "farneback", "rotor")
+
 # sphinx_gallery_thumbnail_number = 5
diff --git a/examples/plot_optical_flow.py b/examples/plot_optical_flow.py
index 2ed26d86..cf196db0 100644
--- a/examples/plot_optical_flow.py
+++ b/examples/plot_optical_flow.py
@@ -130,7 +130,7 @@
 #
 # This module implements the anisotropic diffusion method presented in Proesmans
 # et al. (1994), a robust optical flow technique which employs the notion of
-# inconsitency during the solution of the optical flow equations.
+# inconsistency during the solution of the optical flow equations.
 
 oflow_method = motion.get_method("proesmans")
 R[~np.isfinite(R)] = metadata["zerovalue"]
@@ -141,4 +141,22 @@
 quiver(V4, geodata=metadata, step=25)
 plt.show()
 
+################################################################################
+# Farnebäck smoothed method
+# -------------------------
+#
+# This module implements the pyramidal decomposition method for motion estimation
+# of Farnebäck as implemented in OpenCV, with an option for smoothing and
+# renormalization of the motion fields proposed by Driedger et al.:
+# https://cmosarchives.ca/Congress_P_A/program_abstracts2022.pdf (p. 392).
+
+oflow_method = motion.get_method("farneback")
+R[~np.isfinite(R)] = metadata["zerovalue"]
+V5 = oflow_method(R[-2:, :, :], verbose=True)
+
+# Plot the motion field
+plot_precip_field(R_, geodata=metadata, title="Farneback")
+quiver(V5, geodata=metadata, step=25)
+plt.show()
+
 # sphinx_gallery_thumbnail_number = 1
diff --git a/pysteps/motion/farneback.py b/pysteps/motion/farneback.py
new file mode 100644
index 00000000..79eebeed
--- /dev/null
+++ b/pysteps/motion/farneback.py
@@ -0,0 +1,269 @@
+# -*- coding: utf-8 -*-
+"""
+pysteps.motion.farneback
+========================
+
+The Farneback dense optical flow module.
+
+This module implements the interface to the local `Farneback`_ routine
+available in OpenCV_.
+
+.. _OpenCV: https://opencv.org/
+
+.. _`Farneback`:\
+    https://docs.opencv.org/3.4/dc/d6b/group__video__track.html#ga5d10ebbd59fe09c5f650289ec0ece5af
+
+.. autosummary::
+    :toctree: ../generated/
+
+    farneback
+"""
+
+import numpy as np
+from numpy.ma.core import MaskedArray
+import scipy.ndimage as sndi
+import time
+
+from pysteps.decorators import check_input_frames
+from pysteps.exceptions import MissingOptionalDependency
+from pysteps.utils.images import morph_opening
+
+try:
+    import cv2
+
+    CV2_IMPORTED = True
+except ImportError:
+    CV2_IMPORTED = False
+
+
+@check_input_frames(2)
+def farneback(
+    input_images,
+    pyr_scale=0.5,
+    levels=3,
+    winsize=15,
+    iterations=3,
+    poly_n=5,
+    poly_sigma=1.1,
+    flags=0,
+    size_opening=3,
+    sigma=60.0,
+    verbose=False,
+):
+    """Estimate a dense motion field from a sequence of 2D images using
+    the `Farneback`_ optical flow algorithm.
+
+    This function computes dense optical flow between each pair of consecutive
+    input frames using OpenCV's Farneback method. If more than two frames are
+    provided, the motion fields estimated from all consecutive pairs are
+    averaged to obtain a single representative advection field.
+
+    After the pairwise motion fields are averaged, the resulting motion field
+    can optionally be smoothed with a Gaussian filter. In that case, its
+    amplitude is rescaled so that the mean motion magnitude is preserved.
+
+    .. _OpenCV: https://opencv.org/
+
+    .. _`Farneback`:\
+        https://docs.opencv.org/3.4/dc/d6b/group__video__track.html#ga5d10ebbd59fe09c5f650289ec0ece5af
+
+    .. _MaskedArray:\
+        https://docs.scipy.org/doc/numpy/reference/maskedarray.baseclass.html#numpy.ma.MaskedArray
+
+    .. _ndarray:\
+    https://docs.scipy.org/doc/numpy/reference/generated/numpy.ndarray.html
+
+    Parameters
+    ----------
+    input_images: ndarray_ or MaskedArray_
+        Array of shape (T, m, n) containing a sequence of *T* two-dimensional
+        input images of shape (m, n). The indexing order in **input_images** is
+        assumed to be (time, latitude, longitude).
+
+        *T* = 2 is the minimum required number of images.
+        With *T* > 2, all the resulting motion vectors are averaged together.
+
+        In case of ndarray_, invalid values (Nans or infs) are masked,
+        otherwise the mask of the MaskedArray_ is used. Such mask defines a
+        region where features are not detected for the tracking algorithm.
+
+    pyr_scale : float, optional
+        Parameter specifying the image scale (<1) used to build pyramids for
+        each image; pyr_scale=0.5 means a classical pyramid, where each next
+        layer is twice smaller than the previous one. This and the following
+        parameter descriptions are adapted from the original OpenCV
+        documentation (see https://docs.opencv.org).
+
+    levels : int, optional
+        Number of pyramid layers including the initial image; levels=1 means
+        that no extra layers are created and only the original images are used.
+
+    winsize : int, optional
+        Averaging window size; larger values increase the algorithm robustness
+        to image noise and give more stable motion estimates. Small windows
+        (e.g. 10) lead to unrealistic motion.
+    iterations : int, optional
+        Number of iterations the algorithm does at each pyramid level.
+        
+    poly_n : int
+        Size of the pixel neighborhood used to find polynomial expansion in
+        each pixel; larger values mean that the image will be approximated with
+        smoother surfaces, yielding more robust algorithm and more blurred
+        motion field, typically poly_n = 5 or 7.
+        
+    poly_sigma : float
+        Standard deviation of the Gaussian that is used to smooth derivatives
+        used as a basis for the polynomial expansion; for poly_n=5, you can set
+        poly_sigma=1.1, for poly_n=7, a good value would be poly_sigma=1.5.
+        
+    flags : int, optional
+        Operation flags that can be a combination of the following:
+
+        OPTFLOW_USE_INITIAL_FLOW uses the input 'flow' as an initial flow
+        approximation.
+
+        OPTFLOW_FARNEBACK_GAUSSIAN uses the Gaussian winsize x winsize filter
+        instead of a box filter of the same size for optical flow estimation;
+        usually, this option gives a more accurate flow than with a box filter,
+        at the cost of lower speed; normally, winsize for a Gaussian window
+        should be set to a larger value to achieve the same level of robustness.
+        
+    size_opening : int, optional
+        Non-OpenCV parameter:
+        The structuring element size for the filtering of isolated pixels [px].
+
+    sigma : float, optional
+        Non-OpenCV parameter:
+        The smoothing bandwidth of the motion field. The motion field amplitude
+        is adjusted by multiplying by the ratio of average magnitude before and
+        after smoothing to avoid damping of the motion field.
+
+    verbose: bool, optional
+        If set to True, print some information about the program.
+
+    Returns
+    -------
+    out : ndarray_, shape (2,m,n)
+        Return the advection field having shape
+        (2, m, n), where out[0, :, :] contains the x-components of the motion
+        vectors and out[1, :, :] contains the y-components.
+        The velocities are in units of pixels / timestep, where timestep is the
+        time difference between the two input images.
+        Return a zero motion field of shape (2, m, n) when no motion is
+        detected.
+        
+    References
+    ----------
+    Farnebäck, G.: Two-frame motion estimation based on polynomial expansion, 
+    In Image Analysis, pages 363–370. Springer, 2003.
+    Driedger, N., Mahidjiba, A. and Hortal, A.P. (2022, June 1-8): Evaluation of optical flow
+    methods for radar precipitation extrapolation.
+    Canadian Meteorological and Oceanographic Society Congress, contributed abstract 11801.
+    """
+
+    if len(input_images.shape) != 3:
+        raise ValueError(
+            "input_images has %i dimensions, but a "
+            "three-dimensional array is expected" % len(input_images.shape)
+        )
+
+    input_images = input_images.copy()
+
+    if verbose:
+        print("Computing the motion field with the Farneback method.")
+        t0 = time.time()
+
+    if not CV2_IMPORTED:
+        raise MissingOptionalDependency(
+            "OpenCV (cv2) is required for the Farneback optical flow method, but it is not installed"
+        )
+
+    nr_pairs = input_images.shape[0] - 1
+    domain_size = (input_images.shape[1], input_images.shape[2])
+    u_sum = np.zeros(domain_size)
+    v_sum = np.zeros(domain_size)
+    for n in range(nr_pairs):
+        # extract consecutive images
+        prvs_img = input_images[n, :, :].copy()
+        next_img = input_images[n + 1, :, :].copy()
+
+        # Check if a MaskedArray is used. If not, mask the ndarray
+        if not isinstance(prvs_img, MaskedArray):
+            prvs_img = np.ma.masked_invalid(prvs_img)
+        np.ma.set_fill_value(prvs_img, prvs_img.min())
+
+        if not isinstance(next_img, MaskedArray):
+            next_img = np.ma.masked_invalid(next_img)
+        np.ma.set_fill_value(next_img, next_img.min())
+
+        # scale between 0 and 255
+        im_min = prvs_img.min()
+        im_max = prvs_img.max()
+        if (im_max - im_min) > 1e-8:
+            prvs_img = (prvs_img.filled() - im_min) / (im_max - im_min) * 255
+        else:
+            prvs_img = prvs_img.filled() - im_min
+
+        im_min = next_img.min()
+        im_max = next_img.max()
+        if (im_max - im_min) > 1e-8:
+            next_img = (next_img.filled() - im_min) / (im_max - im_min) * 255
+        else:
+            next_img = next_img.filled() - im_min
+
+        # convert to 8-bit
+        prvs_img = np.ndarray.astype(prvs_img, "uint8")
+        next_img = np.ndarray.astype(next_img, "uint8")
+
+        # remove small noise with a morphological operator (opening)
+        if size_opening > 0:
+            prvs_img = morph_opening(prvs_img, prvs_img.min(), size_opening)
+            next_img = morph_opening(next_img, next_img.min(), size_opening)
+
+        flow = cv2.calcOpticalFlowFarneback(
+            prvs_img,
+            next_img,
+            None,
+            pyr_scale,
+            levels,
+            winsize,
+            iterations,
+            poly_n,
+            poly_sigma,
+            flags,
+        )
+
+        fa, fb = np.dsplit(flow, 2)
+        u_sum += fa.reshape(domain_size)
+        v_sum += fb.reshape(domain_size)
+
+    # Compute the average motion field
+    u = u_sum / nr_pairs
+    v = v_sum / nr_pairs
+
+    # Smoothing
+    if sigma > 0:
+        uv2 = u * u + v * v  # squared magnitude of motion field
+        us = sndi.gaussian_filter(u, sigma, mode="nearest")
+        vs = sndi.gaussian_filter(v, sigma, mode="nearest")
+        uvs2 = us * us + vs * vs  # squared magnitude of smoothed motion field
+
+        mean_uv2 = np.nanmean(uv2)
+        mean_uvs2 = np.nanmean(uvs2)
+        if mean_uvs2 > 0:
+            mult = np.sqrt(mean_uv2 / mean_uvs2)
+        else:
+            mult = 1.0
+    else:
+        mult = 1.0
+        us = u
+        vs = v
+    if verbose:
+        print("mult factor of smoothed motion field=", mult)
+
+    UV = np.stack([us * mult, vs * mult])
+
+    if verbose:
+        print("--- %s seconds ---" % (time.time() - t0))
+
+    return UV
diff --git a/pysteps/motion/interface.py b/pysteps/motion/interface.py
index 52742389..61a8c029 100644
--- a/pysteps/motion/interface.py
+++ b/pysteps/motion/interface.py
@@ -31,6 +31,7 @@
 from pysteps.motion.lucaskanade import dense_lucaskanade
 from pysteps.motion.proesmans import proesmans
 from pysteps.motion.vet import vet
+from pysteps.motion.farneback import farneback
 
 _methods = dict()
 _methods["constant"] = constant
@@ -39,6 +40,7 @@
 _methods["darts"] = DARTS
 _methods["proesmans"] = proesmans
 _methods["vet"] = vet
+_methods["farneback"] = farneback
 _methods[None] = lambda precip, *args, **kw: np.zeros(
     (2, precip.shape[1], precip.shape[2])
 )
@@ -73,6 +75,8 @@ def get_method(name):
     |                   | Laroche and Zawadzki (1995) and                      |
     |                   | Germann and Zawadzki (2002)                          |
     +-------------------+------------------------------------------------------+
+    |  farneback        | OpenCV implementation of the Farneback (2003) method.|
+    +-------------------+------------------------------------------------------+
 
     +--------------------------------------------------------------------------+
     | Methods implemented in C (these require separate compilation and linkage)|
diff --git a/pysteps/tests/test_motion.py b/pysteps/tests/test_motion.py
index 8d198960..ede3d92d 100644
--- a/pysteps/tests/test_motion.py
+++ b/pysteps/tests/test_motion.py
@@ -164,6 +164,8 @@ def _create_observations(input_precip, motion_type, num_times=9):
     (reference_field, "proesmans", "linear_y", 2, 0.45),
     (reference_field, "darts", "linear_x", 9, 20),
     (reference_field, "darts", "linear_y", 9, 20),
+    (reference_field, "farneback", "linear_x", 2, 28),
+    (reference_field, "farneback", "linear_y", 2, 28),
 ]
 
 
@@ -256,6 +258,7 @@ def test_optflow_method_convergence(
     ("vet", 3),
     ("darts", 9),
     ("proesmans", 2),
+    ("farneback", 2),
 ]
 
 
@@ -296,6 +299,7 @@ def test_no_precipitation(optflow_method_name, num_times):
     ("vet", 2, 3),
     ("darts", 9, 9),
     ("proesmans", 2, 2),
+    ("farneback", 2, np.inf),
 ]
 
 
@@ -303,7 +307,7 @@ def test_no_precipitation(optflow_method_name, num_times):
 def test_input_shape_checks(
     optflow_method_name, minimum_input_frames, maximum_input_frames
 ):
-    if optflow_method_name == "lk":
+    if optflow_method_name in ("lk", "farneback"):
         pytest.importorskip("cv2")
     image_size = 100
     motion_method = motion.get_method(optflow_method_name)
@@ -393,26 +397,34 @@ def test_vet_cost_function():
     assert (returned_values[0] - 1548250.87627097) < 0.001
 
 
-def test_lk_masked_array():
+@pytest.mark.parametrize(
+    "method,kwargs",
+    [
+        ("LK", {"fd_kwargs": {"buffer_mask": 20}, "verbose": False}),
+        ("farneback", {"verbose": False}),
+    ],
+)
+def test_motion_masked_array(method, kwargs):
     """
     Passing a ndarray with NaNs or a masked array should produce the same results.
+    Tests for both LK and Farneback motion estimation methods.
     """
     pytest.importorskip("cv2")
 
     __, precip_obs = _create_observations(
         reference_field.copy(), "linear_y", num_times=2
     )
-    motion_method = motion.get_method("LK")
+    motion_method = motion.get_method(method)
 
     # ndarray with nans
     np.ma.set_fill_value(precip_obs, -15)
     ndarray = precip_obs.filled()
     ndarray[ndarray == -15] = np.nan
-    uv_ndarray = motion_method(ndarray, fd_kwargs={"buffer_mask": 20}, verbose=False)
+    uv_ndarray = motion_method(ndarray, **kwargs)
 
     # masked array
     mdarray = np.ma.masked_invalid(ndarray)
     mdarray.data[mdarray.mask] = -15
-    uv_mdarray = motion_method(mdarray, fd_kwargs={"buffer_mask": 20}, verbose=False)
+    uv_mdarray = motion_method(mdarray, **kwargs)
 
     assert np.abs(uv_mdarray - uv_ndarray).max() < 0.01
diff --git a/pysteps/tests/test_motion_farneback.py b/pysteps/tests/test_motion_farneback.py
new file mode 100644
index 00000000..674b80c7
--- /dev/null
+++ b/pysteps/tests/test_motion_farneback.py
@@ -0,0 +1,135 @@
+import pytest
+import numpy as np
+
+from pysteps.motion import farneback
+from pysteps.exceptions import MissingOptionalDependency
+from pysteps.tests.helpers import get_precipitation_fields
+
+fb_arg_names = (
+    "pyr_scale",
+    "levels",
+    "winsize",
+    "iterations",
+    "poly_n",
+    "poly_sigma",
+    "flags",
+    "size_opening",
+    "sigma",
+    "verbose",
+)
+
+fb_arg_values = [
+    (0.5, 3, 15, 3, 5, 1.1, 0, 3, 60.0, False),  # default
+    (0.5, 1, 5, 1, 7, 1.5, 0, 0, 0.0, True),  # minimal settings, sigma=0, verbose
+    (
+        0.3,
+        5,
+        30,
+        10,
+        7,
+        1.5,
+        1,
+        5,
+        10.0,
+        False,
+    ),  # maximal settings, flags=1, big opening
+    (0.5, 3, 15, 3, 5, 1.1, 0, 0, 60.0, True),  # no opening, verbose
+]
+
+
+@pytest.mark.parametrize(fb_arg_names, fb_arg_values)
+def test_farneback_params(
+    pyr_scale,
+    levels,
+    winsize,
+    iterations,
+    poly_n,
+    poly_sigma,
+    flags,
+    size_opening,
+    sigma,
+    verbose,
+):
+    """Test Farneback with various parameters and input types."""
+    pytest.importorskip("cv2")
+    # Input: realistic precipitation fields
+    precip, metadata = get_precipitation_fields(
+        num_prev_files=2,
+        num_next_files=0,
+        return_raw=False,
+        metadata=True,
+        upscale=2000,
+    )
+    precip = precip.filled()
+
+    output = farneback.farneback(
+        precip,
+        pyr_scale=pyr_scale,
+        levels=levels,
+        winsize=winsize,
+        iterations=iterations,
+        poly_n=poly_n,
+        poly_sigma=poly_sigma,
+        flags=flags,
+        size_opening=size_opening,
+        sigma=sigma,
+        verbose=verbose,
+    )
+
+    assert isinstance(output, np.ndarray)
+    assert output.shape[0] == 2
+    assert output.shape[1:] == precip[0].shape
+    assert np.isfinite(output).all() or np.isnan(output).any()
+
+
+def test_farneback_invalid_shape():
+    """Test error when input is wrong shape."""
+    pytest.importorskip("cv2")
+    arr = np.random.rand(64, 64)
+    with pytest.raises(ValueError):
+        farneback.farneback(arr)
+
+
+def test_farneback_nan_input():
+    """Test NaN handling in input."""
+    pytest.importorskip("cv2")
+    arr = np.random.rand(2, 64, 64)
+    arr[0, 0, 0] = np.nan
+    arr[1, 10, 10] = np.inf
+    result = farneback.farneback(arr)
+    assert result.shape == (2, 64, 64)
+
+
+def test_farneback_cv2_missing(monkeypatch):
+    """Test MissingOptionalDependency when cv2 is not injected."""
+    monkeypatch.setattr(farneback, "CV2_IMPORTED", False)
+    arr = np.random.rand(2, 64, 64)
+    with pytest.raises(MissingOptionalDependency):
+        farneback.farneback(arr)
+    monkeypatch.setattr(farneback, "CV2_IMPORTED", True)  # restore
+
+
+def test_farneback_sigma_zero():
+    """Test sigma=0 disables smoothing logic."""
+    pytest.importorskip("cv2")
+    arr = np.random.rand(2, 32, 32)
+    result = farneback.farneback(arr, sigma=0.0)
+    assert isinstance(result, np.ndarray)
+    assert result.shape == (2, 32, 32)
+
+
+def test_farneback_small_window():
+    """Test winsize edge case behavior."""
+    pytest.importorskip("cv2")
+    arr = np.random.rand(2, 16, 16)
+    result = farneback.farneback(arr, winsize=3)
+    assert result.shape == (2, 16, 16)
+
+
+def test_farneback_verbose(capsys):
+    """Test that verbose produces output (side-effect only)."""
+    pytest.importorskip("cv2")
+    arr = np.random.rand(2, 16, 16)
+    farneback.farneback(arr, verbose=True)
+    out = capsys.readouterr().out
+    assert "Farneback method" in out or "mult factor" in out or "---" in out