diff --git a/bin/genEvalSpecializations.py b/bin/genEvalSpecializations.py
index 20e1b1aeee2..bb3af0eb433 100755
--- a/bin/genEvalSpecializations.py
+++ b/bin/genEvalSpecializations.py
@@ -239,6 +239,17 @@ class Evaluation<ValueT, {{ numDerivs }}>
 
         checkDefined_();
     }
+
+    // create a dynamic evaluation which represents a constant function
+    // by using zero number of derivatives.
+    //
+    // i.e., f(x) = c. this implies an evaluation with the given value and all
+    // derivatives being zero.
+    template <class RhsValueType>
+    OPM_HOST_DEVICE Evaluation(const RhsValueType& c)
+        : Evaluation(0, c)
+    {
+    }
 {% else %}\
     // create an evaluation which represents a constant function
     //
@@ -259,7 +270,7 @@ class Evaluation<ValueT, {{ numDerivs }}>
 {% if numDerivs < 0 %}\
     template <class RhsValueType>
     OPM_HOST_DEVICE Evaluation(int nVars, const RhsValueType& c, int varPos)
-     : data_(1 + nVars, 0.0)
+        : data_(1 + nVars, 0.0)
     {
         // The variable position must be in represented by the given variable descriptor
         assert(0 <= varPos && varPos < size());
@@ -400,9 +411,10 @@ class Evaluation<ValueT, {{ numDerivs }}>
     // "evaluate" a constant function (i.e. a function that does not depend on the set of
     // relevant variables, f(x) = c).
     template <class RhsValueType>
-    OPM_HOST_DEVICE static Evaluation createConstant(const RhsValueType&)
+    OPM_HOST_DEVICE static Evaluation createConstant(const RhsValueType& value)
     {
-        throw std::logic_error("Dynamically-sized evaluation objects require to specify the number of derivatives.");
+        // using an evaluation without derivatives
+        return Evaluation(value);
     }
 
     // "evaluate" a constant function (i.e. a function that does not depend on the set of
@@ -457,12 +469,41 @@ class Evaluation<ValueT, {{ numDerivs }}>
     }
 
 
-    // add value and derivatives from other to this values and derivatives
+    // add value and derivatives from other to this value and derivatives
     OPM_HOST_DEVICE Evaluation& operator+=(const Evaluation& other)
     {
+{% if numDerivs < 0 %}\
+        const int thisSize = size();
+        const int otherSize = other.size();
+
+        // Allow operation if sizes match or one is constant (size == 0)
+        assert(thisSize == otherSize || thisSize == 0 || otherSize == 0);
+
+        if (thisSize == otherSize) {
+            for (int i = 0; i < length_(); ++i)
+                data_[i] += other.data_[i];
+            return *this;
+        }
+        if (otherSize == 0) {
+            *this += other.value();
+            return *this;
+        }
+        if (thisSize == 0) {
+            assert(otherSize > 0);
+            this->appendDerivativesToConstant(otherSize);
+            data_[valuepos_()] += other.data_[valuepos_()];
+            for (int i = dstart_(); i < dend_(); ++i)
+                data_[i] = other.data_[i];
+            return *this;
+        }
+        throw std::logic_error(
+                "Cannot operate Evaluations with different number of derivatives "
+                "unless one of them has no derivatives"
+        );
+{% else %}\
         assert(size() == other.size());
 
-{% if numDerivs <= 0 %}\
+{% if numDerivs == 0 %}\
         for (int i = 0; i < length_(); ++i)
             data_[i] += other.data_[i];
 {% else %}\
@@ -472,6 +513,7 @@ class Evaluation<ValueT, {{ numDerivs }}>
 {% endif %}\
 
         return *this;
+{% endif %}\
     }
 
     // add value from other to this values
@@ -487,9 +529,37 @@ class Evaluation<ValueT, {{ numDerivs }}>
     // subtract other's value and derivatives from this values
     OPM_HOST_DEVICE Evaluation& operator-=(const Evaluation& other)
     {
+{% if numDerivs < 0 %}\
+        const int thisSize = size();
+        const int otherSize = other.size();
+
+        // Allow operation if sizes match or one is constant (size == 0)
+        assert(thisSize == otherSize || thisSize == 0 || otherSize == 0);
+
+        if (thisSize == otherSize) {
+            for (int i = 0; i < length_(); ++i)
+                data_[i] -= other.data_[i];
+            return *this;
+        }
+        if (otherSize == 0) {
+            *this -= other.value();
+            return *this;
+        }
+        if (thisSize == 0) {
+            assert(otherSize > 0);
+            this->appendDerivativesToConstant(otherSize);
+            for (int i = 0; i < other.length_(); ++i)
+                data_[i] -= other.data_[i];
+            return *this;
+        }
+        throw std::logic_error(
+                "Cannot operate Evaluations with different number of derivatives "
+                "unless one of them has no derivatives"
+        );
+{% else %}\
         assert(size() == other.size());
 
-{% if numDerivs <= 0 %}\
+{% if numDerivs == 0 %}\
         for (int i = 0; i < length_(); ++i)
             data_[i] -= other.data_[i];
 {% else %}\
@@ -499,6 +569,7 @@ class Evaluation<ValueT, {{ numDerivs }}>
 {% endif %}\
 
         return *this;
+{% endif %}\
     }
 
     // subtract other's value from this values
