diff --git a/.pylintrc b/.pylintrc
index f3194b081..6e8b54072 100644
--- a/.pylintrc
+++ b/.pylintrc
@@ -50,7 +50,7 @@ confidence=
 # --enable=similarities". If you want to run only the classes checker, but have
 # no Warning level messages displayed, use"--disable=all --enable=classes
 # --disable=W"
-disable=print-statement,parameter-unpacking,unpacking-in-except,old-raise-syntax,backtick,long-suffix,old-ne-operator,old-octal-literal,raw-checker-failed,bad-inline-option,locally-disabled,locally-enabled,file-ignored,suppressed-message,deprecated-pragma,apply-builtin,basestring-builtin,buffer-builtin,cmp-builtin,coerce-builtin,execfile-builtin,file-builtin,long-builtin,raw_input-builtin,reduce-builtin,standarderror-builtin,unicode-builtin,xrange-builtin,coerce-method,delslice-method,getslice-method,setslice-method,no-absolute-import,old-division,dict-iter-method,dict-view-method,next-method-called,metaclass-assignment,indexing-exception,raising-string,reload-builtin,oct-method,hex-method,nonzero-method,cmp-method,input-builtin,round-builtin,intern-builtin,unichr-builtin,map-builtin-not-iterating,zip-builtin-not-iterating,range-builtin-not-iterating,filter-builtin-not-iterating,using-cmp-argument,eq-without-hash,div-method,idiv-method,rdiv-method,exception-message-attribute,invalid-str-codec,sys-max-int,bad-python3-import,deprecated-string-function,deprecated-str-translate-call,attribute-defined-outside-init,similarities,bad-continuation,import-error,assignment-from-no-return,unsubscriptable-object
+disable=print-statement,parameter-unpacking,unpacking-in-except,old-raise-syntax,backtick,long-suffix,old-ne-operator,old-octal-literal,raw-checker-failed,bad-inline-option,locally-disabled,locally-enabled,file-ignored,suppressed-message,deprecated-pragma,apply-builtin,basestring-builtin,buffer-builtin,cmp-builtin,coerce-builtin,execfile-builtin,file-builtin,long-builtin,raw_input-builtin,reduce-builtin,standarderror-builtin,unicode-builtin,xrange-builtin,coerce-method,delslice-method,getslice-method,setslice-method,no-absolute-import,old-division,dict-iter-method,dict-view-method,next-method-called,metaclass-assignment,indexing-exception,raising-string,reload-builtin,oct-method,hex-method,nonzero-method,cmp-method,input-builtin,round-builtin,intern-builtin,unichr-builtin,map-builtin-not-iterating,zip-builtin-not-iterating,range-builtin-not-iterating,filter-builtin-not-iterating,using-cmp-argument,eq-without-hash,div-method,idiv-method,rdiv-method,exception-message-attribute,invalid-str-codec,sys-max-int,bad-python3-import,deprecated-string-function,deprecated-str-translate-call,attribute-defined-outside-init,similarities,bad-continuation,import-error,assignment-from-no-return,unsubscriptable-object,abstract-method
 
 # Enable the message, report, category or checker with the given id(s). You can
 # either give multiple identifier separated by comma (,) or put this option
diff --git a/README.rst b/README.rst
index 4b252bef1..c16c3a2de 100644
--- a/README.rst
+++ b/README.rst
@@ -42,7 +42,7 @@ About
 
 *Harmonica* is a Python library for processing and modeling gravity and magnetic data.
 It includes common processing steps, like calculation of Bouguer and terrain
-corrections, reduction to the pole, upward continuation, equivalent layers, and more.
+corrections, reduction to the pole, upward continuation, equivalent sources, and more.
 There are forward modeling functions for basic geometric shapes, like spheres, prisms,
 polygonal prisms, and tesseroids. The inversion methods are implemented as classes with
 an interface inspired by scikit-learn (like `Verde <https://www.fatiando.org/verde>`__).
diff --git a/examples/eql/README.txt b/examples/eql/README.txt
deleted file mode 100644
index 68c396b7e..000000000
--- a/examples/eql/README.txt
+++ /dev/null
@@ -1,2 +0,0 @@
-Equivalent Layers
------------------
diff --git a/examples/equivalent_sources/README.txt b/examples/equivalent_sources/README.txt
new file mode 100644
index 000000000..94d400d30
--- /dev/null
+++ b/examples/equivalent_sources/README.txt
@@ -0,0 +1,2 @@
+Equivalent Sources
+------------------
diff --git a/examples/eql/harmonic.py b/examples/equivalent_sources/cartesian.py
similarity index 74%
rename from examples/eql/harmonic.py
rename to examples/equivalent_sources/cartesian.py
index 7fafc885a..eb0d2d8b8 100644
--- a/examples/eql/harmonic.py
+++ b/examples/equivalent_sources/cartesian.py
@@ -14,16 +14,17 @@
 Upward-continuation is also a routine task for smoothing, noise attenuation,
 source separation, etc.
 
-Both tasks can be done simultaneously through an *equivalent layer*
-[Dampney1969]_. We will use :class:`harmonica.EQLHarmonic` to estimate the
-coefficients of a set of point sources (the equivalent layer) that fit the
-observed data. The fitted layer can then be used to predict data values
-wherever we want, like on a grid at a certain altitude. The sources for
-:class:`~harmonica.EQLHarmonic` in particular are placed one beneath each data
-point at a relative depth from the elevation of the data point
-following [Cooper2000]_.
-
-The advantage of using an equivalent layer is that it takes into account the 3D
+Both tasks can be done simultaneously through an *equivalent sources*
+[Dampney1969]_ (a.k.a *equivalent layer*). We will use
+:class:`harmonica.EquivalentSources` to estimate the coefficients of a set of
+point sources that fit the observed data. The fitted sources can then be used
+to predict data values wherever we want, like on a grid at a certain altitude.
+By default, the sources for :class:`~harmonica.EquivalentSources` are placed
+one beneath each data point at a relative depth from the elevation of the data
+point following [Cooper2000]_. This behaviour can be changed throught the
+`depth_type` optional argument.
+
+The advantage of using equivalent sources is that it takes into account the 3D
 nature of the observations, not just their horizontal positions. It also allows
 data uncertainty to be taken into account and noise to be suppressed though the
 least-squares fitting process. The main disadvantage is the increased
@@ -53,24 +54,25 @@
 easting, northing = projection(data.longitude.values, data.latitude.values)
 coordinates = (easting, northing, data.altitude_m)
 
-# Create the equivalent layer. We'll use the default point source configuration
-# at a relative depth beneath each observation point.
+# Create the equivalent sources.
+# We'll use the default point source configuration at a relative depth beneath
+# each observation point.
 # The damping parameter helps smooth the predicted data and ensure stability.
-eql = hm.EQLHarmonic(depth=1000, damping=1)
+eqs = hm.EquivalentSources(depth=1000, damping=1)
 
