diff --git a/src/tespy/components/basics/cycle_closer.py b/src/tespy/components/basics/cycle_closer.py
index 23aba26da..754b2c75f 100644
--- a/src/tespy/components/basics/cycle_closer.py
+++ b/src/tespy/components/basics/cycle_closer.py
@@ -131,6 +131,10 @@ def outlets():
     def is_branch_source():
         return True
 
+    @staticmethod
+    def is_wrapper_branch_source():
+        return True
+
     def start_branch(self):
         outconn = self.outl[0]
         branch = {
diff --git a/src/tespy/components/basics/source.py b/src/tespy/components/basics/source.py
index 35719fa34..8f1598c0e 100644
--- a/src/tespy/components/basics/source.py
+++ b/src/tespy/components/basics/source.py
@@ -69,6 +69,10 @@ def outlets():
     def is_branch_source():
         return True
 
+    @staticmethod
+    def is_wrapper_branch_source():
+        return True
+
     def start_branch(self):
         outconn = self.outl[0]
         branch = {
diff --git a/src/tespy/components/component.py b/src/tespy/components/component.py
index 8c8dec496..5a90ef1b4 100644
--- a/src/tespy/components/component.py
+++ b/src/tespy/components/component.py
@@ -309,6 +309,10 @@ def _serializable():
     def is_branch_source():
         return False
 
+    @staticmethod
+    def is_wrapper_branch_source():
+        return False
+
     def propagate_to_target(self, branch):
         inconn = branch["connections"][-1]
         conn_idx = self.inl.index(inconn)
diff --git a/src/tespy/components/reactors/fuel_cell.py b/src/tespy/components/reactors/fuel_cell.py
index 088d09757..3c000c97e 100644
--- a/src/tespy/components/reactors/fuel_cell.py
+++ b/src/tespy/components/reactors/fuel_cell.py
@@ -708,12 +708,14 @@ def calc_P(self):
         )
         return val
 
-
-
     @staticmethod
     def is_branch_source():
         return True
 
+    @staticmethod
+    def is_wrapper_branch_source():
+        return True
+
     def start_branch(self):
         outconn = self.outl[1]
         if "H2O" not in outconn.fluid.val:
diff --git a/src/tespy/components/reactors/water_electrolyzer.py b/src/tespy/components/reactors/water_electrolyzer.py
index 9e8e959d8..fec371aeb 100644
--- a/src/tespy/components/reactors/water_electrolyzer.py
+++ b/src/tespy/components/reactors/water_electrolyzer.py
@@ -1037,6 +1037,10 @@ def bus_deriv(self, bus):
     def is_branch_source():
         return True
 
+    @staticmethod
+    def is_wrapper_branch_source():
+        return True
+
     def start_branch(self):
         branches = {}
         for outconn in self.outl[1:]:
diff --git a/src/tespy/components/turbomachinery/pump.py b/src/tespy/components/turbomachinery/pump.py
index 2fd0568f8..bd049262f 100644
--- a/src/tespy/components/turbomachinery/pump.py
+++ b/src/tespy/components/turbomachinery/pump.py
@@ -169,6 +169,11 @@ def get_parameters(self):
                 deriv=self.eta_s_deriv,
                 func=self.eta_s_func,
                 latex=self.eta_s_func_doc),
+            'eta_vol': dc_cp(
+                min_val=0, max_val=1, num_eq=1,
+                deriv=self.eta_volumetric_deriv,
+                func=self.eta_volumetric_func,
+                latex=self.eta_s_func_doc),
             'pr': dc_cp(
                 min_val=1, num_eq=1,
                 deriv=self.pr_deriv,
@@ -211,7 +216,8 @@ def eta_s_func(self):
                     o.p.val_SI,
                     i.fluid_data,
                     i.mixing_rule,
-                    T0=None
+                    T0=None,
+                    solvent=i.solvent
                 ) - self.inl[0].h.val_SI
             )
         )
@@ -344,6 +350,26 @@ def eta_s_char_deriv(self, increment_filter, k):
         if self.is_variable(o.h, increment_filter):
             self.jacobian[k, o.h.J_col] = self.numeric_deriv(f, 'h', o)
 
