diff --git a/doc/api/index.rst b/doc/api/index.rst
index 7e1441026f2..ebed243de36 100644
--- a/doc/api/index.rst
+++ b/doc/api/index.rst
@@ -104,6 +104,7 @@ Operations on grids:
     grdproject
     grdsample
     grdtrack
+    grdvolume
 
 Crossover analysis with x2sys:
 
diff --git a/pygmt/__init__.py b/pygmt/__init__.py
index 5130aa4be0a..0edc03d8659 100644
--- a/pygmt/__init__.py
+++ b/pygmt/__init__.py
@@ -46,6 +46,7 @@
     grdproject,
     grdsample,
     grdtrack,
+    grdvolume,
     info,
     makecpt,
     nearneighbor,
diff --git a/pygmt/src/__init__.py b/pygmt/src/__init__.py
index e733b94badd..01b098a9aa1 100644
--- a/pygmt/src/__init__.py
+++ b/pygmt/src/__init__.py
@@ -24,6 +24,7 @@
 from pygmt.src.grdsample import grdsample
 from pygmt.src.grdtrack import grdtrack
 from pygmt.src.grdview import grdview
+from pygmt.src.grdvolume import grdvolume
 from pygmt.src.histogram import histogram
 from pygmt.src.image import image
 from pygmt.src.info import info