@@ -514,6 +585,48 @@ class Evaluation<ValueT, {{ numDerivs }}>
     // multiply values and apply chain rule to derivatives: (u*v)' = (v'u + u'v)
     OPM_HOST_DEVICE Evaluation& operator*=(const Evaluation& other)
     {
+{% if numDerivs < 0 %}\
+        const int thisSize = size();
+        const int otherSize = other.size();
+
+        // Allow operation if sizes match or one is constant (size == 0)
+        assert(thisSize == otherSize || thisSize == 0 || otherSize == 0);
+
+        if (thisSize == otherSize) {
+            // while the values are multiplied, the derivatives follow the product rule,
+            // i.e., (u*v)' = (v'u + u'v).
+            const ValueType u = this->value();
+            const ValueType v = other.value();
+
+            // value
+            data_[valuepos_()] *= v ;
+
+            //  derivatives
+            for (int i = dstart_(); i < dend_(); ++i)
+                data_[i] = data_[i] * v + other.data_[i] * u;
+
+            return *this;
+        }
+
+        if (otherSize == 0) {
+            *this *= other.value();
+            return *this;
+        }
+
+        if (thisSize == 0) {
+            assert(otherSize > 0);
+            this->appendDerivativesToConstant(otherSize);
+            const ValueType u = this->value();
+            for (int i = 0; i < length_(); ++i) {
+                data_[i] = u * other.data_[i];
+            }
+            return *this;
+        }
+        throw std::logic_error(
+                "Cannot operate Evaluations with different number of derivatives "
+                "unless one of them has no derivatives"
+        );
+{% else %}\
         assert(size() == other.size());
 
         // while the values are multiplied, the derivatives follow the product rule,
@@ -525,7 +638,7 @@ class Evaluation<ValueT, {{ numDerivs }}>
         data_[valuepos_()] *= v ;
 
         //  derivatives
-{% if numDerivs <= 0 %}\
+{% if numDerivs == 0 %}\
         for (int i = dstart_(); i < dend_(); ++i)
             data_[i] = data_[i] * v + other.data_[i] * u;
 {% else %}\
@@ -535,6 +648,7 @@ class Evaluation<ValueT, {{ numDerivs }}>
 {% endif %}\
 
         return *this;
+{% endif %}\
     }
 
     // m(c*u)' = c*u'