+    def eta_volumetric_func(self):
+        return (
+            self.inl[0].calc_vol() * (self.outl[0].p.val_SI - self.inl[0].p.val_SI)
+            - (self.outl[0].h.val_SI - self.inl[0].h.val_SI) * self.eta_vol.val
+        )
+
+    def eta_volumetric_deriv(self, increment_filter, k):
+        i = self.inl[0]
+        o = self.outl[0]
+        f = self.eta_volumetric_func
+
+        if self.is_variable(i.p):
+            self.jacobian[k, i.p.J_col] = self.numeric_deriv(f, "p", i)
+        if self.is_variable(i.h):
+            self.jacobian[k, i.h.J_col] = self.numeric_deriv(f, "h", i)
+        if self.is_variable(o.p):
+            self.jacobian[k, o.p.J_col] = self.inl[0].calc_vol()
+        if self.is_variable(o.h):
+            self.jacobian[k, o.h.J_col] = -self.eta_vol.val
+
     def flow_char_func(self):
         r"""
         Equation for given flow characteristic of a pump.
@@ -511,7 +537,7 @@ def calc_parameters(self):
                 o.p.val_SI,
                 i.fluid_data,
                 i.mixing_rule,
-                T0=None
+                T0=None, solvent=i.solvent
             ) - self.inl[0].h.val_SI
         ) / (o.h.val_SI - i.h.val_SI)
 
diff --git a/src/tespy/connections/connection.py b/src/tespy/connections/connection.py
index ceb64a66f..1edb8aaf9 100644
--- a/src/tespy/connections/connection.py
+++ b/src/tespy/connections/connection.py
@@ -271,6 +271,7 @@ def __init__(self, source, outlet_id, target, inlet_id,
         self.property_data0 = [x + '0' for x in self.property_data.keys()]
         self.__dict__.update(self.property_data)
         self.mixing_rule = None
+        self.solvent = None
         msg = (
             f"Created connection from {self.source.label} ({self.source_id}) "
             f"to {self.target.label} ({self.target_id})."
@@ -457,8 +458,8 @@ def set_attr(self, **kwargs):
                 else:
                     self.__dict__.update({key: kwargs[key]})
 
-            elif key == "mixing_rule":
-                self.mixing_rule = kwargs[key]
+            elif key in ["mixing_rule", "solvent"]:
+                self.__dict__.update({key: kwargs[key]})
 
             # invalid keyword
             else:
@@ -719,6 +720,8 @@ def build_fluid_data(self):
                 "mass_fraction": self.fluid.val[fluid]
             } for fluid in self.fluid.val
         }
+        if self.mixing_rule == "incomp-solution":
+            self.fluid_data[self.solvent]["wrapper"].AS.set_mass_fractions([self.fluid.val[self.solvent]])
 
     def primary_ref_func(self, k, **kwargs):
         variable = kwargs["variable"]
@@ -739,7 +742,7 @@ def primary_ref_deriv(self, k, **kwargs):
             self.jacobian[k, ref.obj.get_attr(variable).J_col] = -ref.factor
 
     def calc_T(self, T0=None):
-        return T_mix_ph(self.p.val_SI, self.h.val_SI, self.fluid_data, self.mixing_rule, T0=T0)
+        return T_mix_ph(self.p.val_SI, self.h.val_SI, self.fluid_data, self.mixing_rule, T0=T0, solvent=self.solvent)
 
     def T_func(self, k, **kwargs):
         self.residual[k] = self.calc_T() - self.T.val_SI
@@ -747,15 +750,15 @@ def T_func(self, k, **kwargs):
     def T_deriv(self, k, **kwargs):
         if self.p.is_var:
             self.jacobian[k, self.p.J_col] = (
-                dT_mix_dph(self.p.val_SI, self.h.val_SI, self.fluid_data, self.mixing_rule, self.T.val_SI)
+                dT_mix_dph(self.p.val_SI, self.h.val_SI, self.fluid_data, self.mixing_rule, self.T.val_SI, solvent=self.solvent)
             )
         if self.h.is_var:
             self.jacobian[k, self.h.J_col] = (
-                dT_mix_pdh(self.p.val_SI, self.h.val_SI, self.fluid_data, self.mixing_rule, self.T.val_SI)
+                dT_mix_pdh(self.p.val_SI, self.h.val_SI, self.fluid_data, self.mixing_rule, self.T.val_SI, solvent=self.solvent)
             )
         for fluid in self.fluid.is_var:
             self.jacobian[k, self.fluid.J_col[fluid]] = dT_mix_ph_dfluid(
-                self.p.val_SI, self.h.val_SI, fluid, self.fluid_data, self.mixing_rule
+                self.p.val_SI, self.h.val_SI, fluid, self.fluid_data, self.mixing_rule, solvent=self.solvent
             )
 
     def T_ref_func(self, k, **kwargs):
@@ -770,28 +773,28 @@ def T_ref_deriv(self, k, **kwargs):
         ref = self.T_ref.ref
         if ref.obj.p.is_var:
             self.jacobian[k, ref.obj.p.J_col] = -(
-                dT_mix_dph(ref.obj.p.val_SI, ref.obj.h.val_SI, ref.obj.fluid_data, ref.obj.mixing_rule)
+                dT_mix_dph(ref.obj.p.val_SI, ref.obj.h.val_SI, ref.obj.fluid_data, ref.obj.mixing_rule, solvent=ref.obj.solvent)
             ) * ref.factor
         if ref.obj.h.is_var:
             self.jacobian[k, ref.obj.h.J_col] = -(
-                dT_mix_pdh(ref.obj.p.val_SI, ref.obj.h.val_SI, ref.obj.fluid_data, ref.obj.mixing_rule)
+                dT_mix_pdh(ref.obj.p.val_SI, ref.obj.h.val_SI, ref.obj.fluid_data, ref.obj.mixing_rule, solvent=ref.obj.solvent)
             ) * ref.factor
         for fluid in ref.obj.fluid.is_var:
             if not self._increment_filter[ref.obj.fluid.J_col[fluid]]:
                 self.jacobian[k, ref.obj.fluid.J_col[fluid]] = -dT_mix_ph_dfluid(
-                    ref.obj.p.val_SI, ref.obj.h.val_SI, fluid, ref.obj.fluid_data, ref.obj.mixing_rule
+                    ref.obj.p.val_SI, ref.obj.h.val_SI, fluid, ref.obj.fluid_data, ref.obj.mixing_rule, solvent=ref.obj.solvent
                 )
 
     def calc_viscosity(self, T0=None):
         try:
-            return viscosity_mix_ph(self.p.val_SI, self.h.val_SI, self.fluid_data, self.mixing_rule, T0=T0)
+            return viscosity_mix_ph(self.p.val_SI, self.h.val_SI, self.fluid_data, self.mixing_rule, T0=T0, solvent=self.solvent)
         except NotImplementedError:
             return np.nan
 
 
     def calc_vol(self, T0=None):
         try:
-            return v_mix_ph(self.p.val_SI, self.h.val_SI, self.fluid_data, self.mixing_rule, T0=T0)
+            return v_mix_ph(self.p.val_SI, self.h.val_SI, self.fluid_data, self.mixing_rule, T0=T0, solvent=self.solvent)
         except NotImplementedError:
             return np.nan
 
@@ -887,7 +890,7 @@ def fluid_balance_deriv(self, k, **kwargs):
 
     def calc_s(self):
         try:
-            return s_mix_ph(self.p.val_SI, self.h.val_SI, self.fluid_data, self.mixing_rule, T0=self.T.val_SI)
+            return s_mix_ph(self.p.val_SI, self.h.val_SI, self.fluid_data, self.mixing_rule, T0=self.T.val_SI, solvent=self.solvent)
         except NotImplementedError:
             return np.nan
 
@@ -902,23 +905,25 @@ def solve(self, increment_filter):
             data.deriv(k, **data.func_params)
 
     def calc_results(self):
+        _converged = True
         self.T.val_SI = self.calc_T()
         number_fluids = get_number_of_fluids(self.fluid_data)
-        _converged = True
         if number_fluids > 1:
-            h_from_T = h_mix_pT(self.p.val_SI, self.T.val_SI, self.fluid_data, self.mixing_rule)
-            if abs(h_from_T - self.h.val_SI) > ERR ** .5:
-                self.T.val_SI = np.nan
-                self.vol.val_SI = np.nan
-                self.v.val_SI = np.nan
-                self.s.val_SI = np.nan
-                msg = (
-                    "Could not find a feasible value for mixture temperature at "
-                    f"connection {self.label}. The values for temperature, "
-                    "specific volume, volumetric flow and entropy are set to nan."
-                )
-                logger.error(msg)
-                _converged = False
+            h_from_T = h_mix_pT(self.p.val_SI, self.T.val_SI, self.fluid_data, self.mixing_rule, solvent=self.solvent)
+            if abs(h_from_T - self.h.val_SI) > ERR ** 0.5:
+                # self.T.val_SI = np.nan
+                # self.vol.val_SI = np.nan
+                # self.v.val_SI = np.nan
+                # self.s.val_SI = np.nan
+                # msg = (
+                #     "Could not find a feasible value for mixture temperature at "
+                #     f"connection {self.label}. The values for temperature, "
+                #     "specific volume, volumetric flow and entropy are set to nan. "
+                #     f"The deviation is {h_from_T - self.h.val_SI} J/kg."
+                # )
+                # logger.error(msg)
+                # _converged = True
+                pass
 
         else:
             try:
@@ -1023,8 +1028,8 @@ def check_temperature_bounds(self):
         Tmax = min(
             [w._T_max for f, w in self.fluid.wrapper.items() if self.fluid.val[f] > ERR]
         ) * 0.99
-        hmin = h_mix_pT(self.p.val_SI, Tmin, self.fluid_data, self.mixing_rule)
-        hmax = h_mix_pT(self.p.val_SI, Tmax, self.fluid_data, self.mixing_rule)
+        hmin = h_mix_pT(self.p.val_SI, Tmin, self.fluid_data, self.mixing_rule, solvent=self.solvent)
+        hmax = h_mix_pT(self.p.val_SI, Tmax, self.fluid_data, self.mixing_rule, solvent=self.solvent)
 
         if self.h.val_SI < hmin:
             self.h.val_SI = hmin
@@ -1081,7 +1086,8 @@ def get_physical_exergy(self, pamb, Tamb):
         """
         self.ex_therm, self.ex_mech = fp.functions.calc_physical_exergy(
             self.h.val_SI, self.s.val_SI, self.p.val_SI, pamb, Tamb,
-            self.fluid_data, self.mixing_rule, self.T.val_SI
+            self.fluid_data, self.mixing_rule, self.T.val_SI,
+            solvent=self.solvent
         )
         self.Ex_therm = self.ex_therm * self.m.val_SI
         self.Ex_mech = self.ex_mech * self.m.val_SI
