diff --git a/doc/interface/binding/generalized_ion_exchange.rst b/doc/interface/binding/generalized_ion_exchange.rst
index c12b2e8ca..98fa74733 100644
--- a/doc/interface/binding/generalized_ion_exchange.rst
+++ b/doc/interface/binding/generalized_ion_exchange.rst
@@ -106,38 +106,70 @@ Generalized Ion Exchange
 
 **Unit:** :math:`m_{MP}^{3} mol^{-1}`
 
-===================  =========================  
+===================  ========================= 
 **Type:** double     **Length:** NCOMP
 ===================  ========================= 
 
 ``GIEX_NU``
    Base value for characteristic charges of the protein; The number of
    sites :math:`\nu` that the protein interacts with on the resin
-   surface
+   surface. If more than NCOMP values are provided, :math:`\nu` is given
+   by a piecewise cubic polynomial. Data is expected in component-major
+   ordering.
 
-===================  =========================  
-**Type:** double     **Length:** NCOMP
-===================  ========================= 
+===================  ===============================  
+**Type:** double     **Length:** multiples of NCOMP
+===================  =============================== 
 
 ``GIEX_NU_LIN``
    Coefficient of linear dependence of characteristic charge on modifier
-   component
+   component. If more than NCOMP values are provided, :math:`\nu` is given
+   by a piecewise cubic polynomial. This parameter is optional, it defaults
+   to all 0. Data is expected in component-major ordering.
 
 **Unit:** :math:`\text{[Mod]}^{-1}`
 
-===================  =========================  
-**Type:** double     **Length:** NCOMP
-===================  ========================= 
+===================  ===============================  
+**Type:** double     **Length:** same as GIEX_NU
+===================  =============================== 
 
 ``GIEX_NU_QUAD``
    Coefficient of quadratic dependence of characteristic charge on
-   modifier component
+   modifier component. If more than NCOMP values are provided,
+   :math:`\nu` is given by a piecewise cubic polynomial. This parameter
+   is optional, it defaults to all 0. Data is expected in component-major
+   ordering.
 
 **Unit:** :math:`\text{[Mod]}^{-2}`
 
-===================  =========================  
-**Type:** double     **Length:** NCOMP
-===================  ========================= 
+===================  ===============================  
+**Type:** double     **Length:** same as GIEX_NU
+===================  =============================== 
+
+``GIEX_NU_CUBE``
+   Coefficient of cubic dependence of characteristic charge on
+   modifier component. If more than NCOMP values are provided,
+   :math:`\nu` is given by a piecewise cubic polynomial. This
+   parameter is optional, it defaults to all 0. Data is expected
+   in component-major ordering.
+
+**Unit:** :math:`\text{[Mod]}^{-3}`
+
+===================  ===============================  
+**Type:** double     **Length:** same as GIEX_NU
+===================  =============================== 
+
+``GIEX_NU_BREAKS``
+   Optional, only required if a piecewise cubic polynomial is
+   used for :math:`\nu`. Contains the breaks of the pieces
+   in component-major ordering. Note that :math:`n` pieces
+   have :math:`n+1` breaks.
+
+**Unit:** :math:`\text{[Mod]}`
+
+===================  ======================================  
+**Type:** double     **Length:** length of GIEX_NU + NCOMP
+===================  ====================================== 
 
 ``GIEX_SIGMA``
    Steric factors of the protein; The number of sites :math:`\sigma` on
diff --git a/src/libcadet/Spline.hpp b/src/libcadet/Spline.hpp
new file mode 100644
index 000000000..469e3cb3a
--- /dev/null
+++ b/src/libcadet/Spline.hpp
@@ -0,0 +1,175 @@
+// =============================================================================
+//  CADET
+//  
+//  Copyright © 2008-2021: The CADET Authors
+//            Please see the AUTHORS and CONTRIBUTORS file.
+//  
+//  All rights reserved. This program and the accompanying materials
+//  are made available under the terms of the GNU Public License v3.0 (or, at
+//  your option, any later version) which accompanies this distribution, and
+//  is available at http://www.gnu.org/licenses/gpl.html
+// =============================================================================
+
+#ifndef LIBCADET_SPLINE_HPP_
+#define LIBCADET_SPLINE_HPP_
+
+#include <algorithm>
+#include <tuple>
+
+#include "cadet/cadetCompilerInfo.hpp"
+#include "AutoDiff.hpp"
+
+namespace cadet 
+{
+
+	/**
+	 * @brief      Evaluates a piecewise cubic polynomial
+	 * @details    If the evaluation point @p x is outside the domain of the piecewise polynomial,
+	 *             constant extrapolation from the first and last value of the polynomial on its
+	 *             domain is applied.
+	 * @param[in]  x           Evaluation position
+	 * @param[in]  breaks      Array with @c nPieces+1 elements (strictly increasing) that defines the domain pieces
+	 * @param[in]  constCoeff  Array with constant coefficients of size @p nPieces
+	 * @param[in]  linCoeff    Array with linear coefficients of size @p nPieces
+	 * @param[in]  quadCoeff   Array with quadratic coefficients of size @p nPieces
+	 * @param[in]  cubeCoeff   Array with cubic coefficients of size @p nPieces
+	 * @param[in]  nPieces     Number of pieces (at least 1)
+	 * @tparam     value_t     Type of the evaluation point
+	 * @tparam     result_t    Type of the result
+	 * @return     Value of the piecewise cubic polynomial at the given point @p x
+	 */
+	template <typename value_t, typename result_t>
+	result_t evaluateCubicPiecewisePolynomial(const value_t& x, active const* breaks, active const* constCoeff,
+		active const* linCoeff, active const* quadCoeff, active const* cubeCoeff, int nPieces) CADET_NOEXCEPT
+	{
+		// Test if outside of domain, apply constant extrapolation
+		if (x < breaks[0])
+		{
+			return constCoeff[0];
+		}
+
+		if (x >= breaks[nPieces])
+		{
+			// Return the value at the right of the domain
+			const result_t y = static_cast<result_t>(breaks[nPieces]) - static_cast<result_t>(breaks[nPieces - 1]);
+			const int idx = nPieces - 1;
+			return static_cast<result_t>(constCoeff[idx]) + y * (
+				static_cast<result_t>(linCoeff[idx]) + y * (
+					static_cast<result_t>(quadCoeff[idx]) + y * static_cast<result_t>(cubeCoeff[idx]) 
+					)
+				);
+		}
+
+		// Find correct piece
+		int idx = 0;
+		for (; idx < nPieces; ++idx)
+		{
+			if (breaks[idx] >= x)
+				break;
+		}
+		
+		// Evaluate polynomial
+		const typename DoubleActivePromoter<value_t, result_t>::type y = x - static_cast<result_t>(breaks[idx]);
+		return static_cast<result_t>(constCoeff[idx]) + y * (
+			static_cast<result_t>(linCoeff[idx]) + y * (
+				static_cast<result_t>(quadCoeff[idx]) + y * static_cast<result_t>(cubeCoeff[idx]) 
+				)
+			);
+	}
+
+	/**
+	 * @brief      Evaluates a piecewise cubic polynomial and returns base and dependent values
+	 * @details    If the evaluation point @p x is outside the domain of the piecewise polynomial,
+	 *             constant extrapolation from the first and last value of the polynomial on its
+	 *             domain is applied.
+	 *             
+	 *             The function returns the constant base value and the variable part depending on @p x.
+	 * @param[in]  x           Evaluation position
+	 * @param[in]  breaks      Array with @c nPieces+1 elements (strictly increasing) that defines the domain pieces
+	 * @param[in]  constCoeff  Array with constant coefficients of size @p nPieces
+	 * @param[in]  linCoeff    Array with linear coefficients of size @p nPieces
+	 * @param[in]  quadCoeff   Array with quadratic coefficients of size @p nPieces
+	 * @param[in]  cubeCoeff   Array with cubic coefficients of size @p nPieces
+	 * @param[in]  nPieces     Number of pieces (at least 1)
+	 * @tparam     value_t     Type of the evaluation point
+	 * @tparam     result_t    Type of the result
+	 * @return     Constant and dynamic value of the piecewise cubic polynomial at the given point @p x
+	 */
+	template <typename value_t, typename result_t>
+	std::tuple<result_t, result_t> evaluateCubicPiecewisePolynomialSplit(const value_t& x, active const* breaks, active const* constCoeff,
+		active const* linCoeff, active const* quadCoeff, active const* cubeCoeff, int nPieces) CADET_NOEXCEPT
+	{
+		// Test if outside of domain, apply constant extrapolation
+		if (x < breaks[0])
+		{
+			return {static_cast<result_t>(constCoeff[0]), result_t(0.0)};
+		}
+
+		if (x >= breaks[nPieces])
+		{
+			// Return the value at the right of the domain
+			const result_t y = static_cast<result_t>(breaks[nPieces]) - static_cast<result_t>(breaks[nPieces - 1]);
+			const int idx = nPieces - 1;
+			return {static_cast<result_t>(constCoeff[idx]) + y * (
+				static_cast<result_t>(linCoeff[idx]) + y * (
+					static_cast<result_t>(quadCoeff[idx]) + y * static_cast<result_t>(cubeCoeff[idx]) 
+					)
+				), result_t(0.0)};
+		}
+
+		// Find correct piece
+		int idx = 0;
+		for (; idx < nPieces; ++idx)
+		{
+			if (breaks[idx] >= x)
+				break;
+		}
+		
+		// Evaluate polynomial
+		const typename DoubleActivePromoter<value_t, result_t>::type y = x - static_cast<result_t>(breaks[idx]);
+		return {static_cast<result_t>(constCoeff[idx]), y * (
+			static_cast<result_t>(linCoeff[idx]) + y * (
+				static_cast<result_t>(quadCoeff[idx]) + y * static_cast<result_t>(cubeCoeff[idx]) 
+				)
+			)};
+	}
+
+	/**
+	 * @brief      Evaluates the derivative of a piecewise cubic polynomial
+	 * @details    If the evaluation point @p x is outside the domain of the piecewise polynomial,
+	 *             constant extrapolation from the first and last value of the polynomial on its
+	 *             domain is applied.
+	 * @param[in]  x           Evaluation position
+	 * @param[in]  breaks      Array with @c nPieces+1 elements (strictly increasing) that defines the domain pieces
+	 * @param[in]  linCoeff    Array with linear coefficients of size @p nPieces
+	 * @param[in]  quadCoeff   Array with quadratic coefficients of size @p nPieces
+	 * @param[in]  cubeCoeff   Array with cubic coefficients of size @p nPieces
+	 * @param[in]  nPieces     Number of pieces (at least 1)
+	 * @return     Derivative of the piecewise cubic polynomial at the given point @p x
+	 */
+	inline double evaluateCubicPiecewisePolynomialDerivative(double x, active const* breaks, active const* linCoeff,
+		active const* quadCoeff, active const* cubeCoeff, int nPieces) CADET_NOEXCEPT
+	{
+		// Test if outside of domain, apply constant extrapolation
+		if (x < breaks[0])
+			return 0.0;
+
+		if (x >= breaks[nPieces])
+			return 0.0;
+
+		// Find correct piece
+		int idx = 0;
+		for (; idx < nPieces; ++idx)
+		{
+			if (breaks[idx] >= x)
+				break;
+		}
+		
+		// Evaluate polynomial
+		const double y = x - static_cast<double>(breaks[idx]);
+		return static_cast<double>(linCoeff[idx]) + y * (2.0 * static_cast<double>(quadCoeff[idx]) + 3.0 * y * static_cast<double>(cubeCoeff[idx]));
+	}
+
+} // namespace cadet
+
+#endif  // LIBCADET_SPLINE_HPP_
diff --git a/src/libcadet/model/Parameters.hpp b/src/libcadet/model/Parameters.hpp
index 808075c5c..ebc69d466 100644
--- a/src/libcadet/model/Parameters.hpp
+++ b/src/libcadet/model/Parameters.hpp
@@ -186,6 +186,42 @@ class ScalarParameter
 };
 
 
