diff --git a/opm/simulators/flow/BlackoilModelEbos.hpp b/opm/simulators/flow/BlackoilModelEbos.hpp
index 9c4511991..e26183ab5 100644
--- a/opm/simulators/flow/BlackoilModelEbos.hpp
+++ b/opm/simulators/flow/BlackoilModelEbos.hpp
@@ -1795,8 +1795,8 @@ namespace Opm {
 
             std::vector<int> domain_order(domains_.size());
             if (param_.local_solve_approach_ == "gauss-seidel") {
-                // TODO: enable flexibility and choice in choosing domain ordering approach.
-                if (true) {
+                switch (param_.local_domain_ordering_) {
+                case DomainOrderingMeasure::AveragePressure: {
                     // Use average pressures to order domains.
                     std::vector<std::pair<double, int>> avgpress_per_domain(domains_.size());
                     for (const auto& domain : domains_) {
@@ -1814,7 +1814,9 @@ namespace Opm {
                     for (size_t ii = 0; ii < domains_.size(); ++ii) {
                         domain_order[ii] = avgpress_per_domain[ii].second;
                     }
-                } else {
+                    break;
+                }
+                case DomainOrderingMeasure::Residual: {
                     // Use maximum residual to order domains.
                     const auto& residual = ebosSimulator().model().linearizer().residual();
                     const int num_vars = residual[0].size();
@@ -1836,6 +1838,7 @@ namespace Opm {
                         domain_order[ii] = maxres_per_domain[ii].second;
                     }
                 }
+                }
             } else {
                 std::iota(domain_order.begin(), domain_order.end(), 0);
             }
diff --git a/opm/simulators/flow/BlackoilModelParametersEbos.hpp b/opm/simulators/flow/BlackoilModelParametersEbos.hpp
index f52b62e71..ba4e3f94d 100644
--- a/opm/simulators/flow/BlackoilModelParametersEbos.hpp
+++ b/opm/simulators/flow/BlackoilModelParametersEbos.hpp
@@ -20,9 +20,15 @@
 #ifndef OPM_BLACKOILMODELPARAMETERS_EBOS_HEADER_INCLUDED
 #define OPM_BLACKOILMODELPARAMETERS_EBOS_HEADER_INCLUDED
 
-#include <opm/models/utils/propertysystem.hh>
-#include <opm/models/utils/parametersystem.hh>
+#include <opm/models/discretization/common/fvbaseproperties.hh>
 
+#include <opm/models/utils/basicproperties.hh>
+#include <opm/models/utils/parametersystem.hh>
+#include <opm/models/utils/propertysystem.hh>
+
+#include <opm/simulators/flow/SubDomain.hpp>
+
+#include <stdexcept>
 #include <string>
 
 namespace Opm::Properties {
@@ -209,6 +215,10 @@ template<class TypeTag, class MyTypeTag>
 struct LocalDomainsPartitioningMethod {
     using type = UndefinedProperty;
 };
+template<class TypeTag, class MyTypeTag>
+struct LocalDomainsOrderingMeasure {
+    using type = UndefinedProperty;
+};
 template<class TypeTag>
 struct DbhpMaxRel<TypeTag, TTag::FlowModelParameters> {
     using type = GetPropType<TypeTag, Scalar>;
@@ -401,6 +411,10 @@ template<class TypeTag>
 struct LocalDomainsPartitioningMethod<TypeTag, TTag::FlowModelParameters> {
     static constexpr auto value = "zoltan";
 };
+template<class TypeTag>
+struct LocalDomainsOrderingMeasure<TypeTag, TTag::FlowModelParameters> {
+    static constexpr auto value = "pressure";
+};
 // if openMP is available, determine the number threads per process automatically.
 #if _OPENMP
 template<class TypeTag>
@@ -535,6 +549,7 @@ namespace Opm
         int num_local_domains_{0};
         double local_domain_partition_imbalance_{1.03};
         std::string local_domain_partition_method_;
+        DomainOrderingMeasure local_domain_ordering_{DomainOrderingMeasure::AveragePressure};
 
         /// Construct from user parameters or defaults.
         BlackoilModelParametersEbos()
@@ -581,7 +596,15 @@ namespace Opm
             local_domain_partition_method_ = EWOMS_GET_PARAM(TypeTag, std::string, LocalDomainsPartitioningMethod);
             deck_file_name_ = EWOMS_GET_PARAM(TypeTag, std::string, EclDeckFileName);
             network_max_strict_iterations_ = EWOMS_GET_PARAM(TypeTag, int, NetworkMaxStrictIterations);
-            network_max_iterations_ = EWOMS_GET_PARAM(TypeTag, int, NetworkMaxIterations);            
+            network_max_iterations_ = EWOMS_GET_PARAM(TypeTag, int, NetworkMaxIterations);
+            std::string measure = EWOMS_GET_PARAM(TypeTag, std::string, LocalDomainsOrderingMeasure);
+            if (measure == "residual") {
+                local_domain_ordering_ = DomainOrderingMeasure::Residual;
+            } else if (measure == "pressure") {
+                local_domain_ordering_ = DomainOrderingMeasure::AveragePressure;
+            } else {
+                throw std::runtime_error("Invalid domain ordering '" + measure + "' specified.");
+            }
         }
 
         static void registerParameters()
@@ -631,6 +654,8 @@ namespace Opm
             EWOMS_REGISTER_PARAM(TypeTag, Scalar, LocalDomainsPartitioningImbalance, "Subdomain partitioning imbalance tolerance. 1.03 is 3 percent imbalance.");
             EWOMS_REGISTER_PARAM(TypeTag, std::string, LocalDomainsPartitioningMethod, "Subdomain partitioning method. "
                                  "Allowed values are 'zoltan', 'simple', and the name of a partition file ending with '.partition'.");
+            EWOMS_REGISTER_PARAM(TypeTag, std::string, LocalDomainsOrderingMeasure, "Subdomain ordering measure. "
+                                 "Allowed values are 'pressure' and  'residual'.");
         }
     };
 } // namespace Opm
diff --git a/opm/simulators/flow/SubDomain.hpp b/opm/simulators/flow/SubDomain.hpp
index 7796ff07e..e79297788 100644
--- a/opm/simulators/flow/SubDomain.hpp
+++ b/opm/simulators/flow/SubDomain.hpp
@@ -28,6 +28,12 @@
 namespace Opm
 {
 
+    //! \brief Measure to use for domain ordering.
+    enum class DomainOrderingMeasure {
+        AveragePressure,
+        Residual
+    };
+
     /// Representing a part of a grid, in a way suitable for performing
     /// local solves.
     template <class Grid>