diff --git a/watertap/costing/units/electroNP.py b/watertap/costing/units/electroNP.py
new file mode 100644
index 0000000000..6a4003d488
--- /dev/null
+++ b/watertap/costing/units/electroNP.py
@@ -0,0 +1,89 @@
+#################################################################################
+# 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/"
+#################################################################################
+
+import pyomo.environ as pyo
+from ..util import (
+    register_costing_parameter_block,
+    make_capital_cost_var,
+)
+
+
+def build_electroNP_cost_param_block(blk):
+
+    blk.HRT = pyo.Var(
+        initialize=1.3333,
+        doc="Hydraulic retention time",
+        units=pyo.units.hr,
+    )
+    blk.sizing_cost = pyo.Var(
+        initialize=1.25,
+        doc="Reactor sizing cost",
+        units=pyo.units.USD_2020 / pyo.units.m**3,
+    )
+
+
+@register_costing_parameter_block(
+    build_rule=build_electroNP_cost_param_block,
+    parameter_block_name="electroNP",
+)
+def cost_electroNP(blk, cost_electricity_flow=True, cost_MgCl2_flow=True):
+    """
+    ElectroNP costing method
+    """
+    cost_electroNP_capital(
+        blk,
+        blk.costing_package.electroNP.HRT,
+        blk.costing_package.electroNP.sizing_cost,
+    )
+
+    t0 = blk.flowsheet().time.first()
+    if cost_electricity_flow:
+        blk.costing_package.cost_flow(
+            pyo.units.convert(
+                blk.unit_model.electricity[t0],
+                to_units=pyo.units.kW,
+            ),
+            "electricity",
+        )
+
+    if cost_MgCl2_flow:
+        blk.costing_package.cost_flow(
+            pyo.units.convert(
+                blk.unit_model.MgCl2_flowrate[t0],
+                to_units=pyo.units.kg / pyo.units.hr,
+            ),
+            "magnesium chloride",
+        )
+
+
+def cost_electroNP_capital(blk, HRT, sizing_cost):
+    """
+    Generic function for costing an ElectroNP system.
+    """
+    make_capital_cost_var(blk)
+
+    blk.HRT = pyo.Expression(expr=HRT)
+    blk.sizing_cost = pyo.Expression(expr=sizing_cost)
+
+    flow_in = pyo.units.convert(
+        blk.unit_model.properties_in[0].flow_vol,
+        to_units=pyo.units.m**3 / pyo.units.hr,
+    )
+
+    print(f"base_currency: {blk.costing_package.base_currency}")
+    blk.capital_cost_constraint = pyo.Constraint(
+        expr=blk.capital_cost
+        == pyo.units.convert(
+            blk.HRT * flow_in * blk.sizing_cost,
+            to_units=blk.costing_package.base_currency,
+        )
+    )
diff --git a/watertap/costing/watertap_costing_package.py b/watertap/costing/watertap_costing_package.py
index fde27cba5b..f568fd3452 100644
--- a/watertap/costing/watertap_costing_package.py
+++ b/watertap/costing/watertap_costing_package.py
@@ -37,6 +37,7 @@
     EnergyRecoveryDevice,
     Electrodialysis0D,
     Electrodialysis1D,
+    ElectroNPZO,
     IonExchange0D,
     GAC,
 )
@@ -55,6 +56,7 @@
 from .units.pump import cost_pump
 from .units.reverse_osmosis import cost_reverse_osmosis
 from .units.uv_aop import cost_uv_aop
+from .units.electroNP import cost_electroNP
 
 
 class _DefinedFlowsDict(MutableMapping, dict):
@@ -95,6 +97,7 @@ class WaterTAPCostingData(FlowsheetCostingBlockData):
         Ultraviolet0D: cost_uv_aop,
         Electrodialysis0D: cost_electrodialysis,
         Electrodialysis1D: cost_electrodialysis,
+        ElectroNPZO: cost_electroNP,
         IonExchange0D: cost_ion_exchange,
         GAC: cost_gac,
     }
@@ -159,6 +162,14 @@ def build_global_params(self):
             units=pyo.units.kg / pyo.units.kWh,
         )
 
+        self.magnesium_chloride_cost = pyo.Param(
+            mutable=True,
+            initialize=0.0786,
+            doc="Magnesium chloride cost",
+            units=pyo.units.USD_2020 / pyo.units.kg,
+        )
+        self.add_defined_flow("magnesium chloride", self.magnesium_chloride_cost)
+
         # fix the parameters
         self.fix_all_vars()
 
