diff --git a/pygmt/datasets/earth_relief.py b/pygmt/datasets/earth_relief.py
index f4cbdf54015..922467e784d 100644
--- a/pygmt/datasets/earth_relief.py
+++ b/pygmt/datasets/earth_relief.py
@@ -76,10 +76,17 @@ def load_earth_relief(
         - **igpp** : IGPP Global Earth Relief [Default option]. See
           :gmt-datasets:`earth-relief.html`.
 
+        - **synbath** : IGPP Global Earth Relief dataset that uses
+          stastical properties of young seafloor to provide more realistic
+          relief of young areas with small seamounts.
+
         - **gebco** : GEBCO Global Earth Relief with only observed relief and
           inferred relief via altimetric gravity. See
           :gmt-datasets:`earth-gebco.html`.
 
+        - **gebcosi** : GEBCO Global Earth Relief that gives sub-ice (si)
+          elevations.
+
     Returns
     -------
     grid : :class:`xarray.DataArray`
@@ -144,12 +151,13 @@ def load_earth_relief(
     earth_relief_sources = {
         "igpp": "earth_relief_",
         "gebco": "earth_gebco_",
+        "gebcosi": "earth_gebcosi_",
         "synbath": "earth_synbath_",
     }
     if data_source not in earth_relief_sources:
         raise GMTInvalidInput(
             f"Invalid earth relief 'data_source' {data_source}, "
-            "valid values are 'igpp', 'gebco', and 'synbath'."
+            "valid values are 'igpp', 'gebco', 'gebcosi' and 'synbath'."
         )
     if data_source != "igpp":
         with Session() as lib:
diff --git a/pygmt/helpers/testing.py b/pygmt/helpers/testing.py
index 1af025e0d8e..57c916f5f83 100644
--- a/pygmt/helpers/testing.py
+++ b/pygmt/helpers/testing.py
@@ -156,6 +156,8 @@ def download_test_data():
     datasets = [
         # Earth relief grids
         "@earth_gebco_01d_g",
+        "@earth_gebcosi_01d_g",
+        "@earth_gebcosi_15m_p",
         "@earth_relief_01d_p",
         "@earth_relief_01d_g",
         "@earth_relief_30m_p",
diff --git a/pygmt/tests/test_datasets_earth_relief.py b/pygmt/tests/test_datasets_earth_relief.py
index 358818d2f0d..23006e80ac9 100644
--- a/pygmt/tests/test_datasets_earth_relief.py
+++ b/pygmt/tests/test_datasets_earth_relief.py
@@ -8,7 +8,7 @@
 from pygmt.exceptions import GMTInvalidInput
 
 
-@pytest.mark.parametrize("data_source", ["igpp", "gebco", "synbath"])
+@pytest.mark.parametrize("data_source", ["igpp", "gebco", "gebcosi", "synbath"])
 def test_earth_relief_fails(data_source):
     """
     Make sure earth relief fails for invalid resolutions.
@@ -37,12 +37,14 @@ def test_earth_relief_01d_igpp_synbath(data_source):
     npt.assert_allclose(data.max(), 5559.0)
 
 
-def test_earth_relief_01d_gebco():
+@pytest.mark.parametrize("data_source", ["gebco", "gebcosi"])
+def test_earth_relief_01d_gebco(data_source):
     """
-    Test some properties of the earth relief 01d data with GEBCO data.
+    Test some properties of the earth relief 01d data with GEBCO and GEBOCSI
+    data.
     """
     data = load_earth_relief(
-        resolution="01d", registration="gridline", data_source="gebco"
+        resolution="01d", registration="gridline", data_source=data_source
     )
     assert data.shape == (181, 361)
     npt.assert_allclose(data.lat, np.arange(-90, 91, 1))
@@ -114,6 +116,23 @@ def test_earth_relief_05m_with_region():
     assert data.sizes["lon"] == 481
 
 
+def test_earth_gebcosi_15m_with_region():
+    """
+    Test loading a subregion of 15 arc-minute resolution earth_gebcosi grid.
+    """
+    data = load_earth_relief(
+        resolution="15m",
+        region=[85, 87, -88, -84],
+        registration="pixel",
+        data_source="gebcosi",
+    )
+    assert data.shape == (16, 8)
+    npt.assert_allclose(data.lat, np.arange(-87.875, -84, 0.25))
+    npt.assert_allclose(data.lon, np.arange(85.125, 87, 0.25))
+    npt.assert_allclose(data.min(), -531)
+    npt.assert_allclose(data.max(), 474)
+
+
 def test_earth_relief_30s_synbath():
     """
     Test some properties of the earth relief 30s data with SYNBATH data.