+/**
+ * @brief Vector parameter
+ * @details A single value per component.
+ */
+class VectorParameter
+{
+public:
+
+	typedef std::vector<active> storage_t;
+
+	VectorParameter(std::vector<active>& p) : _p(&p) { }
+	VectorParameter(std::vector<active>* p) : _p(p) { }
+
+	inline void configure(const std::string& varName, IParameterProvider& paramProvider, unsigned int nComp, unsigned int const* nBoundStates)
+	{
+		readScalarParameterOrArray(*_p, paramProvider, varName, 1);
+	}
+
+	inline void registerParam(const std::string& varName, std::unordered_map<ParameterId, active*>& parameters, UnitOpIdx unitOpIdx, ParticleTypeIdx parTypeIdx, unsigned int nComp, unsigned int const* nBoundStates)
+	{
+		const StringHash hashVar = hashStringRuntime(varName);
+		registerParam1DArray(parameters, *_p, [=](bool multi, unsigned int rid) { return makeParamId(hashVar, unitOpIdx, CompIndep, parTypeIdx, BoundStateIndep, rid, SectionIndep); });
+	}
+
+	inline void reserve(unsigned int numElem, unsigned int numSlices, unsigned int nComp, unsigned int const* nBoundStates) { }
+
+	inline std::size_t size() const CADET_NOEXCEPT { return _p->size(); }
+
+	inline const std::vector<active>& get() const CADET_NOEXCEPT { return *_p; }
+	inline std::vector<active>& get() CADET_NOEXCEPT { return *_p; }
+
+protected:
+	std::vector<active>* _p;
+};
+
+
 /**
  * @brief Component dependent scalar external parameter
  * @details A single value per component.
@@ -304,6 +340,42 @@ class ScalarReactionDependentParameter
 };
 
 
+/**
+ * @brief Component dependent vectorial external parameter
+ * @details A vector per component, the component index changes the slowest.
+ */
+class VectorComponentDependentParameter
+{
+public:
+
+	typedef std::vector<active> storage_t;
+
+	VectorComponentDependentParameter(std::vector<active>& p) : _p(&p) { }
+	VectorComponentDependentParameter(std::vector<active>* p) : _p(p) { }
+
+	inline void configure(const std::string& varName, IParameterProvider& paramProvider, unsigned int nComp, unsigned int const* nBoundStates)
+	{
+		readScalarParameterOrArray(*_p, paramProvider, varName, 1);
+	}
+
+	inline void registerParam(const std::string& varName, std::unordered_map<ParameterId, active*>& parameters, UnitOpIdx unitOpIdx, ParticleTypeIdx parTypeIdx, unsigned int nComp, unsigned int const* nBoundStates)
+	{
+		const StringHash nameHash = hashStringRuntime(varName);
+		registerParam2DArray(parameters, *_p, [=](bool multi, unsigned int comp, unsigned int rid) { return makeParamId(nameHash, unitOpIdx, comp, parTypeIdx, BoundStateIndep, multi ? rid : ReactionIndep, SectionIndep); }, _p->size() / nComp);
+	}
+
+	inline void reserve(unsigned int numElem, unsigned int numSlices, unsigned int nComp, unsigned int const* nBoundStates) { }
+
+	inline std::size_t size() const CADET_NOEXCEPT { return _p->size(); }
+
+	inline const std::vector<active>& get() const CADET_NOEXCEPT { return *_p; }
+	inline std::vector<active>& get() CADET_NOEXCEPT { return *_p; }
+
+protected:
+	std::vector<active>* _p;
+};
+
+
 /**
  * @brief Component and bound state dependent external parameter
  * @details A single value per component and bound state / binding site type.
@@ -749,6 +821,174 @@ class ExternalScalarParameter
 };
 
 
+/**
+ * @brief Vector external parameter
+ * @details A vector.
+ */
+class ExternalVectorParameter
+{
+public:
+
+	/**
+	 * @brief Underlying type
+	 */
+	typedef std::vector<active> storage_t;
+
+	/**
+	 * @brief Reads parameters and verifies them
+	 * @details See IBindingModel::configure() for details.
+	 * @param [in] varName Name of the parameter
+	 * @param [in] paramProvider IParameterProvider used for reading parameters
+	 * @param [in] nComp Number of components
+	 * @param [in] nBoundStates Array with number of bound states for each component
+	 * @return @c true if the parameters were read and validated successfully, otherwise @c false
+	 */
+	inline void configure(const std::string& varName, IParameterProvider& paramProvider, unsigned int nComp, unsigned int const* nBoundStates)
+	{
+		readScalarParameterOrArray(_base, paramProvider, "EXT_" + varName, 1);
+		readScalarParameterOrArray(_linear, paramProvider, "EXT_" + varName + "_T", 1);
+		readScalarParameterOrArray(_quad, paramProvider, "EXT_" + varName + "_TT", 1);
+		readScalarParameterOrArray(_cube, paramProvider, "EXT_" + varName + "_TTT", 1);
+	}
+
+	/**
+	 * @brief Registers the parameters in a map for further use
+	 * @param [in] varName Name of the parameter
+	 * @param [in,out] parameters Map in which the parameters are stored
+	 * @param [in] unitOpIdx Index of the unit operation used for registering the parameters
+	 * @param [in] nComp Number of components
+	 * @param [in] nBoundStates Array with number of bound states for each component
+	 */
+	inline void registerParam(const std::string& varName, std::unordered_map<ParameterId, active*>& parameters, UnitOpIdx unitOpIdx, ParticleTypeIdx parTypeIdx, unsigned int nComp, unsigned int const* nBoundStates)
+	{
+		const StringHash hashConst = hashStringRuntime("EXT_" + varName);
+		registerParam1DArray(parameters, _base, [=](bool multi, unsigned int rid) { return makeParamId(hashConst, unitOpIdx, CompIndep, parTypeIdx, BoundStateIndep, rid, SectionIndep); });
+
+		const StringHash hashLinear = hashStringRuntime("EXT_" + varName + "_T");
+		registerParam1DArray(parameters, _linear, [=](bool multi, unsigned int rid) { return makeParamId(hashLinear, unitOpIdx, CompIndep, parTypeIdx, BoundStateIndep, rid, SectionIndep); });
+
+		const StringHash hashQuad = hashStringRuntime("EXT_" + varName + "_TT");
+		registerParam1DArray(parameters, _quad, [=](bool multi, unsigned int rid) { return makeParamId(hashQuad, unitOpIdx, CompIndep, parTypeIdx, BoundStateIndep, rid, SectionIndep); });
+
+		const StringHash hashCube = hashStringRuntime("EXT_" + varName + "_TTT");
+		registerParam1DArray(parameters, _cube, [=](bool multi, unsigned int rid) { return makeParamId(hashCube, unitOpIdx, CompIndep, parTypeIdx, BoundStateIndep, rid, SectionIndep); });
+	}
+
+	/**
+	 * @brief Reserves space in the storage of the parameters
+	 * @param [in] numElem Total number of components in all slices / binding site types
+	 * @param [in] numSlices Number of slices / binding site types
+	 * @param [in] nComp Number of components
+	 * @param [in] nBoundStates Array with number of bound states for each component
+	 */
+	inline void reserve(unsigned int numElem, unsigned int numSlices, unsigned int nComp, unsigned int const* nBoundStates) { }
+
+	/**
+	 * @brief Calculates a parameter in order to take the external profile into account
+	 * @param [out] result Stores the result of the paramter
+	 * @param [in] extVal Value of the external function
+	 * @param [in] nComp Number of components
+	 * @param [in] nBoundStates Array with number of bound states for each component
+	 */
+	inline void update(std::vector<active>& result, double extVal, unsigned int nComp, unsigned int const* nBoundStates) const
+	{
+		update(result.data(), extVal, nComp, nBoundStates);
+	}
+
+	inline void update(active* result, double extVal, unsigned int nComp, unsigned int const* nBoundStates) const
+	{
+		for (std::size_t i = 0; i < _base.size(); ++i)
+			result[i] = _base[i] + extVal * (_linear[i] + extVal * (_quad[i] + extVal * _cube[i]));
+	}
+
+	/**
+	 * @brief Calculates time derivative of parameter in case of external dependence
+	 * @param [out] result Stores the result of the paramter
+	 * @param [in] extVal Value of the external function
+	 * @param [in] extTimeDiff Time derivative of the external function
+	 * @param [in] nComp Number of components
+	 * @param [in] nBoundStates Array with number of bound states for each component
+	 */
+	inline void updateTimeDerivative(active& result, double extVal, double extTimeDiff, unsigned int nComp, unsigned int const* nBoundStates) const
+	{
+		updateTimeDerivative(&result, extVal, extTimeDiff, nComp, nBoundStates);
+	}
+
+	/**
+	 * @brief Calculates time derivative of parameter in case of external dependence
+	 * @param [out] result Stores the result of the paramter
+	 * @param [in] extVal Value of the external function
+	 * @param [in] extTimeDiff Time derivative of the external function
+	 * @param [in] nComp Number of components
+	 * @param [in] nBoundStates Array with number of bound states for each component
+	 */
+	inline void updateTimeDerivative(std::vector<active>& result, double extVal, double extTimeDiff, unsigned int nComp, unsigned int const* nBoundStates) const
+	{
+		updateTimeDerivative(result.data(), extVal, extTimeDiff, nComp, nBoundStates);
+	}
+
+	inline void updateTimeDerivative(active* result, double extVal, double extTimeDiff, unsigned int nComp, unsigned int const* nBoundStates) const
+	{
+		for (std::size_t i = 0; i < _base.size(); ++i)
+			result[i] = extTimeDiff * (static_cast<double>(_linear[i]) + extVal * (2.0 * static_cast<double>(_quad[i]) + 3.0 * extVal * static_cast<double>(_cube[i])));
+	}
+
+	/**
+	 * @brief Returns the base value that does not depend on an external value
+	 * @return Constant base value
+	 */
+	inline storage_t& base() CADET_NOEXCEPT { return _base; }
+	inline const storage_t& base() const CADET_NOEXCEPT { return _base; }
+
+	/**
+	 * @brief Returns the amount of additional memory (usually dynamically allocated by containers) for storing the final parameters
+	 * @details In a model, externally dependent parameters are stored in a struct, usually called
+	 *          VariableParams. This is sufficient for "static" parameter types that do
+	 *          not require additional memory (which is usually allocated dynamically).
+	 *          For containers using additional dynamic memory, only the container itself
+	 *          is stored in the struct. Memory for the content of the container (i.e., the
+	 *          elements) is still required. This function computes this amount of additional
+	 *          memory.
+	 *
+	 * @param [in] nComp Number of components
+	 * @param [in] totalNumBoundStates Total number of bound states
+	 * @param [in] nBoundStates Array with number of bound states for each component
+	 * @return Amount of additional memory in bytes
+	 */
+	inline std::size_t additionalDynamicMemory(unsigned int nComp, unsigned int totalNumBoundStates, unsigned int const* nBoundStates) const CADET_NOEXCEPT { return 0; }
+
+	/**
+	 * @brief Prepares the cache for the updated values
+	 * @details The cache is a local version of storage_t (e.g., LocalVector).
+	 * @param [in,out] cache Cache object to be prepared
+	 * @param [in] ptr Pointer to cache buffer
+	 */
+	template <typename T>
+	inline void prepareCache(T& cache, LinearBufferAllocator& buffer) const
+	{
+		cache.fromTemplate(buffer, _base);
+	}
+
+	/**
+	 * @brief Returns the number of elements in the parameter
+	 * @return Number of elements in the parameter
+	 */
+	inline std::size_t size() const CADET_NOEXCEPT { return _base.size(); }
+
+	/**
+	 * @brief Returns whether base, linear, quadratic, and cubic arrays have the same size
+	 * @return @c true if the arrays are of the same size, otherwise @c false
+	 */
+	inline bool allSameSize() const CADET_NOEXCEPT { return (_base.size() == _linear.size()) && (_base.size() == _quad.size()) && (_base.size() == _cube.size()); }
+
+protected:
+	std::vector<active> _base;
+	std::vector<active> _linear;
+	std::vector<active> _quad;
+	std::vector<active> _cube;
+};
+
+
 /**
  * @brief Component dependent scalar external parameter
  * @details A single value per component.
@@ -998,6 +1238,95 @@ class ExternalScalarReactionDependentParameter
 };
 
 
+/**
+ * @brief Component dependent vectorial external parameter
+ * @details A vector per component, the component index changes the slowest.
+ */
+class ExternalVectorComponentDependentParameter
+{
+public:
+
+	typedef std::vector<active> storage_t;
+
+	inline void configure(const std::string& varName, IParameterProvider& paramProvider, unsigned int nComp, unsigned int const* nBoundStates)
+	{
+		readScalarParameterOrArray(_base, paramProvider, "EXT_" + varName, 1);
+		readScalarParameterOrArray(_linear, paramProvider, "EXT_" + varName + "_T", 1);
+		readScalarParameterOrArray(_quad, paramProvider, "EXT_" + varName + "_TT", 1);
+		readScalarParameterOrArray(_cube, paramProvider, "EXT_" + varName + "_TTT", 1);
+	}
+
+	inline void registerParam(const std::string& varName, std::unordered_map<ParameterId, active*>& parameters, UnitOpIdx unitOpIdx, ParticleTypeIdx parTypeIdx, unsigned int nComp, unsigned int const* nBoundStates)
+	{
+		const unsigned int stride = _base.size() / nComp;
+		const StringHash hashConst = hashStringRuntime("EXT_" + varName);
+		registerParam2DArray(parameters, _base, [=](bool multi, unsigned int comp, unsigned int rid) { return makeParamId(hashConst, unitOpIdx, comp, parTypeIdx, BoundStateIndep, multi ? rid : ReactionIndep, SectionIndep); }, stride);
+
+		const StringHash hashLinear = hashStringRuntime("EXT_" + varName + "_T");
+		registerParam2DArray(parameters, _linear, [=](bool multi, unsigned int comp, unsigned int rid) { return makeParamId(hashLinear, unitOpIdx, comp, parTypeIdx, BoundStateIndep, multi ? rid : ReactionIndep, SectionIndep); }, stride);
+
+		const StringHash hashQuad = hashStringRuntime("EXT_" + varName + "_TT");
+		registerParam2DArray(parameters, _quad, [=](bool multi, unsigned int comp, unsigned int rid) { return makeParamId(hashQuad, unitOpIdx, comp, parTypeIdx, BoundStateIndep, multi ? rid : ReactionIndep, SectionIndep); }, stride);
+
+		const StringHash hashCube = hashStringRuntime("EXT_" + varName + "_TTT");
+		registerParam2DArray(parameters, _cube, [=](bool multi, unsigned int comp, unsigned int rid) { return makeParamId(hashCube, unitOpIdx, comp, parTypeIdx, BoundStateIndep, multi ? rid : ReactionIndep, SectionIndep); }, stride);
+	}
+
+	inline void reserve(unsigned int numElem, unsigned int numSlices, unsigned int nComp, unsigned int const* nBoundStates) { }
+
+	inline void update(std::vector<active>& result, double extVal, unsigned int nComp, unsigned int const* nBoundStates) const
+	{
+		update(result.data(), extVal, nComp, nBoundStates);
+	}
+
+	inline void update(active* result, double extVal, unsigned int nComp, unsigned int const* nBoundStates) const
+	{
+		for (std::size_t i = 0; i < _base.size(); ++i)
+			result[i] = _base[i] + extVal * (_linear[i] + extVal * (_quad[i] + extVal * _cube[i]));
+	}
+
+	inline void updateTimeDerivative(std::vector<active>& result, double extVal, double extTimeDiff, unsigned int nComp, unsigned int const* nBoundStates) const
+	{
+		updateTimeDerivative(result.data(), extVal, extTimeDiff, nComp, nBoundStates);
+	}
+
+	inline void updateTimeDerivative(active* result, double extVal, double extTimeDiff, unsigned int nComp, unsigned int const* nBoundStates) const
+	{
+		for (std::size_t i = 0; i < _base.size(); ++i)
+			result[i] = extTimeDiff * (static_cast<double>(_linear[i]) + extVal * (2.0 * static_cast<double>(_quad[i]) + 3.0 * extVal * static_cast<double>(_cube[i])));
+	}
+
+	inline std::size_t size() const CADET_NOEXCEPT { return _base.size(); }
+	inline bool allSameSize() const CADET_NOEXCEPT { return (_base.size() == _linear.size()) && (_base.size() == _quad.size()) && (_base.size() == _cube.size()); }
+
+	inline std::size_t additionalDynamicMemory(unsigned int nComp, unsigned int totalNumBoundStates, unsigned int const* nBoundStates) const CADET_NOEXCEPT { return nComp * sizeof(active) + alignof(active); }
+
+	inline storage_t& base() CADET_NOEXCEPT { return _base; }
+	inline const storage_t& base() const CADET_NOEXCEPT { return _base; }
+
+	inline storage_t& linear() CADET_NOEXCEPT { return _linear; }
+	inline const storage_t& linear() const CADET_NOEXCEPT { return _linear; }
+
+	inline storage_t& quadratic() CADET_NOEXCEPT { return _quad; }
+	inline const storage_t& quadratic() const CADET_NOEXCEPT { return _quad; }
+
+	inline storage_t& cubic() CADET_NOEXCEPT { return _cube; }
+	inline const storage_t& cubic() const CADET_NOEXCEPT { return _cube; }
+
+	template <typename T>
+	inline void prepareCache(T& cache, LinearBufferAllocator& buffer) const
+	{
+		cache.fromTemplate(buffer, _base);
+	}
+
+protected:
+	std::vector<active> _base;
+	std::vector<active> _linear;
+	std::vector<active> _quad;
+	std::vector<active> _cube;
+};
+
+
 /**
  * @brief Component and bound state dependent external parameter
  * @details A single value per component and bound state / binding site type.
diff --git a/src/libcadet/model/binding/AntiLangmuirBinding.cpp b/src/libcadet/model/binding/AntiLangmuirBinding.cpp
index b9e3ce6e5..9bb5e2036 100644
--- a/src/libcadet/model/binding/AntiLangmuirBinding.cpp
+++ b/src/libcadet/model/binding/AntiLangmuirBinding.cpp
@@ -52,7 +52,7 @@ namespace model
 
 inline const char* AntiLangmuirParamHandler::identifier() CADET_NOEXCEPT { return "MULTI_COMPONENT_ANTILANGMUIR"; }
 
-inline bool AntiLangmuirParamHandler::validateConfig(unsigned int nComp, unsigned int const* nBoundStates)
+inline bool AntiLangmuirParamHandler::validate(unsigned int nComp, unsigned int const* nBoundStates)
 {
 	if ((_kA.size() != _kD.size()) || (_kA.size() != _qMax.size()) || (_kA.size() != _antiLangmuir.size()) || (_kA.size() < nComp))
 		throw InvalidParameterException("MCAL_KA, MCAL_KD, MCAL_QMAX, and MCAL_ANTILANGMUIR have to have the same size");
@@ -62,7 +62,7 @@ inline bool AntiLangmuirParamHandler::validateConfig(unsigned int nComp, unsigne
 
 inline const char* ExtAntiLangmuirParamHandler::identifier() CADET_NOEXCEPT { return "EXT_MULTI_COMPONENT_ANTILANGMUIR"; }
 
-inline bool ExtAntiLangmuirParamHandler::validateConfig(unsigned int nComp, unsigned int const* nBoundStates)
+inline bool ExtAntiLangmuirParamHandler::validate(unsigned int nComp, unsigned int const* nBoundStates)
 {
 	if ((_kA.size() != _kD.size()) || (_kA.size() != _qMax.size()) || (_kA.size() != _antiLangmuir.size()) || (_kA.size() < nComp))
 		throw InvalidParameterException("MCAL_KA, MCAL_KD, MCAL_QMAX, and MCAL_ANTILANGMUIR have to have the same size");
diff --git a/src/libcadet/model/binding/BiLangmuirBinding.cpp b/src/libcadet/model/binding/BiLangmuirBinding.cpp
index 61cfe4bdb..e7a1ce260 100644
--- a/src/libcadet/model/binding/BiLangmuirBinding.cpp
+++ b/src/libcadet/model/binding/BiLangmuirBinding.cpp
@@ -51,7 +51,7 @@ namespace model
 
 inline const char* BiLangmuirParamHandler::identifier() CADET_NOEXCEPT { return "MULTI_COMPONENT_BILANGMUIR"; }
 
-inline bool BiLangmuirParamHandler::validateConfig(unsigned int nComp, unsigned int const* nBoundStates)
+inline bool BiLangmuirParamHandler::validate(unsigned int nComp, unsigned int const* nBoundStates)
 {
 	if ((_kA.size() != _kD.size()) || (_kA.size() != _qMax.size()) || (_kA.size() < nComp))
 		throw InvalidParameterException("MCBL_KA, MCBL_KD, and MCBL_QMAX have to have the same size");
@@ -61,7 +61,7 @@ inline bool BiLangmuirParamHandler::validateConfig(unsigned int nComp, unsigned
 
 inline const char* ExtBiLangmuirParamHandler::identifier() CADET_NOEXCEPT { return "EXT_MULTI_COMPONENT_BILANGMUIR"; }
 
-inline bool ExtBiLangmuirParamHandler::validateConfig(unsigned int nComp, unsigned int const* nBoundStates)
+inline bool ExtBiLangmuirParamHandler::validate(unsigned int nComp, unsigned int const* nBoundStates)
 {
 	if ((_kA.size() != _kD.size()) || (_kA.size() != _qMax.size()) || (_kA.size() < nComp))
 		throw InvalidParameterException("EXT_MCBL_KA, EXT_MCBL_KD, and EXT_MCBL_QMAX have to have the same size");
diff --git a/src/libcadet/model/binding/BiStericMassActionBinding.cpp b/src/libcadet/model/binding/BiStericMassActionBinding.cpp
index 209b7b663..7ac2c975b 100644
--- a/src/libcadet/model/binding/BiStericMassActionBinding.cpp
+++ b/src/libcadet/model/binding/BiStericMassActionBinding.cpp
@@ -62,7 +62,7 @@ namespace model
 
 inline const char* BiSMAParamHandler::identifier() CADET_NOEXCEPT { return "BI_STERIC_MASS_ACTION"; }
 
-inline bool BiSMAParamHandler::validateConfig(unsigned int nComp, unsigned int const* nBoundStates)
+inline bool BiSMAParamHandler::validate(unsigned int nComp, unsigned int const* nBoundStates)
 {
 	const unsigned int numStates = firstNonEmptyBoundStates(nBoundStates, nComp);
 
@@ -83,7 +83,7 @@ inline bool BiSMAParamHandler::validateConfig(unsigned int nComp, unsigned int c
 
 inline const char* ExtBiSMAParamHandler::identifier() CADET_NOEXCEPT { return "EXT_BI_STERIC_MASS_ACTION"; }
 
-inline bool ExtBiSMAParamHandler::validateConfig(unsigned int nComp, unsigned int const* nBoundStates)
+inline bool ExtBiSMAParamHandler::validate(unsigned int nComp, unsigned int const* nBoundStates)
 {
 	const unsigned int numStates = firstNonEmptyBoundStates(nBoundStates, nComp);
 
diff --git a/src/libcadet/model/binding/BindingModelBase.hpp b/src/libcadet/model/binding/BindingModelBase.hpp
index 85ef09f72..6c2faa996 100644
--- a/src/libcadet/model/binding/BindingModelBase.hpp
+++ b/src/libcadet/model/binding/BindingModelBase.hpp
@@ -133,12 +133,12 @@ class ParamHandlerBindingModelBase : public BindingModelBase
 	virtual bool configureImpl(IParameterProvider& paramProvider, UnitOpIdx unitOpIdx, ParticleTypeIdx parTypeIdx)
 	{
 		// Read parameters
-		_paramHandler.configure(paramProvider, _nComp, _nBoundStates);
+		const bool valid = _paramHandler.configureAndValidate(paramProvider, _nComp, _nBoundStates);
 
 		// Register parameters
 		_paramHandler.registerParameters(_parameters, unitOpIdx, parTypeIdx, _nComp, _nBoundStates);
 
-		return true;
+		return valid;
 	}
 };
 
diff --git a/src/libcadet/model/binding/ColloidalBinding.cpp b/src/libcadet/model/binding/ColloidalBinding.cpp
index 8089050bb..04c582ba4 100644
--- a/src/libcadet/model/binding/ColloidalBinding.cpp
+++ b/src/libcadet/model/binding/ColloidalBinding.cpp
@@ -71,7 +71,7 @@ namespace model
 
 inline const char* ColloidalParamHandler::identifier() CADET_NOEXCEPT { return "MULTI_COMPONENT_COLLOIDAL"; }
 
-inline bool ColloidalParamHandler::validateConfig(unsigned int nComp, unsigned int const* nBoundStates)
+inline bool ColloidalParamHandler::validate(unsigned int nComp, unsigned int const* nBoundStates)
 {
 	if ((_kEqPhExp.size() != _kEqSaltPowerExp.size()) 
 			|| (_kEqPhExp.size() != _kEqSaltPowerFact.size())
@@ -97,7 +97,7 @@ inline bool ColloidalParamHandler::validateConfig(unsigned int nComp, unsigned i
 
 inline const char* ExtColloidalParamHandler::identifier() CADET_NOEXCEPT { return "EXT_MULTI_COMPONENT_COLLOIDAL"; }
 
-inline bool ExtColloidalParamHandler::validateConfig(unsigned int nComp, unsigned int const* nBoundStates)
+inline bool ExtColloidalParamHandler::validate(unsigned int nComp, unsigned int const* nBoundStates)
 {
 	if ((_kEqPhExp.size() != _kEqSaltPowerExp.size()) 
 			|| (_kEqPhExp.size() != _kEqSaltPowerFact.size())
diff --git a/src/libcadet/model/binding/ExtendedMobilePhaseModulatorLangmuirBinding.cpp b/src/libcadet/model/binding/ExtendedMobilePhaseModulatorLangmuirBinding.cpp
index 4d573a01e..14290e643 100644
--- a/src/libcadet/model/binding/ExtendedMobilePhaseModulatorLangmuirBinding.cpp
+++ b/src/libcadet/model/binding/ExtendedMobilePhaseModulatorLangmuirBinding.cpp
@@ -55,7 +55,7 @@ namespace model
 
 inline const char* EMPMLangmuirParamHandler::identifier() CADET_NOEXCEPT { return "EXTENDED_MOBILE_PHASE_MODULATOR"; }
 
-inline bool EMPMLangmuirParamHandler::validateConfig(unsigned int nComp, unsigned int const* nBoundStates)
+inline bool EMPMLangmuirParamHandler::validate(unsigned int nComp, unsigned int const* nBoundStates)
 {
 	if ((_kA.size() != _kD.size()) || (_kA.size() != _qMax.size()) || (_kA.size() != _gamma.size())
 		|| (_kA.size() != _beta.size()) || (_kA.size() < nComp))
@@ -66,7 +66,7 @@ inline bool EMPMLangmuirParamHandler::validateConfig(unsigned int nComp, unsigne
 
 inline const char* ExtEMPMLangmuirParamHandler::identifier() CADET_NOEXCEPT { return "EXT_EXTENDED_MOBILE_PHASE_MODULATOR"; }
 
-inline bool ExtEMPMLangmuirParamHandler::validateConfig(unsigned int nComp, unsigned int const* nBoundStates)
+inline bool ExtEMPMLangmuirParamHandler::validate(unsigned int nComp, unsigned int const* nBoundStates)
 {
 	if ((_kA.size() != _kD.size()) || (_kA.size() != _qMax.size()) || (_kA.size() != _gamma.size())
 		|| (_kA.size() != _beta.size()) || (_kA.size() < nComp))
diff --git a/src/libcadet/model/binding/ExternalFunctionTemplate.cpp b/src/libcadet/model/binding/ExternalFunctionTemplate.cpp
index 7fc91b04f..688b57872 100644
--- a/src/libcadet/model/binding/ExternalFunctionTemplate.cpp
+++ b/src/libcadet/model/binding/ExternalFunctionTemplate.cpp
@@ -83,7 +83,7 @@ class {{ name }} : public cadet::model::ConstParamHandlerBase
 {% endif %}
 	{ }
 
-	inline bool configure(IParameterProvider& paramProvider, unsigned int nComp, unsigned int const* nBoundStates)
+	inline void configure(IParameterProvider& paramProvider, unsigned int nComp, unsigned int const* nBoundStates)
 	{
 {% for p in parameters %}
 	{% if not existsIn(p, "skipConfig") %}
@@ -109,7 +109,12 @@ class {{ name }} : public cadet::model::ConstParamHandlerBase
 		{% endif %}
 	{% endfor %}
 {% endif %}
-		return validateConfig(nComp, nBoundStates);
+	}
+
+	inline bool configureAndValidate(IParameterProvider& paramProvider, unsigned int nComp, unsigned int const* nBoundStates)
+	{
+		configure(paramProvider, nComp, nBoundStates);
+		return validate(nComp, nBoundStates);
 	}
 
 	inline void registerParameters(std::unordered_map<ParameterId, active*>& parameters, UnitOpIdx unitOpIdx, ParticleTypeIdx parTypeIdx, unsigned int nComp, unsigned int const* nBoundStates)
@@ -171,9 +176,9 @@ class {{ name }} : public cadet::model::ConstParamHandlerBase
 {% endif %}
 
 	inline char const* prefixInConfiguration() const CADET_NOEXCEPT { return ""; }
+	inline bool validate(unsigned int nComp, unsigned int const* nBoundStates);
 
 protected:
-	inline bool validateConfig(unsigned int nComp, unsigned int const* nBoundStates);
 	
 	ConstParams _localParams;
 
@@ -262,7 +267,7 @@ class {{ externalName }} : public cadet::model::ExternalParamHandlerBase
 {% endif %}
 	{ }
 
-	inline bool configure(IParameterProvider& paramProvider, unsigned int nComp, unsigned int const* nBoundStates)
+	inline void configure(IParameterProvider& paramProvider, unsigned int nComp, unsigned int const* nBoundStates)
 	{
 {% for p in parameters %}
 	{% if not existsIn(p, "skipConfig") %}
@@ -289,7 +294,12 @@ class {{ externalName }} : public cadet::model::ExternalParamHandlerBase
 	{% endfor %}
 {% endif %}
 		ExternalParamHandlerBase::configure(paramProvider, {{ length(parameters) }});
-		return validateConfig(nComp, nBoundStates);
+	}
+
+	inline bool configureAndValidate(IParameterProvider& paramProvider, unsigned int nComp, unsigned int const* nBoundStates)
+	{
+		configure(paramProvider, nComp, nBoundStates);
+		return validate(nComp, nBoundStates);
 	}
 
 	inline void registerParameters(std::unordered_map<ParameterId, active*>& parameters, UnitOpIdx unitOpIdx, ParticleTypeIdx parTypeIdx, unsigned int nComp, unsigned int const* nBoundStates)
@@ -409,9 +419,9 @@ class {{ externalName }} : public cadet::model::ExternalParamHandlerBase
 {% endif %}
 
 	inline char const* prefixInConfiguration() const CADET_NOEXCEPT { return "EXT_"; }
+	inline bool validate(unsigned int nComp, unsigned int const* nBoundStates);
 
 protected:
-	inline bool validateConfig(unsigned int nComp, unsigned int const* nBoundStates);
 	
 	ConstParams _constParams;
 
diff --git a/src/libcadet/model/binding/GeneralizedIonExchangeBinding.cpp b/src/libcadet/model/binding/GeneralizedIonExchangeBinding.cpp
index e98598bbf..3daec43da 100644
--- a/src/libcadet/model/binding/GeneralizedIonExchangeBinding.cpp
+++ b/src/libcadet/model/binding/GeneralizedIonExchangeBinding.cpp
@@ -19,6 +19,7 @@
 #include "model/Parameters.hpp"
 #include "LocalVector.hpp"
 #include "SimulationTypes.hpp"
+#include "Spline.hpp"
 
 #include <functional>
 #include <unordered_map>
@@ -42,9 +43,11 @@
 			{ "type": "ScalarComponentDependentParameter", "varName": "kDQuad", "confName": "GIEX_KD_QUAD"},
 			{ "type": "ScalarComponentDependentParameter", "varName": "kDSalt", "confName": "GIEX_KD_SALT"},
 			{ "type": "ScalarComponentDependentParameter", "varName": "kDProt", "confName": "GIEX_KD_PROT"},
-			{ "type": "ScalarComponentDependentParameter", "varName": "nu", "confName": "GIEX_NU"},
-			{ "type": "ScalarComponentDependentParameter", "varName": "nuLin", "confName": "GIEX_NU_LIN"},
-			{ "type": "ScalarComponentDependentParameter", "varName": "nuQuad", "confName": "GIEX_NU_QUAD"},
+			{ "type": "VectorComponentDependentParameter", "varName": "nu", "confName": "GIEX_NU"},
+			{ "type": "VectorComponentDependentParameter", "varName": "nuLin", "confName": "GIEX_NU_LIN", "skipConfig": true},
+			{ "type": "VectorComponentDependentParameter", "varName": "nuQuad", "confName": "GIEX_NU_QUAD", "skipConfig": true},
+			{ "type": "VectorComponentDependentParameter", "varName": "nuCube", "confName": "GIEX_NU_CUBE", "skipConfig": true},
+			{ "type": "VectorComponentDependentParameter", "varName": "nuBreaks", "confName": "GIEX_NU_BREAKS", "skipConfig": true},
 			{ "type": "ScalarComponentDependentParameter", "varName": "sigma", "confName": "GIEX_SIGMA"}
 		],
 	"constantParameters":
@@ -66,6 +69,82 @@
  refPhC0,refPhQ = Reference concentrations for pH dependent powers
 */
 