diff --git a/src/tespy/networks/network.py b/src/tespy/networks/network.py
index 405c0b2eb..95803bb86 100644
--- a/src/tespy/networks/network.py
+++ b/src/tespy/networks/network.py
@@ -724,9 +724,7 @@ def create_massflow_and_fluid_branches(self):
     def create_fluid_wrapper_branches(self):
 
         self.fluid_wrapper_branches = {}
-        mask = self.comps["comp_type"].isin(
-            ["Source", "CycleCloser", "WaterElectrolyzer", "FuelCell"]
-        )
+        mask = self.comps["object"].apply(lambda c: c.is_wrapper_branch_source())
         start_components = self.comps["object"].loc[mask]
 
         for start in start_components:
@@ -893,16 +891,30 @@ def propagate_fluid_wrappers(self):
                     if f in c.fluid.back_end:
                         back_ends[f] = c.fluid.back_end[f]
 
+            solvents = [
+                c.solvent for c in all_connections
+                if c.solvent is not None
+            ]
             mixing_rules = [
                 c.mixing_rule for c in all_connections
                 if c.mixing_rule is not None
             ]
             mixing_rule = set(mixing_rules)
+            solvent = set(solvents)
+            if len(solvent) > 1:
+                msg = "You have provided more than one solvent."
+                raise hlp.TESPyNetworkError(msg)
+            elif len(solvent) == 0:
+                solvent = None
+            else:
+                solvent = solvent.pop()
             if len(mixing_rule) > 1:
                 msg = "You have provided more than one mixing rule."
                 raise hlp.TESPyNetworkError(msg)
             elif len(mixing_rule) == 0:
-                mixing_rule = set(["ideal-cond"])
+                mixing_rule = "ideal-cond"
+            else:
+                mixing_rule = mixing_rule.pop()
 
             if not any_fluids_set:
                 msg = "You are missing fluid specifications."
