diff --git a/.github/workflows/run-tests.yml b/.github/workflows/run-tests.yml
index 4611307df..0181bbe1f 100644
--- a/.github/workflows/run-tests.yml
+++ b/.github/workflows/run-tests.yml
@@ -18,7 +18,7 @@ jobs:
           python -m pip install flake8==7.1.0 flake8-docstrings==1.7.0 flake8-annotations==3.1.1
       - name: Run lint with flake8
         run: |
-          flake8 src floodresilience
+          flake8 src floodresilience wrfhydro
 
   linting-pylint:
     name: Lint Python with pylint
@@ -46,7 +46,7 @@ jobs:
           python -m pip install pylint==3.2.6
       - name: Run lint with pylint
         run: |
-          pylint src floodresilience
+          pylint src floodresilience wrfhydro
 
   unit-tests:
     name: Run unit tests with pytest
diff --git a/floodresilience/run_all.py b/floodresilience/run_all.py
index f13c5836d..e5d8f76d9 100644
--- a/floodresilience/run_all.py
+++ b/floodresilience/run_all.py
@@ -29,35 +29,35 @@
 from floodresilience.flood_model import bg_flood_model, process_hydro_dem
 
 DEFAULT_MODULES_TO_PARAMETERS = {
-    retrieve_from_instructions: {
-        "log_level": LogLevel.INFO,
-        "instruction_json_path": pathlib.Path("floodresilience/static_boundary_instructions.json").as_posix()
-    },
-    process_hydro_dem: {
-        "log_level": LogLevel.INFO
-    },
-    main_rainfall: {
-        "rcp": 2.6,
-        "time_period": "2031-2050",
-        "ari": 100,
-        "storm_length_mins": 2880,
-        "time_to_peak_mins": 1440,
-        "increment_mins": 10,
-        "hyeto_method": HyetoMethod.ALT_BLOCK,
-        "input_type": RainInputType.UNIFORM,
-        "log_level": LogLevel.INFO
-    },
-    main_tide_slr: {
-        "tide_length_mins": 2880,
-        "time_to_peak_mins": 1440,
-        "interval_mins": 10,
-        "proj_year": 2030,
-        "confidence_level": "low",
-        "ssp_scenario": "SSP1-2.6",
-        "add_vlm": False,
-        "percentile": 50,
-        "log_level": LogLevel.INFO
-    },
+    # retrieve_from_instructions: {
+    #     "log_level": LogLevel.INFO,
+    #     "instruction_json_path": pathlib.Path("floodresilience/static_boundary_instructions.json").as_posix()
+    # },
+    # process_hydro_dem: {
+    #     "log_level": LogLevel.INFO
+    # },
+    # main_rainfall: {
+    #     "rcp": 2.6,
+    #     "time_period": "2031-2050",
+    #     "ari": 100,
+    #     "storm_length_mins": 2880,
+    #     "time_to_peak_mins": 1440,
+    #     "increment_mins": 10,
+    #     "hyeto_method": HyetoMethod.ALT_BLOCK,
+    #     "input_type": RainInputType.UNIFORM,
+    #     "log_level": LogLevel.INFO
+    # },
+    # main_tide_slr: {
+    #     "tide_length_mins": 2880,
+    #     "time_to_peak_mins": 1440,
+    #     "interval_mins": 10,
+    #     "proj_year": 2030,
+    #     "confidence_level": "low",
+    #     "ssp_scenario": "SSP1-2.6",
+    #     "add_vlm": False,
+    #     "percentile": 50,
+    #     "log_level": LogLevel.INFO
+    # },
     main_river: {
         "flow_length_mins": 2880,
         "time_to_peak_mins": 1440,
diff --git a/src/tasks.py b/src/tasks.py
index 85c8718aa..ebe553c22 100644
--- a/src/tasks.py
+++ b/src/tasks.py
@@ -117,3 +117,4 @@ def wkt_to_gdf(wkt: str) -> gpd.GeoDataFrame:
 
 # These must be imported after app to remove a circular dependency
 import floodresilience.tasks  # pylint: disable=wrong-import-position,unused-import # noqa: E402, F401
+import wrfhydro.tasks  # pylint: disable=wrong-import-position,unused-import # noqa: E402, F401
diff --git a/wflow/.gitignore b/wflow/.gitignore
new file mode 100644
index 000000000..4d2e05d3f
--- /dev/null
+++ b/wflow/.gitignore
@@ -0,0 +1 @@
+**/hydro_mt_files/*
diff --git a/wflow/Convert_nz_to_GLOBCOVER.py b/wflow/Convert_nz_to_GLOBCOVER.py
new file mode 100644
index 000000000..eb8130a43
--- /dev/null
+++ b/wflow/Convert_nz_to_GLOBCOVER.py
@@ -0,0 +1,205 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Mon Nov 24 12:51:55 2025
+
+@author: mng42
+"""
+
+import geopandas as gpd
+
+nz_landcover = gpd.read_file(fr"S:\FloodRiskResearch\Martin\WRF-Hydro\landcover\landv6\lris-lcdb-v60-land-cover-database-version-60-mainland-new-zealand-SHP\lcdb-v60-land-cover-database-version-60-mainland-new-zealand.shp")
+
+nz_to_globcover = {
+    0: 210,   # Not land -> water
+
+    # Urban & infrastructure
+    1: 190,
+    2: 140,
+    5: 190,
+    6: 200,
+    10: 200,
+    12: 200,
+    16: 200,
+
+    # Snow & alpine
+    14: 220,
+    15: 150,
+
+    # Water
+    20: 210,
+    21: 210,
+    22: 210,
+
+    # Croplands
+    30: 14,
+    33: 20,
+
+    # Grasslands
+    40: 140,
+    41: 140,
+    43: 140,
+    44: 140,
+
+    # Wetlands
+    45: 120,
+    46: 180,
+    47: 120,
+
+    # Shrublands
+    50: 130,
+    51: 130,
+    52: 130,
+    54: 130,
+    55: 130,
+    56: 130,
+    58: 130,
+    80: 130,
+    81: 130,
+
+    # Mangrove
+    70: 180,
+
+    # Forests
+    68: 60,
+    69: 40,
+    71: 70,
+    64: 200
+}
+
+# Map NZ LCDB classes to GlobCover
+nz_landcover['GlobCover_2023'] = nz_landcover['Class_2023'].map(nz_to_globcover)
+
+# Check
+nz_landcover[['Class_2023', 'GlobCover_2023']].head(10)
+nz_landcover_4326 = nz_landcover.to_crs(epsg=4326)
+nz_landcover_4326['GlobCover_2023'] = nz_landcover_4326['GlobCover_2023'].astype('int32')
+
+###############################################################################
+
+import xarray as xr
+import rasterio
+from rasterio.features import rasterize
+from affine import Affine
+import numpy as np
+
+# Use the 4326 GeoDataFrame
+gdf = nz_landcover_4326
+
+# Define raster resolution in degrees
+pixel_size = 0.000595  # ~50m, 100m: 0.0009
+minx, miny, maxx, maxy = gdf.total_bounds
+
+# Raster size
+width = int((maxx - minx) / pixel_size)
+height = int((maxy - miny) / pixel_size)
+
+# Affine transform
+transform = Affine.translation(minx, maxy) * Affine.scale(pixel_size, -pixel_size)
+
+shapes = ((geom, value) for geom, value in zip(nz_landcover_4326.geometry, nz_landcover_4326['GlobCover_2023']))
+out_shape = (height, width)  # your raster dimensions
+transform = transform   # rasterio transform
+
+raster = rasterize(
+    shapes,
+    out_shape=out_shape,
+    fill=0,
+    transform=transform,
+    dtype=np.int32
+)
+
+with rasterio.open(
+    r"H:\NZ_landcover_landuse\nz_globcover_4326_50m_003.tif",
+    "w",
+    driver="GTiff",
+    height=raster.shape[0],
+    width=raster.shape[1],
+    count=1,
+    dtype=raster.dtype,  # int32
+    crs="EPSG:4326",
+    transform=transform,
+    nodata=0,            # <-- set nodata explicitly
+) as dst:
+    dst.write(raster, 1)
+
+
+
+# Check
+with rasterio.open(r"H:\NZ_landcover_landuse\nz_globcover_4326_50m_002.tif") as src:
+    print(src.dtypes)       # should be int32
+    print(src.nodata)       # should be 0
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+###############################################################################
+
+
+import xarray as xr
+import rasterio
+from rasterio.features import rasterize
+from affine import Affine
+import numpy as np
+
+# Use the 4326 GeoDataFrame
+gdf = nz_landcover_4326
+
+# Define raster resolution in degrees
+pixel_size = 0.000595  # ~50m, 100m: 0.0009
+minx, miny, maxx, maxy = gdf.total_bounds
+
+# Raster size
+width = int((maxx - minx) / pixel_size)
+height = int((maxy - miny) / pixel_size)
+
+# Affine transform
+transform = Affine.translation(minx, maxy) * Affine.scale(pixel_size, -pixel_size)
+
+# Rasterize
+shapes = ((geom, value) for geom, value in zip(gdf.geometry, gdf['GlobCover_2023']))
+raster = rasterize(
+    shapes=shapes,
+    out_shape=(height, width),
+    fill=np.nan,  # missing areas as NaN
+    transform=transform,
+    dtype='float32'
+)
+
+# Create x and y coordinates (center of pixels)
+x = np.linspace(minx + pixel_size/2, maxx - pixel_size/2, width)
+y = np.linspace(maxy - pixel_size/2, miny + pixel_size/2, height)
+
+# Create DataArray with band dimension
+da = xr.DataArray(
+    raster[np.newaxis, :, :],  # add band dimension
+    dims=('band', 'y', 'x'),
+    coords={'band': [1], 'x': x, 'y': y},
+    name='GlobCover_2023',
+    attrs={'spatial_ref': 4326, 'AREA_OR_POINT': 'Area'}
+)
+
+# Save as GeoTIFF
+da.rio.to_raster(r"H:\NZ_landcover_landuse\nz_globcover_4326_50m_002.tif")
+
+print("GeoTIFF saved successfully!")
+
+
+
+###############################################################################
+
+ex = xr.open_dataarray(fr"C:\Users\mng42\.hydromt_data\artifact_data\v0.0.9\data.tar\globcover.tif")
+
+
+
diff --git a/wflow/Generate_terrain_attributes.py b/wflow/Generate_terrain_attributes.py
new file mode 100644
index 000000000..db2e6701e
--- /dev/null
+++ b/wflow/Generate_terrain_attributes.py
@@ -0,0 +1,977 @@
+from whitebox_workflows import PhotometricInterpretation, RasterDataType, WbEnvironment
+from whitebox.whitebox_tools import WhiteboxTools
+import math
+from pathlib import Path
+
+import xarray as xr
+import rioxarray as rxr
+import numpy as np
+from scipy.ndimage import distance_transform_edt
+
+import geopandas as gpd
+import pandas as pd
+
+from rasterstats import zonal_stats
+
+from typing import Any
+
+wbe = WbEnvironment()
+wbe.verbose = True
+wbe.max_procs = -1
+
+wbt = WhiteboxTools()
+
+class CommonVariable:
+    def __init__(
+            self,
+            path: str
+    ):
+        """
+        A class contains common variables
+
+        Parameters
+        -----------
+        path : str
+            Path to raster that needs manipulating
+        """
+        self.path = path
+
+
+class TerrainAttributesGenerator(CommonVariable):
+    def __init__(
+            self,
+            path: str,
+            raster_name: str = 'dem'
+    ):
+        """
+        A class to generate terrain attributes
+
+        Parameters
+        ----------
+        path : str
+            Path to raster that needs manipulating
+        raster_name : str = 'dem'
+            Name of the raster. Mostly 'dem' and 'roughness'
+        """
+        super().__init__(path)
+        self.raster_name = raster_name
+
+
+    def raster_resampling(
+            self,
+            resolution_crs_4326 : float = 0.00045,
+            resampling_method : str = 'nn'
+    ) -> None:
+        """
+        Resample raster to a specific resolution (good with GeoTiff file)
+
+        Parameters
+        -----------
+        resolution_crs_4326 : float = 0.00045
+            Resolution value in crs 4326. Default is 0.00045 (~100 m)
+        resampling_method : str = 'nn'
+            Resampling methods includes "nn" (nearest neighbor), 'bilinear',
+            and 'cc' (cubic convolution). Default is 'nn'
+        """
+        if not Path(fr"{self.path}\{self.raster_name}_for_wflow_coarser.tif").is_file():
+            # Resample raster
+            wbt.resample(
+                inputs=fr"{self.path}\{self.raster_name}_for_wflow.tif",
+                output=fr"{self.path}\{self.raster_name}_for_wflow_coarser.tif",
+                cell_size=resolution_crs_4326,
+                method=resampling_method
+            )
+        else:
+            print(f"'{self.raster_name}_for_wflow_coarser.tif' already exists!")
+
+
+    def raster_fill_depression(
+            self,
+            flat_increment : float = 0.0001
+    ) -> None:
+        """
+        Fill depressions in raster (specifically in DEM)
+
+        Parameters
+        ----------
+        flat_increment : float = 0.0001
+            If flat surfaces suck as lakes have the slope 0 it will act like a sink.
+            This parameter will set a small slope to flat areas. Default is 0.0001.
+            https://github.com/williamlidberg/Whitebox-tutorial/blob/main/streams.py
+        """
+        if not Path(fr"{self.path}\{self.raster_name}_for_wflow_coarser_nodeps_crs.tif").is_file():
+            # Read the raster using whitebox tool
+            raster_no_deps = wbe.read_raster(
+                fr"{self.path}\{self.raster_name}_for_wflow_coarser.tif"
+            )
+
+            # Fill depressions in the raster
+            raster_no_deps = wbe.fill_depressions(
+                raster_no_deps,
+                flat_increment=flat_increment
+            )
+
+            # Write out
+            wbe.write_raster(
+                raster_no_deps,
+                fr"{self.path}\{self.raster_name}_for_wflow_coarser_nodeps.tif",
+                compress=False
+            )
+
+            # Read data using rioxarray to add crs and then overwrite out
+            raster_no_deps_crs = rxr.open_rasterio(fr"{self.path}\{self.raster_name}_for_wflow_coarser_nodeps.tif")
+            raster_no_deps_crs = raster_no_deps_crs.rio.write_crs('EPSG:4326')
+            raster_no_deps_crs.rio.to_raster(fr"{self.path}\{self.raster_name}_for_wflow_coarser_nodeps_crs.tif")
+
+
+
+        else:
+            print(f"'{self.raster_name}_for_wflow_coarser_nodeps_crs.tif' already exists!")
+
+    def d8_pointer_generator(self) -> None:
+        """
+        Generate D8 pointers based on D8 algorithm (O'Callaghan and Mark, 1984) (mainly from DEM)
+        https://www.whiteboxgeo.com/manual/wbw-user-manual/book/tool_help.html#d8_pointer
+        """
+        if not Path(fr"{self.path}\d8_pointer.tif").is_file():
+            # Read the raster using whitebox tool
+            raster_no_deps = wbe.read_raster(
+                fr"{self.path}\{self.raster_name}_for_wflow_coarser_nodeps_crs.tif"
+            )
+
+            # Generate D8 pointer
+            d8_pointer = wbe.d8_pointer(raster_no_deps)
+
+            # Write out raster
+            wbe.write_raster(
+                d8_pointer,
+                fr"{self.path}\d8_pointer.tif"
+            )
+        else:
+            print(f"'d8_pointer.tif' already exists!")
+
+
+    def d8_stream_generator(
+            self,
+            threshold: int = 25000,
+            catchment_area: bool = True
+    ) -> None:
+        """
+        Generate D8 flow accumulation based on D8 algorithm (O'Callaghan and Mark, 1984) (mainly from DEM)
+
+        Parameters
+        -----------
+        threshold : int = 10000
+            Minimum number of cells/upslope area required to initiate and main a channel.
+        catchment_area : bool = True
+            If True, flow accumulation under catchment are format will be added
+            If False, only flow accumulation under cell format will be generated
+        """
+        if not Path(fr"{self.path}\streams_d8.tif").is_file():
+            # Generate D8 flow accumulation - output is 'cell' type
+            wbt.d8_flow_accumulation(
+                i = fr"{self.path}\{self.raster_name}_for_wflow_coarser_nodeps_crs.tif",
+                out_type = 'cells',
+                output = fr'{self.path}\flow_acc_d8_cells.tif'
+            )
+
+            # Extract streams from the flow accumulation
+            wbt.extract_streams(
+                flow_accum = fr"{self.path}\flow_acc_d8_cells.tif",
+                output = fr"{self.path}\streams_d8.tif",
+                threshold = threshold
+            )
+        else:
+            print(f"'streams_d8.tif' already exists!")
+
+        if catchment_area:
+            if not Path(fr'{self.path}\flow_acc_d8_area_m2.tif').is_file():
+                # Generate D8 flow accumulation - output is 'catchment area' type
+                wbt.d8_flow_accumulation(
+                    i=fr"{self.path}\{self.raster_name}_for_wflow_coarser_nodeps_crs.tif",
+                    out_type="catchment area",
+                    output=fr'{self.path}\flow_acc_d8_area_m2.tif'
+                )
+            else:
+                print(fr"'flow_acc_d8_area_m2.tif' already exists!")
+        else:
+            pass
+
+
+    def strahler_stream_order_generator(self) -> None:
+        """
+        Generate Strahler stream order based on Strahler algorithm (Strahler, A. N., 1957)
+        and convert to vector
+        https://www.whiteboxgeo.com/manual/wbt_book/available_tools/stream_network_analysis.html#strahlerstreamorder
+        https://www.whiteboxgeo.com/manual/wbt_book/available_tools/stream_network_analysis.html?highlight=extract%20stream#RasterStreamsToVector
+        """
+        # Generate stream order
+        wbt.strahler_stream_order(
+            d8_pntr=fr"{self.path}\d8_pointer.tif",
+            streams=fr"{self.path}\streams_d8.tif",
+            output=fr"{self.path}\strahler_d8.tif"
+        )
+
+        # Convert stream raster to vector shapefile
+        wbt.raster_streams_to_vector(
+            d8_pntr=fr"{self.path}\d8_pointer.tif",
+            streams=fr"{self.path}\streams_d8.tif",
+            output=fr"{self.path}\streams_d8.shp"
+        )
+
+
+    def raster_to_points_dataframe(
+            self,
+            file_name: str,
+            output_name: str,
+            column_name: str
+    ) -> gpd.GeoDataFrame:
+        """
+        Convert raster values of pixels in stream network (could be area or strahler) to points
+
+        Parameters
+        -----------
+        file_name : str
+            Name of the file that will be used to convert to points
+        output_name : str
+            Name of output file that contains point shape type
+        column_name : str
+            Name of column that is converted from 'VALUE1' to
+            Options are mostly 'upstream_area_m2' and 'strahler'
+
+        Returns
+        ---------
+        points_df : gpd.GeoDataFrame
+            A GeoDataFrame contains points data of 'upstream_area_m2' or 'strahler'
+        """
+        if file_name != 'strahler_d8':
+            # Convert raster to point shapefile
+            wbt.raster_to_vector_points(
+                i=fr"{self.path}\streams_d8.tif",
+                output=fr"{self.path}\stream_pixels_pts.shp"
+            )
+        else:
+            pass
+
+        # Convert raster to vector of point shape type from flow accumulation
+        wbt.extract_raster_values_at_points(
+            fr"{self.path}\{file_name}.tif",
+            points=fr"{self.path}\stream_pixels_pts.shp"
+        )
+
+        # Convert to geopandas dataframe
+        points_df = gpd.read_file(fr"{self.path}\stream_pixels_pts.shp")
+        points_df = points_df.rename(columns={'VALUE1': f'{column_name}'})
+
+        return points_df
+
+
+    def roughness_to_manning(
+            self,
+            roughness: Any,
+            h: float = 1
+    ):
+        """
+
+        Parameters
+        ----------
+        roughness : Any
+            A raster of roughness data
+        h : float = 1
+            Value of depth. Default is 1
+
+        Returns
+        -------
+
+        """
+        # Convert roughness length to Manning's n
+        manning_n = (0.41 * (h ** (1 / 6)) * ((h / roughness) - 1)) / (np.sqrt(9.80665) * (1 + (h / roughness) * (np.log(h / roughness) - 1)))
+
+        # Write out Manning's n
+        manning_n.rio.to_raster(
+            fr"{self.path}\streams_manning.tif"
+        )
+
+
+class StreamTopologyGenerator(CommonVariable):
+    def __init__(
+            self,
+            path: str
+    ):
+        """
+        A class that generates stream topology: 'upstream_area_m2' and 'strahler'
+
+        Parameters
+        ----------
+        path : str
+            Path to raster that needs manipulating
+        """
+        super().__init__(path)
+        self.stream_topology_data = TerrainAttributesGenerator(self.path, 'dem')
+
+
+    def merge_stream_topology_points(self) -> gpd.GeoDataFrame:
+        """
+        Generate point dataframe that merges upstream catchment area and strahler order.
+        This dataframe already removes the rows with same FIDs.
+
+        Returns
+        --------
+        agg_upstream_area_and_strahler : gpd.GeoDataFrame
+            A geopandas dataframe that contains both 'upstream_area_m2' and 'strahler'.
+            There is no rows with the same FIDs
+        """
+        # Generate points that include upstream area under geopandas dataframe
+        points_area_m2 = self.stream_topology_data.raster_to_points_dataframe(
+            'flow_acc_d8_area_m2',
+            'stream_pixels_pts_from_flow_acc',
+            'upstream_area_m2'
+        )
+
+        # Generate points that include strahler order under geopandas dataframe
+        points_strahler_order = self.stream_topology_data.raster_to_points_dataframe(
+            'strahler_d8',
+            'stream_pixels_pts_from_strahler',
+            'strahler'
+        )
+
+        # Merge two points geopandas dataframe
+        points_merge = points_area_m2.merge(
+            points_strahler_order[['FID', 'strahler']],
+            on='FID',
+            how='left'
+        )
+
+        # After raster-to-point conversions and merging, duplicate FIDs are produced
+        # which leads to multiple rows with the same FID in the merged dataframe.
+        # Hence, the agg function is used to select the max values for both
+        # 'upstream_area_m2' and 'strahler'.
+        # - For 'upstream_area_m2', as flow accumulation
+        # increases downstream, the largest upstream area corresponds to the true
+        # accumulate catchment area at that point. Hence, the maximum upstream area
+        # preserves the most hydrologically meaningful value.
+        # - For 'strahler', as strahler order increases when tributatires merge,
+        # the maximum order ensures the feature keeps the correct stream hierarchy classification.
+        agg_upstream_area_and_strahler = (
+            points_merge.groupby('FID')
+                .agg(
+                    upstream_area = ('upstream_area_m2', 'max'),
+                    strahler = ('strahler', 'max')
+            )
+            .reset_index()
+        )
+
+        return agg_upstream_area_and_strahler
+
+
+    def merge_upstream_area_strahler_stream_geometry(
+            self,
+            agg_upstream_area_and_strahler
+    ):
+        """
+        Merge dataframes of 'upstream_area_m2' and 'strahler' with 'geometry' of stream using FID
+        and write out the merged dataframe
+        """
+        # Read D8 stream dataframe
+        streams = gpd.read_file(fr"{self.path}\streams_d8.shp")
+
+        # Merge the aggregated dataframe with stream dataframe
+        streams = streams.merge(
+            agg_upstream_area_and_strahler,
+            on='FID',
+            how='left'
+        )
+
+        # Convert from km2 to m2
+        streams['upstream_area'] = streams['upstream_area'] * 1e6
+
+        # Rename columns
+        streams_rename = streams.rename(
+            columns={
+                'upstream_area': 'uparea',
+                'strahler': 'strord'
+            }
+        )
+
+        # Add crs
+        streams_rename = streams_rename.set_crs(4326)
+
+        # Write out
+        streams_rename.to_file(
+            fr"{self.path}\streams_d8_area_strahler.shp"
+        )
+
+
+    def dataframe_upstream_area_strahler_geometry_generator(
+            self,
+    ) -> None:
+        """
+        Generate geodataframe of 'upstream_area_m2' and 'strahler'
+        """
+        # Resample raster
+        self.stream_topology_data.raster_resampling(
+            0.00045,
+            'nn'
+        )
+
+        # Fill depression
+        self.stream_topology_data.raster_fill_depression(0.0001)
+
+        # Generate D8 pointer
+        self.stream_topology_data.d8_pointer_generator()
+
+        # Generate D8 stream generator
+        self.stream_topology_data.d8_stream_generator(
+            25000,
+            True
+        )
+
+        # Generate strahler order
+        self.stream_topology_data.strahler_stream_order_generator()
+
+        # Collect dataframe of 'upstream_area_m2' and 'strahler'
+        df_upstream_area_strahler = self.merge_stream_topology_points()
+
+        # Merge dataframe of 'upstream_area_m2' and 'strahler'
+        # with dataframe of stream geometry and write out
+        self.merge_upstream_area_strahler_stream_geometry(df_upstream_area_strahler)
+
+
+class StreamHydraulicsGenerator(CommonVariable):
+    def __init__(
+            self,
+            path: str,
+            outlet_gauge_locations_file: str,
+            streams_bankfull_stage: float = 1.5
+    ):
+        """
+
+        Parameters
+        -----------
+        path : str
+            Path to raster that needs manipulating
+        outlet_gauge_locations_file : str
+            Filename that contains locations of outlet and gauges
+        streams_bankfull_stage : float = 1.5
+            The stage to focus on the area that is considered as stream/river area
+            or bankfull area comparing with HAND.
+            Default is 1.5
+        """
+        super().__init__(path)
+        self.streams_bankfull_stage = streams_bankfull_stage
+        self.outlet_gauge_locations_file = outlet_gauge_locations_file
+        self.stream_topology_data = TerrainAttributesGenerator(self.path, 'dem')
+        self.roughness_data = TerrainAttributesGenerator(self.path, 'roughness')
+
+
+    def watershed_generator(
+            self,
+            snap_dist: float = 5.0,
+            filter_size: int = 5
+
+    ) -> Any:
+        """
+        Generate watershed based on D8 pointer (flow direction) and list of points of outlet and gauges
+        
+        Parameters
+        -----------
+        snap_dist : float = 5.0
+            Measures in map units (e.g. meters, default is meters) the given maximum distance
+            between the pour points to the location coincident with the nearest stream cell.
+            Default is 5 meters.
+            https://www.whiteboxgeo.com/manual/wbt_book/available_tools/hydrological_analysis.html#JensonSnapPourPoints
+        filter_size : int = 5
+            Filter size to smooth a vector coverage of either a Polyline or Polygon base.
+            It can be any integer larger than or equal to 3. Default here is 5.
+            https://www.whiteboxgeo.com/manual/wbw-user-manual/book/tool_help.html#smooth_vectors
+
+        Returns
+        -------
+        watershed_polygon : Any
+            Watershed that leads to the outlet
+        """
+
+        # Read stream raster
+        streams = wbe.read_raster(fr"{self.path}\streams_d8.tif")
+
+        # Read D8 pointer
+        d8_pointer = wbe.read_raster(fr"{self.path}\d8_pointer.tif")
+
+        # Extract watershed for specific points of outlet and gauges
+        outlet_gauge_points = wbe.read_vector(fr"{self.path}\{self.outlet_gauge_locations_file}.shp")
+
+        # Ensure the watershed or streamlines that have points of outlet and gauges
+        outlet_gauge_points_on_streams = wbe.jenson_snap_pour_points(
+            outlet_gauge_points,
+            streams,
+            snap_dist = snap_dist
+        )
+
+        # Extract watershed of the outlet
+        outlet_watershed = wbe.watershed(
+            d8_pointer = d8_pointer,
+            pour_points = outlet_gauge_points_on_streams
+        )
+
+        # Write out watershed raster
+        wbe.write_raster(
+            outlet_watershed,
+            fr"{self.path}\watershed.tif",
+            compress=False
+        )
+
+        # Generate watershed polygon for checking (if necessary)
+        watershed_polygon = wbe.raster_to_vector_polygons(outlet_watershed)
+
+        # Smooth the watershed map
+        watershed_polygon = wbe.smooth_vectors(
+            watershed_polygon,
+            filter_size=filter_size
+        )
+
+        # Write out
+        wbe.write_vector(
+            watershed_polygon,
+            fr"{self.path}\watershed.shp"
+        )
+
+        return outlet_watershed
+
+
+    def stream_watershed_raster_generator(
+            self,
+            watershed_raster: Any
+    ):
+        """
+
+        Returns
+        -------
+
+        """
+
+        # Read stream raster
+        streams = wbe.read_raster(fr"{self.path}\streams_d8.tif")
+
+        # Read D8 pointer raster
+        d8_pointer = wbe.read_raster(fr"{self.path}\d8_pointer.tif")
+
+        # Read DEM that its depressions are filled
+        dem_no_deps = wbe.read_raster(fr"{self.path}\dem_for_wflow_coarser_nodeps_crs.tif")
+
+        # Filter to select only streams inside the watershed
+        streams_watershed = streams * watershed_raster
+
+        # Write out raster with streams within watershed
+        # (this stream data just has 1 and 0 values)
+        wbe.write_raster(
+            streams_watershed,
+            fr"{self.path}\streams_watershed.tif",
+            compress=False
+        )
+
+        # Convert stream raster within watershed to vector (just geometry)
+        streams_watershed_vector = wbe.raster_streams_to_vector(
+            streams_watershed,
+            d8_pointer
+        )
+
+        # Add more information into stream vector such as
+        # reach IDs or FIDs, connectivity, stream orders, flow connectivity information, etc.
+        streams_watershed_vector_more_info, _, _, _ = wbe.vector_stream_network_analysis(
+            streams_watershed_vector,
+            dem_no_deps
+        )
+
+        # Write out vector with streams with watershed
+        # (this stream data has more information like FIDs, connectivity, etc.)
+        wbe.write_vector(
+            streams_watershed_vector_more_info,
+            fr"{self.path}\streams_watershed_more_info.shp"
+        )
+
+        # Convert back to raster once collecting FID
+        streams_watershed_raster = wbe.vector_lines_to_raster(
+            streams_watershed_vector_more_info,
+            'FID',
+            base_raster = dem_no_deps,
+            zero_background = True
+        )
+
+        return streams_watershed, streams_watershed_raster
+
+
+    def hand_generator(self):
+        """
+
+        Returns
+        -------
+
+        """
+        # Read streams within watershed
+        streams_watershed = wbe.read_raster(fr"{self.path}\streams_watershed.tif")
+
+        # Read DEM that its depressions are filled
+        dem_no_deps = wbe.read_raster(fr"{self.path}\dem_for_wflow_coarser_nodeps_crs.tif")
+
+        # Calculate HAND
+        hand = wbe.elevation_above_stream(
+            dem_no_deps,
+            streams_watershed
+        )
+
+        # Write out
+        wbe.write_raster(
+            hand,
+            fr"{self.path}\hand.tif",
+            compress=False
+        )
+
+    def stream_bankfull_width_raster_generator(self):
+        """
+
+        Parameters
+        -----------
+
+        Returns
+        -------
+
+        """
+        # Read streams within watershed using rioxarray
+        streams_watershed = rxr.open_rasterio(fr"{self.path}\streams_watershed.tif").squeeze()
+
+        # Read HAND raster
+        hand = rxr.open_rasterio(fr"{self.path}\hand.tif").squeeze()
+
+        # Set up bankfull
+        bankfull = (hand <= self.streams_bankfull_stage)
+
+        # Get values of bankfull and stream
+        # bankfull_np tells us where the river is
+        # stream_np tells us where the streamline is
+        bankfull_np = bankfull.values  # Raster where river area = 1 (True) and land = 0 (False)
+        stream_np = streams_watershed.values > 0  # A mask of stream centreline pixels (True = this pixel is part of the stream)
+
+        # Calculate distance
+        # Get pixel size: Each pixel represents X meters on the ground
+        pixel_size = abs(hand.rio.resolution()[0])
+
+        # Convert degrees to meters
+        lat = float(bankfull.y.mean())  # Get latitude
+        meters_per_degree = (
+                111300 * np.cos(np.deg2rad(lat))
+        )
+        pixel_size_meter = pixel_size * meters_per_degree
+
+        # Measure distance to the river bank:
+        # For every river pixel, it calculates the distance to the nearest non-river pixel (the river bank)
+        # It here is the function "distance_transform_edt"
+        # So pixels near the bank --> small distahce
+        # Pixels near the center of the river --> larger distance
+        # For example: bank 1m 2m 3m 4m 3m 2m 1m bank (so near the bank --> small distance, near the middle --> larger distance)
+        # ==> Logic is: How far am I from the river edge?
+        distance = distance_transform_edt(bankfull_np) * pixel_size_meter
+
+        # Filter out only the distance from the center line to a bank
+        # So it will be like: bank Nan Nan 4m Nan Nan bank
+        # And then double it for the other bank
+        bankfull_width = 2 * distance[stream_np]
+
+        # Put widths back into a raster
+        streams_bankfull_width = np.full(stream_np.shape, np.nan)
+        streams_bankfull_width[stream_np] = bankfull_width
+
+        # Convert numpy array to xarray data array
+        streams_bankfull_width_da = xr.DataArray(
+            streams_bankfull_width,
+            coords = hand.coords,
+            dims = hand.dims,
+            name = "bankfull_width"
+        )
+
+        # Add crs
+        streams_bankfull_width_da.rio.write_crs(
+            hand.rio.crs,
+            inplace = True
+        )
+
+        # Write out
+        streams_bankfull_width_da.rio.to_raster(
+            fr"{self.path}\streams_bankfull_width.tif"
+        )
+
+
+    def buffer_streams_watershed(self):
+        """
+        Buffer the stream linestrings to capture stream pixels
+
+        Returns
+        -------
+
+        """
+        # Read streams within watershed using geopandas
+        streams_watershed_vector_more_info = gpd.read_file(
+            fr"{self.path}\streams_watershed_more_info.shp"
+        )
+
+        # Read HAND raster
+        hand = rxr.open_rasterio(fr"{self.path}\hand.tif").squeeze()
+
+        # Buffer by half raster pixel
+        # To explain, linestrings have no width and might not be alligned well with the raster.
+        # Hence, buffering helps to ensure the stream intersects well with the raster
+        pixel_size = abs(hand.rio.resolution()[0])
+        streams_watershed_buffer = streams_watershed_vector_more_info.copy()
+        streams_watershed_buffer['geometry'] = streams_watershed_buffer.geometry.buffer(pixel_size / 2)
+
+        return streams_watershed_buffer, streams_watershed_vector_more_info
+
+
+    def assign_streams_hydraulic_values(
+            self,
+            streams_watershed_buffer,
+            streams_watershed_vector_more_info,
+            hydraulic_name
+    ):
+        """
+        Assign stream hydraulic values to each stream segment in the vector network
+
+        Returns
+        -------
+
+        """
+        # Extract stream hydraulic values for each reach
+        streams_hydraulic_values = zonal_stats(
+            vectors = streams_watershed_buffer,
+            raster = fr"{self.path}\streams_{hydraulic_name}.tif",
+            stats = ['mean'],
+            all_touched = True # includes all pixels by touching the buffered ones
+        )
+
+        # Convert list of dicts to dataframe
+        streams_hydraulic_values_df = pd.DataFrame(streams_hydraulic_values)
+
+        # Add stream hydraulic values to streams watershed dataframe
+        streams_hydraulic_values_linestring = streams_watershed_vector_more_info.join(
+            streams_hydraulic_values_df
+        )
+        streams_hydraulic_values_linestring[f'{hydraulic_name}'] = streams_hydraulic_values_linestring['mean']
+
+        # Write out
+        streams_hydraulic_values_linestring.to_file(
+            fr"{self.path}\streams_{hydraulic_name}_linestring.shp"
+        )
+
+
+    def stream_bankfull_width_linestring_generator(self):
+        """
+
+        Returns
+        -------
+
+        """
+        # Buffer stream linestring to capture stream pixels
+        streams_watershed_buffer, streams_watershed_vector_more_info = self.buffer_streams_watershed()
+
+        # Assign stream bankfull width to stream linestring
+        self.assign_streams_hydraulic_values(
+            streams_watershed_buffer,
+            streams_watershed_vector_more_info,
+            'bankfull_width'
+        )
+
+
+    def stream_manning_linestring_generator(self):
+        """
+
+        Returns
+        -------
+
+        """
+        # Resample roughness raster
+        self.roughness_data.raster_resampling(
+            0.00045,
+            'nn'
+        )
+
+        # Read coarse roughness raster
+        roughness_for_wflow_coarser = rxr.open_rasterio(
+            fr"{self.path}\roughness_for_wflow_coarser.tif"
+        )
+
+        # Convert roughness to manning
+        self.roughness_data.roughness_to_manning(
+            roughness_for_wflow_coarser,
+            1
+        )
+
+        # Buffer stream linestring to capture stream pixels
+        streams_watershed_buffer, streams_watershed_vector_more_info = self.buffer_streams_watershed()
+
+        # Assign stream bankfull width to stream linestring
+        self.assign_streams_hydraulic_values(
+            streams_watershed_buffer,
+            streams_watershed_vector_more_info,
+            'manning'
+        )
+
+
+    def stream_slope_linestring_generator(self):
+        """
+
+        Returns
+        -------
+
+        """
+        # Read DEM that its depressions are filled
+        dem_no_deps = wbe.read_raster(fr"{self.path}\dem_for_wflow_coarser_nodeps_crs.tif")
+
+        # Generate slope from DEM
+        slope = wbe.slope(
+            dem_no_deps,
+            units='percent'
+        )
+
+        # Write out slope
+        wbe.write_raster(
+            slope,
+            fr"{self.path}\streams_slope.tif",
+            compress=False
+        )
+
+        # Buffer stream linestring to capture stream pixels
+        streams_watershed_buffer, streams_watershed_vector_more_info = self.buffer_streams_watershed()
+
+        # Assign stream bankfull width to stream linestring
+        self.assign_streams_hydraulic_values(
+            streams_watershed_buffer,
+            streams_watershed_vector_more_info,
+            'slope'
+        )
+
+    def bankfull_discharge_calculation(
+            self,
+            streams_bankfull_width,
+            streams_slope,
+            streams_manning
+    ):
+        """
+
+        Returns
+        -------
+
+        """
+        # Calculate cross-sectional area (rectangular)
+        cross_sectional_area = streams_bankfull_width * self.streams_bankfull_stage
+
+        # Calculate wetted perimeter
+        wetted_perimeter = streams_bankfull_width + 2 * self.streams_bankfull_stage
+
+        # Calculate hydraulic radius
+        hydraulic_radius = cross_sectional_area / wetted_perimeter
+
+        # Calculate bankfull discharge
+        bankfull_discharge = (1 / streams_manning) * wetted_perimeter * (hydraulic_radius ** (2/3)) * np.sqrt(streams_slope)
+
+        return bankfull_discharge
+
+
+    def stream_bankfull_discharge_generator(self):
+        """
+
+        Returns
+        -------
+
+        """
+        # Read stream bankfull width shapefile
+        streams_bankfull_width = gpd.read_file(fr"{self.path}\streams_bankfull_width_linestring.shp").bankfull_w
+        streams_manning = gpd.read_file(fr"{self.path}\streams_manning_linestring.shp").manning
+        streams_slope = gpd.read_file(fr"{self.path}\streams_slope_linestring.shp").slope
+
+        # Generate stream bankfull discharge
+        streams_discharge = self.bankfull_discharge_calculation(
+            streams_bankfull_width,
+            streams_slope,
+            streams_manning
+        )
+
+        return streams_discharge
+
+
+    def stream_bankfull_width_discharge_generator(self):
+        """
+
+        Returns
+        -------
+
+        """
+        # Generate stream bankfull discharge
+        streams_bankfull_discharge = self.stream_bankfull_discharge_generator()
+
+        # Read stream bankfull width
+        streams_bankfull_width = gpd.read_file(
+            fr"{self.path}\streams_bankfull_width_linestring.shp"
+        )
+
+        # Rename stream bankfull width
+        streams_bankfull_width_discharge = streams_bankfull_width.rename(columns={'bankfull_w': 'rivwth'})
+
+        # Add stream bankfull discharge
+        streams_bankfull_width_discharge['qbankfull'] = streams_bankfull_discharge
+
+        # Write out
+        streams_bankfull_width_discharge.to_file(fr"{self.path}\streams_bankfull_width_discharge.gpkg")
+
+
+    def dataframe_stream_bankfull_width_discharge_generator(self):
+        """
+
+        Returns
+        -------
+
+        """
+        # Resample raster
+        self.stream_topology_data.raster_resampling(
+            0.00045,
+            'nn'
+        )
+
+        # Fill depression
+        self.stream_topology_data.raster_fill_depression(0.0001)
+
+        # Generate D8 pointer
+        self.stream_topology_data.d8_pointer_generator()
+
+        # Generate D8 stream generator
+        self.stream_topology_data.d8_stream_generator(
+            10000,
+            True
+        )
+
+        # Generate watershed raster
+        watershed_polygon = self.watershed_generator(
+            5, 5
+        )
+        self.stream_watershed_raster_generator(
+            watershed_polygon
+        )
+
+        # Generate HAND raster
+        self.hand_generator()
+
+        # Generate stream bankfull width raster
+        self.stream_bankfull_width_raster_generator()
+
+        # Buffer stream linestring to capture stream pixels
+        self.buffer_streams_watershed()
+
+        # Generate stream bankfull width linestring
+        self.stream_bankfull_width_linestring_generator()
+
+        # Generate stream manning
+        self.stream_manning_linestring_generator()
+
+        # Generate stream slope
+        self.stream_slope_linestring_generator()
+
+        # Write out stream bankfull width and discharge
+        self.stream_bankfull_width_discharge_generator()
\ No newline at end of file
diff --git a/wflow/Generate_terrain_data_for_wflow.py b/wflow/Generate_terrain_data_for_wflow.py
new file mode 100644
index 000000000..0b6b43d0a
--- /dev/null
+++ b/wflow/Generate_terrain_data_for_wflow.py
@@ -0,0 +1,267 @@
+import os
+import xarray as xr
+import numpy as np
+import geopandas as gpd
+
+import pandas as pd
+
+# pyflwdir
+import pyflwdir
+from PyQt5.QtGui import QVector2D
+
+# hydromt
+from hydromt import DataCatalog, flw
+
+from pathlib import Path
+
+# plot
+import matplotlib.pyplot as plt
+from matplotlib import cm, colors
+
+
+class TerrainDataWflow:
+    def __init__(
+            self,
+            path
+    ):
+        """
+        https://deltares.github.io/hydromt_wflow/stable/_examples/prepare_ldd.html
+        Parameters
+        ----------
+        path
+        """
+        self.path = path
+        self.ds_hydro_org = xr.open_dataarray(fr"{self.path}\dem_for_wflow_coarser_nodeps_crs.tif")
+        self.gdf_riv_org = gpd.read_file(fr"{self.path}\streams_d8_area_strahler.shp")
+
+        # Convert uparea column to numeric
+        self.gdf_riv_org['uparea'] = pd.to_numeric(
+            self.gdf_riv_org['uparea'], errors='coerce'
+        )
+
+
+    def d8_flow_direction_generator(self):
+        """
+
+        Returns
+        -------
+
+        """
+        # Derive flow directions with outlets at the edges
+        da_flwdir = flw.d8_from_dem(
+            da_elv=self.ds_hydro_org.squeeze(),
+            max_depth=-1,  # max depression poir point depth; -1 means no local pits
+            outlets="edge",  # option: "edge" (default), "min", "idxs_pit"
+            idxs_pit=None,
+            gdf_riv=self.gdf_riv_org,  # user supplied river network to aid flow direction derivation
+            riv_burn_method="uparea",  # options: "fixed" (default), "rivdph", "uparea"
+            # riv_depth=5, # fixed river depth in meters, only used if riv_burn_method="fixed"
+            # **kwargs to be passed to pyflwdir.dem.fill_depressions
+        )
+
+        # Convert to vector
+        flwdir = flw.flwdir_from_da(
+            da_flwdir,
+            ftype='infer',
+            check_ftype=True
+        )
+
+        # Create a new ds_hydro dataset with the riverburn flow directions
+        ds_hydro = da_flwdir.to_dataset(name='flwdir')
+        ds_hydro = ds_hydro.raster.gdal_compliant()
+
+        # Extract dimesions
+        dims = ds_hydro.raster.dims
+
+        return flwdir, ds_hydro, dims
+
+
+    def terrain_and_stream_generator(
+            self,
+            flwdir,
+            ds_hydro,
+            dims,
+            long_name,
+            units,
+            variable_name
+    ):
+        """
+
+        Returns
+        -------
+
+        """
+        terrain_data = flwdir.dem_adjust(
+            elevtn = self.ds_hydro_org.squeeze().values
+        )
+        attrs = dict(
+            _FillValue = -9999,
+            long_name = long_name,
+            units = units
+        )
+
+        ds_hydro[f'{variable_name}'] = xr.Variable(
+            dims,
+            terrain_data,
+            attrs=attrs
+        )
+
+        return ds_hydro
+
+
+    def slope_generator(
+            self,
+            ds_hydro,
+            dims
+    ):
+        """
+
+        Returns
+        -------
+
+        """
+        # Generate slope
+        slope = pyflwdir.dem.slope(
+            elevtn = ds_hydro['elevtn'].values,
+            nodata = ds_hydro['elevtn'].raster.nodata,
+            latlon = ds_hydro.raster.crs.is_geographic, # True if geographic crs, False if projected crs
+            transform = ds_hydro['elevtn'].raster.transform
+        )
+
+        # Generate attributes for the slope
+        attrs = dict(
+            _FillValue = -9999,
+            long_name = 'lndslp',
+            units = 'm/m'
+        )
+
+        # Add to hydro dataset
+        ds_hydro['lndslp'] = xr.Variable(
+            dims,
+            slope,
+            attrs=attrs
+        )
+
+        return ds_hydro
+
+
+    def basin_generator(
+            self,
+            flwdir,
+            ds_hydro,
+            dims
+    ):
+        """
+
+        Returns
+        -------
+
+        """
+        # Generate basin
+        basins = flwdir.basins(
+            idxs = flwdir.idxs_pit
+        ).astype(np.int32)
+
+        # Generate attributes
+        attrs = dict(
+            _FillValue = 0,
+            long_name = 'basin ids',
+            units='-'
+        )
+
+        # Add basins to the hydro dataset
+        ds_hydro['basins'] = xr.Variable(
+            dims,
+            basins,
+            attrs = attrs
+        )
+
+        # Add basins to geopandas dataframe
+        gdf_basins = ds_hydro['basins'].raster.vectorize()
+
+        return ds_hydro, gdf_basins
+
+
+    def write_out_hydro_dataset(
+            self,
+            ds_hydro,
+            gdf_basins
+    ):
+        # Write out hydro dataset
+        ds_hydro.raster.to_mapstack(
+            root = os.path.join(
+                self.path,
+                'merit_hydro'
+            ),
+            driver='GTiff'
+        )
+
+        # Rename basins column
+        gdf_basins = gdf_basins.rename(
+            columns = {'value': 'basid'}
+        )
+
+        # Write out basin dataframe
+        gdf_basins.to_file(
+            os.path.join(
+                self.path,
+                'da_hydro_basins.gpkg'
+            )
+        )
+
+
+    def terrain_data_for_wflow_generator(
+            self
+    ):
+
+        # Generate D8 flow direction
+        flwdir, ds_hydro, dims = self.d8_flow_direction_generator()
+
+        # Generate elevation data
+        ds_elevation = self.terrain_and_stream_generator(
+            flwdir,
+            ds_hydro,
+            dims,
+            'corrected elevation',
+            'm',
+            'elevtn'
+        )
+
+        # Generate upstream area
+        ds_upstream_area = self.terrain_and_stream_generator(
+            flwdir,
+            ds_elevation,
+            dims,
+            'upstream area',
+            'km2',
+            'uparea'
+        )
+
+        # Generate stream strahler order
+        ds_strord = self.terrain_and_stream_generator(
+            flwdir,
+            ds_upstream_area,
+            dims,
+            'stream order',
+            '-',
+            'strord'
+        )
+
+        # Generate slope data
+        ds_slope = self.slope_generator(
+            ds_strord,
+            dims
+        )
+
+        # Generate basins
+        ds_basins, gdf_basins = self.basin_generator(
+            flwdir,
+            ds_slope,
+            dims
+        )
+
+        # Write out hydro dataset
+        self.write_out_hydro_dataset(
+            ds_basins,
+            gdf_basins
+        )
\ No newline at end of file
diff --git a/wflow/Generate_wflow_files.py b/wflow/Generate_wflow_files.py
new file mode 100644
index 000000000..1af260646
--- /dev/null
+++ b/wflow/Generate_wflow_files.py
@@ -0,0 +1,37 @@
+import subprocess
+
+
+def generate_wflow_files(
+    folder_path,
+    subbasin,
+    strord,
+    bbox
+):
+    '''
+    Definition:
+        Function to generate necessary files for Wflow model
+    References:
+        https://deltares.github.io/hydromt_wflow/stable/_examples/build_sediment.html
+    Arguments:
+        folder_path (str):
+            Path to folder that contains wflow_build.yml, data_catalog.yml,
+            and a folder that contains all necessary files
+        subbasin (list):
+            Outlet coordinates
+        strord (int):
+            Minimum stream order
+        bbox (list):
+            Given bounding box coordinates that contains the subbasin coordinates
+    '''
+
+    cmd = [
+        "hydromt", "build", "wflow",
+        f"{folder_path}/wflow_test_full",
+        "-r", f"{f'subbasin': {subbasin}, 'strord': {strord}, 'bbox': {bbox}}",
+        "-i", f"{folder_path}/wflow_build.yml",
+        "-d", f"{folder_path}/data_catalog.yml",
+        "-vv"
+    ]
+
+    # Run and print live output
+    subprocess.run(cmd, check=True)
\ No newline at end of file
diff --git a/wflow/Manipulate_terrain_data.py b/wflow/Manipulate_terrain_data.py
new file mode 100644
index 000000000..b864b2923
--- /dev/null
+++ b/wflow/Manipulate_terrain_data.py
@@ -0,0 +1,136 @@
+import rioxarray as rxr
+import geopandas as gpd
+import os
+
+
+def value_change(
+        shapefile_func,
+        file_need_changing_func,
+        value_func,
+        inside=True
+):
+    """
+    A function to change pixel values inside or outside polygons
+
+    Parameters
+    ----------
+    shapefile_func: str
+        Path to shapefile that cover the area that needs changing
+    file_need_changing_func: str
+        Name and path of changed file
+    value_func: float
+        Replaced value
+    inside (boolean):
+        If True, change values inside, else, change values outside
+    """
+    # Set up value changing command
+    if inside:
+        # Change values inside polygons
+        inside_command = fr"gdal_rasterize -burn {value_func} {shapefile_func} {file_need_changing_func}"
+        os.system(inside_command)
+    else:
+        # Change values outside polygons
+        outside_command = fr"gdal_rasterize -i -burn {value_func} {shapefile_func} {file_need_changing_func}"
+        os.system(outside_command)
+
+
+class TerrainFilter:
+    def __init__(
+            self,
+            path: str,
+            origin_crs: int = 2193,
+            converted_crs: int = 4326,
+            roughness: bool = True,
+            sea_value: float = -9999,
+            nodata_value: float = -9999
+    ):
+        """
+        A class to filter terrain data for wflow
+
+        Parameters
+        -----------
+        path: str
+            Path to raster that needs manipulating
+        origin_crs : int = 2193
+            Original CRS (default is 2193)
+        converted_crs : int = 4326
+            Converted CRS (default is 4326)
+        roughness : bool
+            Whether to print out roughness and DEM. Default is True
+        sea_value: float = -9999
+            Value for the sea data. Default is -9999
+        nodata_value: float = -9999
+            Value for the nodata. Default is -9999
+        """
+        self.path = path
+        self.origin_crs = origin_crs
+        self.converted_crs = converted_crs
+        self.roughness = roughness
+        self.sea_value = sea_value
+        self.nodata_value = nodata_value
+
+    def dem_crs_conversion(self):
+        """
+        Convert DEM CRS. Here default is to convert DEM CRS from 2193 to 4326
+        """
+        # Get DEM
+        dem = rxr.open_rasterio(fr"{self.path}\8m_geofabric.nc")
+
+        # Ensure the origin CRS is attached
+        dem_origin_crs = dem.rio.write_crs(f"EPSG:{self.origin_crs}", inplace=True)
+
+        # Reproject the crs to converted crs
+        dem_converted_crs = dem_origin_crs.rio.reproject(f"EPSG:{self.converted_crs}")
+
+        # Save as tif
+        if self.roughness:
+            dem_converted_crs['z'].rio.to_raster(fr"{self.path}\dem_converted_crs.tif")
+            dem_converted_crs['zo'].rio.to_raster(fr"{self.path}\roughness_converted_crs.tif")
+        else:
+            dem_converted_crs.rio.to_raster(fr"{self.path}\dem_converted_crs.tif")
+
+
+    def remove_sea(self):
+        """
+        Clip sea area mainly in DEM and roughness
+        """
+        # Get New Zealand shapefile
+        nz_shapefile = fr"{self.path}\nz_coastline_4326.shp"
+
+        # Files need changing
+        dem = fr"{self.path}\dem_converted_crs.tif"
+        roughness = fr"{self.path}\roughness_converted_crs.tif"
+
+        # Remove sea by changing sea area into -9999
+        value_change(nz_shapefile, dem, self.sea_value, False)
+        value_change(nz_shapefile, roughness, self.sea_value, False)
+
+
+    def nodata_filling(self):
+        """
+        Fill nodata value with -9999
+        """
+        # Fill the nodata value
+        dem_nosea = rxr.open_rasterio(fr"{self.path}\dem_converted_crs.tif")
+        dem_replace_nodata = dem_nosea.fillna(self.nodata_value)
+        dem_write_nodata = dem_replace_nodata.rio.write_nodata(self.nodata_value)
+        dem_write_nodata.rio.to_raster(fr"{self.path}\dem_for_wflow.tif")
+
+        roughness_nosea = rxr.open_rasterio(fr"{self.path}\roughness_converted_crs.tif")
+        roughness_replace_nodata = roughness_nosea.fillna(self.nodata_value)
+        roughness_write_nodata = roughness_replace_nodata.rio.write_nodata(self.nodata_value)
+        roughness_write_nodata.rio.to_raster(fr"{self.path}\roughness_for_wflow.tif")
+
+
+    def filter_dem_for_wflow(self):
+        """
+        Convert DEM into version that can be used by wflow
+        """
+        # Convert DEM and roughness CRS (default: 2193 --> 4326)
+        self.dem_crs_conversion()
+
+        # Remove sea
+        self.remove_sea()
+
+        # Fill nodata
+        self.nodata_filling()
\ No newline at end of file
diff --git a/wflow/Run_wflow.py b/wflow/Run_wflow.py
new file mode 100644
index 000000000..4d89d6487
--- /dev/null
+++ b/wflow/Run_wflow.py
@@ -0,0 +1,31 @@
+import subprocess
+
+def run_wflow(
+    wflow_path,
+    num_threads
+):
+    '''
+    Definition:
+        Function to run wflow model
+    References:
+        https://deltares.github.io/Wflow.jl/v0.8/user_guide/additional_options/
+    Arguments:
+        wflow_path (str):
+            Path to folder that stores all necessary files to run wflow model
+        num_threads (int):
+            Number of threads that controls how fast the wflow model can run
+    '''
+
+    # Build the Julia command
+    cmd = [
+        "julia",
+        "-t", f"{num_threads}",
+        "-e",
+        f'cd("{wflow_path}"); using Wflow; Wflow.run("wflow_sbm.toml")'
+    ]
+
+    # Run the command and write output to log
+    with open(fr"{wflow_path}\wflow_run.log", "w") as f:
+        process = subprocess.run(cmd, stdout=f, stderr=subprocess.STDOUT, text=True)
+
+    print(f"Wflow run completed. Log saved to wflow_run.log")
diff --git a/wrfhydro/__init__.py b/wrfhydro/__init__.py
new file mode 100644
index 000000000..7f03ad5e5
--- /dev/null
+++ b/wrfhydro/__init__.py
@@ -0,0 +1,17 @@
+# Copyright © 2021-2025 Geospatial Research Institute Toi Hangarau
+# LICENSE: https://github.com/GeospatialResearch/Digital-Twins/blob/master/LICENSE
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Affero General Public License as
+# published by the Free Software Foundation, either version 3 of the
+# License, or (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU Affero General Public License for more details.
+#
+# You should have received a copy of the GNU Affero General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+"""Top-level package for the WRF-Hydro® extension for the Digital Twin."""
diff --git a/wrfhydro/blueprint.py b/wrfhydro/blueprint.py
new file mode 100644
index 000000000..d1da1dfc7
--- /dev/null
+++ b/wrfhydro/blueprint.py
@@ -0,0 +1,30 @@
+"""Endpoints and flask configuration for running WRF-Hydro module within the digital twin"""
+
+from flask import Blueprint
+from pywps import Service
+
+from wrfhydro.scenario.flood_scenario_process_service import FloodScenarioProcessService
+from src.check_celery_alive import check_celery_alive
+
+wrf_hydro_blueprint = Blueprint('wrfhydro', __name__)
+processes = [
+    FloodScenarioProcessService()
+]
+
+process_descriptor = {process.identifier: process.abstract for process in processes}
+
+service = Service(processes, ['src/pywps.cfg'])
+
+
+@wrf_hydro_blueprint.route('/wps', methods=['GET', 'POST'])
+@check_celery_alive
+def wps() -> Service:
+    """
+    End point for OGC WebProcessingService spec, allowing clients such as TerriaJS to request processing.
+
+    Returns
+    -------
+    Service
+        The PyWPS WebProcessing Service instance
+    """
+    return service
diff --git a/wrfhydro/forcing_data/__init__.py b/wrfhydro/forcing_data/__init__.py
new file mode 100644
index 000000000..e9d292a89
--- /dev/null
+++ b/wrfhydro/forcing_data/__init__.py
@@ -0,0 +1,18 @@
+# -*- coding: utf-8 -*-
+# Copyright © 2021-2025 Geospatial Research Institute Toi Hangarau
+# LICENSE: https://github.com/GeospatialResearch/Digital-Twins/blob/master/LICENSE
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Affero General Public License as
+# published by the Free Software Foundation, either version 3 of the
+# License, or (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU Affero General Public License for more details.
+#
+# You should have received a copy of the GNU Affero General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+"""Package to handle data and operations specific to meteorological forcing data."""
diff --git a/wrfhydro/forcing_data/read_forcing_data.py b/wrfhydro/forcing_data/read_forcing_data.py
new file mode 100644
index 000000000..824079044
--- /dev/null
+++ b/wrfhydro/forcing_data/read_forcing_data.py
@@ -0,0 +1,39 @@
+# # -*- coding: utf-8 -*-
+# Copyright © 2021-2025 Geospatial Research Institute Toi Hangarau
+# LICENSE: https://github.com/GeospatialResearch/Digital-Twins/blob/master/LICENSE
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Affero General Public License as
+# published by the Free Software Foundation, either version 3 of the
+# License, or (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU Affero General Public License for more details.
+#
+# You should have received a copy of the GNU Affero General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+"""Read meteorological forcing data into a form ready for use within WRF-Hydro modelling."""
+
+import geopandas as gpd
+
+from src.digitaltwin.utils import LogLevel
+
+
+def main(
+    selected_polygon_gdf: gpd.GeoDataFrame,
+    log_level: LogLevel = LogLevel.DEBUG
+) -> None:
+    """
+    Read meteorological forcing data into a form ready for use within WRF-Hydro modelling.
+
+    Parameters
+    ----------
+    selected_polygon_gdf : gpd.GeoDataFrame
+        A GeoDataFrame representing the selected polygon, i.e., the catchment area.
+    log_level : LogLevel = LogLevel.DEBUG
+        The log level to set for the root logger. Defaults to LogLevel.DEBUG.
+    """
+    pass
diff --git a/wrfhydro/land_surface_model/__init__.py b/wrfhydro/land_surface_model/__init__.py
new file mode 100644
index 000000000..c4029ee4f
--- /dev/null
+++ b/wrfhydro/land_surface_model/__init__.py
@@ -0,0 +1,18 @@
+# -*- coding: utf-8 -*-
+# Copyright © 2021-2025 Geospatial Research Institute Toi Hangarau
+# LICENSE: https://github.com/GeospatialResearch/Digital-Twins/blob/master/LICENSE
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Affero General Public License as
+# published by the Free Software Foundation, either version 3 of the
+# License, or (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU Affero General Public License for more details.
+#
+# You should have received a copy of the GNU Affero General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+"""Package to handle data and operations specific to forming land surface models."""
diff --git a/wrfhydro/land_surface_model/create_land_surface_model.py b/wrfhydro/land_surface_model/create_land_surface_model.py
new file mode 100644
index 000000000..8eefa727e
--- /dev/null
+++ b/wrfhydro/land_surface_model/create_land_surface_model.py
@@ -0,0 +1,39 @@
+# # -*- coding: utf-8 -*-
+# Copyright © 2021-2025 Geospatial Research Institute Toi Hangarau
+# LICENSE: https://github.com/GeospatialResearch/Digital-Twins/blob/master/LICENSE
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Affero General Public License as
+# published by the Free Software Foundation, either version 3 of the
+# License, or (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU Affero General Public License for more details.
+#
+# You should have received a copy of the GNU Affero General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+"""Create a land surface model for use within WRF-Hydro modelling."""
+
+import geopandas as gpd
+
+from src.digitaltwin.utils import LogLevel
+
+
+def main(
+    selected_polygon_gdf: gpd.GeoDataFrame,
+    log_level: LogLevel = LogLevel.DEBUG
+) -> None:
+    """
+    Create a land surface model for use within WRF-Hydro modelling.
+
+    Parameters
+    ----------
+    selected_polygon_gdf : gpd.GeoDataFrame
+        A GeoDataFrame representing the selected polygon, i.e., the catchment area.
+    log_level : LogLevel = LogLevel.DEBUG
+        The log level to set for the root logger. Defaults to LogLevel.DEBUG.
+    """
+    pass
diff --git a/wrfhydro/run_all.py b/wrfhydro/run_all.py
new file mode 100644
index 000000000..3b56b0a1d
--- /dev/null
+++ b/wrfhydro/run_all.py
@@ -0,0 +1,48 @@
+# -*- coding: utf-8 -*-
+# Copyright © 2021-2025 Geospatial Research Institute Toi Hangarau
+# LICENSE: https://github.com/GeospatialResearch/Digital-Twins/blob/master/LICENSE
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Affero General Public License as
+# published by the Free Software Foundation, either version 3 of the
+# License, or (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU Affero General Public License for more details.
+#
+# You should have received a copy of the GNU Affero General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+"""This script runs each module in the Digital Twin using a Sample Polygon."""
+import pathlib
+
+from src.digitaltwin import retrieve_from_instructions
+from src.digitaltwin.utils import LogLevel
+from src.run_all import create_sample_polygon, main
+from wrfhydro.forcing_data import read_forcing_data
+from wrfhydro.land_surface_model import create_land_surface_model
+from wrfhydro.scenario import run_wrf_hydro_scenario
+
+DEFAULT_MODULES_TO_PARAMETERS = {
+    retrieve_from_instructions: {
+        "log_level": LogLevel.INFO,
+        "instruction_json_path": pathlib.Path("wrfhydro/static_boundary_instructions.json").as_posix()
+    },
+    create_land_surface_model: {
+        "log_level": LogLevel.INFO,
+    },
+    read_forcing_data: {
+        "log_level": LogLevel.INFO,
+    },
+    run_wrf_hydro_scenario: {
+        "log_level": LogLevel.INFO,
+    },
+}
+
+if __name__ == '__main__':
+    sample_polygon = create_sample_polygon()
+
+    # Run all modules with sample polygon that intentionally contains slight rounding errors.
+    main(sample_polygon, DEFAULT_MODULES_TO_PARAMETERS)
diff --git a/wrfhydro/scenario/__init__.py b/wrfhydro/scenario/__init__.py
new file mode 100644
index 000000000..e9a6fb108
--- /dev/null
+++ b/wrfhydro/scenario/__init__.py
@@ -0,0 +1,18 @@
+# -*- coding: utf-8 -*-
+# Copyright © 2021-2025 Geospatial Research Institute Toi Hangarau
+# LICENSE: https://github.com/GeospatialResearch/Digital-Twins/blob/master/LICENSE
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Affero General Public License as
+# published by the Free Software Foundation, either version 3 of the
+# License, or (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU Affero General Public License for more details.
+#
+# You should have received a copy of the GNU Affero General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+"""Package to handle data and operations specific to running WRF-Hydro scenarios."""
diff --git a/wrfhydro/scenario/flood_scenario_process_service.py b/wrfhydro/scenario/flood_scenario_process_service.py
new file mode 100644
index 000000000..34ecc061d
--- /dev/null
+++ b/wrfhydro/scenario/flood_scenario_process_service.py
@@ -0,0 +1,153 @@
+# # -*- coding: utf-8 -*-
+# Copyright © 2021-2025 Geospatial Research Institute Toi Hangarau
+# LICENSE: https://github.com/GeospatialResearch/Digital-Twins/blob/master/LICENSE
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Affero General Public License as
+# published by the Free Software Foundation, either version 3 of the
+# License, or (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU Affero General Public License for more details.
+#
+# You should have received a copy of the GNU Affero General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+"""Defines PyWPS WebProcessingService process for creating a WRF-Hydro flooding scenario."""
+
+import json
+
+from pywps import BoundingBoxInput, ComplexOutput, Format, Process, WPSRequest
+from pywps.response.execute import ExecuteResponse
+from shapely import box
+
+from wrfhydro import tasks
+from src.config import EnvVariable as EnvVar
+
+
+class FloodScenarioProcessService(Process):
+    """Class representing a WebProcessingService process for creating a flooding scenario"""
+
+    # pylint: disable=too-few-public-methods
+
+    def __init__(self) -> None:
+        """Define inputs and outputs of the WPS process, and assign process handler."""
+        # Create bounding box WPS inputs
+        inputs = [
+            BoundingBoxInput("bboxIn", "Area of Interest", crss=["epsg:4326"]),
+        ]
+        # Create flood inundation WPS outputs
+        outputs = [
+            ComplexOutput("floodDepth", "Maximum Flood Depth",
+                          supported_formats=[Format("application/vnd.terriajs.catalog-member+json")]),
+            ComplexOutput("floodedBuildings", "Flooded Buildings",
+                          supported_formats=[Format("application/vnd.terriajs.catalog-member+json")])
+        ]
+
+        # Initialise the process
+        super().__init__(
+            self._handler,
+            identifier="wrfhydro",
+            title="Model a flood scenario using WRF-Hydro.",
+            inputs=inputs,
+            outputs=outputs,
+        )
+
+    @staticmethod
+    def _handler(request: WPSRequest, response: ExecuteResponse) -> None:
+        """
+        Process handler for modelling a flood scenario
+
+        Parameters
+        ----------
+        request : WPSRequest
+            The WPS request, containing input parameters.
+        response : ExecuteResponse
+            The WPS response, containing output data.
+        """
+        # Get coordinates from bounding box input
+        bounding_box_input = request.inputs['bboxIn'][0]
+        ymin, xmin = bounding_box_input.ll  # lower left
+        ymax, xmax = bounding_box_input.ur  # upper right
+
+        # Form bounding box into standard shapely.box
+        bounding_box = box(xmin, ymin, xmax, ymax)
+
+        modelling_task = tasks.create_scenario_for_area(bounding_box.wkt)
+        scenario_id = modelling_task.get()
+
+        # Add Geoserver JSON Catalog entries to WPS response for use by Terria
+        response.outputs['floodDepth'].data = json.dumps(flood_depth_catalog(scenario_id))
+        response.outputs['floodedBuildings'].data = json.dumps(building_flood_status_catalog(scenario_id))
+
+
+def building_flood_status_catalog(scenario_id: int) -> dict:
+    """
+    Create a dictionary in the format of a terria js catalog json for the building flood status layer.
+
+    Parameters
+    ----------
+    scenario_id : int
+        The ID of the scenario to create the catalog item for.
+
+    Returns
+    ----------
+    dict
+        The TerriaJS catalog item JSON for the building flood status layer.
+    """
+    dataset_name = "Building Flood Status"
+    gs_building_workspace = f"{EnvVar.POSTGRES_DB}-buildings"
+    gs_building_url = f"{EnvVar.GEOSERVER_HOST}:{EnvVar.GEOSERVER_PORT}/geoserver/{gs_building_workspace}/ows"
+    # Open and read HTML/mustache template file for infobox
+    return {
+        "type": "wfs",
+        "name": dataset_name,
+        "url": gs_building_url,
+        "typeNames": f"{gs_building_workspace}:building_flood_status",
+        "parameters": {
+            "viewparams": f"scenario:{scenario_id}",
+        },
+        "maxFeatures": 300000,
+        "heightProperty": "extruded_height",
+        "legends": [{
+            "title": "Building Flood Status",
+            "items": [
+                {
+                    "title": "Non-Flooded",
+                    "color": "darkgreen"
+                },
+                {
+                    "title": "Flooded",
+                    "color": "darkred"
+                }
+            ]
+        }]
+    }
+
+
+def flood_depth_catalog(scenario_id: int) -> dict:
+    """
+    Create a dictionary in the format of a terria js catalog json for the flood depth layer.
+
+    Parameters
+    ----------
+    scenario_id : int
+        The ID of the scenario to create the catalog item for.
+
+    Returns
+    ----------
+    dict
+        The TerriaJS catalog item JSON for the flood depth layer.
+    """
+    gs_flood_model_workspace = f"{EnvVar.POSTGRES_DB}-wh-model-outputs"
+    gs_flood_url = f"{EnvVar.GEOSERVER_HOST}:{EnvVar.GEOSERVER_PORT}/geoserver/{gs_flood_model_workspace}/ows"
+
+    return {
+        "type": "wms",
+        "name": "Flood Depth",
+        "url": gs_flood_url,
+        "layers": f"{gs_flood_model_workspace}:output_{scenario_id}",
+        "styles": "viridis_raster"
+    }
diff --git a/wrfhydro/scenario/run_wrf_hydro_scenario.py b/wrfhydro/scenario/run_wrf_hydro_scenario.py
new file mode 100644
index 000000000..47ca67213
--- /dev/null
+++ b/wrfhydro/scenario/run_wrf_hydro_scenario.py
@@ -0,0 +1,39 @@
+# # -*- coding: utf-8 -*-
+# Copyright © 2021-2025 Geospatial Research Institute Toi Hangarau
+# LICENSE: https://github.com/GeospatialResearch/Digital-Twins/blob/master/LICENSE
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Affero General Public License as
+# published by the Free Software Foundation, either version 3 of the
+# License, or (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU Affero General Public License for more details.
+#
+# You should have received a copy of the GNU Affero General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+"""Run a WRF-Hydro scenario for the selected polygon and serve outputs."""
+
+import geopandas as gpd
+
+from src.digitaltwin.utils import LogLevel
+
+
+def main(
+    selected_polygon_gdf: gpd.GeoDataFrame,
+    log_level: LogLevel = LogLevel.DEBUG
+) -> None:
+    """
+    Run a WRF-Hydro scenario for the selected polygon and serve outputs.
+
+    Parameters
+    ----------
+    selected_polygon_gdf : gpd.GeoDataFrame
+        A GeoDataFrame representing the selected polygon, i.e., the catchment area.
+    log_level : LogLevel = LogLevel.DEBUG
+        The log level to set for the root logger. Defaults to LogLevel.DEBUG.
+    """
+    pass
diff --git a/wrfhydro/static_boundary_instructions.json b/wrfhydro/static_boundary_instructions.json
new file mode 100644
index 000000000..ccc8c6662
--- /dev/null
+++ b/wrfhydro/static_boundary_instructions.json
@@ -0,0 +1,16 @@
+{
+  "instruction_coastlines": {
+    "data_provider": "LINZ",
+    "layer_id": 50258,
+    "table_name": "nz_coastlines",
+    "coverage_area": "New Zealand",
+    "url": "https://data.linz.govt.nz/layer/50258-nz-coastlines-topo-150k/"
+  },
+  "instruction_buildings": {
+    "data_provider": "LINZ",
+    "layer_id": 101292,
+    "table_name": "nz_building_outlines",
+    "unique_column_name": "building_outline_id",
+    "url": "https://data.linz.govt.nz/layer/101292-nz-building-outlines-all-sources/"
+  }
+}
diff --git a/wrfhydro/tasks.py b/wrfhydro/tasks.py
new file mode 100644
index 000000000..41b27d89f
--- /dev/null
+++ b/wrfhydro/tasks.py
@@ -0,0 +1,106 @@
+"""
+Runs backend tasks using Celery. Allowing for multiple long-running tasks to complete in the background.
+Allows the frontend to send tasks and retrieve status later.
+"""
+import logging
+
+from celery import result, signals
+from celery.worker.consumer import Consumer
+import geopandas as gpd
+
+from src.digitaltwin import retrieve_from_instructions
+from src.digitaltwin.utils import setup_logging
+from src.tasks import add_base_data_to_db, app, OnFailureStateTask, wkt_to_gdf  # pylint: disable=cyclic-import
+from wrfhydro.forcing_data import read_forcing_data
+from wrfhydro.land_surface_model import create_land_surface_model
+from wrfhydro.scenario import run_wrf_hydro_scenario
+from wrfhydro.run_all import DEFAULT_MODULES_TO_PARAMETERS
+
+setup_logging()
+log = logging.getLogger(__name__)
+
+
+@signals.worker_ready.connect
+def on_startup(sender: Consumer, **_kwargs: None) -> None:  # pylint: disable=missing-param-doc
+    """
+    Initialise database, runs when Celery instance is ready.
+
+    Parameters
+    ----------
+    sender : Consumer
+        The Celery worker node instance
+    """
+    with sender.app.connection() as conn:
+        # Gather area of interest from file.
+        aoi_wkt = gpd.read_file("selected_polygon.geojson").to_crs(4326).geometry[0].wkt
+        # Send a task to initialise this area of interest.
+        base_data_parameters = DEFAULT_MODULES_TO_PARAMETERS[retrieve_from_instructions]
+        sender.app.send_task("src.tasks.add_base_data_to_db", args=[aoi_wkt, base_data_parameters], connection=conn)
+
+
+def create_scenario_for_area(selected_polygon_wkt: str) -> result.GroupResult:
+    """
+    Create a model for the area using series of chained (sequential) sub-tasks.
+
+    Parameters
+    ----------
+    selected_polygon_wkt : str
+        The polygon defining the selected area to run the model for. Defined in WKT form.
+
+    Returns
+    -------
+    result.GroupResult
+        The task result for the long-running group of tasks. The task ID represents the final task in the group.
+    """
+    base_data_parameters = DEFAULT_MODULES_TO_PARAMETERS[retrieve_from_instructions]
+    return (
+        add_base_data_to_db.si(selected_polygon_wkt, base_data_parameters) |
+        create_land_surface_model_task.si(selected_polygon_wkt) |
+        read_forcing_data_task.si(selected_polygon_wkt) |
+        run_wrf_scenario_task.si(selected_polygon_wkt)
+    )()
+
+
+@app.task(base=OnFailureStateTask)
+def create_land_surface_model_task(selected_polygon_wkt: str) -> None:
+    """
+    Task to ensure rainfall input data for the given area is added to the database and model input files are created.
+
+    Parameters
+    ----------
+    selected_polygon_wkt : str
+        The polygon defining the selected area to add rainfall data for. Defined in WKT form.
+    """
+    parameters = DEFAULT_MODULES_TO_PARAMETERS[create_land_surface_model]
+    selected_polygon = wkt_to_gdf(selected_polygon_wkt)
+    create_land_surface_model.main(selected_polygon, **parameters)
+
+
+@app.task(base=OnFailureStateTask)
+def read_forcing_data_task(selected_polygon_wkt: str) -> None:
+    """
+    Task to ensure meteorological forcing data is processed for the given area and added to the database.
+
+    Parameters
+    ----------
+    selected_polygon_wkt : str
+        The polygon defining the selected area to process the forcing data for. Defined in WKT form.
+    """
+    parameters = DEFAULT_MODULES_TO_PARAMETERS[read_forcing_data]
+    selected_polygon = wkt_to_gdf(selected_polygon_wkt)
+    read_forcing_data.main(selected_polygon, **parameters)
+
+
+@app.task(base=OnFailureStateTask)
+def run_wrf_scenario_task(selected_polygon_wkt: str) -> None:
+    """
+    Task to run a WRF Hydro scenario for a given area.
+
+    Parameters
+    ----------
+    selected_polygon_wkt : str
+        The polygon defining the selected area to run a wrf_hydro scenario for. Defined in WKT form.
+    """
+    parameters = DEFAULT_MODULES_TO_PARAMETERS[run_wrf_hydro_scenario]
+    selected_polygon = wkt_to_gdf(selected_polygon_wkt)
+    run_wrf_hydro_scenario.main(selected_polygon, **parameters)