diff --git a/city_metrix/layers/__init__.py b/city_metrix/layers/__init__.py
index a0bbd82f..0a3817b2 100644
--- a/city_metrix/layers/__init__.py
+++ b/city_metrix/layers/__init__.py
@@ -28,3 +28,4 @@
 from .cams import Cams
 from .vegetation_water_map import VegetationWaterMap
 from .height_above_nearest_drainage import HeightAboveNearestDrainage
+from .riparian_areas import RiparianAreas
diff --git a/city_metrix/layers/height_above_nearest_drainage.py b/city_metrix/layers/height_above_nearest_drainage.py
index fbb2c0ba..c7355152 100644
--- a/city_metrix/layers/height_above_nearest_drainage.py
+++ b/city_metrix/layers/height_above_nearest_drainage.py
@@ -12,12 +12,14 @@ class HeightAboveNearestDrainage(Layer):
                             default is 30, other options - 90
         river_head: number of river head threshold cells
                     default is 1000, other options - 100, 5000
+        thresh: flow accumuation threshold, default is 1
     """
 
-    def __init__(self, spatial_resolution=30, river_head=1000, **kwargs):
+    def __init__(self, spatial_resolution=30, river_head=1000, thresh=1, **kwargs):
         super().__init__(**kwargs)
         self.spatial_resolution = spatial_resolution
         self.river_head = river_head
+        self.thresh = thresh
 
     def get_data(self, bbox):
         if self.spatial_resolution == 30 and self.river_head == 100:
@@ -37,7 +39,7 @@ def get_data(self, bbox):
         swbd = ee.Image('MODIS/MOD44W/MOD44W_005_2000_02_24').select('water_mask')
         swbdMask = swbd.unmask().Not().focal_median(1)
 
-        thresh = 1
+        thresh = self.thresh
         HANDthresh = hand.lte(thresh).focal_max(1).focal_mode(2, 'circle', 'pixels', 5).mask(swbdMask)
         HANDthresh = HANDthresh.mask(HANDthresh)
 
diff --git a/city_metrix/layers/riparian_areas.py b/city_metrix/layers/riparian_areas.py
new file mode 100644
index 00000000..18d1bac1
--- /dev/null
+++ b/city_metrix/layers/riparian_areas.py
@@ -0,0 +1,61 @@
+import xarray as xr
+import numpy as np
+import ee
+from scipy.ndimage import distance_transform_edt
+
+from .layer import Layer, get_image_collection
+from .height_above_nearest_drainage import HeightAboveNearestDrainage
+
+
+class RiparianAreas(Layer):
+    """
+    Attributes:
+        spatial_resolution: raster resolution in meters (see https://github.com/stac-extensions/raster)
+                            default is 30, other options - 90
+        river_head: number of river head threshold cells
+                    default is 1000, other options - 100, 5000
+        thresh: flow accumuation threshold, default is 0
+    """
+
+    def __init__(self, spatial_resolution=30, river_head=1000, thresh=0, **kwargs):
+        super().__init__(**kwargs)
+        self.spatial_resolution = spatial_resolution
+        self.river_head = river_head
+        self.thresh = thresh
+
+    def get_data(self, bbox):
+        # read HAND data to generate drainage paths
+        hand = HeightAboveNearestDrainage(spatial_resolution=self.spatial_resolution, river_head=self.river_head, thresh=self.thresh).get_data(bbox)
+
+        # Read surface water occurance
+        water = ee.Image('JRC/GSW1_3/GlobalSurfaceWater').select(['occurrence']).gte(50)
+        water_da = get_image_collection(
+            ee.ImageCollection(water),
+            bbox,
+            self.spatial_resolution,
+            "water"
+        ).occurrence
+
+        combWater = np.maximum(hand.fillna(0), water_da.fillna(0)) > 0
+        combWater = combWater.fillna(False)
+
+        # Buffer waterways by riparian zone definitions
+        distance_arr = xr.apply_ufunc(
+            distance_transform_edt,
+            ~combWater,  # True where not-water
+            input_core_dims=[('y', 'x')],
+            output_core_dims=[('y', 'x')],
+            dask='parallelized',
+            kwargs={'sampling': self.spatial_resolution}
+        )
+
+        halfpixel = self.spatial_resolution * 0.5
+        # https://doi.org/10.1016/j.jenvman.2019.109391
+        distance_200 = distance_arr.where(distance_arr <= 200)
+        nutrientBuffer = (distance_200 <= (3.0 - halfpixel))
+        floraBuffer = (distance_200 <= (24.0 - halfpixel))
+        birdBuffer = (distance_200 <= (144.0 - halfpixel))
+
+        riparianMask = birdBuffer
+
+        return riparianMask
diff --git a/environment.yml b/environment.yml
index a8d008c4..ce8e11d3 100644
--- a/environment.yml
+++ b/environment.yml
@@ -28,5 +28,7 @@ dependencies:
   - scikit-image=0.24.0
   - exactextract=0.2.0
   - cfgrib=0.9.15.0
+  - scipy=1.13.1
+  - numpy=1.23.0
   - pip:
       - overturemaps==0.6.0
diff --git a/setup.py b/setup.py
index a1ac5824..972be4bc 100644
--- a/setup.py
+++ b/setup.py
@@ -28,6 +28,8 @@
         "timezonefinder",
         "scikit-image>=0.24.0",
         "exactextract>=0.2.0",
-        "cfgrib"
+        "cfgrib",
+        "scipy",
+        "numpy"
     ],
 )
diff --git a/tests/test_layers.py b/tests/test_layers.py
index e16f3b5e..77933b39 100644
--- a/tests/test_layers.py
+++ b/tests/test_layers.py
@@ -25,6 +25,7 @@
     OpenStreetMap,
     OpenStreetMapClass,
     OvertureBuildings,
+    RiparianAreas,
     Sentinel2Level2,
     SmartSurfaceLULC,
     TreeCanopyHeight,
@@ -133,6 +134,10 @@ def test_overture_buildings():
     data = OvertureBuildings().get_data(BBOX)
     assert np.size(data) > 0
 
+def test_riparian_areas():
+    data = RiparianAreas().get_data(BBOX)
+    assert np.size(data) > 0
+
 @pytest.mark.skip(reason="layer is deprecated")
 def test_sentinel_2_level2():
     sentinel_2_bands = ["green"]