+namespace
+{
+	void assignZeros(cadet::model::VectorComponentDependentParameter& p, unsigned int size)
+	{
+		p.get() = std::vector<cadet::active>(size, 0.0);
+	}
+
+	void assignZeros(cadet::model::ExternalVectorComponentDependentParameter& p, unsigned int size)
+	{
+		p.base() = std::vector<cadet::active>(size, 0.0);
+		p.linear() = std::vector<cadet::active>(size, 0.0);
+		p.quadratic() = std::vector<cadet::active>(size, 0.0);
+		p.cubic() = std::vector<cadet::active>(size, 0.0);
+	}
+
+	void assignSinglePiece(cadet::model::VectorComponentDependentParameter& p)
+	{
+		p.get().clear();
+	}
+
+	void assignSinglePiece(cadet::model::ExternalVectorComponentDependentParameter& p)
+	{
+		p.base().clear();
+		p.linear().clear();
+		p.quadratic().clear();
+		p.cubic().clear();
+	}
+
+	template <typename param_t, typename params_t, typename ph_t, typename cp_state_t, typename q_state_t>
+	std::tuple<cp_state_t, q_state_t> cpQNuPowers(int comp, int nPieces, const params_t& p, ph_t pH, cp_state_t cpBase, cp_state_t cpVar, q_state_t qBase, q_state_t qVar)
+	{
+		if (p.nuBreaks.size() == 0)
+		{
+			const cp_state_t nu_i_0_over_nu0 = static_cast<param_t>(p.nu[comp]) / static_cast<param_t>(p.nu[0]);
+			const cp_state_t nu_i_pH_over_nu0 = pH * (static_cast<param_t>(p.nuLin[comp]) + pH * (static_cast<param_t>(p.nuQuad[comp]) + pH * static_cast<param_t>(p.nuCube[comp]))) / static_cast<param_t>(p.nu[0]);
+			return {pow(cpBase, nu_i_0_over_nu0) * pow(cpVar, nu_i_pH_over_nu0), pow(qBase, nu_i_0_over_nu0) * pow(qVar, nu_i_pH_over_nu0)};
+		}
+		else
+		{
+			const int offset = comp * nPieces;
+			const auto [nuConst, nuVar] = cadet::evaluateCubicPiecewisePolynomialSplit<ph_t, typename cadet::DoubleActivePromoter<param_t, cp_state_t>::type>(pH, p.nuBreaks.data() + comp * (nPieces + 1), p.nu.data() + offset, p.nuLin.data() + offset, p.nuQuad.data() + offset, p.nuCube.data() + offset, nPieces);
+			const cp_state_t nu_i_0_over_nu0 = static_cast<param_t>(nuConst) / static_cast<param_t>(p.nu[0]);
+			const cp_state_t nu_i_pH_over_nu0 = nuVar / static_cast<param_t>(p.nu[0]);
+			return {pow(cpBase, nu_i_0_over_nu0) * pow(cpVar, nu_i_pH_over_nu0), pow(qBase, nu_i_0_over_nu0) * pow(qVar, nu_i_pH_over_nu0)};
+		}
+	}
+
+	template <typename params_t, typename ph_t, typename result_t>
+	std::tuple<result_t, result_t> evaluateNu(int comp, int nPieces, const params_t& p, ph_t pH)
+	{
+		if (p.nuBreaks.size() == 0)
+		{
+			return {static_cast<result_t>(p.nu[comp]), pH * (static_cast<result_t>(p.nuLin[comp]) + pH * (static_cast<result_t>(p.nuQuad[comp]) + pH * static_cast<result_t>(p.nuCube[comp])))};
+		}
+		else
+		{
+			const int offset = comp * nPieces;
+			return cadet::evaluateCubicPiecewisePolynomialSplit<ph_t, result_t>(pH, p.nuBreaks.data() + comp * (nPieces + 1), p.nu.data() + offset, p.nuLin.data() + offset, p.nuQuad.data() + offset, p.nuCube.data() + offset, nPieces);
+		}
+	}
+
+	template <typename params_t>
+	double nuDerivative(int comp, int nPieces, const params_t& p, double pH)
+	{
+		if (p.nuBreaks.size() == 0)
+		{
+			return static_cast<double>(p.nuLin[comp]) + pH * (2.0 * static_cast<double>(p.nuQuad[comp]) + 3.0 * pH * static_cast<double>(p.nuCube[comp]));
+		}
+		else
+		{
+			const int offset = comp * nPieces;
+			return cadet::evaluateCubicPiecewisePolynomialDerivative(pH, p.nuBreaks.data() + comp * (nPieces + 1), p.nuLin.data() + offset, p.nuQuad.data() + offset, p.nuCube.data() + offset, nPieces);
+		}
+	}
+}
+
 namespace cadet
 {
 
@@ -74,7 +153,7 @@ namespace model
 
 inline const char* GIEXParamHandler::identifier() CADET_NOEXCEPT { return "GENERALIZED_ION_EXCHANGE"; }
 
-inline bool GIEXParamHandler::validateConfig(unsigned int nComp, unsigned int const* nBoundStates)
+inline bool GIEXParamHandler::validate(unsigned int nComp, unsigned int const* nBoundStates)
 {
 	if (nComp <= 2)
 		throw InvalidParameterException("GENERALIZED_ION_EXCHANGE requires at least 3 components");
@@ -84,20 +163,58 @@ inline bool GIEXParamHandler::validateConfig(unsigned int nComp, unsigned int co
 	if (_kD.size() < nComp)
 		throw InvalidParameterException("GIEX_KD requires NCOMP entries");
 	if (_nu.size() < nComp)
-		throw InvalidParameterException("GIEX_NU requires NCOMP entries");
+		throw InvalidParameterException("GIEX_NU requires at least NCOMP entries");
 	if (_sigma.size() < nComp)
 		throw InvalidParameterException("GIEX_SIGMA requires NCOMP entries");
 
+	if (_nu.size() != _nuLin.size())
+		throw InvalidParameterException("GIEX_NU and GIEX_NU_LIN do not have the same length");
+
+	if (_nu.size() != _nuQuad.size())
+		throw InvalidParameterException("GIEX_NU and GIEX_NU_QUAD do not have the same length");
+
+	if (_nu.size() != _nuCube.size())
+		throw InvalidParameterException("GIEX_NU and GIEX_NU_CUBE do not have the same length");
+
+	if ((_nu.size() % nComp) != 0)
+		throw InvalidParameterException("Length of GIEX_NU must be a multiple of NCOMP");
+
+	if ((_nu.size() + nComp != _nuBreaks.size()) && (_nuBreaks.size() > 0))
+		throw InvalidParameterException("GIEX_NU_BREAKS must have one entry more than polynomial pieces in GIEX_NU");
+
+	if ((_nu.size() != nComp) && (_nuBreaks.size() == 0))
+		throw InvalidParameterException("GIEX_NU is expected to have NCOMP entries");
+
 	// Assume monovalent salt ions by default
-	if (_nu.get()[0] <= 0.0)
-		_nu.get()[0] = 1.0;
+	const int nPieces = _nu.size() / nComp;
+	for (int i = 0; i < nPieces; ++i)
+	{
+		if (_nu.get()[i] <= 0.0)
+			_nu.get()[i] = 1.0;
+	}
+
+	// Check breaks
+	if (_nuBreaks.size() > 1)
+	{
+		for (int i = 0; i < nComp; ++i)
+		{
+			cadet::active const* const b = _nuBreaks.get().data() + (nPieces + 1) * i;
+			for (int j = 0; j < nPieces; ++j)
+			{
+				if (b[j] >= b[j+1])
+				{
+					throw InvalidParameterException("GIEX_NU_BREAKS must be strictly increasing for each component");
+				}
+			}
+		}
+	}
 
 	return true;
 }
 
 inline const char* ExtGIEXParamHandler::identifier() CADET_NOEXCEPT { return "EXT_GENERALIZED_ION_EXCHANGE"; }
 
-inline bool ExtGIEXParamHandler::validateConfig(unsigned int nComp, unsigned int const* nBoundStates)
+inline bool ExtGIEXParamHandler::validate(unsigned int nComp, unsigned int const* nBoundStates)
 {
 	if (nComp <= 2)
 		throw InvalidParameterException("EXT_GENERALIZED_ION_EXCHANGE requires at least 3 components");
@@ -107,13 +224,46 @@ inline bool ExtGIEXParamHandler::validateConfig(unsigned int nComp, unsigned int
 	if (_kD.size() < nComp)
 		throw InvalidParameterException("EXT_GIEX_KD requires NCOMP entries");
 	if (_nu.size() < nComp)
-		throw InvalidParameterException("EXT_GIEX_NU requires NCOMP entries");
+		throw InvalidParameterException("EXT_GIEX_NU requires at least NCOMP entries");
 	if (_sigma.size() < nComp)
 		throw InvalidParameterException("EXT_GIEX_SIGMA requires NCOMP entries");
 
+	if (_nu.size() != _nuLin.size())
+		throw InvalidParameterException("EXT_GIEX_NU and EXT_GIEX_NU_LIN do not have the same length");
+
+	if (_nu.size() != _nuQuad.size())
+		throw InvalidParameterException("EXT_GIEX_NU and EXT_GIEX_NU_QUAD do not have the same length");
+
+	if (_nu.size() != _nuCube.size())
+		throw InvalidParameterException("EXT_GIEX_NU and EXT_GIEX_NU_CUBE do not have the same length");
+
+	if ((_nu.size() % nComp) != 0)
+		throw InvalidParameterException("Length of EXT_GIEX_NU must be a multiple of NCOMP");
+
+	if ((_nu.size() + nComp != _nuBreaks.size()) && (_nuBreaks.size() > 0))
+		throw InvalidParameterException("EXT_GIEX_NU_BREAKS must have one entry more than polynomial pieces in EXT_GIEX_NU");
+
+	if ((_nu.size() != nComp) && (_nuBreaks.size() == 0))
+		throw InvalidParameterException("EXT_GIEX_NU is expected to have NCOMP entries");
+
+	if (!_nu.allSameSize())
+		throw InvalidParameterException("EXT_GIEX_NU, EXT_GIEX_NU_T, EXT_GIEX_NU_TT, and EXT_GIEX_NU_TTT must have the same size");
+	if (!_nuLin.allSameSize())
+		throw InvalidParameterException("EXT_GIEX_NU_LIN, EXT_GIEX_NU_LIN_T, EXT_GIEX_NU_LIN_TT, and EXT_GIEX_NU_LIN_TTT must have the same size");
+	if (!_nuQuad.allSameSize())
+		throw InvalidParameterException("EXT_GIEX_NU_QUAD, EXT_GIEX_NU_QUAD_T, EXT_GIEX_NU_QUAD_TT, and EXT_GIEX_NU_QUAD_TTT must have the same size");
+	if (!_nuCube.allSameSize())
+		throw InvalidParameterException("EXT_GIEX_NU_CUBE, EXT_GIEX_NU_CUBE_T, EXT_GIEX_NU_CUBE_TT, and EXT_GIEX_NU_CUBE_TTT must have the same size");
+	if (!_nuBreaks.allSameSize())
+		throw InvalidParameterException("EXT_GIEX_NU_BREAKS, EXT_GIEX_NU_BREAKS_T, EXT_GIEX_NU_BREAKS_TT, and EXT_GIEX_NU_BREAKS_TTT must have the same size");
+
 	// Assume monovalent salt ions by default
-	if (_nu.base()[0] <= 0.0)
-		_nu.base()[0] = 1.0;
+	const int nPieces = _nu.size() / nComp;
+	for (int i = 0; i < nPieces; ++i)
+	{
+		if ((_nu.base()[i] <= 0.0) && (_nu.linear()[i] <= 0.0) && (_nu.quadratic()[i] <= 0.0) && (_nu.cubic()[i] <= 0.0))
+			_nu.base()[i] = 1.0;
+	}
 
 	return true;
 }
@@ -177,6 +327,7 @@ class GeneralizedIonExchangeBindingBase : public ParamHandlerBindingModelBase<Pa
 
 		// Pseudo component 1 is pH
 		const double pH = yCp[1];
+		const int nPieces = (p->nuBreaks.size() > 0) ? (p->nu.size() / _nComp) : 0;
 
 		// Salt equation: nu_0 * q_0 - Lambda + Sum[nu_j(pH) * q_j, j] == 0
 		//           <=>  q_0 == (Lambda - Sum[nu_j(pH) * q_j, j]) / nu_0
@@ -189,9 +340,8 @@ class GeneralizedIonExchangeBindingBase : public ParamHandlerBindingModelBase<Pa
 			if (_nBoundStates[j] == 0)
 				continue;
 
-			const double nu_j = static_cast<double>(p->nu[j]) + pH * (static_cast<double>(p->nuLin[j]) + pH * static_cast<double>(p->nuQuad[j]));
-
-			y[0] -= nu_j * y[bndIdx];
+			const auto [nuConst, nuVar] = evaluateNu<typename ParamHandler_t::params_t, double, double>(j, nPieces, *p, pH);
+			y[0] -= (nuConst + nuVar) * y[bndIdx];
 
 			// Next bound component
 			++bndIdx;
@@ -215,10 +365,32 @@ class GeneralizedIonExchangeBindingBase : public ParamHandlerBindingModelBase<Pa
 	using ParamHandlerBindingModelBase<ParamHandler_t>::_reactionQuasistationarity;
 	using ParamHandlerBindingModelBase<ParamHandler_t>::_nComp;
 	using ParamHandlerBindingModelBase<ParamHandler_t>::_nBoundStates;
+	using ParamHandlerBindingModelBase<ParamHandler_t>::_parameters;
 
 	virtual bool configureImpl(IParameterProvider& paramProvider, UnitOpIdx unitOpIdx, ParticleTypeIdx parTypeIdx)
 	{
-		const bool valid = ParamHandlerBindingModelBase<ParamHandler_t>::configureImpl(paramProvider, unitOpIdx, parTypeIdx);
+		// Read parameters
+		_paramHandler.configure(paramProvider, _nComp, _nBoundStates);
+
+		if (paramProvider.exists(std::string(_paramHandler.prefixInConfiguration()) + "GIEX_NU_LIN"))
+			_paramHandler.nuLin().configure("GIEX_NU_LIN", paramProvider, _nComp, _nBoundStates);
+		else
+			assignZeros(_paramHandler.nuLin(), _paramHandler.nu().size());
+
+		if (paramProvider.exists(std::string(_paramHandler.prefixInConfiguration()) + "GIEX_NU_QUAD"))
+			_paramHandler.nuQuad().configure("GIEX_NU_QUAD", paramProvider, _nComp, _nBoundStates);
+		else
+			assignZeros(_paramHandler.nuQuad(), _paramHandler.nu().size());
+
+		if (paramProvider.exists(std::string(_paramHandler.prefixInConfiguration()) + "GIEX_NU_CUBE"))
+			_paramHandler.nuCube().configure("GIEX_NU_CUBE", paramProvider, _nComp, _nBoundStates);
+		else
+			assignZeros(_paramHandler.nuCube(), _paramHandler.nu().size());
+
+		if (paramProvider.exists(std::string(_paramHandler.prefixInConfiguration()) + "GIEX_NU_BREAKS"))
+			_paramHandler.nuBreaks().configure("GIEX_NU_BREAKS", paramProvider, _nComp, _nBoundStates);
+		else
+			assignSinglePiece(_paramHandler.nuBreaks());
 
 		if (paramProvider.exists(std::string(_paramHandler.prefixInConfiguration()) + "GIEX_PHREFC0") && paramProvider.exists(std::string(_paramHandler.prefixInConfiguration()) + "GIEX_PHREFQ"))
 		{
@@ -232,7 +404,10 @@ class GeneralizedIonExchangeBindingBase : public ParamHandlerBindingModelBase<Pa
 			_paramHandler.refConcentrationPh().getQ() = _paramHandler.refConcentration().getQ();
 		}
 
-		return valid;
+		// Register parameters
+		_paramHandler.registerParameters(_parameters, unitOpIdx, parTypeIdx, _nComp, _nBoundStates);
+
+		return _paramHandler.validate(_nComp, _nBoundStates);
 	}
 
 	template <typename StateType, typename CpStateType, typename ResidualType, typename ParamType>
@@ -241,11 +416,14 @@ class GeneralizedIonExchangeBindingBase : public ParamHandlerBindingModelBase<Pa
 	{
 		using CpStateParamType = typename DoubleActivePromoter<CpStateType, ParamType>::type;
 		using StateParamType = typename DoubleActivePromoter<StateType, ParamType>::type;
+		using PhType = typename ActiveRefOrDouble<const CpStateType>::type;
 
 		typename ParamHandler_t::ParamsHandle const p = _paramHandler.update(t, secIdx, colPos, _nComp, _nBoundStates, workSpace);
 
+		const int nPieces = (p->nuBreaks.size() > 0) ? (p->nu.size() / _nComp) : 0;
+
 		// Pseudo component 1 is pH
-		const CpStateType pH = yCp[1];
+		const PhType pH = yCp[1];
 
 		// Salt flux: nu_0 * q_0 - Lambda + Sum[nu_j * q_j, j] == 0 
 		//       <=>  nu_0 * q_0 == Lambda - Sum[nu_j * q_j, j] 
@@ -260,9 +438,9 @@ class GeneralizedIonExchangeBindingBase : public ParamHandlerBindingModelBase<Pa
 			if (_nBoundStates[j] == 0)
 				continue;
 
-			const CpStateParamType nu_j = static_cast<ParamType>(p->nu[j]) + pH * (static_cast<ParamType>(p->nuLin[j]) + pH * static_cast<ParamType>(p->nuQuad[j]));
+			const auto [nuConst, nuVar] = evaluateNu<typename ParamHandler_t::params_t,PhType, CpStateParamType>(j, nPieces, *p, pH);
 
-			res[0] += nu_j * y[bndIdx];
+			res[0] += (nuConst + nuVar) * y[bndIdx];
 			q0_bar -= static_cast<ParamType>(p->sigma[j]) * y[bndIdx];
 
 			// Next bound component
@@ -289,11 +467,14 @@ class GeneralizedIonExchangeBindingBase : public ParamHandlerBindingModelBase<Pa
 //			const CpStateParamType nu_i_over_nu0 = (static_cast<ParamType>(p->nu[i]) + pH * (static_cast<ParamType>(p->nuLin[i]) + pH * static_cast<ParamType>(p->nuQuad[i]))) / static_cast<ParamType>(p->nu[0]);
 //			const CpStateParamType c0_pow_nu = pow(yCp0_divRef, nu_i_over_nu0);
 //			const StateParamType q0_bar_pow_nu = pow(q0_bar_divRef, nu_i_over_nu0);
-			const CpStateParamType nu_i_0_over_nu0 = static_cast<ParamType>(p->nu[i]) / static_cast<ParamType>(p->nu[0]);
-			const CpStateParamType nu_i_pH_over_nu0 = pH * (static_cast<ParamType>(p->nuLin[i]) + pH * static_cast<ParamType>(p->nuQuad[i])) / static_cast<ParamType>(p->nu[0]);
-			const CpStateParamType c0_pow_nu = pow(yCp0_divRef, nu_i_0_over_nu0) * pow(yCp0_ph_divRef, nu_i_pH_over_nu0);
-			const StateParamType q0_bar_pow_nu = pow(q0_bar_divRef, nu_i_0_over_nu0) * pow(q0_bar_ph_divRef, nu_i_pH_over_nu0);
+
+//			const CpStateParamType nu_i_0_over_nu0 = static_cast<ParamType>(p->nu[i]) / static_cast<ParamType>(p->nu[0]);
+//			const CpStateParamType nu_i_pH_over_nu0 = pH * (static_cast<ParamType>(p->nuLin[i]) + pH * static_cast<ParamType>(p->nuQuad[i])) / static_cast<ParamType>(p->nu[0]);
+//			const CpStateParamType c0_pow_nu = pow(yCp0_divRef, nu_i_0_over_nu0) * pow(yCp0_ph_divRef, nu_i_pH_over_nu0);
+//			const StateParamType q0_bar_pow_nu = pow(q0_bar_divRef, nu_i_0_over_nu0) * pow(q0_bar_ph_divRef, nu_i_pH_over_nu0);
 			
+			const auto [c0_pow_nu, q0_bar_pow_nu] = cpQNuPowers<ParamType, typename ParamHandler_t::params_t, PhType, CpStateParamType, StateParamType>(i, nPieces, *p, pH, yCp0_divRef, yCp0_ph_divRef, q0_bar_divRef, q0_bar_ph_divRef);
+
 			// k_{a,i}(c_p, q, \mathrm{pH}) = k_{a,i,0} \exp(k_{a,i,1} \mathrm{pH} + k_{a,i,2} \mathrm{pH}^2 + k_{a,i,\mathrm{salt}} c_{p,0} + k_{a,i,\mathrm{prot}} c_{p,i})
 			const CpStateParamType ka_i = static_cast<ParamType>(p->kA[i]) * 
 				exp(pH * (static_cast<ParamType>(p->kALin[i]) + pH * static_cast<ParamType>(p->kAQuad[i])) 
@@ -318,6 +499,7 @@ class GeneralizedIonExchangeBindingBase : public ParamHandlerBindingModelBase<Pa
 	void jacobianImpl(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, double const* yCp, int offsetCp, RowIterator jac, LinearBufferAllocator workSpace) const
 	{
 		typename ParamHandler_t::ParamsHandle const p = _paramHandler.update(t, secIdx, colPos, _nComp, _nBoundStates, workSpace);
+		const int nPieces = (p->nuBreaks.size() > 0) ? (p->nu.size() / _nComp) : 0;
 
 		// Pseudo component 1 is pH
 		const double pH = yCp[1];
@@ -337,10 +519,10 @@ class GeneralizedIonExchangeBindingBase : public ParamHandlerBindingModelBase<Pa
 			if (_nBoundStates[j] == 0)
 				continue;
 
-			const double nu_j = static_cast<double>(p->nu[j]) + pH * (static_cast<double>(p->nuLin[j]) + pH * static_cast<double>(p->nuQuad[j]));
+			const auto [nuConst, nuVar] = evaluateNu<typename ParamHandler_t::params_t, double, double>(j, nPieces, *p, pH);
 
-			jac[bndIdx] = nu_j;
-			jac[1 - offsetCp] += (static_cast<double>(p->nuLin[j]) + 2.0 * pH * static_cast<double>(p->nuQuad[j])) * y[bndIdx];
+			jac[bndIdx] = nuConst + nuVar;
+			jac[1 - offsetCp] += nuDerivative(j, nPieces, *p, pH) * y[bndIdx];
 
 			// Calculate \bar{q}_0 = nu_0 * q_0 - Sum[sigma_j * q_j, j]
 			q0_bar -= static_cast<double>(p->sigma[j]) * y[bndIdx];
@@ -370,12 +552,14 @@ class GeneralizedIonExchangeBindingBase : public ParamHandlerBindingModelBase<Pa
 			if (_nBoundStates[i] == 0)
 				continue;
 
+			const auto [nuConst, nuVar] = evaluateNu<typename ParamHandler_t::params_t, double, double>(i, nPieces, *p, pH);
+
 			const double ka = static_cast<double>(p->kA[i]);
 			const double kd = static_cast<double>(p->kD[i]);
-			const double nu_0 = static_cast<double>(p->nu[i]) / static_cast<double>(p->nu[0]);
-			const double nu_pH = pH * (static_cast<double>(p->nuLin[i]) + pH * static_cast<double>(p->nuQuad[i])) / static_cast<double>(p->nu[0]);
+			const double nu_0 = nuConst / static_cast<double>(p->nu[0]);
+			const double nu_pH = nuVar / static_cast<double>(p->nu[0]);
 			const double nu = nu_0 + nu_pH;
-			const double dNuDpH = (static_cast<double>(p->nuLin[i]) + 2.0 * pH * static_cast<double>(p->nuQuad[i])) / static_cast<double>(p->nu[0]);
+			const double dNuDpH = nuDerivative(i, nPieces, *p, pH) / static_cast<double>(p->nu[0]);
 
 			const double c0_pow_nu     = pow(yCp0_divRef, nu_0) * pow(yCp0_ph_divRef, nu_pH);
 			const double q0_bar_pow_nu = pow(q0_bar_divRef, nu_0) * pow(q0_bar_ph_divRef, nu_pH);
diff --git a/src/libcadet/model/binding/KumarLangmuirBinding.cpp b/src/libcadet/model/binding/KumarLangmuirBinding.cpp
index 93fefaa9c..6f6e4576d 100644
--- a/src/libcadet/model/binding/KumarLangmuirBinding.cpp
+++ b/src/libcadet/model/binding/KumarLangmuirBinding.cpp
@@ -58,7 +58,7 @@ namespace model
 
 inline const char* KumarLangmuirParamHandler::identifier() CADET_NOEXCEPT { return "KUMAR_MULTI_COMPONENT_LANGMUIR"; }
 
-inline bool KumarLangmuirParamHandler::validateConfig(unsigned int nComp, unsigned int const* nBoundStates)
+inline bool KumarLangmuirParamHandler::validate(unsigned int nComp, unsigned int const* nBoundStates)
 {
 	if ((_kA.size() != _kAct.size()) || (_kA.size() != _kD.size()) || (_kA.size() != _qMax.size()) || (_kA.size() != _nu.size()) || (_kA.size() < nComp))
 		throw InvalidParameterException("KMCL_KA, KMCL_KD, KMCL_KACT, KMCL_NU, and KMCL_QMAX have to have the same size");
@@ -68,7 +68,7 @@ inline bool KumarLangmuirParamHandler::validateConfig(unsigned int nComp, unsign
 
 inline const char* ExtKumarLangmuirParamHandler::identifier() CADET_NOEXCEPT { return "EXT_KUMAR_MULTI_COMPONENT_LANGMUIR"; }
 
-inline bool ExtKumarLangmuirParamHandler::validateConfig(unsigned int nComp, unsigned int const* nBoundStates)
+inline bool ExtKumarLangmuirParamHandler::validate(unsigned int nComp, unsigned int const* nBoundStates)
 {
 	if ((_kA.size() != _kAct.size()) || (_kA.size() != _kD.size()) || (_kA.size() != _qMax.size()) || (_kA.size() != _nu.size()) || (_kA.size() < nComp))
 		throw InvalidParameterException("KMCL_KA, KMCL_KD, KMCL_KACT, KMCL_NU, and KMCL_QMAX have to have the same size");
diff --git a/src/libcadet/model/binding/LangmuirBinding.cpp b/src/libcadet/model/binding/LangmuirBinding.cpp
index dcc475a13..650fefc90 100644
--- a/src/libcadet/model/binding/LangmuirBinding.cpp
+++ b/src/libcadet/model/binding/LangmuirBinding.cpp
@@ -51,7 +51,7 @@ namespace model
 
 inline const char* LangmuirParamHandler::identifier() CADET_NOEXCEPT { return "MULTI_COMPONENT_LANGMUIR"; }
 
-inline bool LangmuirParamHandler::validateConfig(unsigned int nComp, unsigned int const* nBoundStates)
+inline bool LangmuirParamHandler::validate(unsigned int nComp, unsigned int const* nBoundStates)
 {
 	if ((_kA.size() != _kD.size()) || (_kA.size() != _qMax.size()) || (_kA.size() < nComp))
 		throw InvalidParameterException("MCL_KA, MCL_KD, and MCL_QMAX have to have the same size");
@@ -61,7 +61,7 @@ inline bool LangmuirParamHandler::validateConfig(unsigned int nComp, unsigned in
 
 inline const char* ExtLangmuirParamHandler::identifier() CADET_NOEXCEPT { return "EXT_MULTI_COMPONENT_LANGMUIR"; }
 
-inline bool ExtLangmuirParamHandler::validateConfig(unsigned int nComp, unsigned int const* nBoundStates)
+inline bool ExtLangmuirParamHandler::validate(unsigned int nComp, unsigned int const* nBoundStates)
 {
 	if ((_kA.size() != _kD.size()) || (_kA.size() != _qMax.size()) || (_kA.size() < nComp))
 		throw InvalidParameterException("EXT_MCL_KA, EXT_MCL_KD, and EXT_MCL_QMAX have to have the same size");
diff --git a/src/libcadet/model/binding/LinearBinding.cpp b/src/libcadet/model/binding/LinearBinding.cpp
index 57dd58d46..6e03833fa 100644
--- a/src/libcadet/model/binding/LinearBinding.cpp
+++ b/src/libcadet/model/binding/LinearBinding.cpp
@@ -99,11 +99,11 @@ class LinearParamHandler : public ConstParamHandlerBase
 	 * @param [in] nBoundStates Array with number of bound states for each component
 	 * @return @c true if the parameters were read and validated successfully, otherwise @c false
 	 */
-	inline bool configure(IParameterProvider& paramProvider, unsigned int nComp, unsigned int const* nBoundStates)
+	inline bool configureAndValidate(IParameterProvider& paramProvider, unsigned int nComp, unsigned int const* nBoundStates)
 	{
 		_kA.configure("LIN_KA", paramProvider, nComp, nBoundStates);
 		_kD.configure("LIN_KD", paramProvider, nComp, nBoundStates);
-		return validateConfig(nComp, nBoundStates);
+		return validate(nComp, nBoundStates);
 	}
 
 	/**
@@ -164,15 +164,13 @@ class LinearParamHandler : public ConstParamHandlerBase
 		return std::make_tuple<ParamsHandle, ParamsHandle>(&_localParams, nullptr);
 	}
 
-protected:
-
 	/**
 	 * @brief Validates recently read parameters
 	 * @param [in] nComp Number of components
 	 * @param [in] nBoundStates Array with number of bound states for each component
 	 * @return @c true if the parameters were validated successfully, otherwise @c false
 	 */
-	inline bool validateConfig(unsigned int nComp, unsigned int const* nBoundStates)
+	inline bool validate(unsigned int nComp, unsigned int const* nBoundStates)
 	{
 		if ((_kA.size() != _kD.size()) || (_kA.size() < nComp))
 			throw InvalidParameterException("LIN_KA and LIN_KD have to have the same size");
@@ -180,6 +178,8 @@ class LinearParamHandler : public ConstParamHandlerBase
 		return true;
 	}
 
+protected:
+
 	ConstParams _localParams; //!< Actual parameter data
 
 	// Handlers provide configure(), reserve(), and registerParam() for parameters
@@ -228,14 +228,14 @@ class ExtLinearParamHandler : public ExternalParamHandlerBase
 	 * @param [in] nBoundStates Array with number of bound states for each component
 	 * @return @c true if the parameters were read and validated successfully, otherwise @c false
 	 */
-	inline bool configure(IParameterProvider& paramProvider, unsigned int nComp, unsigned int const* nBoundStates)
+	inline bool configureAndValidate(IParameterProvider& paramProvider, unsigned int nComp, unsigned int const* nBoundStates)
 	{
 		_kA.configure("LIN_KA", paramProvider, nComp, nBoundStates);
 		_kD.configure("LIN_KD", paramProvider, nComp, nBoundStates);
 		
 		// Number of externally dependent parameters (2) needs to be given to ExternalParamHandlerBase::configure()
 		ExternalParamHandlerBase::configure(paramProvider, 2);
-		return validateConfig(nComp, nBoundStates);
+		return validate(nComp, nBoundStates);
 	}
 
 	/**
@@ -358,15 +358,13 @@ class ExtLinearParamHandler : public ExternalParamHandlerBase
 			+ 2 * (_kA.additionalDynamicMemory(nComp, totalNumBoundStates, nBoundStates) + _kD.additionalDynamicMemory(nComp, totalNumBoundStates, nBoundStates));
 	}
 
-protected:
-
 	/**
 	 * @brief Validates recently read parameters
 	 * @param [in] nComp Number of components
 	 * @param [in] nBoundStates Array with number of bound states for each component
 	 * @return @c true if the parameters were validated successfully, otherwise @c false
 	 */
-	inline bool validateConfig(unsigned int nComp, unsigned int const* nBoundStates)
+	inline bool validate(unsigned int nComp, unsigned int const* nBoundStates)
 	{
 		if ((_kA.size() != _kD.size()) || (_kA.size() < nComp))
 			throw InvalidParameterException("LIN_KA and LIN_KD have to have the same size");
@@ -374,6 +372,8 @@ class ExtLinearParamHandler : public ExternalParamHandlerBase
 		return true;
 	}
 
+protected:
+
 	// Handlers provide configure(), reserve(), and registerParam() for parameters
 	ExternalScalarComponentDependentParameter _kA; //!< Handler for adsorption rate
 	ExternalScalarComponentDependentParameter _kD; //!< Handler for desorption rate
@@ -446,12 +446,12 @@ class LinearBindingBase : public IBindingModel
 		_parameters.clear();
 
 		// Read parameters (k_a and k_d)
-		_paramHandler.configure(paramProvider, _nComp, _nBoundStates);
+		const bool valid = _paramHandler.configureAndValidate(paramProvider, _nComp, _nBoundStates);
 
 		// Register parameters
 		_paramHandler.registerParameters(_parameters, unitOpIdx, parTypeIdx, _nComp, _nBoundStates);
 
-		return true;
+		return valid;
 	}
 
 	virtual void fillBoundPhaseInitialParameters(ParameterId* params, UnitOpIdx unitOpIdx, ParticleTypeIdx parTypeIdx) const CADET_NOEXCEPT
diff --git a/src/libcadet/model/binding/MobilePhaseModulatorLangmuirBinding.cpp b/src/libcadet/model/binding/MobilePhaseModulatorLangmuirBinding.cpp
index 47ecb00a0..9c840c501 100644
--- a/src/libcadet/model/binding/MobilePhaseModulatorLangmuirBinding.cpp
+++ b/src/libcadet/model/binding/MobilePhaseModulatorLangmuirBinding.cpp
@@ -55,7 +55,7 @@ namespace model
 
 inline const char* MPMLangmuirParamHandler::identifier() CADET_NOEXCEPT { return "MOBILE_PHASE_MODULATOR"; }
 
-inline bool MPMLangmuirParamHandler::validateConfig(unsigned int nComp, unsigned int const* nBoundStates)
+inline bool MPMLangmuirParamHandler::validate(unsigned int nComp, unsigned int const* nBoundStates)
 {
 	if ((_kA.size() != _kD.size()) || (_kA.size() != _qMax.size()) || (_kA.size() != _gamma.size())
 		|| (_kA.size() != _beta.size()) || (_kA.size() < nComp))
@@ -66,7 +66,7 @@ inline bool MPMLangmuirParamHandler::validateConfig(unsigned int nComp, unsigned
 
 inline const char* ExtMPMLangmuirParamHandler::identifier() CADET_NOEXCEPT { return "EXT_MOBILE_PHASE_MODULATOR"; }
 
-inline bool ExtMPMLangmuirParamHandler::validateConfig(unsigned int nComp, unsigned int const* nBoundStates)
+inline bool ExtMPMLangmuirParamHandler::validate(unsigned int nComp, unsigned int const* nBoundStates)
 {
 	if ((_kA.size() != _kD.size()) || (_kA.size() != _qMax.size()) || (_kA.size() != _gamma.size())
 		|| (_kA.size() != _beta.size()) || (_kA.size() < nComp))
diff --git a/src/libcadet/model/binding/MultiComponentSpreadingBinding.cpp b/src/libcadet/model/binding/MultiComponentSpreadingBinding.cpp
index dc324a47c..8fc12a84d 100644
--- a/src/libcadet/model/binding/MultiComponentSpreadingBinding.cpp
+++ b/src/libcadet/model/binding/MultiComponentSpreadingBinding.cpp
@@ -55,7 +55,7 @@ namespace model
 
 inline const char* SpreadingParamHandler::identifier() CADET_NOEXCEPT { return "MULTI_COMPONENT_SPREADING"; }
 
-inline bool SpreadingParamHandler::validateConfig(unsigned int nComp, unsigned int const* nBoundStates)
+inline bool SpreadingParamHandler::validate(unsigned int nComp, unsigned int const* nBoundStates)
 {
 	if ((_kA.size() != _kD.size()) || (_kA.size() != _qMax.size()) || (_kA.size() < nComp * 2))
 		throw InvalidParameterException("MCSPR_KA, MCSPR_KD, and MCSPR_QMAX have to have the same size");
@@ -67,7 +67,7 @@ inline bool SpreadingParamHandler::validateConfig(unsigned int nComp, unsigned i
 
 inline const char* ExtSpreadingParamHandler::identifier() CADET_NOEXCEPT { return "EXT_MULTI_COMPONENT_SPREADING"; }
 
-inline bool ExtSpreadingParamHandler::validateConfig(unsigned int nComp, unsigned int const* nBoundStates)
+inline bool ExtSpreadingParamHandler::validate(unsigned int nComp, unsigned int const* nBoundStates)
 {
 	if ((_kA.size() != _kD.size()) || (_kA.size() != _qMax.size()) || (_kA.size() < nComp * 2))
 		throw InvalidParameterException("EXT_MCSPR_KA, EXT_MCSPR_KD, and EXT_MCSPR_QMAX have to have the same size");
diff --git a/src/libcadet/model/binding/MultiStateStericMassActionBinding.cpp b/src/libcadet/model/binding/MultiStateStericMassActionBinding.cpp
index 118273cfe..b921fd739 100644
--- a/src/libcadet/model/binding/MultiStateStericMassActionBinding.cpp
+++ b/src/libcadet/model/binding/MultiStateStericMassActionBinding.cpp
@@ -64,7 +64,7 @@ namespace model
 
 inline const char* MSSMAParamHandler::identifier() CADET_NOEXCEPT { return "MULTISTATE_STERIC_MASS_ACTION"; }
 
-inline bool MSSMAParamHandler::validateConfig(unsigned int nComp, unsigned int const* nBoundStates)
+inline bool MSSMAParamHandler::validate(unsigned int nComp, unsigned int const* nBoundStates)
 {
 	const unsigned int totalBoundStates = numBoundStates(nBoundStates, nComp);
 
@@ -91,7 +91,7 @@ inline bool MSSMAParamHandler::validateConfig(unsigned int nComp, unsigned int c
 
 inline const char* ExtMSSMAParamHandler::identifier() CADET_NOEXCEPT { return "EXT_MULTISTATE_STERIC_MASS_ACTION"; }
 
-inline bool ExtMSSMAParamHandler::validateConfig(unsigned int nComp, unsigned int const* nBoundStates)
+inline bool ExtMSSMAParamHandler::validate(unsigned int nComp, unsigned int const* nBoundStates)
 {
 	const unsigned int totalBoundStates = numBoundStates(nBoundStates, nComp);
 
diff --git a/src/libcadet/model/binding/SaskaBinding.cpp b/src/libcadet/model/binding/SaskaBinding.cpp
index a49c67212..89559e8ec 100644
--- a/src/libcadet/model/binding/SaskaBinding.cpp
+++ b/src/libcadet/model/binding/SaskaBinding.cpp
@@ -50,7 +50,7 @@ namespace model
 
 inline const char* SaskaParamHandler::identifier() CADET_NOEXCEPT { return "SASKA"; }
 
-inline bool SaskaParamHandler::validateConfig(unsigned int nComp, unsigned int const* nBoundStates)
+inline bool SaskaParamHandler::validate(unsigned int nComp, unsigned int const* nBoundStates)
 {
 	if ((_h.size() != _k.slices()) || (_k.size() != nComp * nComp) || (_h.size() < nComp))
 		throw InvalidParameterException("SASKA_K has to have NCOMP^2 and SASKA_H NCOMP many elements");
@@ -60,7 +60,7 @@ inline bool SaskaParamHandler::validateConfig(unsigned int nComp, unsigned int c
 
 inline const char* ExtSaskaParamHandler::identifier() CADET_NOEXCEPT { return "EXT_SASKA"; }
 
-inline bool ExtSaskaParamHandler::validateConfig(unsigned int nComp, unsigned int const* nBoundStates)
+inline bool ExtSaskaParamHandler::validate(unsigned int nComp, unsigned int const* nBoundStates)
 {
 	if ((_h.size() != _k.slices()) || (_k.size() != nComp * nComp) || (_h.size() < nComp))
 		throw InvalidParameterException("EXT_SASKA_K has to have NCOMP^2 and EXT_SASKA_H NCOMP many elements");
diff --git a/src/libcadet/model/binding/SelfAssociationBinding.cpp b/src/libcadet/model/binding/SelfAssociationBinding.cpp
index 90aa8e365..885e21eec 100644
--- a/src/libcadet/model/binding/SelfAssociationBinding.cpp
+++ b/src/libcadet/model/binding/SelfAssociationBinding.cpp
@@ -64,7 +64,7 @@ namespace model
 
 inline const char* SelfAssociationParamHandler::identifier() CADET_NOEXCEPT { return "SELF_ASSOCIATION"; }
 
-inline bool SelfAssociationParamHandler::validateConfig(unsigned int nComp, unsigned int const* nBoundStates)
+inline bool SelfAssociationParamHandler::validate(unsigned int nComp, unsigned int const* nBoundStates)
 {
 	if ((_kA.size() != _kA2.size()) || (_kA.size() != _kD.size()) || (_kA.size() != _nu.size()) 
 		|| (_kA.size() != _sigma.size()) || (_kA.size() < nComp))
@@ -79,7 +79,7 @@ inline bool SelfAssociationParamHandler::validateConfig(unsigned int nComp, unsi
 
 inline const char* ExtSelfAssociationParamHandler::identifier() CADET_NOEXCEPT { return "EXT_SELF_ASSOCIATION"; }
 
-inline bool ExtSelfAssociationParamHandler::validateConfig(unsigned int nComp, unsigned int const* nBoundStates)
+inline bool ExtSelfAssociationParamHandler::validate(unsigned int nComp, unsigned int const* nBoundStates)
 {
 	if ((_kA.size() != _kA2.size()) || (_kA.size() != _kD.size()) || (_kA.size() != _nu.size()) 
 		|| (_kA.size() != _sigma.size()) || (_kA.size() < nComp))
diff --git a/src/libcadet/model/binding/StericMassActionBinding.cpp b/src/libcadet/model/binding/StericMassActionBinding.cpp
index b783cd6b9..0ea5aeaae 100644
--- a/src/libcadet/model/binding/StericMassActionBinding.cpp
+++ b/src/libcadet/model/binding/StericMassActionBinding.cpp
@@ -62,7 +62,7 @@ namespace model
 
 inline const char* SMAParamHandler::identifier() CADET_NOEXCEPT { return "STERIC_MASS_ACTION"; }
 
-inline bool SMAParamHandler::validateConfig(unsigned int nComp, unsigned int const* nBoundStates)
+inline bool SMAParamHandler::validate(unsigned int nComp, unsigned int const* nBoundStates)
 {
 	if ((_kA.size() != _kD.size()) || (_kA.size() != _nu.size()) || (_kA.size() != _sigma.size()) || (_kA.size() < nComp))
 		throw InvalidParameterException("SMA_KA, SMA_KD, SMA_NU, and SMA_SIGMA have to have the same size");
@@ -76,7 +76,7 @@ inline bool SMAParamHandler::validateConfig(unsigned int nComp, unsigned int con
 
 inline const char* ExtSMAParamHandler::identifier() CADET_NOEXCEPT { return "EXT_STERIC_MASS_ACTION"; }
 
-inline bool ExtSMAParamHandler::validateConfig(unsigned int nComp, unsigned int const* nBoundStates)
+inline bool ExtSMAParamHandler::validate(unsigned int nComp, unsigned int const* nBoundStates)
 {
 	if ((_kA.size() != _kD.size()) || (_kA.size() != _nu.size()) || (_kA.size() != _sigma.size()) || (_kA.size() < nComp))
 		throw InvalidParameterException("EXT_SMA_KA, EXT_SMA_KD, EXT_SMA_NU, and EXT_SMA_SIGMA have to have the same size");
diff --git a/src/libcadet/model/reaction/ExternalFunctionTemplate.cpp b/src/libcadet/model/reaction/ExternalFunctionTemplate.cpp
index d0ac4e16b..37f19b530 100644
--- a/src/libcadet/model/reaction/ExternalFunctionTemplate.cpp
+++ b/src/libcadet/model/reaction/ExternalFunctionTemplate.cpp
@@ -83,7 +83,7 @@ class {{ name }} : public cadet::model::ConstParamHandlerBase
 {% endif %}
 	{ }
 
-	inline bool configure(IParameterProvider& paramProvider, unsigned int nReactions, unsigned int nComp, unsigned int const* nBoundStates)
+	inline void configure(IParameterProvider& paramProvider, unsigned int nReactions, unsigned int nComp, unsigned int const* nBoundStates)
 	{
 {% for p in parameters %}
 	{% if not existsIn(p, "skipConfig") %}
@@ -109,7 +109,12 @@ class {{ name }} : public cadet::model::ConstParamHandlerBase
 		{% endif %}
 	{% endfor %}
 {% endif %}
-		return validateConfig(nReactions, nComp, nBoundStates);
+	}
+
+	inline bool configureAndValidate(IParameterProvider& paramProvider, unsigned int nReactions, unsigned int nComp, unsigned int const* nBoundStates)
+	{
+		configure(paramProvider, nReactions, nComp, nBoundStates);
+		return validate(nReactions, nComp, nBoundStates);
 	}
 
 	inline void registerParameters(std::unordered_map<ParameterId, active*>& parameters, UnitOpIdx unitOpIdx, ParticleTypeIdx parTypeIdx, unsigned int nComp, unsigned int const* nBoundStates)
@@ -171,9 +176,9 @@ class {{ name }} : public cadet::model::ConstParamHandlerBase
 {% endif %}
 
 	inline char const* prefixInConfiguration() const CADET_NOEXCEPT { return ""; }
+	inline bool validate(unsigned int nReactions, unsigned int nComp, unsigned int const* nBoundStates);
 
 protected:
-	inline bool validateConfig(unsigned int nReactions, unsigned int nComp, unsigned int const* nBoundStates);
 
 	ConstParams _localParams;
 
@@ -262,7 +267,7 @@ class {{ externalName }} : public cadet::model::ExternalParamHandlerBase
 {% endif %}
 	{ }
 
-	inline bool configure(IParameterProvider& paramProvider, unsigned int nReactions, unsigned int nComp, unsigned int const* nBoundStates)
+	inline void configure(IParameterProvider& paramProvider, unsigned int nReactions, unsigned int nComp, unsigned int const* nBoundStates)
 	{
 {% for p in parameters %}
 	{% if not existsIn(p, "skipConfig") %}
@@ -289,7 +294,12 @@ class {{ externalName }} : public cadet::model::ExternalParamHandlerBase
 	{% endfor %}
 {% endif %}
 		ExternalParamHandlerBase::configure(paramProvider, {{ length(parameters) }});
-		return validateConfig(nReactions, nComp, nBoundStates);
+	}
+
+	inline bool configureAndValidate(IParameterProvider& paramProvider, unsigned int nReactions, unsigned int nComp, unsigned int const* nBoundStates)
+	{
+		configure(paramProvider, nReactions, nComp, nBoundStates);
+		return validate(nReactions, nComp, nBoundStates);
 	}
 
 	inline void registerParameters(std::unordered_map<ParameterId, active*>& parameters, UnitOpIdx unitOpIdx, ParticleTypeIdx parTypeIdx, unsigned int nComp, unsigned int const* nBoundStates)
@@ -409,9 +419,9 @@ class {{ externalName }} : public cadet::model::ExternalParamHandlerBase
 {% endif %}
 
 	inline char const* prefixInConfiguration() const CADET_NOEXCEPT { return "EXT_"; }
+	inline bool validate(unsigned int nReactions, unsigned int nComp, unsigned int const* nBoundStates);
 
 protected:
-	inline bool validateConfig(unsigned int nReactions, unsigned int nComp, unsigned int const* nBoundStates);
 
 	ConstParams _constParams;
 
diff --git a/src/libcadet/model/reaction/MassActionLawReaction.cpp b/src/libcadet/model/reaction/MassActionLawReaction.cpp
index 3a550ae21..e1658802f 100644
--- a/src/libcadet/model/reaction/MassActionLawReaction.cpp
+++ b/src/libcadet/model/reaction/MassActionLawReaction.cpp
@@ -60,14 +60,14 @@ namespace model
 
 inline const char* MassActionLawParamHandler::identifier() CADET_NOEXCEPT { return "MASS_ACTION_LAW"; }
 
-inline bool MassActionLawParamHandler::validateConfig(unsigned int nReactions, unsigned int nComp, unsigned int const* nBoundStates)
+inline bool MassActionLawParamHandler::validate(unsigned int nReactions, unsigned int nComp, unsigned int const* nBoundStates)
 {
 	return true;
 }
 
 inline const char* ExtMassActionLawParamHandler::identifier() CADET_NOEXCEPT { return "EXT_MASS_ACTION_LAW"; }
 
-inline bool ExtMassActionLawParamHandler::validateConfig(unsigned int nReactions, unsigned int nComp, unsigned int const* nBoundStates)
+inline bool ExtMassActionLawParamHandler::validate(unsigned int nReactions, unsigned int nComp, unsigned int const* nBoundStates)
 {
 	return true;
 }
@@ -546,7 +546,7 @@ class MassActionLawReactionBase : public DynamicReactionModelBase
 		readAndRegisterExponents(paramProvider, _parameters, unitOpIdx, parTypeIdx, "MAL_EXPONENTS_SOLID_FWD_MODLIQUID", _expSolidFwdLiquid, _nComp, nullptr);
 		readAndRegisterExponents(paramProvider, _parameters, unitOpIdx, parTypeIdx, "MAL_EXPONENTS_SOLID_BWD_MODLIQUID", _expSolidBwdLiquid, _nComp, nullptr);
 
-		return true;
+		return _paramHandler.validate(maxNumReactions(), _nComp, _nBoundStates);
 	}
 
 	template <typename StateType, typename ResidualType, typename ParamType, typename FactorType>