diff --git a/pygmt/src/grdvolume.py b/pygmt/src/grdvolume.py
new file mode 100644
index 00000000000..bd05bea29eb
--- /dev/null
+++ b/pygmt/src/grdvolume.py
@@ -0,0 +1,103 @@
+"""
+grdvolume - Calculate grid volume and area constrained by a contour.
+"""
+import pandas as pd
+from pygmt.clib import Session
+from pygmt.exceptions import GMTInvalidInput
+from pygmt.helpers import (
+    GMTTempFile,
+    build_arg_string,
+    fmt_docstring,
+    kwargs_to_strings,
+    use_alias,
+)
+
+
+@fmt_docstring
+@use_alias(
+    C="contour",
+    R="region",
+    S="unit",
+    V="verbose",
+)
+@kwargs_to_strings(C="sequence", R="sequence")
+def grdvolume(grid, output_type="pandas", outfile=None, **kwargs):
+    r"""
+    Determine the volume between the surface of a grid and a plane.
+
+    Read a 2-D grid file and calculate the volume contained below the surface
+    and above the plane specified by the given contour (or zero if not given)
+    and return the contour, area, volume, and maximum mean height
+    (volume/area). Alternatively, a range of contours can be specified to
+    return the volume and area inside the contour for all contour values.
+
+    Full option list at :gmt-docs:`grdvolume.html`
+
+    {aliases}
+
+    Parameters
+    ----------
+    grid : str or xarray.DataArray
+        The file name of the input grid or the grid loaded as a DataArray.
+    output_type : str
+        Determine the format the output data will be returned in [Default is
+        ``pandas``]:
+
+            - ``numpy`` - :class:`numpy.ndarray`
+            - ``pandas``-  :class:`pandas.DataFrame`
+            - ``file`` - ASCII file (requires ``outfile``)
+    outfile : str
+        The file name for the output ASCII file.
+    contour : str or int or float or list
+        *cval*\|\ *low/high/delta*\|\ **r**\ *low/high*\|\ **r**\ *cval*.
+        Find area, volume and mean height (volume/area) inside and above the
+        *cval* contour. Alternatively, search using all contours from *low* to
+        *high* in steps of *delta*. [Default returns area, volume and mean
+        height of the entire grid]. The area is measured in the plane of the
+        contour. Adding the **r** prefix computes the volume below the grid
+        surface and above the planes defined by *low* and *high*, or below
+        *cval* and  grid's minimum. Note that this is an *outside* volume
+        whilst the other forms compute an *inside* (below the surface) area
+        volume. Use this form to compute for example the volume of water
+        between two contours. If no *contour* is given then there is no contour
+        and the entire grid area, volume and the mean height is returned and
+        *cval* will be reported as 0.
+    {R}
+    {V}
+
+    Returns
+    -------
+    ret : pandas.DataFrame or numpy.ndarray or None
+        Return type depends on ``outfile`` and ``output_type``:
+
+        - None if ``outfile`` is set (output will be stored in file set by
+          ``outfile``)
+        - :class:`pandas.DataFrame` or :class:`numpy.ndarray` if ``outfile``
+          is not set (depends on ``output_type`` [Default is
+          class:`pandas.DataFrame`])
+    """
+    if output_type not in ["numpy", "pandas", "file"]:
+        raise GMTInvalidInput(
+            """Must specify format as either numpy, pandas, or file."""
+        )
+    if output_type == "file" and outfile is None:
+        raise GMTInvalidInput("""Must specify outfile for ASCII output.""")
+
+    with GMTTempFile() as tmpfile:
+        with Session() as lib:
+            file_context = lib.virtualfile_from_data(check_kind="raster", data=grid)
+            with file_context as infile:
+                if outfile is None:
+                    outfile = tmpfile.name
+                arg_str = " ".join([infile, build_arg_string(kwargs), "->" + outfile])
+                lib.call_module("grdvolume", arg_str)
+
+        # Read temporary csv output to a pandas table
+        if outfile == tmpfile.name:  # if user did not set outfile, return pd.DataFrame
+            result = pd.read_csv(tmpfile.name, sep="\t", header=None, comment=">")
+        elif outfile != tmpfile.name:  # return None if outfile set, output in outfile
+            result = None
+
+        if output_type == "numpy":
+            result = result.to_numpy()
+    return result
diff --git a/pygmt/tests/test_grdvolume.py b/pygmt/tests/test_grdvolume.py
new file mode 100644
index 00000000000..959f8d12f56
--- /dev/null
+++ b/pygmt/tests/test_grdvolume.py
@@ -0,0 +1,112 @@
+"""
+Tests for grdvolume.
+"""
+import os
+
+import numpy as np
+import numpy.testing as npt
+import pandas as pd
+import pytest
+from pygmt import grdvolume
+from pygmt.datasets import load_earth_relief
+from pygmt.exceptions import GMTInvalidInput
+from pygmt.helpers import GMTTempFile
+
+
+@pytest.fixture(scope="module", name="grid")
+def fixture_grid():
+    """
+    Load the grid data from the sample earth_relief file.
+    """
+    return load_earth_relief(resolution="01d", region=[-100, -95, 34, 39])
+
+
+@pytest.fixture(scope="module", name="data")
+def fixture_data():
+    """
+    Load the expected grdvolume data result as a numpy array.
+    """
+    data = np.array(
+        [
+            [
+                2.00000000e02,
+                1.59920815e11,
+                3.16386172e13,
+                1.97839269e02,
+            ],
+            [
+                2.50000000e02,
+                1.44365835e11,
+                2.38676788e13,
+                1.65327751e02,
+            ],
+            [
+                3.00000000e02,
+                1.23788259e11,
+                1.71278707e13,
+                1.38364259e02,
+            ],
+            [
+                3.50000000e02,
+                9.79597525e10,
+                1.15235913e13,
+                1.17635978e02,
+            ],
+            [
+                4.00000000e02,
+                7.26646663e10,
+                7.22303463e12,
+                9.94022955e01,
+            ],
+        ]
+    )
+    return data
+
+
+def test_grdvolume_format(grid):
+    """
+    Test that correct formats are returned.
+    """
+    grdvolume_default = grdvolume(grid=grid)
+    assert isinstance(grdvolume_default, pd.DataFrame)
+    grdvolume_array = grdvolume(grid=grid, output_type="numpy")
+    assert isinstance(grdvolume_array, np.ndarray)
+    grdvolume_df = grdvolume(grid=grid, output_type="pandas")
+    assert isinstance(grdvolume_df, pd.DataFrame)
+
+
+def test_grdvolume_invalid_format(grid):
+    """
+    Test that grdvolume fails with incorrect output_type argument.
+    """
+    with pytest.raises(GMTInvalidInput):
+        grdvolume(grid=grid, output_type=1)
+
+
+def test_grdvolume_no_outfile(grid):
+    """
+    Test that grdvolume fails when output_type set to 'file' but no outfile is
+    specified.
+    """
+    with pytest.raises(GMTInvalidInput):
+        grdvolume(grid=grid, output_type="file")
+
+
+def test_grdvolume_no_outgrid(grid, data):
+    """
+    Test the expected output of grdvolume with no output file set.
+    """
+    test_output = grdvolume(grid=grid, contour=[200, 400, 50], output_type="numpy")
+    npt.assert_allclose(test_output, data)
+
+
+def test_grdvolume_outgrid(grid):
+    """
+    Test the expected output of grdvolume with an output file set.
+    """
+    with GMTTempFile(suffix=".csv") as tmpfile:
+        result = grdvolume(
+            grid=grid, contour=[200, 400, 50], output_type="file", outfile=tmpfile.name
+        )
+        assert result is None  # return value is None
+        assert os.path.exists(path=tmpfile.name)  # check that outfile exists