diff --git a/pygmt/src/grdgradient.py b/pygmt/src/grdgradient.py
index 2ee0c62c7b5..cf4222f58b5 100644
--- a/pygmt/src/grdgradient.py
+++ b/pygmt/src/grdgradient.py
@@ -15,20 +15,74 @@
 __doctest_skip__ = ["grdgradient"]
 
 
+def _alias_option_N(  # noqa: N802
+    normalize=False,
+    norm_amp=None,
+    norm_ambient=None,
+    norm_sigma=None,
+    norm_offset=None,
+):
+    """
+    Helper function to create the alias for the -N option.
+
+    Examples
+    --------
+    >>> def parse(**kwargs):
+    ...     return AliasSystem(N=_alias_option_N(**kwargs)).get("N")
+    >>> parse(normalize=True)
+    ''
+    >>> parse(normalize="laplace")
+    'e'
+    >>> parse(normalize="cauchy")
+    't'
+    >>> parse(
+    ...     normalize="laplace",
+    ...     norm_amp=2,
+    ...     norm_offset=10,
+    ...     norm_sigma=0.5,
+    ...     norm_ambient=0.1,
+    ... )
+    'e2+a0.1+s0.5+o10'
+    >>> # Check for backward compatibility with old syntax
+    >>> parse(normalize="e2+a0.2+s0.5+o10")
+    'e2+a0.2+s0.5+o10'
+    """
+    _normalize_mapping = {"laplace": "e", "cauchy": "t"}
+    # Check for old syntax for normalize
+    if isinstance(normalize, str) and normalize not in _normalize_mapping:
+        if any(
+            v is not None and v is not False
+            for v in [norm_amp, norm_ambient, norm_sigma, norm_offset]
+        ):
+            msg = (
+                "Parameter 'normalize' using old syntax with raw GMT CLI string "
+                "conflicts with parameters 'norm_amp', 'norm_ambient', 'norm_sigma', "
+                "or 'norm_offset'."
+            )
+            raise GMTInvalidInput(msg)
+        _normalize_mapping = None
+
+    return [
+        Alias(normalize, name="normalize", mapping=_normalize_mapping),
+        Alias(norm_amp, name="norm_amp"),
+        Alias(norm_ambient, name="norm_ambient", prefix="+a"),
+        Alias(norm_sigma, name="norm_sigma", prefix="+s"),
+        Alias(norm_offset, name="norm_offset", prefix="+o"),
+    ]
+
+
 @fmt_docstring
