diff --git a/tutorials/RO_flowsheet.png b/tutorials/RO_flowsheet.png
new file mode 100644
index 0000000000..1ce4b3a75a
Binary files /dev/null and b/tutorials/RO_flowsheet.png differ
diff --git a/tutorials/unit_model_customization_example.ipynb b/tutorials/unit_model_customization_example.ipynb
new file mode 100644
index 0000000000..3ff9564a0e
--- /dev/null
+++ b/tutorials/unit_model_customization_example.ipynb
@@ -0,0 +1,886 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Tutorial on customizing unit models in WaterTAP\n",
+    "Demonstration of how to modify existing unit models at flowsheet level\n",
+    "\n",
+    "## Dependencies\n",
+    "* Python - Programming language\n",
+    "* Pyomo - Python package for equation-oriented modeling\n",
+    "* IDAES - Python package extending Pyomo for flowsheet modeling\n",
+    "* WaterTAP - Unit models"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Demonstration structure \n",
+    "* Setting up basic RO flowsheet\n",
+    "* Replace a fixed variable with an equation to study impact of variable performance metrics\n",
+    "* Replace existing constraint with a new one for comparative analysis  \n",
+    "* Replace a fixed cost with variable cost that is a function of operating condition "
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Flowsheet considered in the example\n",
+    "<img src=\"RO_flowsheet.png\" width=\"500\" height=\"200\">"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Import key modules"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "## Import core components \n",
+    "# Pyomo core components\n",
+    "from pyomo.environ import (Param,Var, Constraint, TransformationFactory, Reals,    ConcreteModel,\n",
+    "    value,assert_optimal_termination,\n",
+    "    units as pyunits)\n",
+    "from pyomo.network import Arc\n",
+    "# Ideas core components\n",
+    "from idaes.core import FlowsheetBlock\n",
+    "from idaes.core.util.scaling import calculate_scaling_factors, set_scaling_factor\n",
+    "from idaes.core.util.model_statistics import degrees_of_freedom\n",
+    "from idaes.core.solvers import get_solver\n",
+    "from idaes.core.util.scaling import constraint_scaling_transform\n",
+    "from idaes.core.util.initialization import propagate_state\n",
+    "from idaes.models.unit_models import Feed, Product\n",
+    "# WaterTAP core components \n",
+    "import watertap.property_models.seawater_prop_pack as properties\n",
+    "from watertap.unit_models.reverse_osmosis_1D import (\n",
+    "    ReverseOsmosis1D,\n",
+    "    ConcentrationPolarizationType,\n",
+    "    MassTransferCoefficient,\n",
+    "    PressureChangeType,\n",
+    ")\n",
+    "from watertap.unit_models.pressure_changer import Pump\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Build the flowsheet"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "\n",
+    "m = ConcreteModel()\n",
+    "# create IDAES flowsheet\n",
+    "m.fs = FlowsheetBlock(dynamic=False)\n",
+    "# create seawater property model\n",
+    "m.fs.properties = properties.SeawaterParameterBlock()\n",
+    "\n",
+    "# build feed\n",
+    "m.fs.feed = Feed(property_package=m.fs.properties)\n",
+    "m.fs.product = Product(property_package=m.fs.properties)\n",
+    "#build pump\n",
+    "m.fs.pump = Pump(property_package=m.fs.properties)\n",
+    "m.fs.RO = ReverseOsmosis1D(\n",
+    "    property_package=m.fs.properties,\n",
+    "    has_pressure_change=True,\n",
+    "    pressure_change_type=PressureChangeType.calculated,\n",
+    "    mass_transfer_coefficient=MassTransferCoefficient.calculated,\n",
+    "    concentration_polarization_type=ConcentrationPolarizationType.calculated,\n",
+    "    transformation_scheme=\"BACKWARD\",\n",
+    "    transformation_method=\"dae.finite_difference\",\n",
+    "    finite_elements=10,\n",
+    ")\n",
+    "# connect feed to pump\n",
+    "m.fs.feed_to_pump = Arc(source=m.fs.feed.outlet, destination = m.fs.pump.inlet)\n",
+    "#connect pump to RO unit\n",
+    "m.fs.pump_to_ro = Arc(source=m.fs.pump.outlet, destination = m.fs.RO.inlet)\n",
+    "# connect RO permeate to product\n",
+    "m.fs.ro_to_product = Arc(source=m.fs.RO.permeate, destination = m.fs.product.inlet)\n",
+    "TransformationFactory(\"network.expand_arcs\").apply_to(m) \n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Set default values and calculate scaling factors for flowsheet   "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "\n",
+    "m.fs.feed.properties[0].temperature.fix(273 + 25)                      # temperature (K)\n",
+    "m.fs.feed.properties[0].pressure.fix(101325)                           # pressure (Pa)\n",
+    "m.fs.feed.properties[0].flow_mass_phase_comp['Liq', 'H2O'].fix(0.965)  # mass flowrate of H2O (kg/s)\n",
+    "m.fs.feed.properties[0].flow_mass_phase_comp['Liq', 'TDS'].fix(0.035)  # mass flowrate of TDS (kg/s)\n",
+    "m.fs.feed.properties[0].conc_mass_phase_comp[...] # construct concentration props\n",
+    "m.fs.properties.set_default_scaling(\n",
+    "    \"flow_mass_phase_comp\",\n",
+    "    1/0.965,\n",
+    "    index=(\"Liq\", \"H2O\"),\n",
+    ")\n",
+    "m.fs.properties.set_default_scaling(\n",
+    "    \"flow_mass_phase_comp\",\n",
+    "    1/0.035,\n",
+    "    index=(\"Liq\", \"TDS\"),\n",
+    ")\n",
+    "# to help with initialization, let's build the osmotic pressure variable on the feed block\n",
+    "# which we can use to guess operating pressure for RO unit and set pump pressure during initialization\n",
+    "m.fs.feed.properties[0].pressure_osm_phase[...]\n",
+    "# define pump defaults\n",
+    "m.fs.pump.efficiency_pump[0].fix(0.75)\n",
+    "# scale work and pressures for the pump\n",
+    "set_scaling_factor(m.fs.pump.control_volume.work, 1e-4)\n",
+    "set_scaling_factor(m.fs.pump.control_volume.properties_out[0].pressure, 1e-5)\n",
+    "set_scaling_factor(m.fs.pump.control_volume.properties_in[0].pressure, 1e-5)\n",
+    "\n",
+    "# define RO default values for initialization \n",
+    "# we opt to specify stage area, and inlet velocity\n",
+    "# unfixing width and area\n",
+    "# We also apply variable scaling as we set up each default parameter \n",
+    "\n",
+    "m.fs.RO.feed_side.velocity[0, 0].fix(0.1)\n",
+    "\n",
+    "m.fs.RO.area.fix(100)\n",
+    "set_scaling_factor(m.fs.RO.area,1/50)\n",
+    "m.fs.RO.length.unfix()\n",
+    "set_scaling_factor(m.fs.RO.length, 0.1)\n",
+    "m.fs.RO.width.unfix()\n",
+    "set_scaling_factor(m.fs.RO.width, 0.1)\n",
+    "\n",
+    "# we need to specify RO permeate pressure\n",
+    "m.fs.RO.permeate.pressure[0].fix(101325)\n",
+    "# we need to specify default values for default mass transport\n",
+    "# and friction factor correlations \n",
+    "m.fs.RO.feed_side.channel_height.fix(1e-3)\n",
+    "m.fs.RO.feed_side.spacer_porosity.fix(0.9)\n",
+    "\n",
+    "# Specify default A and B values, these are defined as m/s/Pa and m/s, respectively. \n",
+    "m.fs.RO.A_comp[0, \"H2O\"].fix(3 / (3600 * 1000 * 1e5))\n",
+    "m.fs.RO.B_comp[0, \"TDS\"].fix(0.15 / (3600 * 1000))\n",
+    "\n",
+    "# calculate all the scaling factors \n",
+    "calculate_scaling_factors(m)\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Initialize feed, pump, and RO unit"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "solver = get_solver() # get solver\n",
+    "m.fs.feed.initialize(optarg=solver.options)\n",
+    "propagate_state(m.fs.feed_to_pump)\n",
+    "# get osmotic pressure\n",
+    "osmotic_feed_pressure=value(m.fs.feed.properties[0].pressure_osm_phase['Liq'])\n",
+    "print(\"Osmotic pressure is {} bar\".format(osmotic_feed_pressure/1e5))\n",
+    "m.fs.pump.outlet.pressure[0].fix(osmotic_feed_pressure*1.5) \n",
+    "m.fs.pump.initialize(optarg=solver.options)\n",
+    "propagate_state(m.fs.pump_to_ro)\n",
+    "m.fs.RO.initialize(optarg=solver.options)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Solve box problem at initialized state and solve to operating state of 50% water recovery"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "\n",
+    "# Check degrees of freedom \n",
+    "print('We have {} degrees of freedom and expect 0'.format(degrees_of_freedom(m)))\n",
+    "assert degrees_of_freedom(m) == 0\n",
+    "# Solve box problem at initialized state \n",
+    "result =solver.solve(m, tee=False)\n",
+    "assert_optimal_termination(result)\n",
+    "\n",
+    "# apply operating conditions and solve again\n",
+    "m.fs.pump.outlet.pressure[0].unfix()\n",
+    "m.fs.RO.recovery_vol_phase[0.0, \"Liq\"].fix(0.5)\n",
+    "print('We have {} degrees of freedom and expect 0'.format(degrees_of_freedom(m)))\n",
+    "assert degrees_of_freedom(m) == 0\n",
+    "\n",
+    "result =solver.solve(m, tee=True)\n",
+    "assert_optimal_termination(result)\n",
+    "# lets check current solution \n",
+    "m.fs.RO.report()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Explore relationship between salinity and opertating pressure"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# unfix current concentration \n",
+    "m.fs.feed.properties[0].flow_mass_phase_comp['Liq','TDS'].unfix()\n",
+    "print('Current concentration {}'.format(m.fs.feed.properties[0].conc_mass_phase_comp['Liq','TDS'].value))\n",
+    "\n",
+    "# import numpy \n",
+    "import numpy as np\n",
+    "\n",
+    "# define concentrations we are going to sweep over \n",
+    "concentrations = np.linspace(15,100,10)\n",
+    "pressures = []\n",
+    "for con in concentrations:\n",
+    "    m.fs.feed.properties[0].conc_mass_phase_comp['Liq','TDS'].fix(con)\n",
+    "    assert degrees_of_freedom(m) == 0\n",
+    "    result = solver.solve(m, tee=False)\n",
+    "    assert_optimal_termination(result)\n",
+    "    pressures.append(m.fs.RO.inlet.pressure[0].value/1e5)\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Plot feed salinity vs. operating pressure required to achieve 50% water recovery with 50 $m^2$ of membrane area"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import matplotlib.pyplot as plt\n",
+    "fig, ax = plt.subplots()\n",
+    "ax.plot(concentrations,pressures)\n",
+    "ax.set(xlabel='Concentration (g/L)',ylabel='Pressure (bar)')\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Add relationship between A parameter and inlet pressure to account for compaction effects    \n",
+    "\n",
+    "Here we will relate the A value of our RO unit to RO inlet pressure ($P_{inlet}$) as follows:\n",
+    "$$\n",
+    " A=3 \\frac{LMH}{bar}\\, for \\, P_{inlet}<65\\,bar\n",
+    "$$\n",
+    "$$\n",
+    " A = 3\\frac{LMH}{bar}*65/P_{inlet} \\, for \\, P_{inlet}>=65\\,bar\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Imports smooth min and smooth max functions\n",
+    "import idaes.core.util.math as idaesMath\n",
+    "\n",
+    "# define our default A value \n",
+    "m.fs.A_var_initial=Var(initialize=3.0)\n",
+    "m.fs.A_var_initial.fix()\n",
+    "set_scaling_factor(m.fs.A_var_initial, 1/m.fs.A_var_initial.value)\n",
+    "\n",
+    "# define constraint that relates A value to P\n",
+    "m.fs.RO.A_pressure_constraint=(\n",
+    "    Constraint(expr=m.fs.RO.A_comp[0, \"H2O\"]*(3600 * 1000 * 1e5)==\n",
+    "    idaesMath.smooth_min(m.fs.A_var_initial,(m.fs.A_var_initial*(65*1e5/m.fs.RO.inlet.pressure[0])))))\n",
+    "    \n",
+    "m.fs.RO.A_comp[0,'H2O'].unfix()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Verify that constraint produces expected outcomes "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# import method to calculate variable from constraint, allows evaluation of our constraint \n",
+    "from pyomo.util.calc_var_value import calculate_variable_from_constraint\n",
+    "\n",
+    "A_vals=[]\n",
+    "# define range over which to test the constraint\n",
+    "pressures_for_testing = np.linspace(10,200,20)\n",
+    "for pressure in pressures_for_testing:\n",
+    "    m.fs.RO.inlet.pressure[0]=pressure*1e5 # needs to be in Pa\n",
+    "    # this will estimate the RO.A_comp from our constraint.\n",
+    "    calculate_variable_from_constraint(m.fs.RO.A_comp[0, \"H2O\"],m.fs.RO.A_pressure_constraint)\n",
+    "    A_vals.append(m.fs.RO.A_comp[0, \"H2O\"].value*(3600 * 1000 * 1e5))\n",
+    "\n",
+    "# plot the result \n",
+    "fig, ax = plt.subplots()\n",
+    "ax.plot(pressures_for_testing,A_vals)\n",
+    "ax.set(xlabel='Pressure (bar)',ylabel='A value (LMH/bar)')\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Initalize A value at operating pressure and solve model with new constraint"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Lets initialize the A value to our actual operating pressure \n",
+    "m.fs.RO.inlet.pressure[0]=m.fs.pump.outlet.pressure[0].value \n",
+    "# Calculate the unknown value of A at specified pressure form our constraint \n",
+    "calculate_variable_from_constraint(m.fs.RO.A_comp[0, \"H2O\"],m.fs.RO.A_pressure_constraint)\n",
+    "\n",
+    "\n",
+    "print('We have {} degrees of freedom and expect 0'.format(degrees_of_freedom(m)))\n",
+    "assert degrees_of_freedom(m) == 0\n",
+    "\n",
+    "# Solve model with new constraint\n",
+    "result =solver.solve(m, tee=False)\n",
+    "assert_optimal_termination(result)\n",
+    "\n",
+    "# Explore how performance changes across feed concentrations \n",
+    "concentrations = np.linspace(15,100,10)\n",
+    "pressures_with_a_loss = []\n",
+    "actual_a=[]\n",
+    "for con in concentrations:\n",
+    "    m.fs.feed.properties[0].conc_mass_phase_comp['Liq','TDS'].fix(con)\n",
+    "    result =solver.solve(m, tee=False)\n",
+    "    assert_optimal_termination(result)\n",
+    "    pressures_with_a_loss.append(m.fs.RO.inlet.pressure[0].value/1e5)\n",
+    "    actual_a.append(m.fs.RO.A_comp[0, \"H2O\"].value*(3600 * 1000 * 1e5))\n",
+    "    print(\"Solved con {}\".format(con))\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Plot how operating pressure and A value changed during RO operation"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "\n",
+    "fig, ax = plt.subplots()\n",
+    "ax.plot(concentrations,pressures, color='black', label='Fixed A')\n",
+    "ax.plot(concentrations,pressures_with_a_loss, color='red', label='A as function of pressure')\n",
+    "ax.set(xlabel='Concentration (g/L)',ylabel='Pressure (bar)')\n",
+    "ax.legend()\n",
+    "fig, ax = plt.subplots()\n",
+    "ax.plot(pressures,actual_a,  color='red',label='A as function of pressure')\n",
+    "ax.set(xlabel='Pressure (bar)',ylabel='A value  (LMH/bar)')\n",
+    "ax.legend()\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Replacing an existing constraint with a new one\n",
+    "\n",
+    "A common uncertainty in membrane process is prediction of mass transport rates.\n",
+    "\n",
+    "Current WaterTap model uses a correlation derived by Guillen et al from CFD simulation of 2D flow in RO channel with a circular spacer filaments. The used correlation is shown below, as well as its implementation in RO model. Current correlation is: \n",
+    "$$ Sh = 0.45*(Re*Sc)^{0.36}  $$\n",
+    "\n",
+    "    @self.Constraint(\n",
+    "        self.flowsheet().config.time,\n",
+    "        self.length_domain,\n",
+    "        self.config.property_package.solute_set,\n",
+    "        doc=\"Sherwood number\",\n",
+    "    )\n",
+    "    def eq_N_Sh_comp(b, t, x, j):\n",
+    "        return (\n",
+    "            b.N_Sh_comp[t, x, j]\n",
+    "            == 0.46 * (b.N_Re[t, x] * b.N_Sc_comp[t, x, j]) ** 0.36\n",
+    "        )\n",
+    "\n",
+    "An alternative correlation that could be used has been derived by Schock & Miquel [Desalination, 1987, 64, 339-352] from experiments. This correlation has been shown to be potentially more accurate than the Guillen correlation [Dudchenko et al. ACS ES&T Engineering, https://doi.org/10.1021/acsestengg.1c00496]\n",
+    "\n",
+    "$$ Sh = 0.065 * Re^{0.875} * Sc^{0.33} $$\n",
+    "\n",
+    "Let's implement the Schock & Miquel correlation instead and observe its overall impact on process operation.\n",
+    "\n",
+    "We will further include an Sherwood adjustment factor to enable sensitivity analysis such that our final correlation is: \n",
+    "$$ Sh = 0.065 * Re^{0.875} * Sc^{0.33} * Sh_{multiplier} $$"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Adding a variable so we can adjust magnitude of Sherwood correlation for sensitivity analysis \n",
+    "m.fs.Sh_multiplier=Var(initialize=1)\n",
+    "m.fs.Sh_multiplier.fix()\n",
+    "\n",
+    "# defining new sherwood correlation \n",
+    "@m.fs.RO.feed_side.Constraint(\n",
+    "    [0],\n",
+    "    m.fs.RO.length_domain,\n",
+    "    m.fs.properties.solute_set,\n",
+    "    doc=\"Sherwood number Schock & Miquel\",\n",
+    ")\n",
+    "def eq_N_Sh_comp_S_and_M(b, t, x, j):\n",
+    "    return (\n",
+    "        b.N_Sh_comp[t, x, j]\n",
+    "        == (0.065 * b.N_Re[t, x]**0.875 * b.N_Sc_comp[t, x, j] ** 0.33) *m.fs.Sh_multiplier\n",
+    "    )\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Check output of our new correlation against old one"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "Guillen_sh=[]\n",
+    "Schock_sh=[]\n",
+    "module_position=[]\n",
+    "\n",
+    "# grab Sh values estimated using old correlation \n",
+    "for (t,x,j) in m.fs.RO.feed_side.N_Sh_comp:\n",
+    "    module_position.append(x)    \n",
+    "    Guillen_sh.append(m.fs.RO.feed_side.N_Sh_comp[t,x,j].value)\n",
+    "# grab Sj values using new correlation \n",
+    "for (t,x,j) in m.fs.RO.feed_side.N_Sh_comp:\n",
+    "    calculate_variable_from_constraint(m.fs.RO.feed_side.N_Sh_comp[t,x,j],m.fs.RO.feed_side.eq_N_Sh_comp_S_and_M[t,x,j])\n",
+    "    Schock_sh.append(m.fs.RO.feed_side.N_Sh_comp[t,x,j].value)\n",
+    "\n",
+    "fig, ax = plt.subplots()\n",
+    "ax.plot(module_position,Guillen_sh, label='Guillen correlation')\n",
+    "ax.plot(module_position,Schock_sh, label='Schock & Miquel correlation')\n",
+    "ax.set(xlabel='Module position (-)',ylabel='Sherwood number (-)')\n",
+    "ax.legend()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Solve the model with Sherwood correlation"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Because the old and new Sherwood numbers have similar values (Within 10x of each other)\n",
+    "# it is not necessary to re-initialize them\n",
+    "\n",
+    "# First lets solve with old constraint to get Water Flux and Salt flux, as we did not grab these values earlier on \n",
+    "'''solve with old constraint'''\n",
+    "# make sure old constraint is active\n",
+    "m.fs.RO.feed_side.eq_N_Sh_comp.activate()\n",
+    "# make sure new constraint is NOT active\n",
+    "m.fs.RO.feed_side.eq_N_Sh_comp_S_and_M.deactivate()\n",
+    "print('We have {} degrees of freedom and expect 0'.format(degrees_of_freedom(m)))\n",
+    "assert degrees_of_freedom(m) == 0\n",
+    "result = solver.solve(m, tee=False)\n",
+    "water_flux_Guillen=[]\n",
+    "salt_flux_Guillen=[]\n",
+    "module_position=[]\n",
+    "# get water and salt flux estimated using Guillen correlation \n",
+    "for (t,x,p,j) in m.fs.RO.flux_mass_phase_comp:\n",
+    "    if j=='H2O':\n",
+    "        water_flux_Guillen.append(m.fs.RO.flux_mass_phase_comp[t,x,p,j].value*3600)\n",
+    "    if j=='TDS':\n",
+    "        salt_flux_Guillen.append(m.fs.RO.flux_mass_phase_comp[t,x,p,j].value*3600)\n",
+    "        module_position.append(x)\n",
+    "\n",
+    "# Now let's activate new Sherwood correlation and solve for new water and salt flux \n",
+    "'''solve with new constraint'''\n",
+    "# make sure old constraint is NOT active\n",
+    "m.fs.RO.feed_side.eq_N_Sh_comp.deactivate()\n",
+    "# make sure new constraint is active\n",
+    "m.fs.RO.feed_side.eq_N_Sh_comp_S_and_M.activate()\n",
+    "print('We have {} degrees of freedom and expect 0'.format(degrees_of_freedom(m)))\n",
+    "assert degrees_of_freedom(m) == 0\n",
+    "result = solver.solve(m, tee=False)\n",
+    "assert_optimal_termination(result)\n",
+    "water_flux_Schock=[]\n",
+    "salt_flux_Schock=[]\n",
+    "for (t,x,p,j) in m.fs.RO.flux_mass_phase_comp:\n",
+    "    if j=='H2O':\n",
+    "        water_flux_Schock.append(m.fs.RO.flux_mass_phase_comp[t,x,p,j].value*3600)\n",
+    "    if j=='TDS':\n",
+    "        salt_flux_Schock.append(m.fs.RO.flux_mass_phase_comp[t,x,p,j].value*3600)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Plot comparison of water and salt flux estimated using Guillen and Schock & Miquel correlation"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "\n",
+    "\n",
+    "fig, ax = plt.subplots()\n",
+    "ax.plot(module_position,water_flux_Guillen, label='Guillen correlation')\n",
+    "ax.plot(module_position,water_flux_Schock, label='Schock & Miquel correlation')\n",
+    "ax.set(xlabel='Module position (-)',ylabel='Water flux (LMH)')\n",
+    "ax.legend()\n",
+    "fig, ax = plt.subplots()\n",
+    "ax.plot(module_position,salt_flux_Guillen, label='Guillen correlation')\n",
+    "ax.plot(module_position,salt_flux_Schock, label='Schock & Miquel correlation')\n",
+    "ax.set(xlabel='Module position (-)',ylabel='Salt flux (LMH)')\n",
+    "ax.legend()\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Let's explore how changing Sherwood number impact operating pressure"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# define adjustments for Sh to sweep over\n",
+    "multipliers = np.linspace(1,2,5)\n",
+    "pressure = []\n",
+    "avg_sh =[] \n",
+    "for adj in multipliers:\n",
+    "    m.fs.Sh_multiplier.fix(adj)\n",
+    "    assert degrees_of_freedom(m) == 0\n",
+    "    result =solver.solve(m, tee=False)\n",
+    "    pressure.append(m.fs.RO.inlet.pressure[0].value/1e5)\n",
+    "    sh_avg=np.average([m.fs.RO.feed_side.N_Sh_comp[t,x,j].value for (t,x,j) in m.fs.RO.feed_side.N_Sh_comp])\n",
+    "    avg_sh.append(sh_avg)\n",
+    "    print(\"Solved multiplier {}\".format(adj))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Plot change in operating pressure vs. increase in average Sherwood number "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig, ax = plt.subplots()\n",
+    "ax.plot(multipliers,pressure)\n",
+    "ax.set(xlabel='Sh multiplier (-)',ylabel='Pressure (bar)')\n",
+    "fig, ax = plt.subplots()\n",
+    "ax.plot(multipliers,avg_sh)\n",
+    "ax.set(xlabel='Sh multiplier (-)',ylabel='Average Sh (-)')\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Let's add costing to explore impact of Sh on process cost "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from watertap.costing import WaterTAPCosting\n",
+    "from idaes.core import UnitModelCostingBlock\n",
+    "\n",
+    "# add costing block\n",
+    "m.fs.costing = WaterTAPCosting()\n",
+    "\n",
+    "# add Pump costing\n",
+    "m.fs.pump.costing = UnitModelCostingBlock(flowsheet_costing_block=m.fs.costing)\n",
+    "# add RO costing\n",
+    "m.fs.RO.costing = UnitModelCostingBlock(flowsheet_costing_block=m.fs.costing)\n",
+    "# Cost the process\n",
+    "m.fs.costing.cost_process()\n",
+    "\n",
+    "# track water production rate \n",
+    "m.fs.costing.add_annual_water_production(m.fs.product.properties[0].flow_vol)\n",
+    "\n",
+    "# add LCOW\n",
+    "m.fs.costing.add_LCOW(m.fs.product.properties[0].flow_vol)\n",
+    "# init our costing\n",
+    "m.fs.costing.initialize()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Let's see how adjusting our Sh value impact process cost"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Let's explore how sh impacts cost \n",
+    "multipliers = np.linspace(1,2,5)\n",
+    "pressure = []\n",
+    "avg_sh =[] \n",
+    "lcow = []\n",
+    "for adj in multipliers:\n",
+    "    m.fs.Sh_multiplier.fix(adj)\n",
+    "    assert degrees_of_freedom(m) == 0\n",
+    "    result =solver.solve(m, tee=False)\n",
+    "    pressure.append(m.fs.RO.inlet.pressure[0].value/1e5)\n",
+    "    sh_avg=np.average([m.fs.RO.feed_side.N_Sh_comp[t,x,j].value for (t,x,j) in m.fs.RO.feed_side.N_Sh_comp])\n",
+    "    avg_sh.append(sh_avg)\n",
+    "    lcow.append(m.fs.costing.LCOW.value)\n",
+    "    print(\"Solved multiplier {}\".format(adj))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Plot LCOW result vs. Sh multiplier "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# plot LCOW result \n",
+    "fig, ax = plt.subplots()\n",
+    "ax.plot(multipliers,lcow)\n",
+    "ax.set(xlabel='Sh multiplier',ylabel='LCOW ($\\$$/m$^3$)')\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Lets add a pressure factor to membrane cost\n",
+    "\n",
+    "Typical pressure factors relate capital cost of component to operating pressure \n",
+    "\n",
+    "Here we will relate membrane cost to RO inlet pressure ($P_{inlet}$) such that:\n",
+    "$$\n",
+    "    Membrane\\,cost = 30 \\frac{\\$}{m^2}\\, for\\, P_{inlet} < 60 bar \n",
+    "$$\n",
+    "$$\n",
+    "    Membrane\\,cost = 30\\frac{\\$}{m^2}* \\frac{P_{inlet}}{60}\\, bar\\, for\\, P_{inlet} > 60 bar \n",
+    "$$\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# add a pressure factor to membrane cost \n",
+    "m.fs.base_membrane_cost=Var(initialize=30)\n",
+    "m.fs.base_membrane_cost.fix()\n",
+    "set_scaling_factor(m.fs.base_membrane_cost, 1/m.fs.base_membrane_cost.value)\n",
+    "\n",
+    "# The cost constraint \n",
+    "m.fs.RO_cost_pressure_constraint=(\n",
+    "    Constraint(expr=m.fs.costing.reverse_osmosis.membrane_cost==\n",
+    "    idaesMath.smooth_max(m.fs.base_membrane_cost,(m.fs.base_membrane_cost*(m.fs.RO.inlet.pressure[0]/(60*1e5))))))\n",
+    "\n",
+    "m.fs.costing.reverse_osmosis.membrane_cost.unfix()\n",
+    "\n",
+    "# Lets check that our function works and plot the result \n",
+    "pressures_for_testing=np.linspace(10,350,100)\n",
+    "mem_cost=[]\n",
+    "for pressure in pressures_for_testing:\n",
+    "    m.fs.RO.inlet.pressure[0]=pressure*1e5 # needs to be in kPa\n",
+    "    calculate_variable_from_constraint(m.fs.costing.reverse_osmosis.membrane_cost,m.fs.RO_cost_pressure_constraint)\n",
+    "    mem_cost.append(m.fs.costing.reverse_osmosis.membrane_cost.value)\n",
+    "\n",
+    "# plot the cost as function of pressure \n",
+    "fig, ax = plt.subplots()\n",
+    "ax.plot(pressures_for_testing,mem_cost)\n",
+    "ax.set(xlabel='Pressure (bar)',ylabel='Membrane cost ($\\$$/m$^2$)')\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Lets check how the new cost function impacts LCOW vs. fixed membrane cost"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Define concentrations to solve over \n",
+    "# solving in reverse, as model was previously solved for 100 g/L condition\n",
+    "# this should speed up solving processes\n",
+    "concentrations = np.linspace(100,15,10)\n",
+    "\n",
+    "# arrays to store result\n",
+    "pressures = []\n",
+    "lcow_fixed_mem_cost=[]\n",
+    "lcow_variable_mem_cost=[]\n",
+    "mem_cost_fixed=[]\n",
+    "mem_cost_variable=[]\n",
+    "pressures=[]\n",
+    "for con in concentrations:\n",
+    "\n",
+    "    m.fs.feed.properties[0].conc_mass_phase_comp['Liq','TDS'].fix(con)\n",
+    "\n",
+    "    # first solve for fixed cost\n",
+    "    m.fs.costing.reverse_osmosis.membrane_cost.fix(30)\n",
+    "    m.fs.RO_cost_pressure_constraint.deactivate()\n",
+    "    assert degrees_of_freedom(m) == 0\n",
+    "    result = solver.solve(m, tee=False)\n",
+    "    assert_optimal_termination(result)\n",
+    "    pressures.append(m.fs.RO.inlet.pressure[0].value/1e5)    \n",
+    "    lcow_fixed_mem_cost.append(m.fs.costing.LCOW.value)\n",
+    "    mem_cost_fixed.append(m.fs.costing.reverse_osmosis.membrane_cost.value)\n",
+    "    # second unfix our membrane cost and activate variable cost constraint and solve\n",
+    "    m.fs.costing.reverse_osmosis.membrane_cost.unfix()\n",
+    "    m.fs.RO_cost_pressure_constraint.activate()\n",
+    "    assert degrees_of_freedom(m) == 0\n",
+    "    result = solver.solve(m, tee=False)\n",
+    "    assert_optimal_termination(result)\n",
+    "    lcow_variable_mem_cost.append(m.fs.costing.LCOW.value)\n",
+    "    mem_cost_variable.append(m.fs.costing.reverse_osmosis.membrane_cost.value)\n",
+    "\n",
+    "    print(\"Solved con {}\".format(con))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Let's plot LCOW with variable cost result"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# plot LCOW result \n",
+    "fig, ax = plt.subplots()\n",
+    "ax.plot(pressures,lcow_fixed_mem_cost,  color='black',label='Fixed membrane cost')\n",
+    "ax.plot(pressures,lcow_variable_mem_cost,  color='red',label='Membrane cost as function of pressure')\n",
+    "ax.set(xlabel='Pressure (bar)',ylabel='LCOW ($\\$$/m$^3$)')\n",
+    "plt.legend()\n",
+    "\n",
+    "# plot membrane cost\n",
+    "fig, ax = plt.subplots()\n",
+    "ax.plot(pressures,mem_cost_fixed, color='black', label='Fixed membrane cost')\n",
+    "ax.plot(pressures,mem_cost_variable,  color='red',label='Membrane cost as function of pressure')\n",
+    "ax.set(xlabel='Pressure (bar)',ylabel='Membrane cost ($\\$$/m$^2$)')\n",
+    "plt.legend()\n",
+    "plt.show()"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.9.16"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}