-# Fit the layer coefficients to the observed magnetic anomaly.
-eql.fit(coordinates, data.total_field_anomaly_nt)
+# Fit the sources coefficients to the observed magnetic anomaly.
+eqs.fit(coordinates, data.total_field_anomaly_nt)
 
 # Evaluate the data fit by calculating an R² score against the observed data.
-# This is a measure of how well layer the fits the data NOT how good the
+# This is a measure of how well the sources fit the data, NOT how good the
 # interpolation will be.
-print("R² score:", eql.score(coordinates, data.total_field_anomaly_nt))
+print("R² score:", eqs.score(coordinates, data.total_field_anomaly_nt))
 
 # Interpolate data on a regular grid with 500 m spacing. The interpolation
 # requires the height of the grid points (upward coordinate). By passing in
 # 1500 m, we're effectively upward-continuing the data (mean flight height is
 # 500 m).
-grid = eql.grid(upward=1500, spacing=500, data_names=["magnetic_anomaly"])
+grid = eqs.grid(upward=1500, spacing=500, data_names=["magnetic_anomaly"])
 
 # The grid is a xarray.Dataset with values, coordinates, and metadata
 print("\nGenerated grid:\n", grid)
diff --git a/examples/eql/harmonic_spherical.py b/examples/equivalent_sources/spherical.py
similarity index 86%
rename from examples/eql/harmonic_spherical.py
rename to examples/equivalent_sources/spherical.py
index 378a2768d..a6da62d79 100644
--- a/examples/eql/harmonic_spherical.py
+++ b/examples/equivalent_sources/spherical.py
@@ -13,10 +13,10 @@
 projecting the data may introduce errors due to the distortions caused by the
 projection.
 
-:class:`harmonica.EQLHarmonicSpherical` implements the equivalent layer
+:class:`harmonica.EquivalentSourcesSph` implements the equivalent sources
 technique in spherical coordinates. It has the same advantages as the Cartesian
-equivalent layer (:class:`harmonica.EQLHarmonic`) while taking into account the
-curvature of the Earth.
+equivalent sources (:class:`harmonica.EquivalentSources`) while taking into
+account the curvature of the Earth.
 """
 import numpy as np
 import cartopy.crs as ccrs
@@ -47,23 +47,23 @@
 # spherical (longitude, spherical_latitude, radius).
 coordinates = ellipsoid.geodetic_to_spherical(longitude, latitude, elevation)
 
-# Create the equivalent layer
-eql = hm.EQLHarmonicSpherical(damping=1e-3, relative_depth=10000)
+# Create the equivalent sources
+eqs = hm.EquivalentSourcesSph(damping=1e-3, relative_depth=10000)
 
-# Fit the layer coefficients to the observed magnetic anomaly
-eql.fit(coordinates, gravity_disturbance)
+# Fit the sources coefficients to the observed magnetic anomaly
+eqs.fit(coordinates, gravity_disturbance)
 
 # Evaluate the data fit by calculating an R² score against the observed data.
-# This is a measure of how well layer the fits the data NOT how good the
+# This is a measure of how well the sources fit the data, NOT how good the
 # interpolation will be.
-print("R² score:", eql.score(coordinates, gravity_disturbance))
+print("R² score:", eqs.score(coordinates, gravity_disturbance))
 
 # Interpolate data on a regular grid with 0.2 degrees spacing. The
 # interpolation requires the radius of the grid points (upward coordinate). By
 # passing in the maximum radius of the data, we're effectively
 # upward-continuing the data. The grid will be defined in spherical
 # coordinates.
-grid = eql.grid(
+grid = eqs.grid(
     upward=coordinates[-1].max(),
     spacing=0.2,
     data_names=["gravity_disturbance"],
diff --git a/harmonica/__init__.py b/harmonica/__init__.py
index 50d8ac94a..f0e389a61 100644
--- a/harmonica/__init__.py
+++ b/harmonica/__init__.py
@@ -16,8 +16,11 @@
 from .forward.tesseroid import tesseroid_gravity
 from .forward.prism import prism_gravity
 from .forward.prism_layer import prism_layer, DatasetAccessorPrismLayer
-from .equivalent_layer.harmonic import EQLHarmonic
-from .equivalent_layer.harmonic_spherical import EQLHarmonicSpherical
+from .equivalent_sources.cartesian import EquivalentSources, EQLHarmonic
+from .equivalent_sources.spherical import (
+    EquivalentSourcesSph,
+    EQLHarmonicSpherical,
+)
 
 # This file is generated automatically by setuptools_scm
 from . import _version
diff --git a/harmonica/equivalent_layer/__init__.py b/harmonica/equivalent_sources/__init__.py
similarity index 100%
rename from harmonica/equivalent_layer/__init__.py
rename to harmonica/equivalent_sources/__init__.py
diff --git a/harmonica/equivalent_layer/harmonic.py b/harmonica/equivalent_sources/cartesian.py
similarity index 91%
rename from harmonica/equivalent_layer/harmonic.py
rename to harmonica/equivalent_sources/cartesian.py
index 82b80c8aa..ab238b47d 100644
--- a/harmonica/equivalent_layer/harmonic.py
+++ b/harmonica/equivalent_sources/cartesian.py
@@ -5,7 +5,7 @@
 # This code is part of the Fatiando a Terra project (https://www.fatiando.org)
 #
 """
-Equivalent layer for generic harmonic functions in Cartesian coordinates
+Equivalent sources for generic harmonic functions in Cartesian coordinates
 """
 import warnings
 import numpy as np
@@ -24,11 +24,11 @@
 from ..forward.utils import distance_cartesian
 
 
-class EQLHarmonic(vdb.BaseGridder):
+class EquivalentSources(vdb.BaseGridder):
     r"""
-    Equivalent-layer for generic harmonic functions (gravity, magnetics, etc).
+    Equivalent sources for generic harmonic functions (gravity, magnetics).
 
-    This equivalent layer can be used for:
+    These equivalent sources can be used for:
 
     * Cartesian coordinates (geographic coordinates must be project before use)
     * Gravity and magnetic data (including derivatives)
@@ -37,7 +37,7 @@ class EQLHarmonic(vdb.BaseGridder):
     * Upward continuation
     * Finite-difference based derivative calculations
 
-    It cannot be used for:
+    They cannot be used for:
 
     * Regional or global data where Earth's curvature must be taken into
       account
@@ -69,8 +69,8 @@ class EQLHarmonic(vdb.BaseGridder):
         smoothness is imposed on the estimated coefficients.
         If None, no regularization is used.
     points : None or list of arrays (optional)
-        List containing the coordinates of the point sources used as the
-        equivalent layer. Coordinates are assumed to be in the following order:
+        List containing the coordinates of the equivalent point sources.
+        Coordinates are assumed to be in the following order:
         (``easting``, ``northing``, ``upward``).
         If None, will place one point source below each observation point at
         a fixed relative depth below the observation point [Cooper2000]_.
