diff --git a/plasticparcels/constructors.py b/plasticparcels/constructors.py
index 62dfa7b..04cd5e7 100644
--- a/plasticparcels/constructors.py
+++ b/plasticparcels/constructors.py
@@ -232,6 +232,11 @@ def create_particleset(fieldset, settings, release_locations):
     else:
         plastic_amounts = np.full_like(lons, np.nan)
 
+    if 'depths' in release_locations.keys():
+        depths = release_locations['depths']
+    else:
+        depths = np.full_like(lons, 0.)
+
     # Set particle properties
     plastic_densities = np.full(lons.shape, settings['plastictype']['plastic_density'])
     plastic_diameters = np.full(lons.shape, settings['plastictype']['plastic_diameter'])
@@ -254,6 +259,7 @@ def create_particleset(fieldset, settings, release_locations):
                                  PlasticParticle,
                                  lon=lons,
                                  lat=lats,
+                                 depth=depths,
                                  plastic_diameter=plastic_diameters,
                                  plastic_density=plastic_densities,
                                  wind_coefficient=wind_coefficients,
diff --git a/plasticparcels/kernels.py b/plasticparcels/kernels.py
index 7355fad..a1c59f9 100644
--- a/plasticparcels/kernels.py
+++ b/plasticparcels/kernels.py
@@ -402,7 +402,7 @@ def PolyTEOS10_bsq(particle, fieldset, time):
       Oceanic Technology, 20, 730-741.
 
     """
-    Z = - particle.depth  # note: use negative depths!
+    Z = - math.fabs(particle.depth)  # note: use negative depths!
     SA = fieldset.absolute_salinity[time, particle.depth, particle.lat, particle.lon]
     CT = fieldset.conservative_temperature[time, particle.depth, particle.lat, particle.lon]
 
@@ -410,6 +410,16 @@ def PolyTEOS10_bsq(particle, fieldset, time):
     CTu = 40.
     Zu = 1.0e+04
     deltaS = 32.
+
+    zz = -Z / Zu
+    R00 = 4.6494977072e+01
+    R01 = -5.2099962525e+00
+    R02 = 2.2601900708e-01
+    R03 = 6.4326772569e-02
+    R04 = 1.5616995503e-02
+    R05 = -1.7243708991e-03
+    r0 = (((((R05 * zz + R04) * zz + R03) * zz + R02) * zz + R01) * zz + R00) * zz
+
     R000 = 8.0189615746e+02
     R100 = 8.6672408165e+02
     R200 = -1.7864682637e+03
@@ -469,8 +479,9 @@ def PolyTEOS10_bsq(particle, fieldset, time):
     rz2 = (R022 * tt + R112 * ss + R012) * tt + (R202 * ss + R102) * ss + R002
     rz1 = (((R041 * tt + R131 * ss + R031) * tt + (R221 * ss + R121) * ss + R021) * tt + ((R311 * ss + R211) * ss + R111) * ss + R011) * tt + (((R401 * ss + R301) * ss + R201) * ss + R101) * ss + R001
     rz0 = (((((R060 * tt + R150 * ss + R050) * tt + (R240 * ss + R140) * ss + R040) * tt + ((R330 * ss + R230) * ss + R130) * ss + R030) * tt + (((R420 * ss + R320) * ss + R220) * ss + R120) * ss + R020) * tt + ((((R510 * ss + R410) * ss + R310) * ss + R210) * ss + R110) * ss + R010) * tt + (((((R600 * ss + R500) * ss + R400) * ss + R300) * ss + R200) * ss + R100) * ss + R000
-    particle.seawater_density = ((rz3 * zz + rz2) * zz + rz1) * zz + rz0
-    # particle.sw_surface_density = rz0
+    r = ((rz3 * zz + rz2) * zz + rz1) * zz + rz0
+
+    particle.seawater_density = r0 + r
 
 
 def VerticalMixing(particle, fieldset, time):
diff --git a/tests/test_kernels.py b/tests/test_kernels.py
index 0cf758c..89b644f 100644
--- a/tests/test_kernels.py
+++ b/tests/test_kernels.py
@@ -243,3 +243,50 @@ def test_mixing():
 
     # Assert that the particles move from their initial location
     assert (np.sum(np.abs(pset.lon - pset_mixing.lon)) > 0.) & (np.sum(np.abs(pset.lat - pset_mixing.lat)) > 0.)
+
+@pytest.mark.parametrize("mode", ["scipy", "jit"])
+def test_TEOSdensity_kernels(mode):
+    """ Adapted test from Parcels v3 codebase.
+    """
+    settings_file = 'tests/test_data/test_settings.json'
+    settings = pp.utils.load_settings(settings_file)
+
+    settings['simulation'] = make_standard_simulation_settings()
+    settings['plastictype'] = make_standard_plastictype_settings()
+
+    # Turn off all other processes
+    settings['use_3D'] = False
+    settings['use_biofouling'] = False
+    settings['use_stokes'] = False
+    settings['use_wind'] = False
+    settings['use_mixing'] = False
+
+    def generate_fieldset(xdim=2, ydim=2, zdim=2, tdim=1):
+        lon = np.linspace(0.0, 10.0, xdim, dtype=np.float32)
+        lat = np.linspace(0.0, 10.0, ydim, dtype=np.float32)
+        depth = np.linspace(0, 2000, zdim, dtype=np.float32)
+        time = np.zeros(tdim, dtype=np.float64)
+        U = np.ones((tdim, zdim, ydim, xdim))
+        V = np.ones((tdim, zdim, ydim, xdim))
+        abs_salinity = 30 * np.ones((tdim, zdim, ydim, xdim))
+        cons_temperature = 10 * np.ones((tdim, zdim, ydim, xdim))
+        dimensions = {"lat": lat, "lon": lon, "depth": depth, "time": time}
+        data = {
+            "U": np.array(U, dtype=np.float32),
+            "V": np.array(V, dtype=np.float32),
+            "absolute_salinity": np.array(abs_salinity, dtype=np.float32),
+            "conservative_temperature": np.array(cons_temperature, dtype=np.float32),
+        }
+        return (data, dimensions)
+
+    data, dimensions = generate_fieldset()
+    fieldset = parcels.FieldSet.from_data(data, dimensions)
+
+    release_locations = {'lons': [5], 'lats': [5], 'depths': [1000],
+                         'plastic_amount': [1]}
+
+    pset = pp.constructors.create_particleset(fieldset, settings, release_locations)
+
+    pset.execute(pp.kernels.PolyTEOS10_bsq, runtime=1)
+
+    assert np.allclose(pset[0].seawater_density, 1027.45140)