diff --git a/include/utils/int_range.h b/include/utils/int_range.h
index b967afa3736..de25e823c52 100644
--- a/include/utils/int_range.h
+++ b/include/utils/int_range.h
@@ -61,7 +61,7 @@ class IntRange
     using pointer = T;
     using reference = T&;
 
-    iterator (T i) : _i(i) {}
+    explicit iterator (T i) : _i(i) {}
 
     T operator* () const { return _i; }
 
@@ -78,6 +78,20 @@ class IntRange
       return returnval;
     }
 
+    template <typename T2>
+    void operator+= (T2 n)
+    {
+      _i += cast_int<T>(n);
+    }
+
+    template <typename T2>
+    iterator operator+ (T2 n) const
+    {
+      iterator returnval(_i);
+      returnval += n;
+      return returnval;
+    }
+
     bool operator== (const iterator & j) const
     {
       return ( _i == j._i );
@@ -92,21 +106,40 @@ class IntRange
     T _i;
   };
 
+  typedef typename IntRange<T>::iterator const_iterator;
+
   template <typename U, typename V>
   IntRange(U begin, V end) :
     _begin(cast_int<T>(begin)),
     _end(cast_int<T>(end))
   {}
 
+  // Signature needed for use in threads_pthread.h
+  IntRange(const IntRange &,
+           const const_iterator & begin,
+           const const_iterator & end) :
+    _begin(begin),
+    _end(end)
+  {}
+
   iterator begin() const { return _begin; }
 
   iterator end () const { return _end; }
 
+  std::size_t size () const
+  { return *_end - *_begin; }
+
 private:
   iterator _begin, _end;
 };
 
 
+template <typename T, typename T2>
+typename IntRange<T>::iterator operator+ (T2 n, typename IntRange<T>::iterator i)
+{
+  return i + n;
+}
+
 
 /**
  * Helper function that returns an IntRange<std::size_t> representing
diff --git a/src/numerics/numeric_vector.C b/src/numerics/numeric_vector.C
index a66cda1ec68..36e68236f40 100644
--- a/src/numerics/numeric_vector.C
+++ b/src/numerics/numeric_vector.C
@@ -617,13 +617,34 @@ Real NumericVector<T>::l2_norm_diff (const NumericVector<T> & v) const
 {
   libmesh_assert(this->compatible(v));
 
-  Real norm = 0;
-  for (const auto i : make_range(this->first_local_index(), this->last_local_index()))
-    norm += TensorTools::norm_sq((*this)(i) - v(i));
+  struct Differ {
+    const NumericVector<T> & v1, & v2;
+    Real norm_sq;
 
-  this->comm().sum(norm);
+    Differ (const NumericVector<T> & v1in,
+            const NumericVector<T> & v2in) :
+      v1(v1in), v2(v2in), norm_sq(0) {}
 
-  return std::sqrt(norm);
+    Differ (Differ & other, Threads::split) :
+      v1(other.v1), v2(other.v2), norm_sq(0) {}
+
+    void operator()(const IntRange<numeric_index_type> & index_range) {
+      for (const auto i : index_range)
+        norm_sq += TensorTools::norm_sq(v1(i) - v2(i));
+    }
+
+    void join(const Differ & other) {
+      norm_sq += other.norm_sq;
+    }
+  };
+
+  Differ differ(*this, v);
+
+  Threads::parallel_reduce(make_range(this->first_local_index(), this->last_local_index()), differ);
+
+  this->comm().sum(differ.norm_sq);
+
+  return std::sqrt(differ.norm_sq);
 }