-@use_alias(
-    D="direction",
-    N="normalize",
-    Q="tiles",
-    S="slope_file",
-    f="coltypes",
-    n="interpolation",
-)
-def grdgradient(
+@use_alias(D="direction", Q="tiles", S="slope_file", f="coltypes", n="interpolation")
+def grdgradient(  # noqa: PLR0913
     grid: PathLike | xr.DataArray,
     outgrid: PathLike | None = None,
     azimuth: float | Sequence[float] | None = None,
     radiance: Sequence[float] | str | None = None,
+    normalize: Literal["laplace", "cachy"] | bool = False,
+    norm_amp: float | None = None,
+    norm_ambient: float | None = None,
+    norm_sigma: float | None = None,
+    norm_offset: float | None = None,
     region: Sequence[float | str] | str | None = None,
     verbose: Literal["quiet", "error", "warning", "timing", "info", "compat", "debug"]
     | bool = False,
@@ -44,6 +98,8 @@ def grdgradient(
 
     $aliases
        - A = azimuth
+       - N = normalize, norm_amp, **+a**: norm_ambient, **+s**: norm_sigma,
+         **+o**: norm_offset
        - E = radiance
        - R = region
        - V = verbose
@@ -100,31 +156,37 @@ def grdgradient(
         algorithm; in this case *azim* and *elev* are hardwired to 315
         and 45 degrees. This means that even if you provide other values
         they will be ignored.).
-    normalize : str or bool
-        [**e**\|\ **t**][*amp*][**+a**\ *ambient*][**+s**\ *sigma*]\
-        [**+o**\ *offset*].
-        The actual gradients :math:`g` are offset and scaled to produce
-        normalized gradients :math:`g_n` with a maximum output magnitude of
-        *amp*. If *amp* is not given, default *amp* = 1. If *offset* is not
-        given, it is set to the average of :math:`g`. The following forms are
-        supported:
-
-        - **True**: Normalize using math:`g_n = \mbox{amp}\
-          (\frac{g - \mbox{offset}}{max(|g - \mbox{offset}|)})`
-        - **e**: Normalize using a cumulative Laplace distribution yielding:
-          :math:`g_n = \mbox{amp}(1 - \
-          \exp{(\sqrt{2}\frac{g - \mbox{offset}}{\sigma}))}`, where
-          :math:`\sigma` is estimated using the L1 norm of
-          :math:`(g - \mbox{offset})` if it is not given.
-        - **t**: Normalize using a cumulative Cauchy distribution yielding:
-          :math:`g_n = \
-          \frac{2(\mbox{amp})}{\pi}(\tan^{-1}(\frac{g - \
-          \mbox{offset}}{\sigma}))` where :math:`\sigma` is estimated
-          using the L2 norm of :math:`(g - \mbox{offset})` if it is not
-          given.
-
-        As a final option, you may add **+a**\ *ambient* to add *ambient* to
-        all nodes after gradient calculations are completed.
+    normalize
+        Normalize the output gradients. Valid values are:
+
+        - ``False``: No normalization is done [Default].
+        - ``True``: Normalize using max absolute value.
+        - ``"laplace"``: Normalize using cumulative Laplace distribution.
+        - ``"cauchy"``: Normalize using cumulative Cauchy distribution.
+
+        The normalization process is controlled via the additional parameters
+        ``norm_amp``, ``norm_ambient``, ``norm_sigma``, and ``norm_offset``.
+
+        Let :math:`g` denote the actual gradients, :math:`g_n` the normalized gradients,
+        :math:`a` the maximum output magnitude (``norm_amp``), :math:`o` the offset
+        value (``norm_offset``), and :math:`\sigma` the sigma value (``norm_sigma``).
+        The normalization is computed as follows:
+
+        - ``True``: :math:`g_n = a (\frac{g - o}{\max(|g - o|)})`
+        - ``"laplace"``: :math:`g_n = a(1 - \exp(\sqrt{2}\frac{g - o}{\sigma}))`
+        - ``"cauchy"``: :math:`g_n = \frac{2a}{\pi}\arctan(\frac{g - o}{\sigma})`
+    norm_amp
+        Set the maximum output magnitude [Default is 1].
+    norm_ambient
+        The ambient value to add to all nodes after gradient calculations are completed
+        [Default is 0].
+    norm_offset
+        The offset value used in the normalization. If not given, it is set to the
+        average of :math:`g`.
+    norm_sigma
+        The *sigma* value used in the Laplace or Cauchy normalization. If not given,
+        it is estimated from the L1 norm of :math:`g-o` for Laplace or the L2 norm of
+        :math:`g-o` for Cauchy.
     tiles : str
         **c**\|\ **r**\|\ **R**.
         Control how normalization via ``normalize`` is carried out. When
@@ -181,6 +243,13 @@ def grdgradient(
     aliasdict = AliasSystem(
         A=Alias(azimuth, name="azimuth", sep="/", size=2),
         E=Alias(radiance, name="radiance", sep="/", size=2),
+        N=_alias_option_N(
+            normalize,
+            norm_amp,
+            norm_ambient,
+            norm_sigma,
+            norm_offset,
+        ),
     ).add_common(
         R=region,
         V=verbose,
diff --git a/pygmt/tests/test_grdgradient.py b/pygmt/tests/test_grdgradient.py
index 5749b237aea..49fee09b2f2 100644
--- a/pygmt/tests/test_grdgradient.py
+++ b/pygmt/tests/test_grdgradient.py
@@ -73,6 +73,22 @@ def test_grdgradient_no_outgrid(grid, expected_grid):
     xr.testing.assert_allclose(a=result, b=expected_grid)
 
 
+def test_grdgradient_normalize_mixed_syntax(grid):
+    """
+    Test that mixed syntax for normalize in grdgradient raises GMTInvalidInput.
+    """
+    with pytest.raises(GMTInvalidInput):
+        grdgradient(grid=grid, normalize="t", norm_amp=2)
+    with pytest.raises(GMTInvalidInput):
+        grdgradient(grid=grid, normalize="t1", norm_amp=2)
+    with pytest.raises(GMTInvalidInput):
+        grdgradient(grid=grid, normalize="t1", norm_sigma=2)
+    with pytest.raises(GMTInvalidInput):
+        grdgradient(grid=grid, normalize="e", norm_offset=5)
+    with pytest.raises(GMTInvalidInput):
+        grdgradient(grid=grid, normalize="e", norm_ambient=0.1)
+
+
 def test_grdgradient_fails(grid):
     """
     Check that grdgradient fails correctly.