diff --git a/watertap/core/tests/test_watertap_flowsheet_block.py b/watertap/core/tests/test_watertap_flowsheet_block.py
new file mode 100644
index 0000000000..6a53dcf3ab
--- /dev/null
+++ b/watertap/core/tests/test_watertap_flowsheet_block.py
@@ -0,0 +1,218 @@
+#################################################################################
+# WaterTAP Copyright (c) 2020-2024, 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/"
+#################################################################################
+
+__author__ = "Alexander V. Dudchenko"
+
+from idaes.models.unit_models import Feed, Product
+
+from pyomo.environ import (
+    TransformationFactory,
+    ConcreteModel,
+    Var,
+    Constraint,
+    units as pyunits,
+)
+from idaes.core import (
+    FlowsheetBlock,
+)
+
+from watertap.core.solvers import get_solver
+from pyomo.common.config import ConfigValue
+from watertap.property_models.seawater_prop_pack import SeawaterParameterBlock
+from watertap.core.watertap_flowsheet_block import WaterTapFlowsheetBlockData
+from idaes.core.util.model_statistics import degrees_of_freedom
+from idaes.core import (
+    declare_process_block_class,
+)
+from pyomo.environ import (
+    assert_optimal_termination,
+)
+import pytest
+
+from io import StringIO
+
+
+@declare_process_block_class("FlowsheetFeed")
+class FlowsheetFeedData(WaterTapFlowsheetBlockData):
+    """Test class for a simple feed unit within a WaterTapFlowsheetBlock"""
+
+    CONFIG = WaterTapFlowsheetBlockData.CONFIG()
+    CONFIG.declare(
+        "solute_concentration",
+        ConfigValue(
+            doc="Total Dissolved Solids (TDS) in the feed water (mg/L)",
+            default=35 * pyunits.g / pyunits.L,
+        ),
+    )
+    CONFIG.declare(
+        "feed_flow_rate",
+        ConfigValue(
+            doc="Feed flow rate (m^3/h)",
+            default=0.001 * pyunits.m**3 / pyunits.s,
+        ),
+    )
+
+    def build(self):
+        super().build()
+        self.feed = Feed(property_package=self.config.default_property_package)
+        self.solute_type = list(self.config.default_property_package.solute_set)[0]
+        self.feed.ph = Var(initialize=7.0, units=pyunits.dimensionless)
+        self.register_port("outlet", self.feed.outlet, {"pH": self.feed.ph})
+
+    def set_fixed_operation(self):
+        self.feed.ph.fix(8.5)
+
+        # fix feed flow and concentration based on config arguments
+        self.feed.properties[0].flow_vol_phase["Liq"].fix(self.config.feed_flow_rate)
+        self.feed.properties[0].conc_mass_phase_comp["Liq", self.solute_type].fix(
+            self.config.solute_concentration
+        )
+        self.feed.properties[0].temperature.fix(298.15)
+        self.feed.properties[0].pressure.fix(101325)
+
+        # make sure mass flows are not fixed
+        self.feed.properties[0].flow_mass_phase_comp["Liq", self.solute_type].unfix()
+        self.feed.properties[0].flow_mass_phase_comp["Liq", "H2O"].unfix()
+        assert degrees_of_freedom(self) == 0
+
+        # solve the block so we can get mass flows of solute and water
+        solver = get_solver()
+        results = solver.solve(self.feed, tee=False)
+        # ensure we terminate okay.
+        assert_optimal_termination(results)
+
+    def initialize_unit(self, **kwargs):
+        """custom initialize routine for the feed unit as we
+        are fixing feed concentration and flowrate rather than mass flowrates which will
+        throw an error if we just call self.feed.initialize() since feed mass flow and tds flow as unfixed
+
+        This is probably redundant if the model was already "fixed" as result will not change, but
+        do't want to assume that user has not change flow/concentration since calling set_fixed_operation (or fix_and_scale)
+        """
+        solver = get_solver()
+        result = solver.solve(self.feed, tee=False)
+
+        assert_optimal_termination(result)
+
+    def scale_before_initialization(self):
+        """apply standard scaling factors to feed unit model"""
+        self.config.default_property_package.set_default_scaling(
+            "flow_mass_phase_comp",
+            1 / self.feed.properties[0].flow_mass_phase_comp["Liq", "H2O"].value,
+            index=("Liq", "H2O"),
+        )
+        self.config.default_property_package.set_default_scaling(
+            "flow_mass_phase_comp",
+            1
+            / self.feed.properties[0]
+            .flow_mass_phase_comp["Liq", self.solute_type]
+            .value,
+            index=("Liq", self.solute_type),
+        )
+
+    def get_unit_name(self):
+        return "test feed unit"
+
+    def get_model_state_dict(self):
+        model_state = {
+            "Composition": {},
+            "Physical state": {},
+        }
+        model_state["Composition"]["Mass flow of H2O"] = self.feed.properties[
+            0
+        ].flow_mass_phase_comp["Liq", "H2O"]
+        model_state["Composition"]["Mass flow of TDS"] = self.feed.properties[
+            0
+        ].flow_mass_phase_comp["Liq", "TDS"]
+        for phase, ion in self.feed.properties[0].conc_mass_phase_comp:
+            if ion != "H2O":
+                model_state["Composition"][ion] = self.feed.properties[
+                    0
+                ].conc_mass_phase_comp[phase, ion]
+        model_state["Physical state"]["Temperature"] = self.feed.properties[
+            0
+        ].temperature
+        model_state["Physical state"]["Pressure"] = self.feed.properties[0].pressure
+        model_state["Physical state"]["Volumetric flowrate"] = self.feed.properties[
+            0
+        ].flow_vol_phase["Liq"]
+        return model_state
+
+
+@pytest.fixture
+def feed_flowsheet_fixture():
+    m = ConcreteModel()
+    m.fs = FlowsheetBlock(dynamic=False)
+    m.fs.seawater_props = SeawaterParameterBlock()
+    m.fs.feed = FlowsheetFeed(
+        default_property_package=m.fs.seawater_props,
+        solute_concentration=10 * pyunits.kg / pyunits.m**3,
+        feed_flow_rate=1000 * pyunits.m**3 / pyunits.s,
+    )
+
+    m.fs.product = Product(property_package=m.fs.seawater_props)
+
+    m.fs.feed.outlet.connect_to(m.fs.product.inlet)
+    TransformationFactory("network.expand_arcs").apply_to(m)
+
+    m.fs.feed.fix_and_scale()
+    m.fs.feed.initialize()
+    m.fs.product.initialize()
+    return m
+
+
+def test_feed_flowsheet_initialization(feed_flowsheet_fixture):
+    m = feed_flowsheet_fixture
+
+    # should be zero DOF after initialization
+    assert degrees_of_freedom(m) == 0
+
+    # this is simpel model, so product should have same mass flow as input
+    assert (
+        pytest.approx(
+            m.fs.product.flow_mass_phase_comp[0, "Liq", "H2O"].value, rel=1e-5
+        )
+        == m.fs.feed.feed.properties[0].flow_mass_phase_comp["Liq", "H2O"].value
+    )
+    assert (
+        pytest.approx(
+            m.fs.product.flow_mass_phase_comp[0, "Liq", "TDS"].value, rel=1e-5
+        )
+        == m.fs.feed.feed.properties[0].flow_mass_phase_comp["Liq", "TDS"].value
+    )
+
+
+def test_report(feed_flowsheet_fixture):
+    m = feed_flowsheet_fixture
+    os = StringIO()
+    m.fs.feed.report(ostream=os)
+    result = """
+------------------------------------------------------------------------------------
+    test feed unit state
+
+    Composition:
+    Key              : Value      : Units   : Fixed : Bounds
+    Mass flow of H2O : 9.9448e+05 :    kg/s : False : (0.0, None)
+    Mass flow of TDS :     10000. :    kg/s : False : (0.0, None)
+                 TDS :     10.000 : kg/m**3 :  True : (0.0, 1000000.0)
+
+    Physical state:
+    Key                 : Value      : Units  : Fixed : Bounds
+                Pressure : 1.0132e+05 :     Pa :  True : (1000.0, 50000000.0)
+            Temperature :     298.15 :      K :  True : (273.15, 1000)
+    Volumetric flowrate :     1000.0 : m**3/s :  True : (0.0, None)
+
+------------------------------------------------------------------------------------
+"""
+    print(os.getvalue().replace(" ", ""))
+    # testing with out spaces, as they are hard to control in the report output
+    assert os.getvalue().replace(" ", "") == result.replace(" ", "")
diff --git a/watertap/core/util/connections.py b/watertap/core/util/connections.py
new file mode 100644
index 0000000000..b1299c98eb
--- /dev/null
+++ b/watertap/core/util/connections.py
@@ -0,0 +1,186 @@
+#################################################################################
+# WaterTAP Copyright (c) 2020-2024, 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/"
+#################################################################################
+
+__author__ = "Alexander V. Dudchenko"
+
+from idaes.core.util.initialization import propagate_state
+from pyomo.network import Arc
+from pyomo.environ import (
+    Constraint,
+)
+import idaes.core.util.scaling as iscale
+from pyomo.network import Port
+import idaes.logger as idaeslog
+
+# Set up logger
+_log = idaeslog.getLogger(__name__)
+
+
+class PortContainer:
+    """Port container will aggregate ports and variables
+    user wants to connect between unit models
+    Although, the variable can be registered to the port directly, it can lead
+    to issues when using UI visualization tools"""
+
+    def __init__(self, name, port, var_dict, unit_block_reference):
+        self.name = name
+        self.port = port
+        self.unit_block_reference = unit_block_reference
+        if var_dict != None and isinstance(var_dict, dict) == False:
+            raise TypeError(
+                "Var dict must be a dictionary with structure {'variable name': pyomo.var}"
+            )
+        self.var_dict = var_dict
+
+    def connect_to(self, inlet):
+        """connect current port to provided port, registering
+        the generated connection container using provided function"""
+        connection = ConnectionContainer(self, inlet)
+        if hasattr(self.unit_block_reference, "register_outlet_connection"):
+            self.unit_block_reference.register_outlet_connection(connection)
+        else:
+            self.connection = connection
+
+    def fix(self):
+        """this will fix the port and all variables in the var_dict"""
+        self.port.fix()
+        if self.var_dict != None:
+            for var, obj in self.var_dict.items():
+                obj.fix()
+
+    def unfix(self):
+        """this will unfix the port and all variables in the var_dict"""
+        self.port.unfix()
+        if self.var_dict != None:
+            for var, obj in self.var_dict.items():
+                obj.unfix()
+
+
+class ConnectionContainer:
+    """This serves as container for any ports and equality constraints for use between unit
+    models"""
+
+    def __init__(self, outlet, inlet):
+        self.registered_equality_constraints = []
+        if not isinstance(outlet, (PortContainer, Port)) or not isinstance(
+            inlet, (PortContainer, Port)
+        ):
+            raise TypeError("Provided outlet and inlet must be PortContainer objects")
+        self.build_arc(outlet, inlet)
+        self.build_constraints(outlet, inlet)
+
+    def get_port(self, possible_port_object):
+        """Return port for port container or port object"""
+        if isinstance(possible_port_object, PortContainer):
+            return possible_port_object.port
+        elif isinstance(possible_port_object, Port):
+            return possible_port_object
+        else:
+            raise TypeError("Provided object is not a PortContainer or Port")
+
+    def get_port_unit(self, possible_port_object):
+        """Return unit for port container or port object"""
+        if isinstance(possible_port_object, PortContainer):
+            return possible_port_object.unit_block_reference
+        elif isinstance(possible_port_object, Port):
+            return possible_port_object.parent_block().name
+        else:
+            raise TypeError("Provided object is not a PortContainer or Port")
+
+    def build_arc(self, outlet, inlet):
+        """
+        builds a standard arc while naming it with outlet and inlet name, this should
+        be a unique pair always (e.g. for a unit model you should only havbe single outlet-> inlet conenction)
+        """
+        arc = Arc(source=self.get_port(outlet), destination=self.get_port(inlet))
+
+        def get_safe_name(name):
+            return name.replace(".", "_").replace("-", "_")
+
+        arc_name = f"{get_safe_name(outlet.name)}_to_{get_safe_name(inlet.name)}"
+        outlet.unit_block_reference.add_component(
+            arc_name,
+            arc,
+        )
+
+        # find it, if it was not created correctly, we will get an error
+        self.registered_arc = outlet.unit_block_reference.find_component(arc_name)
+        if self.registered_arc is None:
+            raise ValueError(f"Arc was not created correctly for {arc_name}")
+        self.unit_connection = f"{self.get_port_unit(outlet)}.{get_safe_name(outlet.name)}_to_{self.get_port_unit(inlet)}.{get_safe_name(inlet.name)}"
+        _log.info("Created arc connection: %s", arc_name)
+
+    def build_constraints(self, outlet, inlet):
+        """
+        builds equality constraints for provided variables and scales them.
+        We ensure that outlet and inlet vars are the same.
+        """
+        if isinstance(outlet, PortContainer) and isinstance(inlet, PortContainer):
+            none_test = [d != None for d in [outlet.var_dict, inlet.var_dict]]
+
+            if all(none_test):
+                if set(outlet.var_dict) != set(inlet.var_dict):
+                    raise KeyError(
+                        f"Provided inlet keys: {outlet.var_dict} do not match outlet keys: {inlet.var_dict}"
+                    )
+
+                for outlet_key in outlet.var_dict:
+                    # Do not create constraint if the variable is the same)
+                    if outlet.var_dict[outlet_key] is inlet.var_dict[outlet_key]:
+                        pass
+                    else:
+                        # create equality constraint between outlet and inlet var dicts
+                        outlet.unit_block_reference.add_component(
+                            f"eq_{outlet_key}_{outlet.name}_to_{inlet.name}",
+                            Constraint(
+                                expr=outlet.var_dict[outlet_key]
+                                == inlet.var_dict[outlet_key]
+                            ),
+                        )
+                        constraint = outlet.unit_block_reference.find_component(
+                            f"eq_{outlet_key}_{outlet.name}_to_{inlet.name}"
+                        )
+                        # ensure we register it for propagation later on
+                        self.registered_equality_constraints.append(
+                            (
+                                constraint,
+                                outlet.var_dict[outlet_key],
+                                inlet.var_dict[outlet_key],
+                            )
+                        )
+                        _log.info(
+                            "Created arc constraint: %s",
+                            f"eq_{outlet_key}_{outlet.name}_to_{inlet.name}",
+                        )
+                        # scale the constraint
+                        sf = iscale.get_scaling_factor(outlet.var_dict[outlet_key])
+                        if sf != None:
+                            iscale.constraint_scaling_transform(
+                                constraint,
+                                sf,
+                            )
+
+    def propagate(self):
+        """this should prop any arcs and also ensure all equality constraints are satisfied"""
+        _log.info("Propagating connection: %s", self.unit_connection)
+        propagate_state(self.registered_arc)
+        self.propagate_equality_constraints()
+
+    def propagate_equality_constraints(self):
+        """this will ensure that all equality constraints are satisfied by setting the inlet var to the outlet var"""
+        for (
+            constraint,
+            outlet_var,
+            inlet_var,
+        ) in self.registered_equality_constraints:
+            # set the inlet var to the outlet var value
+            inlet_var.value = outlet_var.value
diff --git a/watertap/core/util/report.py b/watertap/core/util/report.py
new file mode 100644
index 0000000000..ce6fa8c498
--- /dev/null
+++ b/watertap/core/util/report.py
@@ -0,0 +1,137 @@
+#################################################################################
+# WaterTAP Copyright (c) 2020-2024, 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/"
+#################################################################################
+
+__author__ = "Alexander V. Dudchenko"
+
+from pyomo.common.formatting import tabular_writer
+from pyomo.environ import (
+    value,
+    units as pyunits,
+)
+import sys
+from idaes.core.util.units_of_measurement import report_quantity
+
+
+def build_report_table(
+    unit_name, data_dict, ostream=None, prefix="", use_default_units=False
+):
+    """Builds a report table for the unit model using supplied unit name and data dict
+    that contains variables to display
+
+    # Example implementation:
+    #
+    # model_state = {
+    #     "Composition": {},
+    #     "Physical state": {},
+    # }
+    # for phase, ion in self.feed.properties[0].conc_mass_phase_comp:
+    #     model_state["Composition"][ion] = self.feed.properties[
+    #         0
+    #     ].conc_mass_phase_comp[phase, ion]
+    # model_state["Physical state"]["pH"] = self.feed.pH
+    # model_state["Physical state"]["Temperature"] = self.feed.properties[
+    #     0
+    # ].temperature
+    # model_state["Physical state"]["Pressure"] = self.feed.properties[0].pressure
+    #
+    #
+    #Output should look like:
+    #
+    # ------------------------------------------------------------------------------------
+    #     fs.feed state
+    #
+    #     Composition:
+    #     Key    : Value     : Units                 : Fixed : Bounds
+    #      IONA :   0.25800 : kilogram / meter ** 3 :  True : (0, 2000.0)
+    #      IONB :   0.87000 : kilogram / meter ** 3 :  True : (0, 2000.0)
+    #
+    #     Physical state:
+    #     Key         : Value      : Units         : Fixed : Bounds
+    #        Pressure : 1.0000e+05 :        pascal :  True : (100000.0, None)
+    #     Temperature :     293.15 :        kelvin :  True : (273.15, 373.15)
+    #              pH :     7.0700 : dimensionless :  True : (None, None)
+    #
+    # ------------------------------------------------------------------------------------
+
+    Args:
+        unit_name: Name of the unit
+        data_dict: Dictionary containing data to report
+        ostream: Output stream for the report (default: sys.stdout)
+        prefix: String prefix for formatting the report
+        use_default_units: Boolean to indicate if default units should be used
+    """
+
+    def _get_fixed_state(v):
+        """Get the fixed state of a variable, if none exists then return N/A"""
+        try:
+            return v.fixed
+        except AttributeError:
+            return "N/A"
+
+    def _get_bounds(v):
+        """Get the bounds of a variable, if none exists then return N/A"""
+        try:
+            return v.bounds
+        except AttributeError:
+            return "N/A"
+
+    def get_values(k, v):
+        """Get the values of a variable,  dimensions, fixed state and bounds"""
+        if isinstance(v, int):
+            return [
+                v,
+                "dimensionless",
+                "N/A",
+                "N/A",
+            ]
+        elif isinstance(v, float):
+            return [
+                "{:#.5g}".format(v),
+                "dimensionless",
+                "N/A",
+                "N/A",
+            ]
+        elif use_default_units:
+            return [
+                "{:#.5g}".format(report_quantity(v).m),
+                report_quantity(v).u,
+                _get_fixed_state(v),
+                _get_bounds(v),
+            ]
+        else:
+            return [
+                "{:#.5g}".format(value(v)),
+                pyunits.get_units(v),
+                _get_fixed_state(v),
+                _get_bounds(v),
+            ]
+
+    if ostream is None:
+        ostream = sys.stdout
+    max_str_length = 84
+    tab = " " * 4
+    ostream.write("\n" + "-" * max_str_length + "\n")
+    ostream.write(f"{prefix}{tab}{unit_name} state")
+    ostream.write("\n" * 2)
+    for key, sub_data in data_dict.items():
+
+        ostream.write(f"{prefix}{tab}{key}: \n")
+        tabular_writer(
+            ostream,
+            prefix + tab,
+            ((k, v) for k, v in sub_data.items()),
+            ("Value", "Units", "Fixed", "Bounds"),
+            get_values,
+        )
+        ostream.write(f"\n")
+
+    ostream.write("-" * max_str_length + "\n")
diff --git a/watertap/core/util/tests/test_connections.py b/watertap/core/util/tests/test_connections.py
new file mode 100644
index 0000000000..502804904a
--- /dev/null
+++ b/watertap/core/util/tests/test_connections.py
@@ -0,0 +1,168 @@
+#################################################################################
+# WaterTAP Copyright (c) 2020-2024, 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/"
+#################################################################################
+
+__author__ = "Alexander V. Dudchenko"
+
+from watertap.core.util.connections import ConnectionContainer, PortContainer
+
+from idaes.models.unit_models import Feed, Product
+
+from pyomo.environ import (
+    TransformationFactory,
+    ConcreteModel,
+    Var,
+    Constraint,
+    units as pyunits,
+)
+from idaes.core import (
+    FlowsheetBlock,
+)
+from pyomo.network import Arc
+
+from watertap.property_models.seawater_prop_pack import SeawaterParameterBlock
+
+import pytest
+
+from idaes.core.util.model_statistics import degrees_of_freedom
+
+
+@pytest.fixture
+def build_test_model():
+    """Test that ConnectionContainer can be created without error."""
+    m = ConcreteModel()
+    m.fs = FlowsheetBlock(dynamic=False)
+    m.fs.seawater_props = SeawaterParameterBlock()
+    m.fs.feed = Feed(property_package=m.fs.seawater_props)
+    m.fs.feed.ph = Var(initialize=7.0, units=pyunits.dimensionless)
+    m.fs.product = Product(property_package=m.fs.seawater_props)
+    m.fs.product.ph = Var(initialize=7.0, units=pyunits.dimensionless)
+    # Create PortContainers for the outlet of unit1 and inlet of unit2
+    m.fs.outlet_port_container = PortContainer(
+        name="outlet1",
+        port=m.fs.feed.outlet,
+        var_dict={"pH": m.fs.feed.ph},
+        unit_block_reference=m.fs.feed,
+    )
+    m.fs.inlet_port_container = PortContainer(
+        name="inlet2",
+        port=m.fs.product.inlet,
+        var_dict={"pH": m.fs.product.ph},
+        unit_block_reference=m.fs.product,
+    )
+    m.fs.feed.ph.fix(8.5)
+    m.fs.feed.flow_mass_phase_comp[0, "Liq", "H2O"].fix(1000)
+    m.fs.feed.flow_mass_phase_comp[0, "Liq", "TDS"].fix(10)
+    m.fs.feed.properties[0].temperature.fix(298.15)
+    m.fs.feed.properties[0].pressure.fix(101325)
+    return m
+
+
+def test_connect_with_port_container(build_test_model):
+    m = build_test_model
+    # Connect the two PortContainers
+    m.fs.outlet_port_container.connect_to(m.fs.inlet_port_container)
+
+    assert m.fs.feed.find_component("outlet1_to_inlet2") is not None
+    assert isinstance(m.fs.feed.find_component("outlet1_to_inlet2"), Arc)
+    assert m.fs.feed.find_component("eq_pH_outlet1_to_inlet2") is not None
+    assert isinstance(m.fs.feed.find_component("eq_pH_outlet1_to_inlet2"), Constraint)
+
+
+def test_connect_with_port_registration(build_test_model):
+    m = build_test_model
+
+    # Connect the two PortContainers
+    def register_outlet_connection(connection):
+        m.fs.feed.registered_connections = []
+        m.fs.feed.registered_connections.append(connection)
+        print(m.fs.feed.registered_connections)
+
+    m.fs.feed.register_outlet_connection = register_outlet_connection
+    m.fs.outlet_port_container.connect_to(m.fs.inlet_port_container)
+
+    assert isinstance(m.fs.feed.registered_connections[0], ConnectionContainer)
+    assert isinstance(m.fs.feed.registered_connections[0].registered_arc, Arc)
+    assert isinstance(
+        m.fs.feed.registered_connections[0].registered_equality_constraints[0][0],
+        Constraint,
+    )
+    assert (
+        m.fs.feed.registered_connections[0].registered_equality_constraints[0][1]
+        == m.fs.feed.ph
+    )
+
+    assert (
+        m.fs.feed.registered_connections[0].registered_equality_constraints[0][2]
+        == m.fs.product.ph
+    )
+
+
+def test_connect_with_port(build_test_model):
+    m = build_test_model
+    # Connect the two PortContainers
+    m.fs.outlet_port_container.connect_to(m.fs.product.inlet)
+
+    assert m.fs.feed.find_component("outlet1_to_fs_product_inlet") is not None
+    assert isinstance(m.fs.feed.find_component("outlet1_to_fs_product_inlet"), Arc)
+    assert m.fs.feed.find_component("eq_pH_outlet1_to_fs_product_inlet") is None
+
+
+def test_network_propagation(build_test_model):
+    m = build_test_model
+    # Connect the two PortContainers
+    m.fs.outlet_port_container.connect_to(m.fs.inlet_port_container)
+    TransformationFactory("network.expand_arcs").apply_to(m)
+    # Set pH on feed unit
+
+    # Propagate values through the network
+    m.fs.outlet_port_container.connection.propagate()
+
+    assert pytest.approx(m.fs.product.ph.value, rel=1e-5) == 8.5
+    assert (
+        pytest.approx(
+            m.fs.product.flow_mass_phase_comp[0, "Liq", "H2O"].value, rel=1e-5
+        )
+        == 1000
+    )
+    assert (
+        pytest.approx(
+            m.fs.product.flow_mass_phase_comp[0, "Liq", "TDS"].value, rel=1e-5
+        )
+        == 10
+    )
+    m.fs.inlet_port_container.fix()
+    m.fs.product.properties[0].display()
+    assert degrees_of_freedom(m) == -5
+    assert degrees_of_freedom(m.fs.product) == 0
+
+
+def test_network_propagation_with_normal_port_destination(build_test_model):
+    m = build_test_model
+    # Connect the two PortContainers
+    m.fs.outlet_port_container.connect_to(m.fs.product.inlet)
+    TransformationFactory("network.expand_arcs").apply_to(m)
+
+    m.fs.outlet_port_container.connection.propagate()
+    assert degrees_of_freedom(m) == 0
+    assert pytest.approx(m.fs.product.ph.value, rel=1e-5) == 7
+    assert (
+        pytest.approx(
+            m.fs.product.flow_mass_phase_comp[0, "Liq", "H2O"].value, rel=1e-5
+        )
+        == 1000
+    )
+    assert (
+        pytest.approx(
+            m.fs.product.flow_mass_phase_comp[0, "Liq", "TDS"].value, rel=1e-5
+        )
+        == 10
+    )
diff --git a/watertap/core/util/tests/test_report.py b/watertap/core/util/tests/test_report.py
new file mode 100644
index 0000000000..02e5b8a1e0
--- /dev/null
+++ b/watertap/core/util/tests/test_report.py
@@ -0,0 +1,63 @@
+#################################################################################
+# WaterTAP Copyright (c) 2020-2024, 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/"
+#################################################################################
+
+__author__ = "Alexander V. Dudchenko"
+
+
+from watertap.core.util import report
+
+from pyomo.environ import ConcreteModel, Var, units as pyunits
+from io import StringIO
+
+
+def test_report_table():
+    """Test that the report table can be built without error."""
+    # Create a dummy model with a variable
+
+    m = ConcreteModel()
+    m.x = Var(["t1", "t2", "t3"], initialize=1.0, units=pyunits.m)
+    os = StringIO()
+    # Build the report table
+    # should look like this:
+    # ------------------------------------------------------------------------------------
+    #     test_var state
+
+    #     test_value:
+    #     Key : Value  : Units : Fixed : Bounds
+    #      t1 : 1.0000 :     m : False : (None, None)
+    #      t2 : 1.0000 :     m : False : (None, None)
+    #      t3 : 1.0000 :     m : False : (None, None)
+
+    # ------------------------------------------------------------------------------------
+    report.build_report_table("test_var", {"test_value": m.x}, os)
+    result = "\n------------------------------------------------------------------------------------\n    test_var state\n\n    test_value: \n    Key : Value  : Units : Fixed : Bounds\n     t1 : 1.0000 :     m : False : (None, None)\n     t2 : 1.0000 :     m : False : (None, None)\n     t3 : 1.0000 :     m : False : (None, None)\n\n------------------------------------------------------------------------------------\n"
+    print(os.getvalue())
+    assert os.getvalue() == result, "Report table does not match expected output"
+    os = StringIO()
+    # Build the report table
+    report.build_report_table(
+        "test_var", {"test_value": m.x}, os, use_default_units=True
+    )
+    # should look like this:
+    # ------------------------------------------------------------------------------------
+    #     test_var state
+
+    #     test_value:
+    #     Key : Value  : Units : Fixed : Bounds
+    #      t1 : 1.0000 :     meter : False : (None, None)
+    #      t2 : 1.0000 :     meter : False : (None, None)
+    #      t3 : 1.0000 :     meter : False : (None, None)
+
+    # ------------------------------------------------------------------------------------
+    result = "\n------------------------------------------------------------------------------------\n    test_var state\n\n    test_value: \n    Key : Value  : Units : Fixed : Bounds\n     t1 : 1.0000 : meter : False : (None, None)\n     t2 : 1.0000 : meter : False : (None, None)\n     t3 : 1.0000 : meter : False : (None, None)\n\n------------------------------------------------------------------------------------\n"
+    print(os.getvalue())
+    assert os.getvalue() == result, "Report table does not match expected output"
diff --git a/watertap/core/watertap_flowsheet_block.py b/watertap/core/watertap_flowsheet_block.py
new file mode 100644
index 0000000000..0f617d760d
--- /dev/null
+++ b/watertap/core/watertap_flowsheet_block.py
@@ -0,0 +1,177 @@
+#################################################################################
+# WaterTAP Copyright (c) 2020-2024, 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/"
+#################################################################################
+
+__author__ = "Alexander V. Dudchenko"
+
+from idaes.core.base.flowsheet_model import FlowsheetBlockData
+
+from idaes.core import (
+    declare_process_block_class,
+)
+from pyomo.common.config import ConfigValue
+
+from watertap.core.util.connections import (
+    PortContainer,
+    ConnectionContainer,
+)
+from watertap.core.util.report import (
+    build_report_table,
+)
+
+
+@declare_process_block_class("WaterTapFlowsheetBlock")
+class WaterTapFlowsheetBlockData(FlowsheetBlockData):
+    CONFIG = FlowsheetBlockData.CONFIG()
+    CONFIG.declare(
+        "default_costing_package",
+        ConfigValue(
+            default=None,
+            description="defines default costing package",
+            doc="""
+                defines default costing package
+            """,
+        ),
+    )
+    CONFIG.declare(
+        "default_costing_package_kwargs",
+        ConfigValue(
+            default={},
+            description="kwargs to pass into the Costing unit block",
+            doc="""
+                kwargs to pass into the Costing unit block,
+            """,
+        ),
+    )
+
+    def build(self):
+        self.outlet_connections = []
+        super().build()
+
+    def fix_and_scale(self):
+        self.set_fixed_operation()
+        self.scale_before_initialization()
+
+    def initialize_unit(self, **kwargs):
+        """Developer should implement an initialize routine for their flowsheet model"""
+
+    def initialize(self, **kwargs):
+        """routine to initialize a unit and propagate its connections"""
+        self.initialize_unit()
+        self.propagate_outlets()
+
+    def propagate_outlets(self):
+        """propagates registered outlet connections"""
+        for outlet in self.outlet_connections:
+            outlet.propagate()
+
+    def set_fixed_operation(self, **kwargs):
+        """Developer should implement a routine to fix unit operation for initialization and 0DOF solving"""
+
+    def set_optimization_operation(self, **kwargs):
+        """Developer should implement a routine to unfix variables for unit to perform optimization"""
+
+    def scale_before_initialization(self, **kwargs):
+        """Developer should implement scaling function to scale unit using
+        default values before initialization routine is ran"""
+
+    def scale_post_initialization(self, **kwargs):
+        """Developer should implement scaling function to scale unit after initialization routine is ran"""
+
+    def register_port(self, name, port=None, var_list=None):
+        """Registers a port for the flowsheet unit, including variables that should
+        be connected through equality constraints
+
+        Args:
+        name - Name of the port
+        port - port from a unit model, can be none
+        var_list - list of variables that should be connected through equality constraints
+        """
+        # create our var on the block so we can reference port and variables
+        setattr(
+            self,
+            name,
+            PortContainer(name, port, var_list, self),
+        )
+
+    def register_outlet_connection(self, connection):
+        """registers outlet connections to enable automatic propagation"""
+        if isinstance(connection, ConnectionContainer):
+            self.outlet_connections.append(connection)
+        else:
+            raise TypeError("Outlet connection must be a ConnectionContainer")
+
+    def get_unit_name(self):
+        """returns the name of the unit block, developer can overwrite this
+        to provide more descriptive name"""
+        return self.name
+
+    def get_model_state_dict(self):
+        """
+        Developer should over write this function so it generates a dict of model
+        variables to display.
+
+        # Example implementation:
+        # def get_model_state_dict(self):
+        #     model_state = {
+        #         "Composition": {},
+        #         "Physical state": {},
+        #     }
+        #     for phase, ion in self.feed.properties[0].conc_mass_phase_comp:
+        #         model_state["Composition"][ion] = self.feed.properties[
+        #             0
+        #         ].conc_mass_phase_comp[phase, ion]
+        #     model_state["Physical state"]["pH"] = self.feed.pH
+        #     model_state["Physical state"]["Temperature"] = self.feed.properties[
+        #         0
+        #     ].temperature
+        #     model_state["Physical state"]["Pressure"] = self.feed.properties[0].pressure
+        #
+        #     return model_state
+        #
+        #Output should look like:
+        #
+        # ------------------------------------------------------------------------------------
+        #     fs.feed state
+        #
+        #     Composition:
+        #     Key    : Value     : Units                 : Fixed : Bounds
+        #       IONA :   0.25800 : kilogram / meter ** 3 :  True : (0, 2000.0)
+        #       IONB :   0.87000 : kilogram / meter ** 3 :  True : (0, 2000.0)
+        #     Physical state:
+        #     Key         : Value      : Units         : Fixed : Bounds
+        #        Pressure : 1.0000e+05 :        pascal :  True : (100000.0, None)
+        #     Temperature :     293.15 :        kelvin :  True : (273.15, 373.15)
+        #              pH :     7.0700 : dimensionless :  True : (None, None)
+        #
+        # ------------------------------------------------------------------------------------
+        """
+        return None, None
+
+    def _get_stream_table_contents(self, time_point=0):
+        """override default as developer should manually define model state"""
+        if self.get_model_state_dict() is None:
+            return super()._get_stream_table_contents(time_point)
+        else:
+            return None
+
+    def report(
+        self, time_point=0, dof=False, ostream=None, prefix="", use_default_units=False
+    ):
+        defined_variables = self.get_model_state_dict()
+        unit_name = self.get_unit_name()
+
+        if unit_name is not None:
+            build_report_table(
+                unit_name, defined_variables, ostream, prefix, use_default_units
+            )
+        if unit_name is None:
+            super().report(time_point, dof, ostream, prefix)