diff --git a/watertap/unit_models/__init__.py b/watertap/unit_models/__init__.py
index 3bd828e717..ed3d867e8d 100644
--- a/watertap/unit_models/__init__.py
+++ b/watertap/unit_models/__init__.py
@@ -26,3 +26,4 @@
 from .ion_exchange_0D import IonExchange0D
 from .thickener import Thickener
 from .dewatering import DewateringUnit
+from .electroNP_ZO import ElectroNPZO
diff --git a/watertap/unit_models/electroNP_ZO.py b/watertap/unit_models/electroNP_ZO.py
new file mode 100644
index 0000000000..02fc712828
--- /dev/null
+++ b/watertap/unit_models/electroNP_ZO.py
@@ -0,0 +1,508 @@
+#################################################################################
+# 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/"
+#################################################################################
+
+# Import Pyomo libraries
+from pyomo.environ import (
+    Var,
+    check_optimal_termination,
+    Param,
+    Suffix,
+    NonNegativeReals,
+    Reference,
+    units as pyunits,
+)
+from pyomo.common.config import Bool, ConfigBlock, ConfigValue, In
+
+# Import IDAES cores
+from idaes.core import (
+    declare_process_block_class,
+    MaterialBalanceType,
+    EnergyBalanceType,
+    MomentumBalanceType,
+    UnitModelBlockData,
+    useDefault,
+    MaterialFlowBasis,
+)
+from idaes.core.solvers import get_solver
+from idaes.core.util.tables import create_stream_table_dataframe
+from idaes.core.util.config import is_physical_parameter_block
+from idaes.core.util.exceptions import ConfigurationError, InitializationError
+import idaes.core.util.scaling as iscale
+import idaes.logger as idaeslog
+
+from watertap.core import ControlVolume0DBlock, InitializationMixin
+
+__author__ = "Chenyu Wang"
+
+_log = idaeslog.getLogger(__name__)
+
+
+@declare_process_block_class("ElectroNPZO")
+class ElectroNPZOData(InitializationMixin, UnitModelBlockData):
+    """
+    Zero order electrochemical nutrient removal (ElectroNP) model based on specified water flux and ion rejection.
+    """
+
+    CONFIG = ConfigBlock()
+
+    CONFIG.declare(
+        "dynamic",
+        ConfigValue(
+            domain=In([False]),
+            default=False,
+            description="Dynamic model flag - must be False",
+            doc="""Indicates whether this model will be dynamic or not,
+    **default** = False. NF units do not support dynamic
+    behavior.""",
+        ),
+    )
+    CONFIG.declare(
+        "has_holdup",
+        ConfigValue(
+            default=False,
+            domain=In([False]),
+            description="Holdup construction flag - must be False",
+            doc="""Indicates whether holdup terms should be constructed or not.
+    **default** - False. NF units do not have defined volume, thus
+    this must be False.""",
+        ),
+    )
+    CONFIG.declare(
+        "property_package",
+        ConfigValue(
+            default=useDefault,
+            domain=is_physical_parameter_block,
+            description="Property package to use for control volume",
+            doc="""Property parameter object used to define property calculations,
+    **default** - useDefault.
+    **Valid values:** {
+    **useDefault** - use default package from parent model or flowsheet,
+    **PhysicalParameterObject** - a PhysicalParameterBlock object.}""",
+        ),
+    )
+    CONFIG.declare(
+        "property_package_args",
+        ConfigBlock(
+            implicit=True,
+            description="Arguments to use for constructing property packages",
+            doc="""A ConfigBlock with arguments to be passed to a property block(s)
+    and used when constructing these,
+    **default** - None.
+    **Valid values:** {
+    see property package for documentation.}""",
+        ),
+    )
+
+    def _process_config(self):
+        if len(self.config.property_package.solvent_set) > 1:
+            raise ConfigurationError(
+                "ElectroNP model only supports one solvent component,"
+                "the provided property package has specified {} solvent components".format(
+                    len(self.config.property_package.solvent_set)
+                )
+            )
+
+        if len(self.config.property_package.solvent_set) == 0:
+            raise ConfigurationError(
+                "The ElectroNP model was expecting a solvent and did not receive it."
+            )
+
+        if (
+            len(self.config.property_package.solute_set) == 0
+            and len(self.config.property_package.ion_set) == 0
+        ):
+            raise ConfigurationError(
+                "The ElectroNP model was expecting at least one solute or ion and did not receive any."
+            )
+
+    def build(self):
+        # Call UnitModel.build to setup dynamics
+        super().build()
+
+        self.scaling_factor = Suffix(direction=Suffix.EXPORT)
+
+        units_meta = self.config.property_package.get_metadata().get_derived_units
+
+        # Check configs for errors
+        self._process_config()
+
+        # Create state blocks for inlet and outlets
+        tmp_dict = dict(**self.config.property_package_args)
+        tmp_dict["has_phase_equilibrium"] = False
+        tmp_dict["defined_state"] = True
+
+        self.properties_in = self.config.property_package.build_state_block(
+            self.flowsheet().time, doc="Material properties at inlet", **tmp_dict
+        )
+
+        tmp_dict_2 = dict(**tmp_dict)
+        tmp_dict_2["defined_state"] = False
+
+        self.properties_treated = self.config.property_package.build_state_block(
+            self.flowsheet().time,
+            doc="Material properties of treated water",
+            **tmp_dict_2,
+        )
+        self.properties_byproduct = self.config.property_package.build_state_block(
+            self.flowsheet().time,
+            doc="Material properties of byproduct stream",
+            **tmp_dict_2,
+        )
+
+        # Create Ports
+        self.add_port("inlet", self.properties_in, doc="Inlet port")
+        self.add_port(
+            "treated", self.properties_treated, doc="Treated water outlet port"
+        )
+        self.add_port(
+            "byproduct", self.properties_byproduct, doc="Byproduct outlet port"
+        )
+
+        # Add isothermal constraints
+        @self.Constraint(
+            self.flowsheet().config.time,
+            doc="Isothermal assumption for treated flow",
+        )
+        def eq_treated_isothermal(b, t):
+            return b.properties_in[t].temperature == b.properties_treated[t].temperature
+
+        @self.Constraint(
+            self.flowsheet().config.time,
+            doc="Isothermal assumption for byproduct flow",
+        )
+        def eq_byproduct_isothermal(b, t):
+            return (
+                b.properties_in[t].temperature == b.properties_byproduct[t].temperature
+            )
+
+        # Add performance variables
+        self.recovery_frac_mass_H2O = Var(
+            self.flowsheet().time,
+            initialize=1,
+            domain=NonNegativeReals,
+            units=pyunits.dimensionless,
+            bounds=(0.0, 1.0000001),
+            doc="Mass recovery fraction of water in the treated stream",
+        )
+        self.recovery_frac_mass_H2O.fix()
+
+        self.removal_frac_mass_comp = Var(
+            self.flowsheet().time,
+            self.config.property_package.solute_set,
+            domain=NonNegativeReals,
+            initialize=0.01,
+            units=pyunits.dimensionless,
+            doc="Solute removal fraction on a mass basis",
+        )
+
+        # Add performance constraints
+        # Water recovery
+        @self.Constraint(self.flowsheet().time, doc="Water recovery equation")
+        def water_recovery_equation(b, t):
+            return b.recovery_frac_mass_H2O[t] * b.properties_in[
+                t
+            ].get_material_flow_terms("Liq", "H2O") == b.properties_treated[
+                t
+            ].get_material_flow_terms(
+                "Liq", "H2O"
+            )
+
+        # Flow balance
+        @self.Constraint(self.flowsheet().time, doc="Overall flow balance")
+        def water_balance(b, t):
+            return b.properties_in[t].get_material_flow_terms(
+                "Liq", "H2O"
+            ) == b.properties_treated[t].get_material_flow_terms(
+                "Liq", "H2O"
+            ) + b.properties_byproduct[
+                t
+            ].get_material_flow_terms(
+                "Liq", "H2O"
+            )
+
+        # Solute removal
+        @self.Constraint(
+            self.flowsheet().time,
+            self.config.property_package.phase_list,
+            self.config.property_package.solute_set,
+            doc="Solute removal equations",
+        )
+        def solute_removal_equation(b, t, p, j):
+            return b.removal_frac_mass_comp[t, j] * b.properties_in[
+                t
+            ].get_material_flow_terms(p, j) == b.properties_byproduct[
+                t
+            ].get_material_flow_terms(
+                p, j
+            )
+
+        # Solute concentration of treated stream
+        @self.Constraint(
+            self.flowsheet().time,
+            self.config.property_package.phase_list,
+            self.config.property_package.solute_set,
+            doc="Constraint for solute concentration in treated stream.",
+        )
+        def solute_treated_equation(b, t, p, j):
+            return (1 - b.removal_frac_mass_comp[t, j]) * b.properties_in[
+                t
+            ].get_material_flow_terms(p, j) == b.properties_treated[
+                t
+            ].get_material_flow_terms(
+                p, j
+            )
+
+        # Default solute concentration
+        self.P_removal = Param(
+            within=NonNegativeReals,
+            mutable=True,
+            default=0.98,
+            doc="Reference phosphorus removal fraction on a mass basis",
+            units=pyunits.dimensionless,
+        )
+
+        self.N_removal = Param(
+            within=NonNegativeReals,
+            mutable=True,
+            default=0.3,
+            doc="Reference ammonia removal fraction on a mass basis",
+            units=pyunits.dimensionless,
+        )
+
+        # NOTE: revisit if we should sum soluble nitrogen as total nitrogen
+        @self.Constraint(
+            self.flowsheet().time,
+            self.config.property_package.solute_set,
+            doc="Default solute removal equations",
+        )
+        def default_solute_removal_equation(b, t, j):
+            if j == "S_PO4":
+                return b.removal_frac_mass_comp[t, j] == b.P_removal
+            elif j == "S_NH4":
+                return b.removal_frac_mass_comp[t, j] == b.N_removal
+            else:
+                return b.removal_frac_mass_comp[t, j] == 0
+
+        self._stream_table_dict = {
+            "Inlet": self.inlet,
+            "Treated": self.treated,
+            "Byproduct": self.byproduct,
+        }
+
+        self.electricity = Var(
+            self.flowsheet().time,
+            units=pyunits.kW,
+            bounds=(0, None),
+            doc="Electricity consumption of unit",
+        )
+
+        self.energy_electric_flow_mass = Var(
+            units=pyunits.kWh / pyunits.kg,
+            doc="Electricity intensity with respect to phosphorus removal",
+        )
+
+        @self.Constraint(
+            self.flowsheet().time,
+            doc="Constraint for electricity consumption based on phosphorus removal",
+        )
+        def electricity_consumption(b, t):
+            return b.electricity[t] == (
+                b.energy_electric_flow_mass
+                * pyunits.convert(
+                    b.properties_treated[t].get_material_flow_terms("Liq", "S_PO4"),
+                    to_units=pyunits.kg / pyunits.hour,
+                )
+            )
+
+        self.magnesium_chloride_dosage = Var(
+            units=pyunits.dimensionless,
+            bounds=(0, None),
+            doc="Dosage of magnesium chloride per treated phosphorus",
+        )
+
+        self.MgCl2_flowrate = Var(
+            self.flowsheet().time,
+            units=pyunits.kg / pyunits.hr,
+            bounds=(0, None),
+            doc="Magnesium chloride flowrate",
+        )
+
+        @self.Constraint(
+            self.flowsheet().time,
+            doc="Constraint for magnesium chloride demand based on phosphorus removal.",
+        )
+        def MgCl2_demand(b, t):
+            return b.MgCl2_flowrate[t] == (
+                b.magnesium_chloride_dosage
+                * pyunits.convert(
+                    b.properties_treated[t].get_material_flow_terms("Liq", "S_PO4"),
+                    to_units=pyunits.kg / pyunits.hour,
+                )
+            )
+
+    def initialize_build(
+        self, state_args=None, outlvl=idaeslog.NOTSET, solver=None, optarg=None
+    ):
+        """
+        Initialization routine for single inlet-double outlet unit models.
+
+        Keyword Arguments:
+            state_args : a dict of arguments to be passed to the property
+                           package(s) to provide an initial state for
+                           initialization (see documentation of the specific
+                           property package) (default = {}).
+            outlvl : sets output level of initialization routine
+            optarg : solver options dictionary object (default=None, use
+                     default solver options)
+            solver : str indicating which solver to use during
+                     initialization (default = None, use default IDAES solver)
+
+        Returns:
+            None
+        """
+        if optarg is None:
+            optarg = {}
+
+        # Set solver options
+        init_log = idaeslog.getInitLogger(self.name, outlvl, tag="unit")
+        solve_log = idaeslog.getSolveLogger(self.name, outlvl, tag="unit")
+
+        solver_obj = get_solver(solver, optarg)
+
+        # Get initial guesses for inlet if none provided
+        if state_args is None:
+            state_args = {}
+            state_dict = self.properties_in[
+                self.flowsheet().time.first()
+            ].define_port_members()
+
+            for k in state_dict.keys():
+                if state_dict[k].is_indexed():
+                    state_args[k] = {}
+                    for m in state_dict[k].keys():
+                        state_args[k][m] = state_dict[k][m].value
+                else:
+                    state_args[k] = state_dict[k].value
+
+        # ---------------------------------------------------------------------
+        # Initialize control volume block
+        flags = self.properties_in.initialize(
+            outlvl=outlvl,
+            optarg=optarg,
+            solver=solver,
+            state_args=state_args,
+            hold_state=True,
+        )
+        self.properties_treated.initialize(
+            outlvl=outlvl,
+            optarg=optarg,
+            solver=solver,
+            state_args=state_args,
+            hold_state=False,
+        )
+        self.properties_byproduct.initialize(
+            outlvl=outlvl,
+            optarg=optarg,
+            solver=solver,
+            state_args=state_args,
+            hold_state=False,
+        )
+
+        init_log.info_high("Initialization Step 1 Complete.")
+
+        # ---------------------------------------------------------------------
+        # Solve unit
+        with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc:
+            results = solver_obj.solve(self, tee=slc.tee)
+
+        init_log.info_high(
+            "Initialization Step 2 {}.".format(idaeslog.condition(results))
+        )
+
+        # ---------------------------------------------------------------------
+        # Release Inlet state
+        self.properties_in.release_state(flags, outlvl)
+
+        init_log.info("Initialization Complete: {}".format(idaeslog.condition(results)))
+
+        if not check_optimal_termination(results):
+            raise InitializationError(
+                f"{self.name} failed to initialize successfully. Please check "
+                f"the output logs for more information."
+            )
+
+    def _get_performance_contents(self, time_point=0):
+        var_dict = {}
+        var_dict["Water Recovery"] = self.recovery_frac_mass_H2O[time_point]
+        for j in self.config.property_package.solute_set:
+            var_dict[f"Solute Removal {j}"] = self.removal_frac_mass_comp[time_point, j]
+        var_dict["Electricity Demand"] = self.electricity[time_point]
+        var_dict["Electricity Intensity"] = self.energy_electric_flow_mass
+        var_dict[
+            "Dosage of magnesium chloride per treated phosphorus"
+        ] = self.magnesium_chloride_dosage
+        var_dict["Magnesium Chloride Demand"] = self.MgCl2_flowrate[time_point]
+        return {"vars": var_dict}
+
+    def _get_stream_table_contents(self, time_point=0):
+        return create_stream_table_dataframe(
+            {
+                "Inlet": self.inlet,
+                "Treated": self.treated,
+                "Byproduct": self.byproduct,
+            },
+            time_point=time_point,
+        )
+
+    def calculate_scaling_factors(self):
+        # Get default scale factors and do calculations from base classes
+        for t, v in self.water_recovery_equation.items():
+            iscale.constraint_scaling_transform(
+                v,
+                iscale.get_scaling_factor(
+                    self.properties_in[t].get_material_flow_terms("Liq", "H2O"),
+                    default=1,
+                    warning=True,
+                    hint=" for water recovery",
+                ),
+            )
+
+        for t, v in self.water_balance.items():
+            iscale.constraint_scaling_transform(
+                v,
+                iscale.get_scaling_factor(
+                    self.properties_in[t].get_material_flow_terms("Liq", "H2O"),
+                    default=1,
+                    warning=False,
+                ),
+            )  # would just be a duplicate of above
+
+        for (t, p, j), v in self.solute_removal_equation.items():
+            iscale.constraint_scaling_transform(
+                v,
+                iscale.get_scaling_factor(
+                    self.properties_in[t].get_material_flow_terms(p, j),
+                    default=1,
+                    warning=True,
+                    hint=" for solute removal",
+                ),
+            )
+
+        for (t, p, j), v in self.solute_treated_equation.items():
+            iscale.constraint_scaling_transform(
+                v,
+                iscale.get_scaling_factor(
+                    self.properties_in[t].get_material_flow_terms(p, j),
+                    default=1,
+                    warning=False,
+                ),
+            )  # would just be a duplicate of above
diff --git a/watertap/unit_models/tests/test_electroNP_ZO.py b/watertap/unit_models/tests/test_electroNP_ZO.py
new file mode 100644
index 0000000000..fc9d39eac1
--- /dev/null
+++ b/watertap/unit_models/tests/test_electroNP_ZO.py
@@ -0,0 +1,269 @@
+#################################################################################
+# 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/"
+#################################################################################
+
+import pytest
+from pyomo.environ import (
+    ConcreteModel,
+    assert_optimal_termination,
+    value,
+    Var,
+    log10,
+    units,
+)
+from pyomo.network import Port
+from idaes.core import (
+    FlowsheetBlock,
+    MaterialBalanceType,
+    EnergyBalanceType,
+    MomentumBalanceType,
+)
+from watertap.unit_models.electroNP_ZO import ElectroNPZO
+from watertap.property_models.activated_sludge.modified_asm2d_properties import (
+    ModifiedASM2dParameterBlock,
+)
+from watertap.property_models.activated_sludge.modified_asm2d_reactions import (
+    ModifiedASM2dReactionParameterBlock,
+    ModifiedASM2dReactionBlock,
+)
+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,
+)
+from idaes.core.util.testing import initialization_tester
+from idaes.core.util.scaling import (
+    calculate_scaling_factors,
+    unscaled_variables_generator,
+    unscaled_constraints_generator,
+    badly_scaled_var_generator,
+)
+from pyomo.util.check_units import assert_units_consistent
+from idaes.core import UnitModelCostingBlock
+from watertap.costing import WaterTAPCosting
+
+# -----------------------------------------------------------------------------
+# Get default solver for testing
+solver = get_solver()
+
+
+class TestElectroNP:
+    @pytest.fixture(scope="class")
+    def ElectroNP_frame(self):
+        m = ConcreteModel()
+        m.fs = FlowsheetBlock(dynamic=False)
+
+        m.fs.properties = ModifiedASM2dParameterBlock(
+            additional_solute_list=["S_K", "S_Mg"]
+        )
+
+        m.fs.unit = ElectroNPZO(property_package=m.fs.properties)
+
+        EPS = 1e-8
+
+        m.fs.unit.inlet.temperature.fix(298.15 * units.K)
+        m.fs.unit.inlet.pressure.fix(1 * units.atm)
+
+        m.fs.unit.inlet.flow_vol.fix(18446 * units.m**3 / units.day)
+        m.fs.unit.inlet.conc_mass_comp[0, "S_O2"].fix(10 * units.mg / units.liter)
+        m.fs.unit.inlet.conc_mass_comp[0, "S_N2"].fix(EPS * units.mg / units.liter)
+        m.fs.unit.inlet.conc_mass_comp[0, "S_NH4"].fix(16 * units.mg / units.liter)
+        m.fs.unit.inlet.conc_mass_comp[0, "S_NO3"].fix(EPS * units.mg / units.liter)
+        m.fs.unit.inlet.conc_mass_comp[0, "S_PO4"].fix(3.6 * units.mg / units.liter)
+        m.fs.unit.inlet.conc_mass_comp[0, "S_F"].fix(30 * units.mg / units.liter)
+        m.fs.unit.inlet.conc_mass_comp[0, "S_A"].fix(20 * units.mg / units.liter)
+        m.fs.unit.inlet.conc_mass_comp[0, "S_I"].fix(30 * units.mg / units.liter)
+        m.fs.unit.inlet.conc_mass_comp[0, "X_I"].fix(25 * units.mg / units.liter)
+        m.fs.unit.inlet.conc_mass_comp[0, "X_S"].fix(125 * units.mg / units.liter)
+        m.fs.unit.inlet.conc_mass_comp[0, "X_H"].fix(30 * units.mg / units.liter)
+        m.fs.unit.inlet.conc_mass_comp[0, "X_PAO"].fix(EPS * units.mg / units.liter)
+        m.fs.unit.inlet.conc_mass_comp[0, "X_PP"].fix(EPS * units.mg / units.liter)
+        m.fs.unit.inlet.conc_mass_comp[0, "X_PHA"].fix(EPS * units.mg / units.liter)
+        m.fs.unit.inlet.conc_mass_comp[0, "X_AUT"].fix(EPS * units.mg / units.liter)
+        m.fs.unit.inlet.conc_mass_comp[0, "X_MeOH"].fix(EPS * units.mg / units.liter)
+        m.fs.unit.inlet.conc_mass_comp[0, "X_MeP"].fix(EPS * units.mg / units.liter)
+        m.fs.unit.inlet.conc_mass_comp[0, "X_TSS"].fix(EPS * units.mg / units.liter)
+        m.fs.unit.inlet.conc_mass_comp[0, "S_K"].fix(EPS * units.mg / units.liter)
+        m.fs.unit.inlet.conc_mass_comp[0, "S_Mg"].fix(EPS * units.mg / units.liter)
+
+        # Alkalinity was givien in mg/L based on C
+        m.fs.unit.inlet.alkalinity[0].fix(61 / 12 * units.mmol / units.liter)
+
+        # Unit option
+        m.fs.unit.energy_electric_flow_mass.fix(0.044 * units.kWh / units.kg)
+        m.fs.unit.magnesium_chloride_dosage.fix(0.388)
+
+        return m
+
+    @pytest.mark.unit
+    def test_dof(self, ElectroNP_frame):
+        m = ElectroNP_frame
+        assert degrees_of_freedom(m) == 0
+
+    @pytest.mark.unit
+    def test_units(self, ElectroNP_frame):
+        assert_units_consistent(ElectroNP_frame)
+
+    @pytest.mark.unit
+    def test_calculate_scaling(self, ElectroNP_frame):
+        m = ElectroNP_frame
+        calculate_scaling_factors(m)
+
+    @pytest.mark.solver
+    @pytest.mark.skipif(solver is None, reason="Solver not available")
+    @pytest.mark.component
+    def test_initialize(self, ElectroNP_frame):
+        initialization_tester(ElectroNP_frame)
+
+    @pytest.mark.solver
+    @pytest.mark.skipif(solver is None, reason="Solver not available")
+    @pytest.mark.component
+    def test_solve(self, ElectroNP_frame):
+        m = ElectroNP_frame
+        results = solver.solve(m)
+
+        # Check for optimal solution
+        assert_optimal_termination(results)
+
+    @pytest.mark.solver
+    @pytest.mark.skipif(solver is None, reason="Solver not available")
+    @pytest.mark.component
+    def test_conservation(self, ElectroNP_frame):
+        m = ElectroNP_frame
+        assert (
+            abs(
+                value(
+                    m.fs.unit.inlet.flow_vol[0] * m.fs.properties.dens_mass
+                    - m.fs.unit.treated.flow_vol[0] * m.fs.properties.dens_mass
+                    - m.fs.unit.byproduct.flow_vol[0] * m.fs.properties.dens_mass
+                )
+            )
+            <= 1e-6
+        )
+        for j in m.fs.properties.solute_set:
+            assert 1e-6 >= abs(
+                value(
+                    m.fs.unit.inlet.flow_vol[0] * m.fs.unit.inlet.conc_mass_comp[0, j]
+                    - m.fs.unit.treated.flow_vol[0]
+                    * m.fs.unit.treated.conc_mass_comp[0, j]
+                    - m.fs.unit.byproduct.flow_vol[0]
+                    * m.fs.unit.byproduct.conc_mass_comp[0, j]
+                )
+            )
+
+    @pytest.mark.solver
+    @pytest.mark.skipif(solver is None, reason="Solver not available")
+    @pytest.mark.component
+    def test_solution(self, ElectroNP_frame):
+        m = ElectroNP_frame
+        assert value(m.fs.unit.treated.flow_vol[0]) == pytest.approx(0.213495, rel=1e-3)
+
+        assert value(m.fs.unit.treated.temperature[0]) == pytest.approx(
+            298.15, rel=1e-4
+        )
+        assert value(m.fs.unit.treated.pressure[0]) == pytest.approx(101325, rel=1e-4)
+        assert value(m.fs.unit.treated.conc_mass_comp[0, "S_A"]) == pytest.approx(
+            0.02, rel=1e-4
+        )
+        assert value(m.fs.unit.treated.conc_mass_comp[0, "S_F"]) == pytest.approx(
+            0.03, rel=1e-2
+        )
+        assert value(m.fs.unit.treated.conc_mass_comp[0, "S_I"]) == pytest.approx(
+            0.03, rel=1e-4
+        )
+        assert value(m.fs.unit.treated.conc_mass_comp[0, "S_N2"]) == pytest.approx(
+            0, abs=1e-4
+        )
+        assert value(m.fs.unit.treated.conc_mass_comp[0, "S_NH4"]) == pytest.approx(
+            0.0112, rel=1e-4
+        )
+        assert value(m.fs.unit.treated.conc_mass_comp[0, "S_NO3"]) == pytest.approx(
+            0, abs=1e-4
+        )
+        assert value(m.fs.unit.treated.conc_mass_comp[0, "S_O2"]) == pytest.approx(
+            0.01, rel=1e-4
+        )
+        assert value(m.fs.unit.treated.conc_mass_comp[0, "S_PO4"]) == pytest.approx(
+            7.2e-5, rel=1e-4
+        )
+        assert value(m.fs.unit.treated.conc_mass_comp[0, "X_AUT"]) == pytest.approx(
+            0, abs=1e-4
+        )
+        assert value(m.fs.unit.treated.conc_mass_comp[0, "X_H"]) == pytest.approx(
+            0.03, rel=1e-4
+        )
+        assert value(m.fs.unit.treated.conc_mass_comp[0, "X_I"]) == pytest.approx(
+            0.025, rel=1e-4
+        )
+        assert value(m.fs.unit.treated.conc_mass_comp[0, "X_MeOH"]) == pytest.approx(
+            0, abs=1e-4
+        )
+        assert value(m.fs.unit.treated.conc_mass_comp[0, "X_MeP"]) == pytest.approx(
+            0, abs=1e-4
+        )
+        assert value(m.fs.unit.treated.conc_mass_comp[0, "X_PAO"]) == pytest.approx(
+            0, abs=1e-4
+        )
+        assert value(m.fs.unit.treated.conc_mass_comp[0, "X_PHA"]) == pytest.approx(
+            0, abs=1e-4
+        )
+        assert value(m.fs.unit.treated.conc_mass_comp[0, "X_PP"]) == pytest.approx(
+            0, abs=1e-4
+        )
+        assert value(m.fs.unit.treated.conc_mass_comp[0, "X_S"]) == pytest.approx(
+            0.125, rel=1e-4
+        )
+        assert value(m.fs.unit.treated.conc_mass_comp[0, "X_TSS"]) == pytest.approx(
+            0, abs=1e-4
+        )
+        assert value(m.fs.unit.treated.conc_mass_comp[0, "S_Mg"]) == pytest.approx(
+            0, abs=1e-4
+        )
+        assert value(m.fs.unit.treated.conc_mass_comp[0, "S_K"]) == pytest.approx(
+            0, abs=1e-4
+        )
+        assert value(m.fs.unit.byproduct.conc_mass_comp[0, "S_NH4"]) == pytest.approx(
+            4.5289e14, rel=1e-4
+        )
+        assert value(m.fs.unit.byproduct.conc_mass_comp[0, "S_PO4"]) == pytest.approx(
+            3.328756e14, rel=1e-4
+        )
+        assert value(m.fs.unit.treated.alkalinity[0]) == pytest.approx(
+            0.005083, rel=1e-4
+        )
+
+    @pytest.mark.solver
+    @pytest.mark.skipif(solver is None, reason="Solver not available")
+    @pytest.mark.component
+    def test_costing(self, ElectroNP_frame):
+        m = ElectroNP_frame
+
+        m.fs.costing = WaterTAPCosting()
+
+        m.fs.unit.costing = UnitModelCostingBlock(flowsheet_costing_block=m.fs.costing)
+        m.fs.costing.cost_process()
+        m.fs.costing.add_LCOW(m.fs.unit.properties_treated[0].flow_vol)
+        results = solver.solve(m)
+
+        assert_optimal_termination(results)
+
+        # Check solutions
+        assert pytest.approx(1295.765, rel=1e-5) == value(
+            m.fs.unit.costing.capital_cost
+        )
+        assert pytest.approx(5.800325e-5, rel=1e-5) == value(m.fs.costing.LCOW)
+
+    @pytest.mark.unit
+    def test_report(self, ElectroNP_frame):
+        m = ElectroNP_frame
+        m.fs.unit.report()