diff --git a/watertap/unit_models/__init__.py b/watertap/unit_models/__init__.py
index 7d0c4d7b12..3bd828e717 100644
--- a/watertap/unit_models/__init__.py
+++ b/watertap/unit_models/__init__.py
@@ -24,3 +24,5 @@
 from .electrodialysis_1D import Electrodialysis1D
 from .gac import GAC
 from .ion_exchange_0D import IonExchange0D
+from .thickener import Thickener
+from .dewatering import DewateringUnit
diff --git a/watertap/unit_models/dewatering.py b/watertap/unit_models/dewatering.py
new file mode 100644
index 0000000000..9ed7205b19
--- /dev/null
+++ b/watertap/unit_models/dewatering.py
@@ -0,0 +1,163 @@
+###############################################################################
+# WaterTAP Copyright (c) 2021, The Regents of the University of California,
+# through Lawrence Berkeley National Laboratory, Oak Ridge National
+# Laboratory, National Renewable Energy Laboratory, and National Energy
+# Technology Laboratory (subject to receipt of any required approvals from
+# the U.S. Dept. of Energy). All rights reserved.
+#
+# Please see the files COPYRIGHT.md and LICENSE.md for full copyright and license
+# information, respectively. These files are also available online at the URL
+# "https://github.com/watertap-org/watertap/"
+#
+###############################################################################
+"""
+Dewatering unit model for BSM2. Based on IDAES separator unit
+
+Model based on 
+
+J. Alex, L. Benedetti, J.B. Copp, K.V. Gernaey, U. Jeppsson,
+I. Nopens, M.N. Pons, C. Rosen, J.P. Steyer and
+P. A. Vanrolleghem
+Benchmark Simulation Model no. 2 (BSM2)
+"""
+# Import IDAES cores
+from idaes.core import (
+    declare_process_block_class,
+)
+from idaes.models.unit_models.separator import SeparatorData, SplittingType
+
+from idaes.core.util.tables import create_stream_table_dataframe
+import idaes.logger as idaeslog
+
+from pyomo.environ import (
+    Param,
+    units as pyunits,
+    Set,
+)
+
+from idaes.core.util.exceptions import (
+    ConfigurationError,
+)
+
+__author__ = "Alejandro Garciadiego"
+
+
+# Set up logger
+_log = idaeslog.getLogger(__name__)
+
+
+@declare_process_block_class("DewateringUnit")
+class DewateringData(SeparatorData):
+    """
+    Dewatering unit block for BSM2
+    """
+
+    CONFIG = SeparatorData.CONFIG()
+    CONFIG.outlet_list = ["underflow", "overflow"]
+    CONFIG.split_basis = SplittingType.componentFlow
+
+    def build(self):
+        """
+        Begin building model.
+        Args:
+            None
+        Returns:
+            None
+        """
+
+        # Call UnitModel.build to set up dynamics
+        super(DewateringData, self).build()
+
+        if "underflow" and "overflow" not in self.config.outlet_list:
+            raise ConfigurationError(
+                "{} encountered unrecognised "
+                "outlet_list. This should not "
+                "occur - please use overflow "
+                "and underflow as outlets.".format(self.name)
+            )
+
+        self.p_dewat = Param(
+            initialize=0.28,
+            units=pyunits.dimensionless,
+            mutable=True,
+            doc="Percentage of suspended solids in the underflow",
+        )
+
+        self.TSS_rem = Param(
+            initialize=0.98,
+            units=pyunits.dimensionless,
+            mutable=True,
+            doc="Percentage of suspended solids removed",
+        )
+
+        @self.Expression(self.flowsheet().time, doc="Suspended solid concentration")
+        def TSS(blk, t):
+            return 0.75 * (
+                blk.inlet.conc_mass_comp[t, "X_I"]
+                + blk.inlet.conc_mass_comp[t, "X_P"]
+                + blk.inlet.conc_mass_comp[t, "X_BH"]
+                + blk.inlet.conc_mass_comp[t, "X_BA"]
+                + blk.inlet.conc_mass_comp[t, "X_S"]
+            )
+
+        @self.Expression(self.flowsheet().time, doc="Dewatering factor")
+        def f_dewat(blk, t):
+            return blk.p_dewat * (10 / (blk.TSS[t]))
+
+        @self.Expression(self.flowsheet().time, doc="Remove factor")
+        def f_q_du(blk, t):
+            return blk.TSS_rem / (pyunits.kg / pyunits.m**3) / 100 / blk.f_dewat[t]
+
+        self.non_particulate_components = Set(
+            initialize=[
+                "S_I",
+                "S_S",
+                "S_O",
+                "S_NO",
+                "S_NH",
+                "S_ND",
+                "H2O",
+                "S_ALK",
+            ]
+        )
+
+        self.particulate_components = Set(
+            initialize=["X_I", "X_S", "X_P", "X_BH", "X_BA", "X_ND"]
+        )
+
+        @self.Constraint(
+            self.flowsheet().time,
+            self.particulate_components,
+            doc="particulate fraction",
+        )
+        def overflow_particulate_fraction(blk, t, i):
+            return blk.split_fraction[t, "overflow", i] == 1 - blk.TSS_rem
+
+        @self.Constraint(
+            self.flowsheet().time,
+            self.non_particulate_components,
+            doc="soluble fraction",
+        )
+        def non_particulate_components(blk, t, i):
+            return blk.split_fraction[t, "overflow", i] == 1 - blk.f_q_du[t]
+
+    def _get_performance_contents(self, time_point=0):
+        var_dict = {}
+        for k in self.split_fraction.keys():
+            if k[0] == time_point:
+                var_dict[f"Split Fraction [{str(k[1:])}]"] = self.split_fraction[k]
+        return {"vars": var_dict}
+
+    def _get_stream_table_contents(self, time_point=0):
+        outlet_list = self.create_outlet_list()
+
+        io_dict = {}
+        if self.config.mixed_state_block is None:
+            io_dict["Inlet"] = self.mixed_state
+        else:
+            io_dict["Inlet"] = self.config.mixed_state_block
+
+        for o in outlet_list:
+            io_dict[o] = getattr(self, o + "_state")
+
+        return create_stream_table_dataframe(io_dict, time_point=time_point)
diff --git a/watertap/unit_models/tests/test_dewatering_unit.py b/watertap/unit_models/tests/test_dewatering_unit.py
new file mode 100644
index 0000000000..fa97744c64
--- /dev/null
+++ b/watertap/unit_models/tests/test_dewatering_unit.py
@@ -0,0 +1,272 @@
+#################################################################################
+# WaterTAP Copyright (c) 2020-2023, The Regents of the University of California,
+# through Lawrence Berkeley National Laboratory, Oak Ridge National Laboratory,
+# National Renewable Energy Laboratory, and National Energy Technology
+# Laboratory (subject to receipt of any required approvals from the U.S. Dept.
+# of Energy). All rights reserved.
+#
+# Please see the files COPYRIGHT.md and LICENSE.md for full copyright and license
+# information, respectively. These files are also available online at the URL
+# "https://github.com/watertap-org/watertap/"
+#################################################################################
+"""
+Tests for dewatering unit example.
+"""
+
+import pytest
+from pyomo.environ import (
+    ConcreteModel,
+    value,
+    assert_optimal_termination,
+)
+
+from idaes.core import (
+    FlowsheetBlock,
+    MaterialBalanceType,
+    MomentumBalanceType,
+)
+
+from idaes.models.unit_models.separator import SplittingType
+
+from pyomo.environ import (
+    units,
+)
+
+from idaes.core.solvers import get_solver
+from idaes.core.util.model_statistics import (
+    degrees_of_freedom,
+    number_variables,
+    number_total_constraints,
+    number_unused_variables,
+)
+import idaes.core.util.scaling as iscale
+from idaes.core.util.testing import (
+    initialization_tester,
+)
+
+from idaes.core.util.exceptions import (
+    ConfigurationError,
+)
+
+from watertap.unit_models.dewatering import DewateringUnit
+from watertap.property_models.activated_sludge.asm1_properties import (
+    ASM1ParameterBlock,
+)
+
+from pyomo.util.check_units import assert_units_consistent
+
+# -----------------------------------------------------------------------------
+# Get default solver for testing
+solver = get_solver()
+
+# -----------------------------------------------------------------------------
+@pytest.mark.unit
+def test_config():
+    m = ConcreteModel()
+
+    m.fs = FlowsheetBlock(dynamic=False)
+
+    m.fs.props = ASM1ParameterBlock()
+
+    m.fs.unit = DewateringUnit(property_package=m.fs.props)
+
+    assert len(m.fs.unit.config) == 15
+
+    assert not m.fs.unit.config.dynamic
+    assert not m.fs.unit.config.has_holdup
+    assert m.fs.unit.config.material_balance_type == MaterialBalanceType.useDefault
+    assert m.fs.unit.config.momentum_balance_type == MomentumBalanceType.pressureTotal
+    assert "underflow" in m.fs.unit.config.outlet_list
+    assert "overflow" in m.fs.unit.config.outlet_list
+    assert SplittingType.componentFlow is m.fs.unit.config.split_basis
+
+
+@pytest.mark.unit
+def test_list_error():
+    m = ConcreteModel()
+
+    m.fs = FlowsheetBlock(dynamic=False)
+
+    m.fs.props = ASM1ParameterBlock()
+
+    with pytest.raises(
+        ConfigurationError,
+        match="fs.unit encountered unrecognised "
+        "outlet_list. This should not "
+        "occur - please use overflow "
+        "and underflow as outlets.",
+    ):
+        m.fs.unit = DewateringUnit(
+            property_package=m.fs.props,
+            outlet_list=["outlet1", "outlet2"],
+        )
+
+
+# -----------------------------------------------------------------------------
+class TestDu(object):
+    @pytest.fixture(scope="class")
+    def du(self):
+        m = ConcreteModel()
+        m.fs = FlowsheetBlock(dynamic=False)
+
+        m.fs.props = ASM1ParameterBlock()
+
+        m.fs.unit = DewateringUnit(property_package=m.fs.props)
+
+        m.fs.unit.inlet.flow_vol.fix(178.4674 * units.m**3 / units.day)
+        m.fs.unit.inlet.temperature.fix(308.15 * units.K)
+        m.fs.unit.inlet.pressure.fix(1 * units.atm)
+
+        m.fs.unit.inlet.conc_mass_comp[0, "S_I"].fix(130.867 * units.mg / units.liter)
+        m.fs.unit.inlet.conc_mass_comp[0, "S_S"].fix(258.5789 * units.mg / units.liter)
+        m.fs.unit.inlet.conc_mass_comp[0, "X_I"].fix(
+            17216.2434 * units.mg / units.liter
+        )
+        m.fs.unit.inlet.conc_mass_comp[0, "X_S"].fix(2611.4843 * units.mg / units.liter)
+        m.fs.unit.inlet.conc_mass_comp[0, "X_BH"].fix(1e-6 * units.mg / units.liter)
+        m.fs.unit.inlet.conc_mass_comp[0, "X_BA"].fix(1e-6 * units.mg / units.liter)
+        m.fs.unit.inlet.conc_mass_comp[0, "X_P"].fix(626.0652 * units.mg / units.liter)
+        m.fs.unit.inlet.conc_mass_comp[0, "S_O"].fix(1e-6 * units.mg / units.liter)
+        m.fs.unit.inlet.conc_mass_comp[0, "S_NO"].fix(1e-6 * units.mg / units.liter)
+        m.fs.unit.inlet.conc_mass_comp[0, "S_NH"].fix(
+            1442.7882 * units.mg / units.liter
+        )
+        m.fs.unit.inlet.conc_mass_comp[0, "S_ND"].fix(0.54323 * units.mg / units.liter)
+        m.fs.unit.inlet.conc_mass_comp[0, "X_ND"].fix(100.8668 * units.mg / units.liter)
+        m.fs.unit.inlet.alkalinity.fix(97.8459 * units.mol / units.m**3)
+
+        return m
+
+    @pytest.mark.build
+    @pytest.mark.unit
+    def test_build(self, du):
+
+        assert hasattr(du.fs.unit, "inlet")
+        assert len(du.fs.unit.inlet.vars) == 5
+        assert hasattr(du.fs.unit.inlet, "flow_vol")
+        assert hasattr(du.fs.unit.inlet, "conc_mass_comp")
+        assert hasattr(du.fs.unit.inlet, "temperature")
+        assert hasattr(du.fs.unit.inlet, "pressure")
+        assert hasattr(du.fs.unit.inlet, "alkalinity")
+
+        assert hasattr(du.fs.unit, "underflow")
+        assert len(du.fs.unit.underflow.vars) == 5
+        assert hasattr(du.fs.unit.underflow, "flow_vol")
+        assert hasattr(du.fs.unit.underflow, "conc_mass_comp")
+        assert hasattr(du.fs.unit.underflow, "temperature")
+        assert hasattr(du.fs.unit.underflow, "pressure")
+        assert hasattr(du.fs.unit.underflow, "alkalinity")
+
+        assert hasattr(du.fs.unit, "overflow")
+        assert len(du.fs.unit.overflow.vars) == 5
+        assert hasattr(du.fs.unit.overflow, "flow_vol")
+        assert hasattr(du.fs.unit.overflow, "conc_mass_comp")
+        assert hasattr(du.fs.unit.overflow, "temperature")
+        assert hasattr(du.fs.unit.overflow, "pressure")
+        assert hasattr(du.fs.unit.overflow, "alkalinity")
+
+        assert number_variables(du) == 76
+        assert number_total_constraints(du) == 60
+        assert number_unused_variables(du) == 0
+
+    @pytest.mark.unit
+    def test_dof(self, du):
+        assert degrees_of_freedom(du) == 0
+
+    @pytest.mark.unit
+    def test_units(self, du):
+        assert_units_consistent(du)
+
+    @pytest.mark.solver
+    @pytest.mark.skipif(solver is None, reason="Solver not available")
+    @pytest.mark.component
+    def test_initialize(self, du):
+
+        iscale.calculate_scaling_factors(du)
+        initialization_tester(du)
+
+    @pytest.mark.solver
+    @pytest.mark.skipif(solver is None, reason="Solver not available")
+    @pytest.mark.component
+    def test_solve(self, du):
+        solver = get_solver()
+        results = solver.solve(du)
+        assert_optimal_termination(results)
+
+    @pytest.mark.solver
+    @pytest.mark.skipif(solver is None, reason="Solver not available")
+    @pytest.mark.component
+    def test_solution(self, du):
+        assert pytest.approx(101325.0, rel=1e-3) == value(
+            du.fs.unit.overflow.pressure[0]
+        )
+        assert pytest.approx(308.15, rel=1e-3) == value(
+            du.fs.unit.overflow.temperature[0]
+        )
+        assert pytest.approx(0.001954, rel=1e-3) == value(
+            du.fs.unit.overflow.flow_vol[0]
+        )
+        assert pytest.approx(0.1308, rel=1e-3) == value(
+            du.fs.unit.overflow.conc_mass_comp[0, "S_I"]
+        )
+        assert pytest.approx(0.2585, rel=1e-3) == value(
+            du.fs.unit.overflow.conc_mass_comp[0, "S_S"]
+        )
+        assert pytest.approx(0.3638, rel=1e-3) == value(
+            du.fs.unit.overflow.conc_mass_comp[0, "X_I"]
+        )
+        assert pytest.approx(0.0552, rel=1e-3) == value(
+            du.fs.unit.overflow.conc_mass_comp[0, "X_S"]
+        )
+        assert value(du.fs.unit.overflow.conc_mass_comp[0, "X_BH"]) <= 1e-6
+        assert value(du.fs.unit.overflow.conc_mass_comp[0, "X_BA"]) <= 1e-6
+        assert pytest.approx(0.01323, rel=1e-3) == value(
+            du.fs.unit.overflow.conc_mass_comp[0, "X_P"]
+        )
+        assert value(du.fs.unit.overflow.conc_mass_comp[0, "S_O"]) <= 1e-6
+        assert value(du.fs.unit.overflow.conc_mass_comp[0, "S_NO"]) <= 1e-6
+        assert pytest.approx(1.4427, rel=1e-3) == value(
+            du.fs.unit.overflow.conc_mass_comp[0, "S_NH"]
+        )
+        assert pytest.approx(0.000543, rel=1e-3) == value(
+            du.fs.unit.overflow.conc_mass_comp[0, "S_ND"]
+        )
+        assert pytest.approx(0.00213, rel=1e-3) == value(
+            du.fs.unit.overflow.conc_mass_comp[0, "X_ND"]
+        )
+        assert pytest.approx(0.09784, rel=1e-3) == value(
+            du.fs.unit.overflow.alkalinity[0]
+        )
+
+    @pytest.mark.solver
+    @pytest.mark.skipif(solver is None, reason="Solver not available")
+    @pytest.mark.component
+    def test_conservation(self, du):
+        assert (
+            abs(
+                value(
+                    du.fs.unit.inlet.flow_vol[0] * du.fs.props.dens_mass
+                    - du.fs.unit.overflow.flow_vol[0] * du.fs.props.dens_mass
+                    - du.fs.unit.underflow.flow_vol[0] * du.fs.props.dens_mass
+                )
+            )
+            <= 1e-6
+        )
+        for i in du.fs.props.solute_set:
+            assert (
+                abs(
+                    value(
+                        du.fs.unit.inlet.flow_vol[0]
+                        * du.fs.unit.inlet.conc_mass_comp[0, i]
+                        - du.fs.unit.overflow.flow_vol[0]
+                        * du.fs.unit.overflow.conc_mass_comp[0, i]
+                        - du.fs.unit.underflow.flow_vol[0]
+                        * du.fs.unit.underflow.conc_mass_comp[0, i]
+                    )
+                )
+                <= 1e-6
+            )
+
+    @pytest.mark.unit
+    def test_report(self, du):
+        du.fs.unit.report()
diff --git a/watertap/unit_models/tests/test_thickener_unit.py b/watertap/unit_models/tests/test_thickener_unit.py
new file mode 100644
index 0000000000..356f806740
--- /dev/null
+++ b/watertap/unit_models/tests/test_thickener_unit.py
@@ -0,0 +1,277 @@
+#################################################################################
+# WaterTAP Copyright (c) 2020-2023, The Regents of the University of California,
+# through Lawrence Berkeley National Laboratory, Oak Ridge National Laboratory,
+# National Renewable Energy Laboratory, and National Energy Technology
+# Laboratory (subject to receipt of any required approvals from the U.S. Dept.
+# of Energy). All rights reserved.
+#
+# Please see the files COPYRIGHT.md and LICENSE.md for full copyright and license
+# information, respectively. These files are also available online at the URL
+# "https://github.com/watertap-org/watertap/"
+#################################################################################
+"""
+Tests for thickener unit example.
+"""
+
+import pytest
+from pyomo.environ import (
+    ConcreteModel,
+    value,
+    assert_optimal_termination,
+)
+
+from idaes.core import (
+    FlowsheetBlock,
+    MaterialBalanceType,
+    MomentumBalanceType,
+    PhysicalParameterBlock,
+)
+
+from idaes.models.unit_models.separator import SplittingType
+
+from pyomo.environ import (
+    units,
+)
+
+from idaes.core.solvers import get_solver
+from idaes.core.util.model_statistics import (
+    degrees_of_freedom,
+    number_variables,
+    number_total_constraints,
+    number_unused_variables,
+)
+import idaes.core.util.scaling as iscale
+from idaes.core.util.testing import (
+    initialization_tester,
+)
+
+from idaes.core.util.exceptions import (
+    ConfigurationError,
+)
+from watertap.unit_models.thickener import Thickener
+from watertap.property_models.activated_sludge.asm1_properties import (
+    ASM1ParameterBlock,
+)
+
+from pyomo.util.check_units import assert_units_consistent, assert_units_equivalent
+
+# -----------------------------------------------------------------------------
+# Get default solver for testing
+solver = get_solver()
+
+# -----------------------------------------------------------------------------
+@pytest.mark.unit
+def test_config():
+    m = ConcreteModel()
+
+    m.fs = FlowsheetBlock(dynamic=False)
+
+    m.fs.props = ASM1ParameterBlock()
+
+    m.fs.unit = Thickener(property_package=m.fs.props)
+
+    assert len(m.fs.unit.config) == 15
+
+    assert not m.fs.unit.config.dynamic
+    assert not m.fs.unit.config.has_holdup
+    assert m.fs.unit.config.material_balance_type == MaterialBalanceType.useDefault
+    assert m.fs.unit.config.momentum_balance_type == MomentumBalanceType.pressureTotal
+    assert "underflow" in m.fs.unit.config.outlet_list
+    assert "overflow" in m.fs.unit.config.outlet_list
+    assert SplittingType.componentFlow is m.fs.unit.config.split_basis
+
+
+@pytest.mark.unit
+def test_list_error():
+    m = ConcreteModel()
+
+    m.fs = FlowsheetBlock(dynamic=False)
+
+    m.fs.props = ASM1ParameterBlock()
+
+    with pytest.raises(
+        ConfigurationError,
+        match="fs.unit encountered unrecognised "
+        "outlet_list. This should not "
+        "occur - please use overflow "
+        "and underflow as outlets.",
+    ):
+        m.fs.unit = Thickener(
+            property_package=m.fs.props, outlet_list=["outlet1", "outlet2"]
+        )
+
+
+# -----------------------------------------------------------------------------
+class TestThick(object):
+    @pytest.fixture(scope="class")
+    def tu(self):
+        m = ConcreteModel()
+        m.fs = FlowsheetBlock(dynamic=False)
+
+        m.fs.props = ASM1ParameterBlock()
+
+        m.fs.unit = Thickener(property_package=m.fs.props)
+
+        m.fs.unit.inlet.flow_vol.fix(300 * units.m**3 / units.day)
+        m.fs.unit.inlet.temperature.fix(308.15 * units.K)
+        m.fs.unit.inlet.pressure.fix(1 * units.atm)
+
+        m.fs.unit.inlet.conc_mass_comp[0, "S_I"].fix(28.0643 * units.mg / units.liter)
+        m.fs.unit.inlet.conc_mass_comp[0, "S_S"].fix(0.67336 * units.mg / units.liter)
+        m.fs.unit.inlet.conc_mass_comp[0, "X_I"].fix(3036.2175 * units.mg / units.liter)
+        m.fs.unit.inlet.conc_mass_comp[0, "X_S"].fix(63.2392 * units.mg / units.liter)
+        m.fs.unit.inlet.conc_mass_comp[0, "X_BH"].fix(
+            4442.8377 * units.mg / units.liter
+        )
+        m.fs.unit.inlet.conc_mass_comp[0, "X_BA"].fix(332.5958 * units.mg / units.liter)
+        m.fs.unit.inlet.conc_mass_comp[0, "X_P"].fix(1922.8108 * units.mg / units.liter)
+        m.fs.unit.inlet.conc_mass_comp[0, "S_O"].fix(1.3748 * units.mg / units.liter)
+        m.fs.unit.inlet.conc_mass_comp[0, "S_NO"].fix(9.1948 * units.mg / units.liter)
+        m.fs.unit.inlet.conc_mass_comp[0, "S_NH"].fix(0.15845 * units.mg / units.liter)
+        m.fs.unit.inlet.conc_mass_comp[0, "S_ND"].fix(0.55943 * units.mg / units.liter)
+        m.fs.unit.inlet.conc_mass_comp[0, "X_ND"].fix(4.7411 * units.mg / units.liter)
+        m.fs.unit.inlet.alkalinity.fix(4.5646 * units.mol / units.m**3)
+
+        return m
+
+    @pytest.mark.build
+    @pytest.mark.unit
+    def test_build(self, tu):
+
+        assert hasattr(tu.fs.unit, "inlet")
+        assert len(tu.fs.unit.inlet.vars) == 5
+        assert hasattr(tu.fs.unit.inlet, "flow_vol")
+        assert hasattr(tu.fs.unit.inlet, "conc_mass_comp")
+        assert hasattr(tu.fs.unit.inlet, "temperature")
+        assert hasattr(tu.fs.unit.inlet, "pressure")
+        assert hasattr(tu.fs.unit.inlet, "alkalinity")
+
+        assert hasattr(tu.fs.unit, "underflow")
+        assert len(tu.fs.unit.underflow.vars) == 5
+        assert hasattr(tu.fs.unit.underflow, "flow_vol")
+        assert hasattr(tu.fs.unit.underflow, "conc_mass_comp")
+        assert hasattr(tu.fs.unit.underflow, "temperature")
+        assert hasattr(tu.fs.unit.underflow, "pressure")
+        assert hasattr(tu.fs.unit.underflow, "alkalinity")
+
+        assert hasattr(tu.fs.unit, "overflow")
+        assert len(tu.fs.unit.overflow.vars) == 5
+        assert hasattr(tu.fs.unit.overflow, "flow_vol")
+        assert hasattr(tu.fs.unit.overflow, "conc_mass_comp")
+        assert hasattr(tu.fs.unit.overflow, "temperature")
+        assert hasattr(tu.fs.unit.overflow, "pressure")
+        assert hasattr(tu.fs.unit.overflow, "alkalinity")
+
+        assert number_variables(tu) == 76
+        assert number_total_constraints(tu) == 60
+        assert number_unused_variables(tu) == 0
+
+    @pytest.mark.unit
+    def test_dof(self, tu):
+        assert degrees_of_freedom(tu) == 0
+
+    @pytest.mark.unit
+    def test_units(self, tu):
+        assert_units_consistent(tu)
+
+    @pytest.mark.solver
+    @pytest.mark.skipif(solver is None, reason="Solver not available")
+    @pytest.mark.component
+    def test_initialize(self, tu):
+
+        iscale.calculate_scaling_factors(tu)
+        initialization_tester(tu)
+
+    @pytest.mark.solver
+    @pytest.mark.skipif(solver is None, reason="Solver not available")
+    @pytest.mark.component
+    def test_solve(self, tu):
+        solver = get_solver()
+        results = solver.solve(tu)
+        assert_optimal_termination(results)
+
+    @pytest.mark.solver
+    @pytest.mark.skipif(solver is None, reason="Solver not available")
+    @pytest.mark.component
+    def test_solution(self, tu):
+        assert pytest.approx(101325.0, rel=1e-3) == value(
+            tu.fs.unit.overflow.pressure[0]
+        )
+        assert pytest.approx(308.15, rel=1e-3) == value(
+            tu.fs.unit.overflow.temperature[0]
+        )
+        assert pytest.approx(0.003115, rel=1e-3) == value(
+            tu.fs.unit.overflow.flow_vol[0]
+        )
+        assert pytest.approx(0.02806, rel=1e-3) == value(
+            tu.fs.unit.overflow.conc_mass_comp[0, "S_I"]
+        )
+        assert pytest.approx(0.000673, rel=1e-3) == value(
+            tu.fs.unit.overflow.conc_mass_comp[0, "S_S"]
+        )
+        assert pytest.approx(0.06768, rel=1e-3) == value(
+            tu.fs.unit.overflow.conc_mass_comp[0, "X_I"]
+        )
+        assert pytest.approx(0.001409, rel=1e-3) == value(
+            tu.fs.unit.overflow.conc_mass_comp[0, "X_S"]
+        )
+        assert pytest.approx(0.04286, rel=1e-3) == value(
+            tu.fs.unit.overflow.conc_mass_comp[0, "X_P"]
+        )
+        assert pytest.approx(0.099046, rel=1e-3) == value(
+            tu.fs.unit.overflow.conc_mass_comp[0, "X_BH"]
+        )
+        assert pytest.approx(0.007414, rel=1e-3) == value(
+            tu.fs.unit.overflow.conc_mass_comp[0, "X_BA"]
+        )
+        assert pytest.approx(0.001374, rel=1e-3) == value(
+            tu.fs.unit.overflow.conc_mass_comp[0, "S_O"]
+        )
+        assert pytest.approx(0.009194, rel=1e-3) == value(
+            tu.fs.unit.overflow.conc_mass_comp[0, "S_NO"]
+        )
+        assert pytest.approx(0.0001584, rel=1e-3) == value(
+            tu.fs.unit.overflow.conc_mass_comp[0, "S_NH"]
+        )
+        assert pytest.approx(0.0005594, rel=1e-3) == value(
+            tu.fs.unit.overflow.conc_mass_comp[0, "S_ND"]
+        )
+        assert pytest.approx(0.0001056, rel=1e-3) == value(
+            tu.fs.unit.overflow.conc_mass_comp[0, "X_ND"]
+        )
+        assert pytest.approx(0.004564, rel=1e-3) == value(
+            tu.fs.unit.overflow.alkalinity[0]
+        )
+
+    @pytest.mark.solver
+    @pytest.mark.skipif(solver is None, reason="Solver not available")
+    @pytest.mark.component
+    def test_conservation(self, tu):
+        assert (
+            abs(
+                value(
+                    tu.fs.unit.inlet.flow_vol[0] * tu.fs.props.dens_mass
+                    - tu.fs.unit.overflow.flow_vol[0] * tu.fs.props.dens_mass
+                    - tu.fs.unit.underflow.flow_vol[0] * tu.fs.props.dens_mass
+                )
+            )
+            <= 1e-6
+        )
+        for i in tu.fs.props.solute_set:
+            assert (
+                abs(
+                    value(
+                        tu.fs.unit.inlet.flow_vol[0]
+                        * tu.fs.unit.inlet.conc_mass_comp[0, i]
+                        - tu.fs.unit.overflow.flow_vol[0]
+                        * tu.fs.unit.overflow.conc_mass_comp[0, i]
+                        - tu.fs.unit.underflow.flow_vol[0]
+                        * tu.fs.unit.underflow.conc_mass_comp[0, i]
+                    )
+                )
+                <= 1e-6
+            )
+
+    @pytest.mark.unit
+    def test_report(self, tu):
+        tu.fs.unit.report()
diff --git a/watertap/unit_models/thickener.py b/watertap/unit_models/thickener.py
new file mode 100644
index 0000000000..721c88f3bf
--- /dev/null
+++ b/watertap/unit_models/thickener.py
@@ -0,0 +1,167 @@
+###############################################################################
+# WaterTAP Copyright (c) 2021, The Regents of the University of California,
+# through Lawrence Berkeley National Laboratory, Oak Ridge National
+# Laboratory, National Renewable Energy Laboratory, and National Energy
+# Technology Laboratory (subject to receipt of any required approvals from
+# the U.S. Dept. of Energy). All rights reserved.
+#
+# Please see the files COPYRIGHT.md and LICENSE.md for full copyright and license
+# information, respectively. These files are also available online at the URL
+# "https://github.com/watertap-org/watertap/"
+#
+###############################################################################
+"""
+Thickener unit model for BSM2. Based on IDAES separator unit
+
+Model based on 
+
+J. Alex, L. Benedetti, J.B. Copp, K.V. Gernaey, U. Jeppsson,
+I. Nopens, M.N. Pons, C. Rosen, J.P. Steyer and
+P. A. Vanrolleghem
+Benchmark Simulation Model no. 2 (BSM2)
+"""
+
+# Import IDAES cores
+from idaes.core import (
+    declare_process_block_class,
+)
+from idaes.models.unit_models.separator import SeparatorData, SplittingType
+
+from idaes.core.util.tables import create_stream_table_dataframe
+import idaes.logger as idaeslog
+
+from pyomo.environ import (
+    Param,
+    units as pyunits,
+    Set,
+)
+
+from idaes.core.util.exceptions import (
+    ConfigurationError,
+)
+
+__author__ = "Alejandro Garciadiego"
+
+
+# Set up logger
+_log = idaeslog.getLogger(__name__)
+
+
+@declare_process_block_class("Thickener")
+class ThickenerData(SeparatorData):
+    """
+    Thickener unit block for BSM2
+    """
+
+    CONFIG = SeparatorData.CONFIG()
+    CONFIG.outlet_list = ["underflow", "overflow"]
+    CONFIG.split_basis = SplittingType.componentFlow
+
+    def build(self):
+        """
+        Begin building model.
+        Args:
+            None
+        Returns:
+            None
+        """
+
+        # Call UnitModel.build to set up dynamics
+        super(ThickenerData, self).build()
+
+        if "underflow" and "overflow" not in self.config.outlet_list:
+            raise ConfigurationError(
+                "{} encountered unrecognised "
+                "outlet_list. This should not "
+                "occur - please use overflow "
+                "and underflow as outlets.".format(self.name)
+            )
+
+        self.p_thick = Param(
+            initialize=0.07,
+            units=pyunits.dimensionless,
+            mutable=True,
+            doc="Percentage of suspended solids in the underflow",
+        )
+
+        self.TSS_rem = Param(
+            initialize=0.98,
+            units=pyunits.dimensionless,
+            mutable=True,
+            doc="Percentage of suspended solids removed",
+        )
+
+        @self.Expression(self.flowsheet().time, doc="Suspended solid concentration")
+        def TSS(blk, t):
+            return 0.75 * (
+                blk.inlet.conc_mass_comp[t, "X_I"]
+                + blk.inlet.conc_mass_comp[t, "X_P"]
+                + blk.inlet.conc_mass_comp[t, "X_BH"]
+                + blk.inlet.conc_mass_comp[t, "X_BA"]
+                + blk.inlet.conc_mass_comp[t, "X_S"]
+            )
+
+        @self.Expression(self.flowsheet().time, doc="Thickening factor")
+        def f_thick(blk, t):
+            return blk.p_thick * (10 / (blk.TSS[t]))
+
+        @self.Expression(self.flowsheet().time, doc="Remove factor")
+        def f_q_du(blk, t):
+            return blk.TSS_rem / (pyunits.kg / pyunits.m**3) / 100 / blk.f_thick[t]
+
+        self.non_particulate_components = Set(
+            initialize=[
+                "S_I",
+                "S_S",
+                "S_O",
+                "S_NO",
+                "S_NH",
+                "S_ND",
+                "H2O",
+                "S_ALK",
+            ]
+        )
+
+        self.particulate_components = Set(
+            initialize=["X_I", "X_S", "X_P", "X_BH", "X_BA", "X_ND"]
+        )
+
+        @self.Constraint(
+            self.flowsheet().time,
+            self.particulate_components,
+            doc="particulate fraction",
+        )
+        def overflow_particulate_fraction(blk, t, i):
+            return blk.split_fraction[t, "overflow", i] == 1 - blk.TSS_rem
+
+        @self.Constraint(
+            self.flowsheet().time,
+            self.non_particulate_components,
+            doc="soluble fraction",
+        )
+        def non_particulate_components(blk, t, i):
+            return blk.split_fraction[t, "overflow", i] == 1 - blk.f_q_du[t]
+
+    def _get_performance_contents(self, time_point=0):
+        if hasattr(self, "split_fraction"):
+            var_dict = {}
+            for k in self.split_fraction.keys():
+                if k[0] == time_point:
+                    var_dict[f"Split Fraction [{str(k[1:])}]"] = self.split_fraction[k]
+            return {"vars": var_dict}
+        else:
+            return None
+
+    def _get_stream_table_contents(self, time_point=0):
+        outlet_list = self.create_outlet_list()
+
+        io_dict = {}
+        if self.config.mixed_state_block is None:
+            io_dict["Inlet"] = self.mixed_state
+        else:
+            io_dict["Inlet"] = self.config.mixed_state_block
+
+        for o in outlet_list:
+            io_dict[o] = getattr(self, o + "_state")
+
+        return create_stream_table_dataframe(io_dict, time_point=time_point)