@@ -556,13 +670,67 @@ class Evaluation<ValueT, {{ numDerivs }}>
     // m(u*v)' = (vu' - uv')/v^2
     OPM_HOST_DEVICE Evaluation& operator/=(const Evaluation& other)
     {
+{% if numDerivs < 0 %}\
+        const int thisSize = size();
+        const int otherSize = other.size();
+
+        // Allow operation if sizes match or one is constant (size == 0)
+        assert(thisSize == otherSize || thisSize == 0 || otherSize == 0);
+
+        if (thisSize == otherSize) {
+            // while the values are divided, the derivatives follow the rule for division,
+            // i.e., (u/v)' = (v'u - u'v)/v^2.
+            ValueType& u = data_[valuepos_()];
+            const ValueType& v = other.value();
+
+            // derivatives
+            for (int idx = dstart_(); idx < dend_(); ++idx) {
+                const ValueType& uPrime = data_[idx];
+                const ValueType& vPrime = other.data_[idx];
+
+                data_[idx] = (v*uPrime - u*vPrime)/(v*v);
+            }
+
+            // value
+            u /= v;
+
+            return *this;
+        }
+
+        if (otherSize == 0) {
+            *this /= other.value();
+            return *this;
+        }
+
+        if (thisSize == 0) {
+            assert(otherSize > 0);
+            // if *this is a constant, extend to hold derivatives before division
+            this->appendDerivativesToConstant(otherSize);
+
+            // values are divided, derivatives follow the rule for division, i.e., (u/v)' = (v'u -
+            // u'v)/v^2.
+            ValueType& u = data_[valuepos_()];
+            const ValueType& v = other.value();
+            for (int idx = dstart_(); idx < dend_(); ++idx) {
+                const ValueType& vPrime = other.data_[idx];
+                data_[idx] =  -u * vPrime / (v * v);
+            }
+            u /= v;
+
+            return *this;
+        }
+        throw std::logic_error(
+                "Cannot operate Evaluations with different number of derivatives "
+                "unless one of them has no derivatives"
+        );
+{% else %}\
         assert(size() == other.size());
 
         // values are divided, derivatives follow the rule for division, i.e., (u/v)' = (v'u -
         // u'v)/v^2.
         ValueType& u = data_[valuepos_()];
         const ValueType& v = other.value();
-{% if numDerivs <= 0 %}\
+{% if numDerivs == 0 %}\
         for (int idx = dstart_(); idx < dend_(); ++idx) {
             const ValueType& uPrime = data_[idx];
             const ValueType& vPrime = other.data_[idx];
@@ -577,6 +745,7 @@ class Evaluation<ValueT, {{ numDerivs }}>
         u /= v;
 
         return *this;
+{% endif %}\
     }
 
     // divide value and derivatives by value of other
@@ -600,7 +769,15 @@ class Evaluation<ValueT, {{ numDerivs }}>
     // add two evaluation objects
     OPM_HOST_DEVICE Evaluation operator+(const Evaluation& other) const
     {
+{% if numDerivs < 0 %}\
+        const int thisSize = size();
+        const int otherSize = other.size();
+
+        // Allow operation if sizes match or one is constant (size == 0)
+        assert(thisSize == otherSize || thisSize == 0 || otherSize == 0);
+{% else %}\
         assert(size() == other.size());
+{% endif %}\
 
         Evaluation result(*this);
 
@@ -623,7 +800,15 @@ class Evaluation<ValueT, {{ numDerivs }}>
     // subtract two evaluation objects
     OPM_HOST_DEVICE Evaluation operator-(const Evaluation& other) const
     {
+{% if numDerivs < 0 %}\
+        const int thisSize = size();
+        const int otherSize = other.size();
+
+        // Allow operation if sizes match or one is constant (size == 0)
+        assert(thisSize == otherSize || thisSize == 0 || otherSize == 0);
+{% else %}\
         assert(size() == other.size());
+{% endif %}\
 
         Evaluation result(*this);
 
@@ -667,7 +852,15 @@ class Evaluation<ValueT, {{ numDerivs }}>
 
     OPM_HOST_DEVICE Evaluation operator*(const Evaluation& other) const
     {
+{% if numDerivs < 0 %}\
+        const int thisSize = size();
+        const int otherSize = other.size();
+
+        // Allow operation if sizes match or one is constant (size == 0)
+        assert(thisSize == otherSize || thisSize == 0 || otherSize == 0);
+{% else %}\
         assert(size() == other.size());
+{% endif %}\
 
         Evaluation result(*this);
 
@@ -688,7 +881,15 @@ class Evaluation<ValueT, {{ numDerivs }}>
 
     OPM_HOST_DEVICE Evaluation operator/(const Evaluation& other) const
     {
+{% if numDerivs < 0 %}\
+        const int thisSize = size();
+        const int otherSize = other.size();
+
+        // Allow operation if sizes match or one is constant (size == 0)
+        assert(thisSize == otherSize || thisSize == 0 || otherSize == 0);
+{% else %}\
         assert(size() == other.size());
+{% endif %}\
 
         Evaluation result(*this);
 
@@ -725,6 +926,23 @@ class Evaluation<ValueT, {{ numDerivs }}>
 
     OPM_HOST_DEVICE bool operator==(const Evaluation& other) const
     {
+{% if numDerivs < 0 %}\
+        const int thisSize = size();
+        const int otherSize = other.size();
+
+        // Allow operation if sizes match or one is constant (size == 0)
+        assert(thisSize == otherSize || thisSize == 0 || otherSize == 0);
+
+        if (thisSize == otherSize) {
+            for (int idx = 0; idx < length_(); ++idx) {
+                if (data_[idx] != other.data_[idx]) {
+                    return false;
+                }
+            }
+        } else {
+            return value() == other.value();
+        }
+{% else %}\
         assert(size() == other.size());
 
         for (int idx = 0; idx < length_(); ++idx) {
@@ -732,6 +950,7 @@ class Evaluation<ValueT, {{ numDerivs }}>
                 return false;
             }
         }
+{% endif %}\
         return true;
     }
 
@@ -748,7 +967,11 @@ class Evaluation<ValueT, {{ numDerivs }}>
 
     OPM_HOST_DEVICE bool operator>(const Evaluation& other) const
     {
+{% if numDerivs < 0 %}\
+        assert(size() == other.size() || size() == 0 || other.size() == 0);
+{% else %}\
         assert(size() == other.size());
+{% endif %}\
 
         return value() > other.value();
     }
@@ -759,7 +982,11 @@ class Evaluation<ValueT, {{ numDerivs }}>
 
     OPM_HOST_DEVICE bool operator<(const Evaluation& other) const
     {
+{% if numDerivs < 0 %}\
+        assert(size() == other.size() || size() == 0 || other.size() == 0);
+{% else %}\
         assert(size() == other.size());
+{% endif %}\
 
         return value() < other.value();
     }
@@ -770,7 +997,11 @@ class Evaluation<ValueT, {{ numDerivs }}>
 
     OPM_HOST_DEVICE bool operator>=(const Evaluation& other) const
     {
+{% if numDerivs < 0 %}\
+        assert(size() == other.size() || size() == 0 || other.size() == 0);
+{% else %}\
         assert(size() == other.size());
+{% endif %}\
 
         return value() >= other.value();
     }
@@ -781,7 +1012,11 @@ class Evaluation<ValueT, {{ numDerivs }}>
 
     OPM_HOST_DEVICE bool operator<=(const Evaluation& other) const
     {
+{% if numDerivs < 0 %}\
+        assert(size() == other.size() || size() == 0 || other.size() == 0);
+{% else %}\
         assert(size() == other.size());
+{% endif %}\
 
         return value() <= other.value();
     }
@@ -798,14 +1033,43 @@ class Evaluation<ValueT, {{ numDerivs }}>
     // return varIdx'th derivative
     OPM_HOST_DEVICE const ValueType& derivative(int varIdx) const
     {
+{% if numDerivs < 0 %}\
+        assert(size() == 0 || (0 <= varIdx && varIdx < size()) );
+
+        if (size() == 0) {
+            static const ValueType zero {0.0};
+            return zero;
+        }
+{% else %}\
         assert(0 <= varIdx && varIdx < size());
+{% endif %}\
 
         return data_[dstart_() + varIdx];
     }
 
     // set derivative at position varIdx
+{% if numDerivs < 0 %}\
+    OPM_HOST_DEVICE void setDerivative(int varIdx, const ValueType& derVal, const int nVars = -1)
+{% else %}\
     OPM_HOST_DEVICE void setDerivative(int varIdx, const ValueType& derVal)
+{% endif %}\
     {
+{% if numDerivs < 0 %}\
+        // if size() == 0, we need nVars to be positive to extend the number of derivatives
+        // if size() > 0, nVars must be either negative (i.e. ignore) or match size()
+        assert( (size() == 0 && nVars > 0) ||
+                ((size() > 0 && (nVars < 0 || nVars == size()))) );
+
+        if (size() == 0) {
+            if (nVars < 0) {
+                throw std::logic_error("Cannot set derivative for a DynamicEvaluation initialized from a scalar "
+                                       "without specifying a positive number of derivatives");
+            }
+
+            this->appendDerivativesToConstant(nVars);
+        }
+
+{% endif %}\
         assert(0 <= varIdx && varIdx < size());
 
         data_[dstart_() + varIdx] = derVal;
@@ -820,6 +1084,13 @@ class Evaluation<ValueT, {{ numDerivs }}>
 private:
 {% if numDerivs < 0 %}\
     FastSmallVector<ValueT, staticSize> data_;
+
+    void appendDerivativesToConstant(size_t numDer) {
+        assert(size() == 0); // we only append derivatives to a constant
+        if (numDer > 0) {
+            data_.resize(1 + numDer);
+        }
+    }
 {% elif numDerivs == 0 %}\
     std::array<ValueT, numDerivs + 1> data_;
 {% else %}\
diff --git a/opm/material/common/FastSmallVector.hpp b/opm/material/common/FastSmallVector.hpp
index 1ad71178f10..75e2fb0f839 100644
--- a/opm/material/common/FastSmallVector.hpp
+++ b/opm/material/common/FastSmallVector.hpp
@@ -99,7 +99,7 @@ class FastSmallVector
     FastSmallVector& operator=(FastSmallVector&& other)
     {
         size_ = other.size_;
-        if (size_ <= N) {
+        if (other.usingSmallBuf()) {
             smallBuf_ = std::move(other.smallBuf_);
             dataPtr_ = smallBuf_.data();
         }
@@ -119,7 +119,7 @@ class FastSmallVector
     {
         size_ = other.size_;
 
-        if (size_ <= N) {
+        if (other.usingSmallBuf()) {
             smallBuf_ = other.smallBuf_;
             dataPtr_ = smallBuf_.data();
         }
@@ -176,20 +176,23 @@ class FastSmallVector
     { return dataPtr_ + size_; }
 
     size_type capacity()
-    { return size_ <= N ? N : data_.capacity(); }
+    { return this->usingSmallBuf() ? N : data_.capacity(); }
 
     void push_back(const ValueType& value)
     {
-        if (size_ < N) {
-            // Data is contained in smallBuf_
-            smallBuf_[size_++] = value;
-        } else if (size_ == N) {
-            // Must switch from using smallBuf_ to using data_
-            data_.reserve(N + 1);
-            data_.assign(smallBuf_.begin(), smallBuf_.end());
-            data_.push_back(value);
-            ++size_;
-            dataPtr_ = data_.data();
+        if (this->usingSmallBuf()) {
+            if (size_ < N) {
+                // Data is contained in smallBuf_
+                smallBuf_[size_++] = value;
+            } else if (size_ == N) {
+                // Must switch from using smallBuf_ to using data_
+                data_.reserve(N + 1); // Since we will push_back to it later
+                data_.resize(N);      // So we can move the existing data to it
+                std::move(smallBuf_.begin(), smallBuf_.end(), data_.begin());
+                data_.push_back(value);
+                ++size_;
+                dataPtr_ = data_.data();
+            }
         } else {
             // Data is contained in data_
             data_.push_back(value);
@@ -198,6 +201,31 @@ class FastSmallVector
         }
     }
 
+    void resize(size_t numElem)
+    {
+        if (numElem == size_) return; // nothing to do
+
+        if (this->usingSmallBuf()) {
+            if (numElem > N) {
+                data_.resize(numElem);
+                if (size_ > 0 && N > 0) {
+                    std::copy(smallBuf_.begin(), smallBuf_.begin() + size_, data_.begin());
+                }
+                dataPtr_ = data_.data();
+            } else if (numElem < size_) {
+                // when shrinking, remove the values after numElem so that the space
+                // is ready to use in the potentional future resize
+                for (auto ii = numElem; ii < size_; ++ii) {
+                    smallBuf_[ii].~ValueType();
+                }
+            }
+        } else {
+            // when shriking to numElem < N, we do not switch back to use smallBuf_
+            data_.resize(numElem);
+        }
+        size_ = numElem;
+    }
+
 private:
     void init_(size_t numElem)
     {
@@ -210,6 +238,11 @@ class FastSmallVector
             dataPtr_ = smallBuf_.data();
     }
 
+    bool usingSmallBuf() const
+    {
+        return dataPtr_ == smallBuf_.data();
+    }
+
     std::array<ValueType, N> smallBuf_{};
     std::vector<ValueType> data_;
     size_type size_;
diff --git a/opm/material/densead/DynamicEvaluation.hpp b/opm/material/densead/DynamicEvaluation.hpp
index 4073bf358a2..20f2fea885c 100644
--- a/opm/material/densead/DynamicEvaluation.hpp
+++ b/opm/material/densead/DynamicEvaluation.hpp
@@ -131,11 +131,22 @@ class Evaluation<ValueT, DynamicSize, staticSize>
         checkDefined_();
     }
 
+    // create a dynamic evaluation which represents a constant function
+    // by using zero number of derivatives.
+    //
+    // i.e., f(x) = c. this implies an evaluation with the given value and all
+    // derivatives being zero.
+    template <class RhsValueType>
+    OPM_HOST_DEVICE Evaluation(const RhsValueType& c)
+        : Evaluation(0, c)
+    {
+    }
+
     // create an evaluation representing a variable with the variable position of varPos
     // The value is set to c, all derivatives are zero except for the one at varPos, which is set to 1.
     template <class RhsValueType>
     OPM_HOST_DEVICE Evaluation(int nVars, const RhsValueType& c, int varPos)
-     : data_(1 + nVars, 0.0)
+        : data_(1 + nVars, 0.0)
     {
         // The variable position must be in represented by the given variable descriptor
         assert(0 <= varPos && varPos < size());
@@ -209,9 +220,10 @@ class Evaluation<ValueT, DynamicSize, staticSize>
     // "evaluate" a constant function (i.e. a function that does not depend on the set of
     // relevant variables, f(x) = c).
     template <class RhsValueType>
-    OPM_HOST_DEVICE static Evaluation createConstant(const RhsValueType&)
+    OPM_HOST_DEVICE static Evaluation createConstant(const RhsValueType& value)
     {
-        throw std::logic_error("Dynamically-sized evaluation objects require to specify the number of derivatives.");
+        // using an evaluation without derivatives
+        return Evaluation(value);
     }
 
     // "evaluate" a constant function (i.e. a function that does not depend on the set of
@@ -232,15 +244,36 @@ class Evaluation<ValueT, DynamicSize, staticSize>
     }
 
 
-    // add value and derivatives from other to this values and derivatives
+    // add value and derivatives from other to this value and derivatives
     OPM_HOST_DEVICE Evaluation& operator+=(const Evaluation& other)
     {
-        assert(size() == other.size());
+        const int thisSize = size();
+        const int otherSize = other.size();
 
-        for (int i = 0; i < length_(); ++i)
-            data_[i] += other.data_[i];
+        // Allow operation if sizes match or one is constant (size == 0)
+        assert(thisSize == otherSize || thisSize == 0 || otherSize == 0);
 
-        return *this;
+        if (thisSize == otherSize) {
+            for (int i = 0; i < length_(); ++i)
+                data_[i] += other.data_[i];
+            return *this;
+        }
+        if (otherSize == 0) {
+            *this += other.value();
+            return *this;
+        }
+        if (thisSize == 0) {
+            assert(otherSize > 0);
+            this->appendDerivativesToConstant(otherSize);
+            data_[valuepos_()] += other.data_[valuepos_()];
+            for (int i = dstart_(); i < dend_(); ++i)
+                data_[i] = other.data_[i];
+            return *this;
+        }
+        throw std::logic_error(
+                "Cannot operate Evaluations with different number of derivatives "
+                "unless one of them has no derivatives"
+        );
     }
 
     // add value from other to this values
@@ -256,12 +289,32 @@ class Evaluation<ValueT, DynamicSize, staticSize>
     // subtract other's value and derivatives from this values
     OPM_HOST_DEVICE Evaluation& operator-=(const Evaluation& other)
     {
-        assert(size() == other.size());
+        const int thisSize = size();
+        const int otherSize = other.size();
 
-        for (int i = 0; i < length_(); ++i)
-            data_[i] -= other.data_[i];
+        // Allow operation if sizes match or one is constant (size == 0)
+        assert(thisSize == otherSize || thisSize == 0 || otherSize == 0);
 
-        return *this;
+        if (thisSize == otherSize) {
+            for (int i = 0; i < length_(); ++i)
+                data_[i] -= other.data_[i];
+            return *this;
+        }
+        if (otherSize == 0) {
+            *this -= other.value();
+            return *this;
+        }
+        if (thisSize == 0) {
+            assert(otherSize > 0);
+            this->appendDerivativesToConstant(otherSize);
+            for (int i = 0; i < other.length_(); ++i)
+                data_[i] -= other.data_[i];
+            return *this;
+        }
+        throw std::logic_error(
+                "Cannot operate Evaluations with different number of derivatives "
+                "unless one of them has no derivatives"
+        );
     }
 
     // subtract other's value from this values
@@ -277,21 +330,46 @@ class Evaluation<ValueT, DynamicSize, staticSize>
     // multiply values and apply chain rule to derivatives: (u*v)' = (v'u + u'v)
     OPM_HOST_DEVICE Evaluation& operator*=(const Evaluation& other)
     {
-        assert(size() == other.size());
+        const int thisSize = size();
+        const int otherSize = other.size();
 
-        // while the values are multiplied, the derivatives follow the product rule,
-        // i.e., (u*v)' = (v'u + u'v).
-        const ValueType u = this->value();
-        const ValueType v = other.value();
+        // Allow operation if sizes match or one is constant (size == 0)
+        assert(thisSize == otherSize || thisSize == 0 || otherSize == 0);
 
-        // value
-        data_[valuepos_()] *= v ;
+        if (thisSize == otherSize) {
+            // while the values are multiplied, the derivatives follow the product rule,
+            // i.e., (u*v)' = (v'u + u'v).
+            const ValueType u = this->value();
+            const ValueType v = other.value();
 
-        //  derivatives
-        for (int i = dstart_(); i < dend_(); ++i)
-            data_[i] = data_[i] * v + other.data_[i] * u;
+            // value
+            data_[valuepos_()] *= v ;
 
-        return *this;
+            //  derivatives
+            for (int i = dstart_(); i < dend_(); ++i)
+                data_[i] = data_[i] * v + other.data_[i] * u;
+
+            return *this;
+        }
+
+        if (otherSize == 0) {
+            *this *= other.value();
+            return *this;
+        }
+
+        if (thisSize == 0) {
+            assert(otherSize > 0);
+            this->appendDerivativesToConstant(otherSize);
+            const ValueType u = this->value();
+            for (int i = 0; i < length_(); ++i) {
+                data_[i] = u * other.data_[i];
+            }
+            return *this;
+        }
+        throw std::logic_error(
+                "Cannot operate Evaluations with different number of derivatives "
+                "unless one of them has no derivatives"
+        );
     }
 
     // m(c*u)' = c*u'
@@ -307,21 +385,58 @@ class Evaluation<ValueT, DynamicSize, staticSize>
     // m(u*v)' = (vu' - uv')/v^2
     OPM_HOST_DEVICE Evaluation& operator/=(const Evaluation& other)
     {
-        assert(size() == other.size());
+        const int thisSize = size();
+        const int otherSize = other.size();
+
+        // Allow operation if sizes match or one is constant (size == 0)
+        assert(thisSize == otherSize || thisSize == 0 || otherSize == 0);
 
-        // values are divided, derivatives follow the rule for division, i.e., (u/v)' = (v'u -
-        // u'v)/v^2.
-        ValueType& u = data_[valuepos_()];
-        const ValueType& v = other.value();
-        for (int idx = dstart_(); idx < dend_(); ++idx) {
-            const ValueType& uPrime = data_[idx];
-            const ValueType& vPrime = other.data_[idx];
+        if (thisSize == otherSize) {
+            // while the values are divided, the derivatives follow the rule for division,
+            // i.e., (u/v)' = (v'u - u'v)/v^2.
+            ValueType& u = data_[valuepos_()];
+            const ValueType& v = other.value();
 
-            data_[idx] = (v*uPrime - u*vPrime)/(v*v);
+            // derivatives
+            for (int idx = dstart_(); idx < dend_(); ++idx) {
+                const ValueType& uPrime = data_[idx];
+                const ValueType& vPrime = other.data_[idx];
+
+                data_[idx] = (v*uPrime - u*vPrime)/(v*v);
+            }
+
+            // value
+            u /= v;
+
+            return *this;
         }
-        u /= v;
 
-        return *this;
+        if (otherSize == 0) {
+            *this /= other.value();
+            return *this;
+        }
+
+        if (thisSize == 0) {
+            assert(otherSize > 0);
+            // if *this is a constant, extend to hold derivatives before division
+            this->appendDerivativesToConstant(otherSize);
+
+            // values are divided, derivatives follow the rule for division, i.e., (u/v)' = (v'u -
+            // u'v)/v^2.
+            ValueType& u = data_[valuepos_()];
+            const ValueType& v = other.value();
+            for (int idx = dstart_(); idx < dend_(); ++idx) {
+                const ValueType& vPrime = other.data_[idx];
+                data_[idx] =  -u * vPrime / (v * v);
+            }
+            u /= v;
+
+            return *this;
+        }
+        throw std::logic_error(
+                "Cannot operate Evaluations with different number of derivatives "
+                "unless one of them has no derivatives"
+        );
     }
 
     // divide value and derivatives by value of other
@@ -339,7 +454,11 @@ class Evaluation<ValueT, DynamicSize, staticSize>
     // add two evaluation objects
     OPM_HOST_DEVICE Evaluation operator+(const Evaluation& other) const
     {
-        assert(size() == other.size());
+        const int thisSize = size();
+        const int otherSize = other.size();
+
+        // Allow operation if sizes match or one is constant (size == 0)
+        assert(thisSize == otherSize || thisSize == 0 || otherSize == 0);
 
         Evaluation result(*this);
 
@@ -362,7 +481,11 @@ class Evaluation<ValueT, DynamicSize, staticSize>
     // subtract two evaluation objects
     OPM_HOST_DEVICE Evaluation operator-(const Evaluation& other) const
     {
-        assert(size() == other.size());
+        const int thisSize = size();
+        const int otherSize = other.size();
+
+        // Allow operation if sizes match or one is constant (size == 0)
+        assert(thisSize == otherSize || thisSize == 0 || otherSize == 0);
 
         Evaluation result(*this);
 
@@ -396,7 +519,11 @@ class Evaluation<ValueT, DynamicSize, staticSize>
 
     OPM_HOST_DEVICE Evaluation operator*(const Evaluation& other) const
     {
-        assert(size() == other.size());
+        const int thisSize = size();
+        const int otherSize = other.size();
+
+        // Allow operation if sizes match or one is constant (size == 0)
+        assert(thisSize == otherSize || thisSize == 0 || otherSize == 0);
 
         Evaluation result(*this);
 
@@ -417,7 +544,11 @@ class Evaluation<ValueT, DynamicSize, staticSize>
 
     OPM_HOST_DEVICE Evaluation operator/(const Evaluation& other) const
     {
-        assert(size() == other.size());
+        const int thisSize = size();
+        const int otherSize = other.size();
+
+        // Allow operation if sizes match or one is constant (size == 0)
+        assert(thisSize == otherSize || thisSize == 0 || otherSize == 0);
 
         Evaluation result(*this);
 
@@ -454,12 +585,20 @@ class Evaluation<ValueT, DynamicSize, staticSize>
 
     OPM_HOST_DEVICE bool operator==(const Evaluation& other) const
     {
-        assert(size() == other.size());
+        const int thisSize = size();
+        const int otherSize = other.size();
 
-        for (int idx = 0; idx < length_(); ++idx) {
-            if (data_[idx] != other.data_[idx]) {
-                return false;
+        // Allow operation if sizes match or one is constant (size == 0)
+        assert(thisSize == otherSize || thisSize == 0 || otherSize == 0);
+
+        if (thisSize == otherSize) {
+            for (int idx = 0; idx < length_(); ++idx) {
+                if (data_[idx] != other.data_[idx]) {
+                    return false;
+                }
             }
+        } else {
+            return value() == other.value();
         }
         return true;
     }
@@ -477,7 +616,7 @@ class Evaluation<ValueT, DynamicSize, staticSize>
 
     OPM_HOST_DEVICE bool operator>(const Evaluation& other) const
     {
-        assert(size() == other.size());
+        assert(size() == other.size() || size() == 0 || other.size() == 0);
 
         return value() > other.value();
     }
@@ -488,7 +627,7 @@ class Evaluation<ValueT, DynamicSize, staticSize>
 
     OPM_HOST_DEVICE bool operator<(const Evaluation& other) const
     {
-        assert(size() == other.size());
+        assert(size() == other.size() || size() == 0 || other.size() == 0);
 
         return value() < other.value();
     }
@@ -499,7 +638,7 @@ class Evaluation<ValueT, DynamicSize, staticSize>
 
     OPM_HOST_DEVICE bool operator>=(const Evaluation& other) const
     {
-        assert(size() == other.size());
+        assert(size() == other.size() || size() == 0 || other.size() == 0);
 
         return value() >= other.value();
     }
@@ -510,7 +649,7 @@ class Evaluation<ValueT, DynamicSize, staticSize>
 
     OPM_HOST_DEVICE bool operator<=(const Evaluation& other) const
     {
-        assert(size() == other.size());
+        assert(size() == other.size() || size() == 0 || other.size() == 0);
 
         return value() <= other.value();
     }
@@ -527,14 +666,33 @@ class Evaluation<ValueT, DynamicSize, staticSize>
     // return varIdx'th derivative
     OPM_HOST_DEVICE const ValueType& derivative(int varIdx) const
     {
-        assert(0 <= varIdx && varIdx < size());
+        assert(size() == 0 || (0 <= varIdx && varIdx < size()) );
+
+        if (size() == 0) {
+            static const ValueType zero {0.0};
+            return zero;
+        }
 
         return data_[dstart_() + varIdx];
     }
 
     // set derivative at position varIdx
-    OPM_HOST_DEVICE void setDerivative(int varIdx, const ValueType& derVal)
+    OPM_HOST_DEVICE void setDerivative(int varIdx, const ValueType& derVal, const int nVars = -1)
     {
+        // if size() == 0, we need nVars to be positive to extend the number of derivatives
+        // if size() > 0, nVars must be either negative (i.e. ignore) or match size()
+        assert( (size() == 0 && nVars > 0) ||
+                ((size() > 0 && (nVars < 0 || nVars == size()))) );
+
+        if (size() == 0) {
+            if (nVars < 0) {
+                throw std::logic_error("Cannot set derivative for a DynamicEvaluation initialized from a scalar "
+                                       "without specifying a positive number of derivatives");
+            }
+
+            this->appendDerivativesToConstant(nVars);
+        }
+
         assert(0 <= varIdx && varIdx < size());
 
         data_[dstart_() + varIdx] = derVal;
@@ -548,6 +706,13 @@ class Evaluation<ValueT, DynamicSize, staticSize>
 
 private:
     FastSmallVector<ValueT, staticSize> data_;
+
+    void appendDerivativesToConstant(size_t numDer) {
+        assert(size() == 0); // we only append derivatives to a constant
+        if (numDer > 0) {
+            data_.resize(1 + numDer);
+        }
+    }
 };
 
 template <class Scalar, unsigned staticSize = 0>
diff --git a/opm/material/densead/Evaluation.hpp b/opm/material/densead/Evaluation.hpp
index ca75448d2d8..dce3d197541 100644
--- a/opm/material/densead/Evaluation.hpp
+++ b/opm/material/densead/Evaluation.hpp
@@ -230,7 +230,7 @@ class Evaluation
     }
 
 
-    // add value and derivatives from other to this values and derivatives
+    // add value and derivatives from other to this value and derivatives
     OPM_HOST_DEVICE Evaluation& operator+=(const Evaluation& other)
     {
         assert(size() == other.size());
diff --git a/opm/material/densead/Evaluation1.hpp b/opm/material/densead/Evaluation1.hpp
index dc497b74314..90d8eceee3c 100644
--- a/opm/material/densead/Evaluation1.hpp
+++ b/opm/material/densead/Evaluation1.hpp
@@ -215,7 +215,7 @@ class Evaluation<ValueT, 1>
     }
 
 
-    // add value and derivatives from other to this values and derivatives
+    // add value and derivatives from other to this value and derivatives
     OPM_HOST_DEVICE Evaluation& operator+=(const Evaluation& other)
     {
         assert(size() == other.size());
diff --git a/opm/material/densead/Evaluation10.hpp b/opm/material/densead/Evaluation10.hpp
index a10bb0286af..b733dbdb22f 100644
--- a/opm/material/densead/Evaluation10.hpp
+++ b/opm/material/densead/Evaluation10.hpp
@@ -233,7 +233,7 @@ class Evaluation<ValueT, 10>
     }
 
 
-    // add value and derivatives from other to this values and derivatives
+    // add value and derivatives from other to this value and derivatives
     OPM_HOST_DEVICE Evaluation& operator+=(const Evaluation& other)
     {
         assert(size() == other.size());
diff --git a/opm/material/densead/Evaluation11.hpp b/opm/material/densead/Evaluation11.hpp
index a36278a7fd7..b564b475247 100644
--- a/opm/material/densead/Evaluation11.hpp
+++ b/opm/material/densead/Evaluation11.hpp
@@ -235,7 +235,7 @@ class Evaluation<ValueT, 11>
     }
 
 
-    // add value and derivatives from other to this values and derivatives
+    // add value and derivatives from other to this value and derivatives
     OPM_HOST_DEVICE Evaluation& operator+=(const Evaluation& other)
     {
         assert(size() == other.size());
diff --git a/opm/material/densead/Evaluation12.hpp b/opm/material/densead/Evaluation12.hpp
index 80de74e23cb..af9517af4ef 100644
--- a/opm/material/densead/Evaluation12.hpp
+++ b/opm/material/densead/Evaluation12.hpp
@@ -237,7 +237,7 @@ class Evaluation<ValueT, 12>
     }
 
 
-    // add value and derivatives from other to this values and derivatives
+    // add value and derivatives from other to this value and derivatives
     OPM_HOST_DEVICE Evaluation& operator+=(const Evaluation& other)
     {
         assert(size() == other.size());
diff --git a/opm/material/densead/Evaluation2.hpp b/opm/material/densead/Evaluation2.hpp
index 77dd7df49cd..d468b73c4d8 100644
--- a/opm/material/densead/Evaluation2.hpp
+++ b/opm/material/densead/Evaluation2.hpp
@@ -217,7 +217,7 @@ class Evaluation<ValueT, 2>
     }
 
 
-    // add value and derivatives from other to this values and derivatives
+    // add value and derivatives from other to this value and derivatives
     OPM_HOST_DEVICE Evaluation& operator+=(const Evaluation& other)
     {
         assert(size() == other.size());
diff --git a/opm/material/densead/Evaluation3.hpp b/opm/material/densead/Evaluation3.hpp
index dbb255f2ca0..9017127e1e8 100644
--- a/opm/material/densead/Evaluation3.hpp
+++ b/opm/material/densead/Evaluation3.hpp
@@ -219,7 +219,7 @@ class Evaluation<ValueT, 3>
     }
 
 
-    // add value and derivatives from other to this values and derivatives
+    // add value and derivatives from other to this value and derivatives
     OPM_HOST_DEVICE Evaluation& operator+=(const Evaluation& other)
     {
         assert(size() == other.size());
diff --git a/opm/material/densead/Evaluation4.hpp b/opm/material/densead/Evaluation4.hpp
index bd618343099..be9888e2323 100644
--- a/opm/material/densead/Evaluation4.hpp
+++ b/opm/material/densead/Evaluation4.hpp
@@ -221,7 +221,7 @@ class Evaluation<ValueT, 4>
     }
 
 
-    // add value and derivatives from other to this values and derivatives
+    // add value and derivatives from other to this value and derivatives
     OPM_HOST_DEVICE Evaluation& operator+=(const Evaluation& other)
     {
         assert(size() == other.size());
diff --git a/opm/material/densead/Evaluation5.hpp b/opm/material/densead/Evaluation5.hpp
index 9a27d8258fa..8b1b656e812 100644
--- a/opm/material/densead/Evaluation5.hpp
+++ b/opm/material/densead/Evaluation5.hpp
@@ -223,7 +223,7 @@ class Evaluation<ValueT, 5>
     }
 
 
-    // add value and derivatives from other to this values and derivatives
+    // add value and derivatives from other to this value and derivatives
     OPM_HOST_DEVICE Evaluation& operator+=(const Evaluation& other)
     {
         assert(size() == other.size());
diff --git a/opm/material/densead/Evaluation6.hpp b/opm/material/densead/Evaluation6.hpp
index d7c6a5aacbf..58b0441d175 100644
--- a/opm/material/densead/Evaluation6.hpp
+++ b/opm/material/densead/Evaluation6.hpp
@@ -225,7 +225,7 @@ class Evaluation<ValueT, 6>
     }
 
 
-    // add value and derivatives from other to this values and derivatives
+    // add value and derivatives from other to this value and derivatives
     OPM_HOST_DEVICE Evaluation& operator+=(const Evaluation& other)
     {
         assert(size() == other.size());
diff --git a/opm/material/densead/Evaluation7.hpp b/opm/material/densead/Evaluation7.hpp
index 7b7257bf007..720a8c90904 100644
--- a/opm/material/densead/Evaluation7.hpp
+++ b/opm/material/densead/Evaluation7.hpp
@@ -227,7 +227,7 @@ class Evaluation<ValueT, 7>
     }
 
 
-    // add value and derivatives from other to this values and derivatives
+    // add value and derivatives from other to this value and derivatives
     OPM_HOST_DEVICE Evaluation& operator+=(const Evaluation& other)
     {
         assert(size() == other.size());
diff --git a/opm/material/densead/Evaluation8.hpp b/opm/material/densead/Evaluation8.hpp
index 82849a21caf..140c979a196 100644
--- a/opm/material/densead/Evaluation8.hpp
+++ b/opm/material/densead/Evaluation8.hpp
@@ -229,7 +229,7 @@ class Evaluation<ValueT, 8>
     }
 
 
-    // add value and derivatives from other to this values and derivatives
+    // add value and derivatives from other to this value and derivatives
     OPM_HOST_DEVICE Evaluation& operator+=(const Evaluation& other)
     {
         assert(size() == other.size());
diff --git a/opm/material/densead/Evaluation9.hpp b/opm/material/densead/Evaluation9.hpp
index 66df07155cc..80db819d722 100644
--- a/opm/material/densead/Evaluation9.hpp
+++ b/opm/material/densead/Evaluation9.hpp
@@ -231,7 +231,7 @@ class Evaluation<ValueT, 9>
     }
 
 
-    // add value and derivatives from other to this values and derivatives
+    // add value and derivatives from other to this value and derivatives
     OPM_HOST_DEVICE Evaluation& operator+=(const Evaluation& other)
     {
         assert(size() == other.size());