@@ -104,7 +104,7 @@ class EQLHarmonic(vdb.BaseGridder):
     Attributes
     ----------
     points_ : 2d-array
-        Coordinates of the point sources used to build the equivalent layer.
+        Coordinates of the equivalent point sources.
     coefs_ : array
         Estimated coefficients of every point source.
     region_ : tuple
@@ -157,7 +157,7 @@ def __init__(
 
     def fit(self, coordinates, data, weights=None):
         """
-        Fit the coefficients of the equivalent layer.
+        Fit the coefficients of the equivalent sources.
 
         The data region is captured and used as default for the
         :meth:`~harmonica.EQLHarmonic.grid` method.
@@ -217,9 +217,8 @@ def _build_points(self, coordinates):
         Returns
         -------
         points : tuple of arrays
-            Tuple containing the coordinates of the point sources used as the
-            equivalent layer, in the following order:
-            (``easting``, ``northing``, ``upward``).
+            Tuple containing the coordinates of the equivalent point sources,
+            in the following order: (``easting``, ``northing``, ``upward``).
         """
         if self.depth_type == "relative":
             return (
@@ -237,7 +236,7 @@ def _build_points(self, coordinates):
 
     def predict(self, coordinates):
         """
-        Evaluate the estimated equivalent layer on the given set of points.
+        Evaluate the estimated equivalent sources on the given set of points.
 
         Requires a fitted estimator (see :meth:`~harmonica.EQLHarmonic.fit`).
 
@@ -270,7 +269,7 @@ def jacobian(
         self, coordinates, points, dtype="float64"
     ):  # pylint: disable=no-self-use
         """
-        Make the Jacobian matrix for the equivalent layer.
+        Make the Jacobian matrix for the equivalent sources.
 
         Each column of the Jacobian is the Green's function for a single point
         source evaluated on all observation points.
@@ -283,8 +282,8 @@ def jacobian(
             Only ``easting``, ``northing`` and ``upward`` will be used, all
             subsequent coordinates will be ignored.
         points : tuple of arrays
-            Tuple of arrays containing the coordinates of the point sources
-            used as equivalent layer in the following order:
+            Tuple of arrays containing the coordinates of the equivalent point
+            sources in the following order:
             (``easting``, ``northing``, ``upward``).
         dtype : str or numpy dtype
             The type of the Jacobian array.
@@ -507,10 +506,42 @@ def profile(
         return table
 
 
+class EQLHarmonic(EquivalentSources):
+    """
+    DEPRECATED, use ``harmonica.EquivalentSources`` instead.
+
+    This class exists to support backward compatibility until next release.
+    """
+
+    def __init__(
+        self,
+        damping=None,
+        points=None,
+        depth=500,
+        depth_type="relative",
+        parallel=True,
+        **kwargs,
+    ):
+        warnings.warn(
+            "The 'EQLHarmonic' class has been renamed to 'EquivalentSources' "
+            + "and will be deprecated on the next release, "
+            + "please use 'EquivalentSources' instead.",
+            FutureWarning,
+        )
+        super().__init__(
+            damping=damping,
+            points=points,
+            depth=depth,
+            depth_type=depth_type,
+            parallel=parallel,
+            **kwargs,
+        )
+
+
 @jit(nopython=True)
 def greens_func_cartesian(east, north, upward, point_east, point_north, point_upward):
     """
-    Green's function for the equivalent layer in Cartesian coordinates
+    Green's function for the equivalent sources in Cartesian coordinates
 
     Uses Numba to speed up things.
     """
diff --git a/harmonica/equivalent_layer/harmonic_spherical.py b/harmonica/equivalent_sources/spherical.py
similarity index 87%
rename from harmonica/equivalent_layer/harmonic_spherical.py
rename to harmonica/equivalent_sources/spherical.py
index 416c62193..444afed53 100644
--- a/harmonica/equivalent_layer/harmonic_spherical.py
+++ b/harmonica/equivalent_sources/spherical.py
@@ -5,8 +5,9 @@
 # This code is part of the Fatiando a Terra project (https://www.fatiando.org)
 #
 """
-Equivalent layer for generic harmonic functions in spherical coordinates
+Equivalent sources for generic harmonic functions in spherical coordinates
 """
+import warnings
 import numpy as np
 from numba import jit
 from sklearn.utils.validation import check_is_fitted
@@ -23,11 +24,11 @@
 from ..forward.utils import distance_spherical
 
 
-class EQLHarmonicSpherical(vdb.BaseGridder):
+class EquivalentSourcesSph(vdb.BaseGridder):
     r"""
-    Equivalent-layer for generic harmonic functions in spherical coordinates
+    Equivalent sources for generic harmonic functions in spherical coordinates
 
-    This equivalent layer can be used for:
+    These equivalent sources can be used for:
 
     * Spherical coordinates (geographic coordinates must be converted before
       use)
@@ -39,7 +40,7 @@ class EQLHarmonicSpherical(vdb.BaseGridder):
     * Upward continuation
     * Finite-difference based derivative calculations
 
-    It cannot be used for:
+    They cannot be used for:
 
     * Joint inversion of multiple data types (e.g., gravity + gravity
       gradients)
@@ -69,8 +70,8 @@ class EQLHarmonicSpherical(vdb.BaseGridder):
         smoothness is imposed on the estimated coefficients.
         If None, no regularization is used.
     points : None or list of arrays (optional)
-        List containing the coordinates of the point sources used as the
-        equivalent layer. Coordinates are assumed to be in the following order:
+        List containing the coordinates of the equivalent point sources.
+        Coordinates are assumed to be in the following order:
         (``longitude``, ``latitude``, ``radius``). Both ``longitude`` and
         ``latitude`` must be in degrees and ``radius`` in meters.
         If None, will place one point source below each observation point at
@@ -91,13 +92,13 @@ class EQLHarmonicSpherical(vdb.BaseGridder):
     Attributes
     ----------
     points_ : 2d-array
-        Coordinates of the point sources used to build the equivalent layer.
+        Coordinates of the equivalent point sources.
     coefs_ : array
         Estimated coefficients of every point source.
     region_ : tuple
         The boundaries (``[W, E, S, N]``) of the data used to fit the
         interpolator. Used as the default region for the
-        :meth:`~harmonica.EQLHarmonicSpherical.grid` method.
+        :meth:`~harmonica.EQLHarmonicSph.grid` method.
     """
 
     # Set the default dimension names for generated outputs
@@ -127,7 +128,7 @@ def __init__(
 
     def fit(self, coordinates, data, weights=None):
         """
-        Fit the coefficients of the equivalent layer.
+        Fit the coefficients of the equivalent sources.
 
         The data region is captured and used as default for the
         :meth:`~harmonica.EQLHarmonicSpherical.grid` method.
@@ -170,7 +171,7 @@ def fit(self, coordinates, data, weights=None):
 
     def predict(self, coordinates):
         """
-        Evaluate the estimated equivalent layer on the given set of points.
+        Evaluate the estimated equivalent sources on the given set of points.
 
         Requires a fitted estimator
         (see :meth:`~harmonica.EQLHarmonicSpherical.fit`).
@@ -204,7 +205,7 @@ def jacobian(
         self, coordinates, points, dtype="float64"
     ):  # pylint: disable=no-self-use
         """
-        Make the Jacobian matrix for the equivalent layer.
+        Make the Jacobian matrix for the equivalent sources.
 
         Each column of the Jacobian is the Green's function for a single point
         source evaluated on all observation points.
@@ -217,8 +218,8 @@ def jacobian(
             Only ``longitude``, ``latitude`` and ``radius`` will be used, all
             subsequent coordinates will be ignored.
         points : tuple of arrays
-            Tuple of arrays containing the coordinates of the point sources
-            used as equivalent layer in the following order:
+            Tuple of arrays containing the coordinates of the equivalent point
+            sources in the following order:
             (``longitude``, ``latitude``, ``radius``).
         dtype : str or numpy dtype
             The type of the Jacobian array.
@@ -245,7 +246,7 @@ def grid(
         spacing=None,
         dims=None,
         data_names=None,
-        **kwargs
+        **kwargs,
     ):  # pylint: disable=arguments-differ
         """
         Interpolate the data onto a regular grid.
@@ -323,7 +324,7 @@ def scatter(
         dims=None,
         data_names=None,
         projection=None,
-        **kwargs
+        **kwargs,
     ):
         """
         .. warning ::
@@ -342,7 +343,7 @@ def profile(
         dims=None,
         data_names=None,
         projection=None,
-        **kwargs
+        **kwargs,
     ):
         """
         .. warning ::
@@ -354,12 +355,40 @@ def profile(
         raise NotImplementedError
 
 
+class EQLHarmonicSpherical(EquivalentSourcesSph):
+    """
+    DEPRECATED, use ``harmonica.EquivalentSourcesSph`` instead.
+
+    This class exists to support backward compatibility until next release.
+    """
+
+    def __init__(
+        self,
+        damping=None,
+        points=None,
+        relative_depth=500,
+        parallel=True,
+    ):
+        warnings.warn(
+            "The 'EQLHarmonic' class has been renamed to 'EquivalentSources' "
+            + "and will be removed on the next release, "
+            + "please use 'EquivalentSources' instead.",
+            FutureWarning,
+        )
+        super().__init__(
+            damping=damping,
+            points=points,
+            relative_depth=relative_depth,
+            parallel=parallel,
+        )
+
+
 @jit(nopython=True)
 def greens_func_spherical(
     longitude, latitude, radius, point_longitude, point_latitude, point_radius
 ):
     """
-    Green's function for the equivalent layer in spherical coordinates
+    Green's function for the equivalent sources in spherical coordinates
 
     Uses Numba to speed up things.
     """
diff --git a/harmonica/equivalent_layer/utils.py b/harmonica/equivalent_sources/utils.py
similarity index 97%
rename from harmonica/equivalent_layer/utils.py
rename to harmonica/equivalent_sources/utils.py
index 59de53422..1fdd62b73 100644
--- a/harmonica/equivalent_layer/utils.py
+++ b/harmonica/equivalent_sources/utils.py
@@ -5,7 +5,7 @@
 # This code is part of the Fatiando a Terra project (https://www.fatiando.org)
 #
 """
-Utility functions for equivalent layer gridders
+Utility functions for equivalent sources gridders
 """
 from warnings import warn
 from numba import jit, prange
diff --git a/harmonica/tests/test_eql_harmonic.py b/harmonica/tests/test_eq_sources_cartesian.py
similarity index 52%
rename from harmonica/tests/test_eql_harmonic.py
rename to harmonica/tests/test_eq_sources_cartesian.py
index 4498b16bd..3bdd74a76 100644
--- a/harmonica/tests/test_eql_harmonic.py
+++ b/harmonica/tests/test_eq_sources_cartesian.py
@@ -6,18 +6,19 @@
 #
 # pylint: disable=protected-access
 """
-Test the EQLHarmonic gridder
+Test the EquivalentSources gridder
 """
 import warnings
 import pytest
 import numpy as np
 import numpy.testing as npt
+import xarray.testing as xrt
 import verde as vd
 import verde.base as vdb
 
-from .. import EQLHarmonic, EQLHarmonicSpherical, point_mass_gravity
-from ..equivalent_layer.harmonic import greens_func_cartesian
-from ..equivalent_layer.utils import (
+from .. import EquivalentSources, EQLHarmonic, point_mass_gravity
+from ..equivalent_sources.cartesian import greens_func_cartesian
+from ..equivalent_sources.utils import (
     jacobian_numba_serial,
     pop_extra_coords,
 )
@@ -43,7 +44,7 @@ def test_pop_extra_coords():
 
 
 @require_numba
-def test_eql_harmonic_cartesian():
+def test_equivalent_sources_cartesian():
     """
     Check that predictions are reasonable when interpolating from one grid to
     a denser grid. Use Cartesian coordinates.
@@ -58,9 +59,9 @@ def test_eql_harmonic_cartesian():
     data = point_mass_gravity(coordinates, points, masses, field="g_z")
 
     # The interpolation should be perfect on the data points
-    eql = EQLHarmonic()
-    eql.fit(coordinates, data)
-    npt.assert_allclose(data, eql.predict(coordinates), rtol=1e-5)
+    eqs = EquivalentSources()
+    eqs.fit(coordinates, data)
+    npt.assert_allclose(data, eqs.predict(coordinates), rtol=1e-5)
 
     # Gridding onto a denser grid should be reasonably accurate when compared
     # to synthetic values
@@ -68,23 +69,23 @@ def test_eql_harmonic_cartesian():
     shape = (60, 60)
     grid = vd.grid_coordinates(region=region, shape=shape, extra_coords=upward)
     true = point_mass_gravity(grid, points, masses, field="g_z")
-    npt.assert_allclose(true, eql.predict(grid), rtol=1e-3)
+    npt.assert_allclose(true, eqs.predict(grid), rtol=1e-3)
 
     # Test grid method
-    grid = eql.grid(upward, shape=shape, region=region)
+    grid = eqs.grid(upward, shape=shape, region=region)
     npt.assert_allclose(true, grid.scalars, rtol=1e-3)
 
     # Test profile method
     point1 = (region[0], region[2])
     point2 = (region[0], region[3])
-    profile = eql.profile(point1, point2, upward, shape[0])
+    profile = eqs.profile(point1, point2, upward, shape[0])
     true = point_mass_gravity(
         (profile.easting, profile.northing, profile.upward), points, masses, field="g_z"
     )
     npt.assert_allclose(true, profile.scalars, rtol=1e-3)
 
 
-def test_eql_harmonic_small_data_cartesian():
+def test_equivalent_sources_small_data_cartesian():
     """
     Check predictions against synthetic data using few data points for speed
     Use Cartesian coordinates.
@@ -99,14 +100,14 @@ def test_eql_harmonic_small_data_cartesian():
     data = point_mass_gravity(coordinates, points, masses, field="g_z")
 
     # The interpolation should be perfect on the data points
-    eql = EQLHarmonic(depth=500)
-    eql.fit(coordinates, data)
-    npt.assert_allclose(data, eql.predict(coordinates), rtol=1e-5)
+    eqs = EquivalentSources(depth=500)
+    eqs.fit(coordinates, data)
+    npt.assert_allclose(data, eqs.predict(coordinates), rtol=1e-5)
 
     # Check that the proper source locations were set
     tmp = [i.ravel() for i in coordinates]
-    npt.assert_allclose(tmp[:2], eql.points_[:2], rtol=1e-5)
-    npt.assert_allclose(tmp[2] - 500, eql.points_[2], rtol=1e-5)
+    npt.assert_allclose(tmp[:2], eqs.points_[:2], rtol=1e-5)
+    npt.assert_allclose(tmp[2] - 500, eqs.points_[2], rtol=1e-5)
 
     # Gridding at higher altitude should be reasonably accurate when compared
     # to synthetic values
@@ -114,16 +115,16 @@ def test_eql_harmonic_small_data_cartesian():
     shape = (8, 8)
     grid = vd.grid_coordinates(region=region, shape=shape, extra_coords=upward)
     true = point_mass_gravity(grid, points, masses, field="g_z")
-    npt.assert_allclose(true, eql.predict(grid), rtol=0.08)
+    npt.assert_allclose(true, eqs.predict(grid), rtol=0.08)
 
     # Test grid method
-    grid = eql.grid(upward, shape=shape, region=region)
+    grid = eqs.grid(upward, shape=shape, region=region)
     npt.assert_allclose(true, grid.scalars, rtol=0.08)
 
     # Test profile method
     point1 = (region[0], region[2])
     point2 = (region[0], region[3])
-    profile = eql.profile(point1, point2, upward, 10)
+    profile = eqs.profile(point1, point2, upward, 10)
     true = point_mass_gravity(
         (profile.easting, profile.northing, profile.upward), points, masses, field="g_z"
     )
@@ -151,7 +152,7 @@ def fixture_coordinates():
     ],
     ids=["relative", "constant"],
 )
-def test_eql_harmonic_build_points(
+def test_equivalent_sources_build_points(
     coordinates,
     depth_type,
     upward_expected,
@@ -159,13 +160,13 @@ def test_eql_harmonic_build_points(
     """
     Check if build_points method works as expected
     """
-    eql = EQLHarmonic(depth=1.5e3, depth_type=depth_type)
-    points = eql._build_points(coordinates)
+    eqs = EquivalentSources(depth=1.5e3, depth_type=depth_type)
+    points = eqs._build_points(coordinates)
     expected = (*coordinates[:2], upward_expected)
     npt.assert_allclose(points, expected)
 
 
-def test_eql_harmonic_build_points_bacwkards(coordinates):
+def test_equivalent_sources_build_points_bacwkards(coordinates):
     """
     Check if the old relative_depth argument is well supported
 
@@ -177,27 +178,27 @@ def test_eql_harmonic_build_points_bacwkards(coordinates):
     expected_upward = coordinates[2] - depth
     # Check if FutureWarning is raised after passing relative_depth
     with warnings.catch_warnings(record=True) as warn:
-        eql = EQLHarmonic(relative_depth=depth)
+        eqs = EquivalentSources(relative_depth=depth)
         assert len(warn) == 1
         assert issubclass(warn[-1].category, FutureWarning)
     # Check if the `depth` and `depth_type` attributes are well fixed
-    npt.assert_allclose(eql.depth, depth)
-    assert eql.depth_type == "relative"
+    npt.assert_allclose(eqs.depth, depth)
+    assert eqs.depth_type == "relative"
     # Check if location of sources are correct
-    points = eql._build_points(coordinates)
+    points = eqs._build_points(coordinates)
     expected = (*coordinates[:2], expected_upward)
     npt.assert_allclose(points, expected)
 
 
-def test_eql_harmonic_invalid_depth_type():
+def test_equivalent_sources_invalid_depth_type():
     """
     Check if ValueError is raised if invalid depth_type is passed
     """
     with pytest.raises(ValueError):
-        EQLHarmonic(depth=300, depth_type="blabla")
+        EquivalentSources(depth=300, depth_type="blabla")
 
 
-def test_eql_harmonic_points_depth():
+def test_equivalent_sources_points_depth():
     """
     Check if the points coordinates are properly defined by the fit method
     """
@@ -213,29 +214,31 @@ def test_eql_harmonic_points_depth():
     data = point_mass_gravity(coordinates, points, masses, field="g_z")
 
     # Test with constant depth
-    eql = EQLHarmonic(depth=1.3e3, depth_type="constant")
-    eql.fit(coordinates, data)
+    eqs = EquivalentSources(depth=1.3e3, depth_type="constant")
+    eqs.fit(coordinates, data)
     expected_points = vdb.n_1d_arrays(
         (easting, northing, -1.3e3 * np.ones_like(easting)), n=3
     )
-    npt.assert_allclose(expected_points, eql.points_)
+    npt.assert_allclose(expected_points, eqs.points_)
 
     # Test with relative depth
-    eql = EQLHarmonic(depth=1.3e3, depth_type="relative")
-    eql.fit(coordinates, data)
+    eqs = EquivalentSources(depth=1.3e3, depth_type="relative")
+    eqs.fit(coordinates, data)
     expected_points = vdb.n_1d_arrays((easting, northing, upward - 1.3e3), n=3)
-    npt.assert_allclose(expected_points, eql.points_)
+    npt.assert_allclose(expected_points, eqs.points_)
 
     # Test with invalid depth_type
-    eql = EQLHarmonic(depth=300, depth_type="constant")  # init with valid depth_type
-    eql.depth_type = "blabla"  # change depth_type afterwards
-    points = eql._build_points(
+    eqs = EquivalentSources(
+        depth=300, depth_type="constant"
+    )  # init with valid depth_type
+    eqs.depth_type = "blabla"  # change depth_type afterwards
+    points = eqs._build_points(
         vd.grid_coordinates(region=(-1, 1, -1, 1), spacing=0.25, extra_coords=1)
     )
     assert points is None
 
 
-def test_eql_harmonic_custom_points_cartesian():
+def test_equivalent_sources_custom_points_cartesian():
     """
     Check that passing in custom points works and actually uses the points
     Use Cartesian coordinates.
@@ -254,24 +257,24 @@ def test_eql_harmonic_custom_points_cartesian():
         i.ravel()
         for i in vd.grid_coordinates(region=region, shape=(3, 3), extra_coords=-550)
     )
-    eql = EQLHarmonic(points=points_custom)
-    eql.fit(coordinates, data)
+    eqs = EquivalentSources(points=points_custom)
+    eqs.fit(coordinates, data)
 
     # Check that the proper source locations were set
-    npt.assert_allclose(points_custom, eql.points_, rtol=1e-5)
+    npt.assert_allclose(points_custom, eqs.points_, rtol=1e-5)
 
 
-def test_eql_harmonic_scatter_not_implemented():
+def test_equivalent_sources_scatter_not_implemented():
     """
     Check if scatter method raises a NotImplementedError
     """
-    eql = EQLHarmonic()
+    eqs = EquivalentSources()
     with pytest.raises(NotImplementedError):
-        eql.scatter()
+        eqs.scatter()
 
 
 @pytest.mark.use_numba
-def test_eql_harmonic_jacobian_cartesian():
+def test_equivalent_sources_jacobian_cartesian():
     """
     Test Jacobian matrix under symmetric system of point sources.
     Use Cartesian coordinates.
@@ -299,7 +302,7 @@ def test_eql_harmonic_jacobian_cartesian():
 
 
 @require_numba
-def test_eql_harmonic_cartesian_parallel():
+def test_equivalent_sources_cartesian_parallel():
     """
     Check predictions when parallel is enabled and disabled
     """
@@ -313,194 +316,52 @@ def test_eql_harmonic_cartesian_parallel():
     data = point_mass_gravity(coordinates, points, masses, field="g_z")
 
     # The predictions should be equal whether are run in parallel or in serial
-    eql_serial = EQLHarmonic(parallel=False)
-    eql_serial.fit(coordinates, data)
-    eql_parallel = EQLHarmonic(parallel=True)
-    eql_parallel.fit(coordinates, data)
+    eqs_serial = EquivalentSources(parallel=False)
+    eqs_serial.fit(coordinates, data)
+    eqs_parallel = EquivalentSources(parallel=True)
+    eqs_parallel.fit(coordinates, data)
 
     upward = 0
     shape = (60, 60)
-    grid_serial = eql_serial.grid(upward, shape=shape, region=region)
-    grid_parallel = eql_parallel.grid(upward, shape=shape, region=region)
+    grid_serial = eqs_serial.grid(upward, shape=shape, region=region)
+    grid_parallel = eqs_parallel.grid(upward, shape=shape, region=region)
     npt.assert_allclose(grid_serial.scalars, grid_parallel.scalars, rtol=1e-7)
 
 
-@require_numba
-def test_eql_harmonic_spherical():
-    """
-    Check that predictions are reasonable when interpolating from one grid to
-    a denser grid. Use spherical coordinates.
+@pytest.mark.parametrize("depth_type", ("constant", "relative"))
+def test_backward_eqlharmonic(depth_type):
     """
-    region = (-70, -60, -40, -30)
-    radius = 6400e3
-    # Build synthetic point masses
-    points = vd.grid_coordinates(
-        region=region, shape=(6, 6), extra_coords=radius - 500e3
-    )
-    masses = vd.datasets.CheckerBoard(amplitude=1e13, region=region).predict(points)
-    # Define a set of observation points
-    coordinates = vd.grid_coordinates(
-        region=region, shape=(40, 40), extra_coords=radius
-    )
-    # Get synthetic data
-    data = point_mass_gravity(
-        coordinates, points, masses, field="g_z", coordinate_system="spherical"
-    )
+    Check backward compatibility with to-be-deprecated EQLHarmonic class
 
-    # The interpolation should be perfect on the data points
-    eql = EQLHarmonicSpherical(relative_depth=500e3)
-    eql.fit(coordinates, data)
-    npt.assert_allclose(data, eql.predict(coordinates), rtol=1e-5)
-
-    # Gridding onto a denser grid should be reasonably accurate when compared
-    # to synthetic values
-    upward = radius
-    shape = (60, 60)
-    grid = vd.grid_coordinates(region=region, shape=shape, extra_coords=upward)
-    true = point_mass_gravity(
-        grid, points, masses, field="g_z", coordinate_system="spherical"
-    )
-    npt.assert_allclose(true, eql.predict(grid), rtol=1e-3)
-
-    # Test grid method
-    grid = eql.grid(upward, shape=shape, region=region)
-    npt.assert_allclose(true, grid.scalars, rtol=1e-3)
-
-
-def test_eql_harmonic_small_data_spherical():
-    """
-    Check predictions against synthetic data using few data points for speed
-    Use spherical coordinates.
+    Check if FutureWarning is raised on initialization
     """
-    region = (-70, -60, -40, -30)
-    radius = 6400e3
+    region = (-3e3, -1e3, 5e3, 7e3)
     # Build synthetic point masses
-    points = vd.grid_coordinates(
-        region=region, shape=(6, 6), extra_coords=radius - 500e3
-    )
+    points = vd.grid_coordinates(region=region, shape=(6, 6), extra_coords=-1e3)
     masses = vd.datasets.CheckerBoard(amplitude=1e13, region=region).predict(points)
-    # Define a set of observation points
-    coordinates = vd.grid_coordinates(region=region, shape=(8, 8), extra_coords=radius)
+    # Define a set of observation points with variable elevation coordinates
+    easting, northing = vd.grid_coordinates(region=region, shape=(5, 5))
+    upward = np.arange(25, dtype=float).reshape((5, 5))
+    coordinates = (easting, northing, upward)
     # Get synthetic data
-    data = point_mass_gravity(
-        coordinates, points, masses, field="g_z", coordinate_system="spherical"
-    )
+    data = point_mass_gravity(coordinates, points, masses, field="g_z")
 
-    # The interpolation should be perfect on the data points
-    eql = EQLHarmonicSpherical(relative_depth=500e3)
-    eql.fit(coordinates, data)
-    npt.assert_allclose(data, eql.predict(coordinates), rtol=1e-5)
+    # Fit EquivalentSources instance
+    eqs = EquivalentSources(depth=1.3e3, depth_type=depth_type)
+    eqs.fit(coordinates, data)
 
-    # Check that the proper source locations were set
-    tmp = [i.ravel() for i in coordinates]
-    npt.assert_allclose(tmp[:2], eql.points_[:2], rtol=1e-5)
-    npt.assert_allclose(tmp[2] - 500e3, eql.points_[2], rtol=1e-5)
+    # Fit deprecated EQLHarmonic instance
+    # (check if FutureWarning is raised)
+    with warnings.catch_warnings(record=True) as warn:
+        eql_harmonic = EQLHarmonic(depth=1.3e3, depth_type=depth_type)
+        assert len(warn) == 1
+        assert issubclass(warn[-1].category, FutureWarning)
+    eql_harmonic.fit(coordinates, data)
 
-    # Gridding at higher altitude should be reasonably accurate when compared
-    # to synthetic values
-    upward = radius + 2e3
+    # Check if both gridders are equivalent
+    npt.assert_allclose(eqs.points_, eql_harmonic.points_)
     shape = (8, 8)
-    grid = vd.grid_coordinates(region=region, shape=shape, extra_coords=upward)
-    true = point_mass_gravity(
-        grid, points, masses, field="g_z", coordinate_system="spherical"
+    xrt.assert_allclose(
+        eqs.grid(upward=2e3, shape=shape, region=region),
+        eql_harmonic.grid(upward=2e3, shape=shape, region=region),
     )
-    npt.assert_allclose(true, eql.predict(grid), rtol=0.05)
-
-    # Test grid method
-    grid = eql.grid(upward, shape=shape, region=region)
-    npt.assert_allclose(true, grid.scalars, rtol=0.05)
-
-
-def test_eql_harmonic_custom_points_spherical():
-    """
-    Check that passing in custom points works and actually uses the points
-    Use spherical coordinates.
-    """
-    region = (-70, -60, -40, -30)
-    radius = 6400e3
-    # Build synthetic point masses
-    points = vd.grid_coordinates(
-        region=region, shape=(6, 6), extra_coords=radius - 500e3
-    )
-    masses = vd.datasets.CheckerBoard(amplitude=1e13, region=region).predict(points)
-    # Define a set of observation points
-    coordinates = vd.grid_coordinates(region=region, shape=(5, 5), extra_coords=radius)
-    # Get synthetic data
-    data = point_mass_gravity(
-        coordinates, points, masses, field="g_z", coordinate_system="spherical"
-    )
-
-    # Pass a custom set of point sources
-    points_custom = tuple(
-        i.ravel()
-        for i in vd.grid_coordinates(
-            region=region, shape=(3, 3), extra_coords=radius - 500e3
-        )
-    )
-    eql = EQLHarmonicSpherical(points=points_custom)
-    eql.fit(coordinates, data)
-
-    # Check that the proper source locations were set
-    npt.assert_allclose(points_custom, eql.points_, rtol=1e-5)
-
-
-def test_eql_harmonic_spherical_scatter_not_implemented():
-    """
-    Check if scatter method raises a NotImplementedError
-    """
-    eql = EQLHarmonicSpherical()
-    with pytest.raises(NotImplementedError):
-        eql.scatter()
-
-
-def test_eql_harmonic_spherical_profile_not_implemented():
-    """
-    Check if scatter method raises a NotImplementedError
-    """
-    eql = EQLHarmonicSpherical()
-    with pytest.raises(NotImplementedError):
-        eql.profile(point1=(1, 1), point2=(2, 2), size=3)
-
-
-def test_eql_harmonic_spherical_no_projection():
-    """
-    Check if projection is not a valid argument of grid method
-    """
-    eql = EQLHarmonicSpherical()
-    with pytest.raises(TypeError):
-        eql.grid(upward=10, projection=lambda a, b: (a * 2, b * 2))
-
-
-@require_numba
-def test_eql_harmonic_spherical_parallel():
-    """
-    Check predictions when parallel is enabled and disabled
-    """
-    region = (-70, -60, -40, -30)
-    radius = 6400e3
-    # Build synthetic point masses
-    points = vd.grid_coordinates(
-        region=region, shape=(6, 6), extra_coords=radius - 500e3
-    )
-    masses = vd.datasets.CheckerBoard(amplitude=1e13, region=region).predict(points)
-    # Define a set of observation points
-    coordinates = vd.grid_coordinates(
-        region=region, shape=(40, 40), extra_coords=radius
-    )
-    # Get synthetic data
-    data = point_mass_gravity(
-        coordinates, points, masses, field="g_z", coordinate_system="spherical"
-    )
-
-    # The predictions should be equal whether are run in parallel or in serial
-    relative_depth = 500e3
-    eql_serial = EQLHarmonicSpherical(relative_depth=relative_depth, parallel=False)
-    eql_serial.fit(coordinates, data)
-    eql_parallel = EQLHarmonicSpherical(relative_depth=relative_depth, parallel=True)
-    eql_parallel.fit(coordinates, data)
-
-    upward = radius
-    shape = (60, 60)
-    grid_serial = eql_serial.grid(upward, shape=shape, region=region)
-    grid_parallel = eql_parallel.grid(upward, shape=shape, region=region)
-    npt.assert_allclose(grid_serial.scalars, grid_parallel.scalars, rtol=1e-7)
diff --git a/harmonica/tests/test_eq_sources_spherical.py b/harmonica/tests/test_eq_sources_spherical.py
new file mode 100644
index 000000000..f40bd5249
--- /dev/null
+++ b/harmonica/tests/test_eq_sources_spherical.py
@@ -0,0 +1,240 @@
+# Copyright (c) 2018 The Harmonica Developers.
+# Distributed under the terms of the BSD 3-Clause License.
+# SPDX-License-Identifier: BSD-3-Clause
+#
+# This code is part of the Fatiando a Terra project (https://www.fatiando.org)
+#
+# pylint: disable=protected-access
+"""
+Test the EquivalentSourcesSph gridder
+"""
+import warnings
+import pytest
+import numpy.testing as npt
+import xarray.testing as xrt
+import verde as vd
+
+from .. import EquivalentSourcesSph, EQLHarmonicSpherical, point_mass_gravity
+from .utils import require_numba
+
+
+@require_numba
+def test_equivalent_sources_spherical():
+    """
+    Check that predictions are reasonable when interpolating from one grid to
+    a denser grid.
+    """
+    region = (-70, -60, -40, -30)
+    radius = 6400e3
+    # Build synthetic point masses
+    points = vd.grid_coordinates(
+        region=region, shape=(6, 6), extra_coords=radius - 500e3
+    )
+    masses = vd.datasets.CheckerBoard(amplitude=1e13, region=region).predict(points)
+    # Define a set of observation points
+    coordinates = vd.grid_coordinates(
+        region=region, shape=(40, 40), extra_coords=radius
+    )
+    # Get synthetic data
+    data = point_mass_gravity(
+        coordinates, points, masses, field="g_z", coordinate_system="spherical"
+    )
+
+    # The interpolation should be perfect on the data points
+    eqs = EquivalentSourcesSph(relative_depth=500e3)
+    eqs.fit(coordinates, data)
+    npt.assert_allclose(data, eqs.predict(coordinates), rtol=1e-5)
+
+    # Gridding onto a denser grid should be reasonably accurate when compared
+    # to synthetic values
+    upward = radius
+    shape = (60, 60)
+    grid = vd.grid_coordinates(region=region, shape=shape, extra_coords=upward)
+    true = point_mass_gravity(
+        grid, points, masses, field="g_z", coordinate_system="spherical"
+    )
+    npt.assert_allclose(true, eqs.predict(grid), rtol=1e-3)
+
+    # Test grid method
+    grid = eqs.grid(upward, shape=shape, region=region)
+    npt.assert_allclose(true, grid.scalars, rtol=1e-3)
+
+
+def test_equivalent_sources_small_data_spherical():
+    """
+    Check predictions against synthetic data using few data points for speed
+    Use spherical coordinates.
+    """
+    region = (-70, -60, -40, -30)
+    radius = 6400e3
+    # Build synthetic point masses
+    points = vd.grid_coordinates(
+        region=region, shape=(6, 6), extra_coords=radius - 500e3
+    )
+    masses = vd.datasets.CheckerBoard(amplitude=1e13, region=region).predict(points)
+    # Define a set of observation points
+    coordinates = vd.grid_coordinates(region=region, shape=(8, 8), extra_coords=radius)
+    # Get synthetic data
+    data = point_mass_gravity(
+        coordinates, points, masses, field="g_z", coordinate_system="spherical"
+    )
+
+    # The interpolation should be perfect on the data points
+    eqs = EquivalentSourcesSph(relative_depth=500e3)
+    eqs.fit(coordinates, data)
+    npt.assert_allclose(data, eqs.predict(coordinates), rtol=1e-5)
+
+    # Check that the proper source locations were set
+    tmp = [i.ravel() for i in coordinates]
+    npt.assert_allclose(tmp[:2], eqs.points_[:2], rtol=1e-5)
+    npt.assert_allclose(tmp[2] - 500e3, eqs.points_[2], rtol=1e-5)
+
+    # Gridding at higher altitude should be reasonably accurate when compared
+    # to synthetic values
+    upward = radius + 2e3
+    shape = (8, 8)
+    grid = vd.grid_coordinates(region=region, shape=shape, extra_coords=upward)
+    true = point_mass_gravity(
+        grid, points, masses, field="g_z", coordinate_system="spherical"
+    )
+    npt.assert_allclose(true, eqs.predict(grid), rtol=0.05)
+
+    # Test grid method
+    grid = eqs.grid(upward, shape=shape, region=region)
+    npt.assert_allclose(true, grid.scalars, rtol=0.05)
+
+
+def test_equivalent_sources_custom_points_spherical():
+    """
+    Check that passing in custom points works and actually uses the points
+    Use spherical coordinates.
+    """
+    region = (-70, -60, -40, -30)
+    radius = 6400e3
+    # Build synthetic point masses
+    points = vd.grid_coordinates(
+        region=region, shape=(6, 6), extra_coords=radius - 500e3
+    )
+    masses = vd.datasets.CheckerBoard(amplitude=1e13, region=region).predict(points)
+    # Define a set of observation points
+    coordinates = vd.grid_coordinates(region=region, shape=(5, 5), extra_coords=radius)
+    # Get synthetic data
+    data = point_mass_gravity(
+        coordinates, points, masses, field="g_z", coordinate_system="spherical"
+    )
+
+    # Pass a custom set of point sources
+    points_custom = tuple(
+        i.ravel()
+        for i in vd.grid_coordinates(
+            region=region, shape=(3, 3), extra_coords=radius - 500e3
+        )
+    )
+    eqs = EquivalentSourcesSph(points=points_custom)
+    eqs.fit(coordinates, data)
+
+    # Check that the proper source locations were set
+    npt.assert_allclose(points_custom, eqs.points_, rtol=1e-5)
+
+
+def test_equivalent_sources_scatter_not_implemented():
+    """
+    Check if scatter method raises a NotImplementedError
+    """
+    eqs = EquivalentSourcesSph()
+    with pytest.raises(NotImplementedError):
+        eqs.scatter()
+
+
+def test_equivalent_sources_profile_not_implemented():
+    """
+    Check if scatter method raises a NotImplementedError
+    """
+    eqs = EquivalentSourcesSph()
+    with pytest.raises(NotImplementedError):
+        eqs.profile(point1=(1, 1), point2=(2, 2), size=3)
+
+
+def test_equivalent_sources_spherical_no_projection():
+    """
+    Check if projection is not a valid argument of grid method
+    """
+    eqs = EquivalentSourcesSph()
+    with pytest.raises(TypeError):
+        eqs.grid(upward=10, projection=lambda a, b: (a * 2, b * 2))
+
+
+@require_numba
+def test_equivalent_sources_spherical_parallel():
+    """
+    Check predictions when parallel is enabled and disabled
+    """
+    region = (-70, -60, -40, -30)
+    radius = 6400e3
+    # Build synthetic point masses
+    points = vd.grid_coordinates(
+        region=region, shape=(6, 6), extra_coords=radius - 500e3
+    )
+    masses = vd.datasets.CheckerBoard(amplitude=1e13, region=region).predict(points)
+    # Define a set of observation points
+    coordinates = vd.grid_coordinates(
+        region=region, shape=(40, 40), extra_coords=radius
+    )
+    # Get synthetic data
+    data = point_mass_gravity(
+        coordinates, points, masses, field="g_z", coordinate_system="spherical"
+    )
+
+    # The predictions should be equal whether are run in parallel or in serial
+    relative_depth = 500e3
+    eqs_serial = EquivalentSourcesSph(relative_depth=relative_depth, parallel=False)
+    eqs_serial.fit(coordinates, data)
+    eqs_parallel = EquivalentSourcesSph(relative_depth=relative_depth, parallel=True)
+    eqs_parallel.fit(coordinates, data)
+
+    upward = radius
+    shape = (60, 60)
+    grid_serial = eqs_serial.grid(upward, shape=shape, region=region)
+    grid_parallel = eqs_parallel.grid(upward, shape=shape, region=region)
+    npt.assert_allclose(grid_serial.scalars, grid_parallel.scalars, rtol=1e-7)
+
+
+def test_backward_eqlharmonicspherical():
+    """
+    Check backward compatibility with to-be-deprecated EQLHarmonicSph
+
+    Check if FutureWarning is raised on initialization
+    """
+    region = (-70, -60, -40, -30)
+    radius = 6400e3
+    # Build synthetic point masses
+    points = vd.grid_coordinates(
+        region=region, shape=(6, 6), extra_coords=radius - 500e3
+    )
+    masses = vd.datasets.CheckerBoard(amplitude=1e13, region=region).predict(points)
+    # Define a set of observation points
+    coordinates = vd.grid_coordinates(region=region, shape=(8, 8), extra_coords=radius)
+    # Get synthetic data
+    data = point_mass_gravity(
+        coordinates, points, masses, field="g_z", coordinate_system="spherical"
+    )
+
+    # Fit EquivalentSourcesSph instance
+    eqs = EquivalentSourcesSph(relative_depth=1.3e3)
+    eqs.fit(coordinates, data)
+
+    # Fit deprecated EQLHarmonicSpherical instance
+    # (check if FutureWarning is raised)
+    with warnings.catch_warnings(record=True) as warn:
+        eql_harmonic = EQLHarmonicSpherical(relative_depth=1.3e3)
+        assert len(warn) == 1
+        assert issubclass(warn[-1].category, FutureWarning)
+    eql_harmonic.fit(coordinates, data)
+
+    # Check if both gridders are equivalent
+    npt.assert_allclose(eqs.points_, eql_harmonic.points_)
+    shape = (8, 8)
+    xrt.assert_allclose(
+        eqs.grid(upward=6405e3, shape=shape, region=region),
+        eql_harmonic.grid(upward=6405e3, shape=shape, region=region),
+    )