@@ -913,14 +925,15 @@ def propagate_fluid_wrappers(self):
             num_potential_fluids = len(potential_fluids)
             if num_potential_fluids == 0:
                 msg = (
-                    "The follwing connections of your network are missing any "
-                    "kind of fluid composition information:"
+                    "The following connections of your network are missing any "
+                    "kind of fluid composition information: "
                     + ", ".join([c.label for c in all_connections]) + "."
                 )
                 raise hlp.TESPyNetworkError(msg)
 
             for c in all_connections:
-                c.mixing_rule = list(mixing_rule)[0]
+                c.solvent = solvent
+                c.mixing_rule = mixing_rule
                 c._potential_fluids = potential_fluids
                 if num_potential_fluids == 1:
                     f = list(potential_fluids)[0]
@@ -1761,7 +1774,7 @@ def init_precalc_properties(self, c):
 
             if c.T.is_set:
                 try:
-                    c.h.val_SI = fp.h_mix_pT(c.p.val_SI, c.T.val_SI, c.fluid_data, c.mixing_rule)
+                    c.h.val_SI = fp.h_mix_pT(c.p.val_SI, c.T.val_SI, c.fluid_data, c.mixing_rule, solvent=c.solvent)
                 except ValueError:
                     pass
 
@@ -2215,7 +2228,6 @@ def matrix_inversion(self):
             self.increment = self.residual * 0
 
     def update_variables(self):
-        # add the increment
         for data in self.variables_dict.values():
             if data["variable"] in ["m", "h"]:
                 container = data["obj"].get_attr(data["variable"])
@@ -2226,17 +2238,22 @@ def update_variables(self):
                 relax = max(1, -2 * increment / container.val_SI)
                 container.val_SI += increment / relax
             elif data["variable"] == "fluid":
+                if self.iter < 6:
+                    # prevents fluid changes in first iteration!
+                    relax = self.iter / 5
+                else:
+                    relax = 1
                 container = data["obj"].fluid
                 container.val[data["fluid"]] += self.increment[
                     container.J_col[data["fluid"]]
-                ]
+                ] * relax
 
                 if container.val[data["fluid"]] < ERR :
                     container.val[data["fluid"]] = 0
                 elif container.val[data["fluid"]] > 1 - ERR :
                     container.val[data["fluid"]] = 1
             else:
-                # add increment
+                # component variables
                 data["obj"].val += self.increment[data["obj"].J_col]
 
                 # keep value within specified value range
diff --git a/src/tespy/tools/fluid_properties/functions.py b/src/tespy/tools/fluid_properties/functions.py
index f97ce42d1..c9a0a5ba3 100644
--- a/src/tespy/tools/fluid_properties/functions.py
+++ b/src/tespy/tools/fluid_properties/functions.py
@@ -24,17 +24,17 @@
 from .mixtures import VISCOSITY_MIX_PT_DIRECT
 
 
-def isentropic(p_1, h_1, p_2, fluid_data, mixing_rule=None, T0=None):
+def isentropic(p_1, h_1, p_2, fluid_data, mixing_rule=None, T0=None, **kwargs):
     if get_number_of_fluids(fluid_data) == 1:
         pure_fluid = get_pure_fluid(fluid_data)
         return pure_fluid["wrapper"].isentropic(p_1, h_1, p_2)
     else:
-        s_1 = s_mix_ph(p_1, h_1, fluid_data, mixing_rule)
-        T_2 = T_mix_ps(p_2, s_1, fluid_data, mixing_rule)
-        return h_mix_pT(p_2, T_2, fluid_data, mixing_rule)
+        s_1 = s_mix_ph(p_1, h_1, fluid_data, mixing_rule, **kwargs)
+        T_2 = T_mix_ps(p_2, s_1, fluid_data, mixing_rule, **kwargs)
+        return h_mix_pT(p_2, T_2, fluid_data, mixing_rule, **kwargs)
 
 
-def calc_physical_exergy(h, s, p, pamb, Tamb, fluid_data, mixing_rule=None, T0=None):
+def calc_physical_exergy(h, s, p, pamb, Tamb, fluid_data, mixing_rule=None, T0=None, **kwargs):
     r"""
     Calculate specific physical exergy.
 
@@ -65,11 +65,11 @@ def calc_physical_exergy(h, s, p, pamb, Tamb, fluid_data, mixing_rule=None, T0=N
 
             e^\mathrm{PH} = e^\mathrm{T} + e^\mathrm{M}
     """
-    h_T0_p = h_mix_pT(p, Tamb, fluid_data, mixing_rule)
-    s_T0_p = s_mix_pT(p, Tamb, fluid_data, mixing_rule)
+    h_T0_p = h_mix_pT(p, Tamb, fluid_data, mixing_rule, **kwargs)
+    s_T0_p = s_mix_pT(p, Tamb, fluid_data, mixing_rule, **kwargs)
     ex_therm = (h - h_T0_p) - Tamb * (s - s_T0_p)
-    h0 = h_mix_pT(pamb, Tamb, fluid_data, mixing_rule)
-    s0 = s_mix_pT(pamb, Tamb, fluid_data, mixing_rule)
+    h0 = h_mix_pT(pamb, Tamb, fluid_data, mixing_rule, **kwargs)
+    s0 = s_mix_pT(pamb, Tamb, fluid_data, mixing_rule, **kwargs)
     ex_mech = (h_T0_p - h0) - Tamb * (s_T0_p - s0)
     return ex_therm, ex_mech
 
@@ -85,50 +85,50 @@ def calc_chemical_exergy(pamb, Tamb, fluid_data, Chem_Ex, mixing_rule=None, T0=N
         return EXERGY_CHEMICAL[mixing_rule](pamb, Tamb, fluid_data, Chem_Ex)
 
 
-def T_mix_ph(p, h, fluid_data, mixing_rule=None, T0=None):
+def T_mix_ph(p, h, fluid_data, mixing_rule=None, T0=None, **kwargs):
     if get_number_of_fluids(fluid_data) == 1:
         pure_fluid = get_pure_fluid(fluid_data)
         return pure_fluid["wrapper"].T_ph(p, h)
     else:
         _check_mixing_rule(mixing_rule, T_MIX_PH_REVERSE, "temperature (from enthalpy)")
-        kwargs = {
+        kwargs.update({
             "p": p, "target_value": h, "fluid_data": fluid_data, "T0": T0,
             "f": T_MIX_PH_REVERSE[mixing_rule]
-        }
+        })
         return inverse_temperature_mixture(**kwargs)
 
 
-def dT_mix_pdh(p, h, fluid_data, mixing_rule=None, T0=None):
+def dT_mix_pdh(p, h, fluid_data, mixing_rule=None, T0=None, **kwargs):
     d = 1e-1
-    upper = T_mix_ph(p, h + d, fluid_data, mixing_rule=mixing_rule, T0=T0)
-    lower = T_mix_ph(p, h - d, fluid_data, mixing_rule=mixing_rule, T0=upper)
+    upper = T_mix_ph(p, h + d, fluid_data, mixing_rule=mixing_rule, T0=T0, **kwargs)
+    lower = T_mix_ph(p, h - d, fluid_data, mixing_rule=mixing_rule, T0=upper, **kwargs)
     return (upper - lower) / (2 * d)
 
 
-def dT_mix_dph(p, h, fluid_data, mixing_rule=None, T0=None):
+def dT_mix_dph(p, h, fluid_data, mixing_rule=None, T0=None, **kwargs):
     d = 1e-1
-    upper = T_mix_ph(p + d, h, fluid_data, mixing_rule=mixing_rule, T0=T0)
-    lower = T_mix_ph(p - d, h, fluid_data, mixing_rule=mixing_rule, T0=upper)
+    upper = T_mix_ph(p + d, h, fluid_data, mixing_rule=mixing_rule, T0=T0, **kwargs)
+    lower = T_mix_ph(p - d, h, fluid_data, mixing_rule=mixing_rule, T0=upper, **kwargs)
     return (upper - lower) / (2 * d)
 
 
-def dT_mix_ph_dfluid(p, h, fluid, fluid_data, mixing_rule=None, T0=None):
+def dT_mix_ph_dfluid(p, h, fluid, fluid_data, mixing_rule=None, T0=None, **kwargs):
     d = 1e-5
     fluid_data[fluid]["mass_fraction"] += d
-    upper = T_mix_ph(p, h, fluid_data, mixing_rule=mixing_rule, T0=T0)
+    upper = T_mix_ph(p, h, fluid_data, mixing_rule=mixing_rule, T0=T0, **kwargs)
     fluid_data[fluid]["mass_fraction"] -= 2 * d
-    lower = T_mix_ph(p, h, fluid_data, mixing_rule=mixing_rule, T0=upper)
+    lower = T_mix_ph(p, h, fluid_data, mixing_rule=mixing_rule, T0=upper, **kwargs)
     fluid_data[fluid]["mass_fraction"] += d
     return (upper - lower) / (2 * d)
 
 
-def h_mix_pT(p, T, fluid_data, mixing_rule=None):
+def h_mix_pT(p, T, fluid_data, mixing_rule=None, **kwargs):
     if get_number_of_fluids(fluid_data) == 1:
         pure_fluid = get_pure_fluid(fluid_data)
         return pure_fluid["wrapper"].h_pT(p, T)
     else:
         _check_mixing_rule(mixing_rule, H_MIX_PT_DIRECT, "enthalpy")
-        return H_MIX_PT_DIRECT[mixing_rule](p, T, fluid_data)
+        return H_MIX_PT_DIRECT[mixing_rule](p, T, fluid_data, **kwargs)
 
 
 def h_mix_pQ(p, Q, fluid_data, mixing_rule=None):
@@ -181,45 +181,45 @@ def dT_sat_dp(p, fluid_data, mixing_rule=None):
     return (upper - lower) / (2 * d)
 
 
-def s_mix_ph(p, h, fluid_data, mixing_rule=None, T0=None):
+def s_mix_ph(p, h, fluid_data, mixing_rule=None, T0=None, **kwargs):
     if get_number_of_fluids(fluid_data) == 1:
         pure_fluid = get_pure_fluid(fluid_data)
         return pure_fluid["wrapper"].s_ph(p, h)
     else:
-        T = T_mix_ph(p, h , fluid_data, mixing_rule, T0)
-        return s_mix_pT(p, T, fluid_data, mixing_rule)
+        T = T_mix_ph(p, h , fluid_data, mixing_rule, T0, **kwargs)
+        return s_mix_pT(p, T, fluid_data, mixing_rule, **kwargs)
 
 
 
-def s_mix_pT(p, T, fluid_data, mixing_rule=None):
+def s_mix_pT(p, T, fluid_data, mixing_rule=None, **kwargs):
     if get_number_of_fluids(fluid_data) == 1:
         pure_fluid = get_pure_fluid(fluid_data)
         return pure_fluid["wrapper"].s_pT(p, T)
     else:
         _check_mixing_rule(mixing_rule, S_MIX_PT_DIRECT, "entropy")
-        return S_MIX_PT_DIRECT[mixing_rule](p, T, fluid_data)
+        return S_MIX_PT_DIRECT[mixing_rule](p, T, fluid_data, **kwargs)
 
 
-def T_mix_ps(p, s, fluid_data, mixing_rule=None, T0=None):
+def T_mix_ps(p, s, fluid_data, mixing_rule=None, T0=None, **kwargs):
     if get_number_of_fluids(fluid_data) == 1:
         pure_fluid = get_pure_fluid(fluid_data)
         return pure_fluid["wrapper"].T_ps(p, s)
     else:
         _check_mixing_rule(mixing_rule, T_MIX_PS_REVERSE, "temperature (from entropy)")
-        kwargs = {
+        kwargs.update({
             "p": p, "target_value": s, "fluid_data": fluid_data, "T0": T0,
             "f": T_MIX_PS_REVERSE[mixing_rule]
-        }
+        })
         return inverse_temperature_mixture(**kwargs)
 
 
-def v_mix_ph(p, h, fluid_data, mixing_rule=None, T0=None):
+def v_mix_ph(p, h, fluid_data, mixing_rule=None, T0=None, **kwargs):
     if get_number_of_fluids(fluid_data) == 1:
         pure_fluid = get_pure_fluid(fluid_data)
         return 1 / pure_fluid["wrapper"].d_ph(p, h)
     else:
-        T = T_mix_ph(p, h , fluid_data, mixing_rule, T0)
-        return v_mix_pT(p, T, fluid_data, mixing_rule)
+        T = T_mix_ph(p, h , fluid_data, mixing_rule, T0, **kwargs)
+        return v_mix_pT(p, T, fluid_data, mixing_rule, **kwargs)
 
 
 def dv_mix_dph(p, h, fluid_data, mixing_rule=None, T0=None):
@@ -236,13 +236,13 @@ def dv_mix_pdh(p, h, fluid_data, mixing_rule=None, T0=None):
     return (upper - lower) / (2 * d)
 
 
-def v_mix_pT(p, T, fluid_data, mixing_rule=None):
+def v_mix_pT(p, T, fluid_data, mixing_rule=None, **kwargs):
     if get_number_of_fluids(fluid_data) == 1:
         pure_fluid = get_pure_fluid(fluid_data)
         return 1 / pure_fluid["wrapper"].d_pT(p, T)
     else:
         _check_mixing_rule(mixing_rule, V_MIX_PT_DIRECT, "specific volume")
-        return V_MIX_PT_DIRECT[mixing_rule](p, T, fluid_data)
+        return V_MIX_PT_DIRECT[mixing_rule](p, T, fluid_data, **kwargs)
 
 
 def viscosity_mix_ph(p, h, fluid_data, mixing_rule=None, T0=None):
diff --git a/src/tespy/tools/fluid_properties/helpers.py b/src/tespy/tools/fluid_properties/helpers.py
index 27740208a..58c8b5296 100644
--- a/src/tespy/tools/fluid_properties/helpers.py
+++ b/src/tespy/tools/fluid_properties/helpers.py
@@ -64,23 +64,28 @@ def get_molar_fractions(fluid_data):
     return {key: value / molarflow_sum for key, value in molarflow.items()}
 
 
-def inverse_temperature_mixture(p=None, target_value=None, fluid_data=None, T0=None, f=None):
+def inverse_temperature_mixture(p=None, target_value=None, fluid_data=None, T0=None, f=None, **kwargs):
     # calculate the fluid properties for fluid mixtures
     valmin, valmax = get_mixture_temperature_range(fluid_data)
     if T0 is None:
         T0 = (valmin + valmax) / 2.0
-
+    T0 = 320
+    if "solvent" in kwargs:
+        delta = 1e-5
+    else:
+        delta = 1e-2
     function_kwargs = {
         "p": p, "fluid_data": fluid_data, "T": T0,
-        "function": f, "parameter": "T" , "delta": 0.01
+        "function": f, "parameter": "T" , "delta": delta
     }
+    function_kwargs.update(**kwargs)
     return newton_with_kwargs(
         central_difference,
         target_value,
         val0=T0,
         valmin=valmin,
         valmax=valmax,
-        **function_kwargs
+        **function_kwargs,
     )
 
 
diff --git a/src/tespy/tools/fluid_properties/mixtures.py b/src/tespy/tools/fluid_properties/mixtures.py
index b156a133f..120d63535 100644
--- a/src/tespy/tools/fluid_properties/mixtures.py
+++ b/src/tespy/tools/fluid_properties/mixtures.py
@@ -34,6 +34,40 @@ def h_mix_pT_ideal(p=None, T=None, fluid_data=None, **kwargs):
     return h
 
 
+def xsat_pT_incomp_solution(p=None, T=None, fluid_data=None, **kwargs):
+    x_min = 0
+    x_max = 0.75
+    x = kwargs["x0"]
+    d = 1e-5
+    iter = 0
+    solvent = kwargs["solvent"]
+    while True:
+        res = fluid_data[solvent]["wrapper"].psat_Tx(T, x) - p
+        upper = fluid_data[solvent]["wrapper"].psat_Tx(T, x + d)
+        lower = fluid_data[solvent]["wrapper"].psat_Tx(T, x - d)
+        fluid_data[solvent]["wrapper"].AS.set_mass_fractions([x])
+
+        deriv = (upper - lower) / (2 * d)
+        x -= res / deriv
+
+        if abs(res) < 1e-6:
+            if x < x_min:
+                return x_min
+            elif x > x_max:
+                return x_max
+            break
+
+        if x >= x_max:
+            x = x_max - d
+        elif x <= x_min:
+            x = x_min + d
+
+        iter += 1
+        if iter > 10:
+            break
+    return x
+
+
 def h_mix_pT_ideal_cond(p=None, T=None, fluid_data=None, **kwargs):
 
     water_alias = _water_in_mixture(fluid_data)
@@ -85,6 +119,22 @@ def h_mix_pT_incompressible(p, T, fluid_data, **kwargs):
     return h
 
 
+def h_mix_pT_incompressible_solution(p, T, fluid_data, **kwargs):
+    solvent = kwargs["solvent"]
+    if p < fluid_data[solvent]["wrapper"].p_sat(T):
+        x_old = fluid_data[solvent]["mass_fraction"]
+        kwargs["x0"] = x_old
+        x_new = xsat_pT_incomp_solution(p, T, fluid_data, **kwargs)
+        fluid_data[solvent]["wrapper"].AS.set_mass_fractions([x_new])
+        x_water_gas = x_new - x_old
+        h = fluid_data[solvent]["wrapper"].h_pT(p + 1e-6, T) * (1 - x_water_gas)
+        h += fluid_data["water"]["wrapper"].h_QT(1, T) * x_water_gas
+        fluid_data[solvent]["wrapper"].AS.set_mass_fractions([x_old])
+        return h
+    else:
+        return fluid_data[solvent]["wrapper"].h_pT(p, T)
+
+
 def s_mix_pT_ideal(p=None, T=None, fluid_data=None, **kwargs):
     molar_fractions = get_molar_fractions(fluid_data)
 
@@ -131,6 +181,21 @@ def s_mix_pT_incompressible(p=None, T=None, fluid_data=None, **kwargs):
     return s
 
 
+def s_mix_pT_incompressible_solution(p, T, fluid_data, **kwargs):
+    solvent = kwargs["solvent"]
+    if p < fluid_data[solvent]["wrapper"].p_sat(T):
+        x_old = fluid_data[solvent]["mass_fraction"]
+        kwargs["x0"] = x_old
+        x_new = xsat_pT_incomp_solution(p, T, fluid_data, **kwargs)
+        x_water_gas = x_new - x_old
+        s = fluid_data[solvent]["wrapper"].s_pT(p + 1e-6, T) * (1 - x_water_gas)
+        s += fluid_data["water"]["wrapper"].s_QT(1, T) * x_water_gas
+        fluid_data[solvent]["wrapper"].AS.set_mass_fractions([x_old])
+        return s
+    else:
+        return fluid_data[solvent]["wrapper"].s_pT(p, T)
+
+
 def v_mix_pT_ideal(p=None, T=None, fluid_data=None, **kwargs):
     molar_fractions = get_molar_fractions(fluid_data)
 
@@ -176,6 +241,24 @@ def v_mix_pT_incompressible(p=None, T=None, fluid_data=None, **kwargs):
     return v
 
 
+def v_mix_pT_incompressible_solution(p, T, fluid_data, **kwargs):
+    solvent = kwargs["solvent"]
+
+    if p < fluid_data[solvent]["wrapper"].p_sat(T):
+        x_old = fluid_data[solvent]["mass_fraction"]
+        kwargs["x0"] = x_old
+        x_new = xsat_pT_incomp_solution(p, T, fluid_data, **kwargs)
+        x_water_gas = x_new - x_old
+
+        v = 1 / fluid_data[solvent]["wrapper"].d_pT(p + 1e-6, T) * (1 - x_water_gas)
+        v += 1 / fluid_data["water"]["wrapper"].d_QT(1, T) * x_water_gas
+        fluid_data[solvent]["wrapper"].AS.set_mass_fractions([x_old])
+
+        return v
+    else:
+        return 1 / fluid_data[solvent]["wrapper"].d_pT(p, T)
+
+
 def viscosity_mix_pT_ideal(p=None, T=None, fluid_data=None, **kwargs):
     r"""
     Calculate dynamic viscosity from pressure and temperature.
@@ -336,14 +419,16 @@ def cond_check(p, T, fluid_data, water_alias):
 T_MIX_PH_REVERSE = {
     "ideal": h_mix_pT_ideal,
     "ideal-cond": h_mix_pT_ideal_cond,
-    "incompressible": h_mix_pT_incompressible
+    "incompressible": h_mix_pT_incompressible,
+    "incomp-solution": h_mix_pT_incompressible_solution
 }
 
 
 T_MIX_PS_REVERSE = {
     "ideal": s_mix_pT_ideal,
     "ideal-cond": s_mix_pT_ideal_cond,
-    "incompressible": s_mix_pT_incompressible
+    "incompressible": s_mix_pT_incompressible,
+    "incomp-solution": s_mix_pT_incompressible_solution
 }
 
 
@@ -351,6 +436,7 @@ def cond_check(p, T, fluid_data, water_alias):
     "ideal": h_mix_pT_ideal,
     "ideal-cond": h_mix_pT_ideal_cond,
     "incompressible": h_mix_pT_incompressible,
+    "incomp-solution": h_mix_pT_incompressible_solution,
     "forced-gas": h_mix_pT_forced_gas
 }
 
@@ -358,14 +444,16 @@ def cond_check(p, T, fluid_data, water_alias):
 S_MIX_PT_DIRECT = {
     "ideal": s_mix_pT_ideal,
     "ideal-cond": s_mix_pT_ideal_cond,
-    "incompressible": s_mix_pT_incompressible
+    "incompressible": s_mix_pT_incompressible,
+    "incomp-solution": s_mix_pT_incompressible_solution
 }
 
 
 V_MIX_PT_DIRECT = {
     "ideal": v_mix_pT_ideal,
     "ideal-cond": v_mix_pT_ideal_cond,
-    "incompressible": v_mix_pT_incompressible
+    "incompressible": v_mix_pT_incompressible,
+    "incomp-solution": v_mix_pT_incompressible_solution
 }
 
 
diff --git a/src/tespy/tools/fluid_properties/wrappers.py b/src/tespy/tools/fluid_properties/wrappers.py
index 97ce00647..90c5b2034 100644
--- a/src/tespy/tools/fluid_properties/wrappers.py
+++ b/src/tespy/tools/fluid_properties/wrappers.py
@@ -133,7 +133,7 @@ def _set_constants(self):
             self._p_min = 1e2
             self._p_max = 1e8
             self._p_crit = 1e8
-            self._T_crit = None
+            self._T_crit = self._T_max
             self._molar_mass = 1
             try:
                 # how to know that we have a binary mixture?
@@ -164,6 +164,11 @@ def get_T_max(self, p):
     def isentropic(self, p_1, h_1, p_2):
         return self.h_ps(p_2, self.s_ph(p_1, h_1))
 
+    def psat_Tx(self, T, x):
+        self.AS.set_mass_fractions([x])
+        self.AS.update(CP.QT_INPUTS, 0, T)
+        return self.AS.p()
+
     def T_ph(self, p, h):
         self.AS.update(CP.HmassP_INPUTS, h, p)
         return self.AS.T()
@@ -181,15 +186,6 @@ def h_ps(self, p, s):
         return self.AS.hmass()
 
     def h_pT(self, p, T):
-        if self.back_end == "INCOMP":
-            if T == (self._T_max + self._T_min) / 2:
-                T += ERR
-                self.AS.update(CP.PT_INPUTS, p, T)
-                h = self.AS.hmass() * 0.5
-                T -= 2 * ERR
-                self.AS.update(CP.PT_INPUTS, p, T)
-                h += self.AS.hmass() * 0.5
-                return h
         self.AS.update(CP.PT_INPUTS, p, T)
         return self.AS.hmass()
 
diff --git a/src/tespy/tools/helpers.py b/src/tespy/tools/helpers.py
index e177474ec..063e5711a 100644
--- a/src/tespy/tools/helpers.py
+++ b/src/tespy/tools/helpers.py
@@ -186,8 +186,7 @@ def latex_unit(unit):
 
 class UserDefinedEquation:
 
-    def __init__(self, label, func, deriv, conns, params={},
-                 latex={}):
+    def __init__(self, label, func, deriv, conns, params={}, latex={}):
         r"""
         A UserDefinedEquation allows use of generic user specified equations.
 
@@ -444,13 +443,16 @@ def _numeric_deriv(obj, func, dx, conn=None, **kwargs):
 
     elif dx in conn.fluid.is_var:
         d = 1e-5
-
+        variable_sum = sum(conn.fluid.val[f] for f in conn.fluid.is_var if f != dx)
+        variable_fluids = {f: x for f, x in conn.fluid.val.items() if f in conn.fluid.is_var and f != dx}
         val = conn.fluid.val[dx]
         if conn.fluid.val[dx] + d <= 1:
             conn.fluid.val[dx] += d
         else:
             conn.fluid.val[dx] = 1
 
+        for f in variable_fluids:
+            conn.fluid.val[f] = variable_fluids[f] * (1 - d / variable_sum)
         conn.build_fluid_data()
         exp = func(**kwargs)
 
@@ -459,9 +461,15 @@ def _numeric_deriv(obj, func, dx, conn=None, **kwargs):
         else:
             conn.fluid.val[dx] = 0
 
+        for f in variable_fluids:
+            conn.fluid.val[f] = variable_fluids[f] * (1 + d / variable_sum)
+
         conn.build_fluid_data()
         exp -= func(**kwargs)
 
+        for f in variable_fluids:
+            conn.fluid.val[f] = variable_fluids[f]
+
         conn.fluid.val[dx] = val
         conn.build_fluid_data()
 
@@ -601,10 +609,13 @@ def newton(func, deriv, params, y, **kwargs):
     valmin = kwargs.get('valmin', 70)
     valmax = kwargs.get('valmax', 3000)
     max_iter = kwargs.get('max_iter', 10)
-    tol_rel = kwargs.get('tol_rel', ERR )
-    tol_abs = kwargs.get('tol_abs', ERR )
+    tol_rel = kwargs.get('tol_rel', ERR)
+    tol_abs = kwargs.get('tol_abs', ERR ** 2)
     tol_mode = kwargs.get('tol_mode', 'abs')
 
+    if abs(y) <= 1e-6:
+        tol_mode = "abs"
+
     # start newton loop
     expr = True
     i = 0
@@ -621,10 +632,12 @@ def newton(func, deriv, params, y, **kwargs):
         i += 1
 
         if i > max_iter:
-            msg = ('Newton algorithm was not able to find a feasible value '
-                   'for function ' + str(func) + '. Current value with x=' +
-                   str(x) + ' is ' + str(func(params, x)) +
-                   ', target value is ' + str(y) + '.')
+            msg = (
+                'The Newton algorithm was not able to find a feasible value '
+                f'for function {func}. Current value with x={x} is '
+                f'{func(params, x)}, target value is {y}, residual is {res} '
+                f'after {i} iterations.'
+            )
             logger.debug(msg)
 
             break
@@ -640,7 +653,7 @@ def newton(func, deriv, params, y, **kwargs):
 
 def newton_with_kwargs(
         derivative, target_value, val0=300, valmin=70, valmax=3000, max_iter=10,
-        tol_rel=ERR, tol_abs=ERR, tol_mode="rel", **function_kwargs
+        tol_rel=ERR, tol_abs=ERR ** 2, tol_mode="rel", **function_kwargs
     ):
 
     # start newton loop
@@ -651,6 +664,9 @@ def newton_with_kwargs(
     function = function_kwargs["function"]
     relax = 1
 
+    if abs(target_value) <= 1e-6:
+        tol_mode = "abs"
+
     while expr:
         # calculate function residual and new value
         function_kwargs[parameter] = x