diff --git a/CHANGELOG.md b/CHANGELOG.md
index 28d3dcc6e..164a07279 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -5,6 +5,7 @@
 ### Added
 
 - Recommendation on how to cite in publications (#943)
+- Store distance matrices (#956)
 
 ### Changed
 
diff --git a/docs/API.md b/docs/API.md
index 0c5f0ebdd..dcbe0c9ee 100644
--- a/docs/API.md
+++ b/docs/API.md
@@ -281,11 +281,12 @@ integers filed under the `profile` key, then under:
 
 - `durations` for a custom travel-time matrix that will be used for
   all checks against timing constraints;
+- `distances` for a custom distance matrix;
 - `costs` for a custom cost matrix that will be used within all route
   cost evaluations.
 
-If only the `durations` value is provided, it's implied that it should
-also be used for costs evaluations.
+If only the `durations` matrix is provided, internal costs are derived from
+durations based on vehicles `costs` properties.
 
 Example of describing different matrices for different vehicle
 profiles:
diff --git a/src/routing/http_wrapper.cpp b/src/routing/http_wrapper.cpp
index d79798e4d..99a87f06a 100644
--- a/src/routing/http_wrapper.cpp
+++ b/src/routing/http_wrapper.cpp
@@ -22,14 +22,16 @@ HttpWrapper::HttpWrapper(const std::string& profile,
                          Server server,
                          std::string matrix_service,
                          std::string matrix_durations_key,
+                         std::string matrix_distances_key,
                          std::string route_service,
-                         std::string extra_args)
+                         std::string routing_args)
   : Wrapper(profile),
     _server(std::move(server)),
     _matrix_service(std::move(matrix_service)),
     _matrix_durations_key(std::move(matrix_durations_key)),
+    _matrix_distances_key(std::move(matrix_distances_key)),
     _route_service(std::move(route_service)),
-    _extra_args(std::move(extra_args)) {
+    _routing_args(std::move(routing_args)) {
 }
 
 std::string HttpWrapper::send_then_receive(const std::string& query) const {
@@ -144,8 +146,7 @@ void HttpWrapper::parse_response(rapidjson::Document& json_result,
 #endif
 }
 
-Matrix<UserCost>
-HttpWrapper::get_matrix(const std::vector<Location>& locs) const {
+Matrices HttpWrapper::get_matrices(const std::vector<Location>& locs) const {
   std::string query = this->build_query(locs, _matrix_service);
   std::string json_string = this->run_query(query);
 
@@ -161,25 +162,34 @@ HttpWrapper::get_matrix(const std::vector<Location>& locs) const {
   }
   assert(json_result[_matrix_durations_key.c_str()].Size() == m_size);
 
-  // Build matrix while checking for unfound routes ('null' values) to
-  // avoid unexpected behavior.
-  Matrix<UserCost> m(m_size);
+  if (!json_result.HasMember(_matrix_distances_key.c_str())) {
+    throw RoutingException("Missing " + _matrix_distances_key + ".");
+  }
+  assert(json_result[_matrix_distances_key.c_str()].Size() == m_size);
+
+  // Build matrices while checking for unfound routes ('null' values)
+  // to avoid unexpected behavior.
+  Matrices m(m_size);
 
   std::vector<unsigned> nb_unfound_from_loc(m_size, 0);
   std::vector<unsigned> nb_unfound_to_loc(m_size, 0);
 
   for (rapidjson::SizeType i = 0; i < m_size; ++i) {
-    const auto& line = json_result[_matrix_durations_key.c_str()][i];
-    assert(line.Size() == m_size);
-    for (rapidjson::SizeType j = 0; j < line.Size(); ++j) {
-      if (duration_value_is_null(line[j])) {
+    const auto& duration_line = json_result[_matrix_durations_key.c_str()][i];
+    const auto& distance_line = json_result[_matrix_distances_key.c_str()][i];
+    assert(duration_line.Size() == m_size);
+    assert(distance_line.Size() == m_size);
+    for (rapidjson::SizeType j = 0; j < m_size; ++j) {
+      if (duration_value_is_null(duration_line[j]) or
+          distance_value_is_null(distance_line[j])) {
         // No route found between i and j. Just storing info as we
         // don't know yet which location is responsible between i
         // and j.
         ++nb_unfound_from_loc[i];
         ++nb_unfound_to_loc[j];
       } else {
-        m[i][j] = get_duration_value(line[j]);
+        m.durations[i][j] = get_duration_value(duration_line[j]);
+        m.distances[i][j] = get_distance_value(distance_line[j]);
       }
     }
   }
@@ -208,8 +218,7 @@ void HttpWrapper::add_route_info(Route& route) const {
   }
   assert(!non_break_locations.empty());
 
-  std::string query =
-    build_query(non_break_locations, _route_service, _extra_args);
+  std::string query = build_query(non_break_locations, _route_service);
 
   std::string json_string = this->run_query(query);
 
@@ -218,7 +227,7 @@ void HttpWrapper::add_route_info(Route& route) const {
   this->check_response(json_result, _route_service);
 
   // Total distance and route geometry.
-  route.distance = round_cost(get_total_distance(json_result));
+  route.distance = round_cost<UserDistance>(get_total_distance(json_result));
   route.geometry = get_geometry(json_result);
 
   auto nb_legs = get_legs_number(json_result);
@@ -245,17 +254,17 @@ void HttpWrapper::add_route_info(Route& route) const {
     for (unsigned b = 1; b <= number_breaks_after[i]; ++b) {
       auto& break_step = route.steps[steps_rank + b];
       if (next_duration == 0) {
-        break_step.distance = round_cost(sum_distance);
+        break_step.distance = round_cost<UserDistance>(sum_distance);
       } else {
-        break_step.distance =
-          round_cost(sum_distance +
-                     ((break_step.duration - step.duration) * next_distance) /
-                       next_duration);
+        break_step.distance = round_cost<UserDistance>(
+          sum_distance +
+          ((break_step.duration - step.duration) * next_distance) /
+            next_duration);
       }
     }
 
     sum_distance += next_distance;
-    next_step.distance = round_cost(sum_distance);
+    next_step.distance = round_cost<UserDistance>(sum_distance);
 
     steps_rank += number_breaks_after[i] + 1;
   }
diff --git a/src/routing/http_wrapper.h b/src/routing/http_wrapper.h
index 183cddcca..f6177ca95 100644
--- a/src/routing/http_wrapper.h
+++ b/src/routing/http_wrapper.h
@@ -28,15 +28,17 @@ class HttpWrapper : public Wrapper {
   const Server _server;
   const std::string _matrix_service;
   const std::string _matrix_durations_key;
+  const std::string _matrix_distances_key;
   const std::string _route_service;
-  const std::string _extra_args;
+  const std::string _routing_args;
 
   HttpWrapper(const std::string& profile,
               Server server,
               std::string matrix_service,
               std::string matrix_durations_key,
+              std::string matrix_distances_key,
               std::string route_service,
-              std::string extra_args);
+              std::string routing_args);
 
   std::string run_query(const std::string& query) const;
 
@@ -44,20 +46,25 @@ class HttpWrapper : public Wrapper {
                              const std::string& json_content);
 
   virtual std::string build_query(const std::vector<Location>& locations,
-                                  const std::string& service,
-                                  const std::string& extra_args = "") const = 0;
+                                  const std::string& service) const = 0;
 
   virtual void check_response(const rapidjson::Document& json_result,
                               const std::string& service) const = 0;
 
-  Matrix<UserCost> get_matrix(const std::vector<Location>& locs) const override;
+  Matrices get_matrices(const std::vector<Location>& locs) const override;
 
   virtual bool
   duration_value_is_null(const rapidjson::Value& matrix_entry) const = 0;
 
+  virtual bool
+  distance_value_is_null(const rapidjson::Value& matrix_entry) const = 0;
+
   virtual UserDuration
   get_duration_value(const rapidjson::Value& matrix_entry) const = 0;
 
+  virtual UserDistance
+  get_distance_value(const rapidjson::Value& matrix_entry) const = 0;
+
   virtual double get_total_distance(const rapidjson::Value& result) const = 0;
 
   virtual unsigned get_legs_number(const rapidjson::Value& result) const = 0;
diff --git a/src/routing/libosrm_wrapper.cpp b/src/routing/libosrm_wrapper.cpp
index cb8bbf551..0a0669024 100644
--- a/src/routing/libosrm_wrapper.cpp
+++ b/src/routing/libosrm_wrapper.cpp
@@ -32,9 +32,10 @@ LibosrmWrapper::LibosrmWrapper(const std::string& profile)
   : Wrapper(profile), _config(get_config(profile)), _osrm(_config) {
 }
 
-Matrix<UserCost>
-LibosrmWrapper::get_matrix(const std::vector<Location>& locs) const {
+Matrices LibosrmWrapper::get_matrices(const std::vector<Location>& locs) const {
   osrm::TableParameters params;
+  params.annotations = osrm::engine::api::TableParameters::AnnotationsType::All;
+
   for (auto const& location : locs) {
     assert(location.has_coordinates());
     params.coordinates
@@ -51,34 +52,43 @@ LibosrmWrapper::get_matrix(const std::vector<Location>& locs) const {
       ": " + result.values["message"].get<osrm::json::String>().value);
   }
 
-  auto& table = result.values["durations"].get<osrm::json::Array>();
+  const auto& durations = result.values["durations"].get<osrm::json::Array>();
+  const auto& distances = result.values["distances"].get<osrm::json::Array>();
 
   // Expected matrix size.
   std::size_t m_size = locs.size();
-  assert(table.values.size() == m_size);
+  assert(durations.values.size() == m_size);
+  assert(distances.values.size() == m_size);
 
   // Build matrix while checking for unfound routes to avoid
   // unexpected behavior (OSRM raises 'null').
-  Matrix<UserCost> m(m_size);
+  Matrices m(m_size);
 
   std::vector<unsigned> nb_unfound_from_loc(m_size, 0);
   std::vector<unsigned> nb_unfound_to_loc(m_size, 0);
 
   std::string reason;
   for (std::size_t i = 0; i < m_size; ++i) {
-    const auto& line = table.values.at(i).get<osrm::json::Array>();
-    assert(line.values.size() == m_size);
+    const auto& duration_line = durations.values.at(i).get<osrm::json::Array>();
+    const auto& distance_line = distances.values.at(i).get<osrm::json::Array>();
+    assert(duration_line.values.size() == m_size);
+    assert(distance_line.values.size() == m_size);
+
     for (std::size_t j = 0; j < m_size; ++j) {
-      const auto& el = line.values.at(j);
-      if (el.is<osrm::json::Null>()) {
+      const auto& duration_el = duration_line.values.at(j);
+      const auto& distance_el = distance_line.values.at(j);
+      if (duration_el.is<osrm::json::Null>() or
+          distance_el.is<osrm::json::Null>()) {
         // No route found between i and j. Just storing info as we
         // don't know yet which location is responsible between i
         // and j.
         ++nb_unfound_from_loc[i];
         ++nb_unfound_to_loc[j];
       } else {
-        auto cost = round_cost(el.get<osrm::json::Number>().value);
-        m[i][j] = cost;
+        m.durations[i][j] =
+          round_cost<UserDuration>(duration_el.get<osrm::json::Number>().value);
+        m.distances[i][j] =
+          round_cost<UserDistance>(distance_el.get<osrm::json::Number>().value);
       }
     }
   }
@@ -130,8 +140,8 @@ void LibosrmWrapper::add_route_info(Route& route) const {
   auto& json_route = result_routes.values.at(0).get<osrm::json::Object>();
 
   // Total distance and route geometry.
-  route.distance =
-    round_cost(json_route.values["distance"].get<osrm::json::Number>().value);
+  route.distance = round_cost<UserDistance>(
+    json_route.values["distance"].get<osrm::json::Number>().value);
   route.geometry =
     std::move(json_route.values["geometry"].get<osrm::json::String>().value);
 
@@ -170,13 +180,13 @@ void LibosrmWrapper::add_route_info(Route& route) const {
     // non-breaks steps.
     for (unsigned b = 1; b <= number_breaks_after[i]; ++b) {
       auto& break_step = route.steps[steps_rank + b];
-      break_step.distance = round_cost(
+      break_step.distance = round_cost<UserDistance>(
         sum_distance + ((break_step.duration - step.duration) * next_distance) /
                          next_duration);
     }
 
     sum_distance += next_distance;
-    next_step.distance = round_cost(sum_distance);
+    next_step.distance = round_cost<UserDistance>(sum_distance);
 
     steps_rank += number_breaks_after[i] + 1;
   }
diff --git a/src/routing/libosrm_wrapper.h b/src/routing/libosrm_wrapper.h
index d0797064d..8c031c4d0 100644
--- a/src/routing/libosrm_wrapper.h
+++ b/src/routing/libosrm_wrapper.h
@@ -29,8 +29,8 @@ class LibosrmWrapper : public Wrapper {
 public:
   LibosrmWrapper(const std::string& profile);
 
-  virtual Matrix<UserCost>
-  get_matrix(const std::vector<Location>& locs) const override;
+  virtual Matrices
+  get_matrices(const std::vector<Location>& locs) const override;
 
   virtual void add_route_info(Route& route) const override;
 };
diff --git a/src/routing/ors_wrapper.cpp b/src/routing/ors_wrapper.cpp
index 16282890e..e0f653d10 100644
--- a/src/routing/ors_wrapper.cpp
+++ b/src/routing/ors_wrapper.cpp
@@ -16,14 +16,14 @@ OrsWrapper::OrsWrapper(const std::string& profile, const Server& server)
                 server,
                 "matrix",
                 "durations",
+                "distances",
                 "directions",
                 "\"geometry_simplify\":\"false\",\"continue_straight\":"
                 "\"false\"") {
 }
 
 std::string OrsWrapper::build_query(const std::vector<Location>& locations,
-                                    const std::string& service,
-                                    const std::string& extra_args) const {
+                                    const std::string& service) const {
   // Adding locations.
   std::string body = "{\"";
   if (service == "directions") {
@@ -38,8 +38,11 @@ std::string OrsWrapper::build_query(const std::vector<Location>& locations,
   }
   body.pop_back(); // Remove trailing ','.
   body += "]";
-  if (!extra_args.empty()) {
-    body += "," + extra_args;
+  if (service == _route_service) {
+    body += "," + _routing_args;
+  } else {
+    assert(service == _matrix_service);
+    body += ",\"metrics\":[\"duration\",\"distance\"]";
   }
   body += "}";
 
@@ -70,9 +73,19 @@ bool OrsWrapper::duration_value_is_null(
   return matrix_entry.IsNull();
 }
 
+bool OrsWrapper::distance_value_is_null(
+  const rapidjson::Value& matrix_entry) const {
+  return matrix_entry.IsNull();
+}
+
 UserDuration
 OrsWrapper::get_duration_value(const rapidjson::Value& matrix_entry) const {
-  return round_cost(matrix_entry.GetDouble());
+  return round_cost<UserDuration>(matrix_entry.GetDouble());
+}
+
+UserDistance
+OrsWrapper::get_distance_value(const rapidjson::Value& matrix_entry) const {
+  return round_cost<UserDistance>(matrix_entry.GetDouble());
 }
 
 double OrsWrapper::get_total_distance(const rapidjson::Value& result) const {
diff --git a/src/routing/ors_wrapper.h b/src/routing/ors_wrapper.h
index 923072eca..ae6656e8a 100644
--- a/src/routing/ors_wrapper.h
+++ b/src/routing/ors_wrapper.h
@@ -17,8 +17,7 @@ namespace vroom::routing {
 class OrsWrapper : public HttpWrapper {
 private:
   std::string build_query(const std::vector<Location>& locations,
-                          const std::string& service,
-                          const std::string& extra_args) const override;
+                          const std::string& service) const override;
 
   void check_response(const rapidjson::Document& json_result,
                       const std::string& service) const override;
@@ -26,9 +25,15 @@ class OrsWrapper : public HttpWrapper {
   bool
   duration_value_is_null(const rapidjson::Value& matrix_entry) const override;
 
+  bool
+  distance_value_is_null(const rapidjson::Value& matrix_entry) const override;
+
   UserDuration
   get_duration_value(const rapidjson::Value& matrix_entry) const override;
 
+  UserDistance
+  get_distance_value(const rapidjson::Value& matrix_entry) const override;
+
   double get_total_distance(const rapidjson::Value& result) const override;
 
   unsigned get_legs_number(const rapidjson::Value& result) const override;
diff --git a/src/routing/osrm_routed_wrapper.cpp b/src/routing/osrm_routed_wrapper.cpp
index b20fd2ca3..4dc0a08f9 100644
--- a/src/routing/osrm_routed_wrapper.cpp
+++ b/src/routing/osrm_routed_wrapper.cpp
@@ -17,6 +17,7 @@ OsrmRoutedWrapper::OsrmRoutedWrapper(const std::string& profile,
                 server,
                 "table",
                 "durations",
+                "distances",
                 "route",
                 "alternatives=false&steps=false&overview=full&continue_"
                 "straight=false") {
@@ -24,8 +25,7 @@ OsrmRoutedWrapper::OsrmRoutedWrapper(const std::string& profile,
 
 std::string
 OsrmRoutedWrapper::build_query(const std::vector<Location>& locations,
-                               const std::string& service,
-                               const std::string& extra_args) const {
+                               const std::string& service) const {
   // Building query for osrm-routed
   std::string query = "GET /" + service;
 
@@ -38,8 +38,11 @@ OsrmRoutedWrapper::build_query(const std::vector<Location>& locations,
   }
   query.pop_back(); // Remove trailing ';'.
 
-  if (!extra_args.empty()) {
-    query += "?" + extra_args;
+  if (service == _route_service) {
+    query += "?" + _routing_args;
+  } else {
+    assert(service == _matrix_service);
+    query += "?annotations=duration,distance";
   }
 
   query += " HTTP/1.1\r\n";
@@ -63,9 +66,19 @@ bool OsrmRoutedWrapper::duration_value_is_null(
   return matrix_entry.IsNull();
 }
 
+bool OsrmRoutedWrapper::distance_value_is_null(
+  const rapidjson::Value& matrix_entry) const {
+  return matrix_entry.IsNull();
+}
+
 UserDuration OsrmRoutedWrapper::get_duration_value(
   const rapidjson::Value& matrix_entry) const {
-  return round_cost(matrix_entry.GetDouble());
+  return round_cost<UserDuration>(matrix_entry.GetDouble());
+}
+
+UserDistance OsrmRoutedWrapper::get_distance_value(
+  const rapidjson::Value& matrix_entry) const {
+  return round_cost<UserDistance>(matrix_entry.GetDouble());
 }
 
 double
diff --git a/src/routing/osrm_routed_wrapper.h b/src/routing/osrm_routed_wrapper.h
index 86f33795a..f57cd159a 100644
--- a/src/routing/osrm_routed_wrapper.h
+++ b/src/routing/osrm_routed_wrapper.h
@@ -17,8 +17,7 @@ namespace vroom::routing {
 class OsrmRoutedWrapper : public HttpWrapper {
 private:
   std::string build_query(const std::vector<Location>& locations,
-                          const std::string& service,
-                          const std::string& extra_args) const override;
+                          const std::string& service) const override;
 
   void check_response(const rapidjson::Document& json_result,
                       const std::string& service) const override;
@@ -26,9 +25,15 @@ class OsrmRoutedWrapper : public HttpWrapper {
   bool
   duration_value_is_null(const rapidjson::Value& matrix_entry) const override;
 
+  bool
+  distance_value_is_null(const rapidjson::Value& matrix_entry) const override;
+
   UserDuration
   get_duration_value(const rapidjson::Value& matrix_entry) const override;
 
+  UserDistance
+  get_distance_value(const rapidjson::Value& matrix_entry) const override;
+
   double get_total_distance(const rapidjson::Value& result) const override;
 
   unsigned get_legs_number(const rapidjson::Value& result) const override;
diff --git a/src/routing/valhalla_wrapper.cpp b/src/routing/valhalla_wrapper.cpp
index a82977e96..335820375 100644
--- a/src/routing/valhalla_wrapper.cpp
+++ b/src/routing/valhalla_wrapper.cpp
@@ -23,6 +23,7 @@ ValhallaWrapper::ValhallaWrapper(const std::string& profile,
                 server,
                 "sources_to_targets",
                 "sources_to_targets",
+                "sources_to_targets",
                 "route",
                 R"("directions_type":"none")") {
 }
@@ -53,8 +54,7 @@ std::string ValhallaWrapper::get_matrix_query(
 }
 
 std::string
-ValhallaWrapper::get_route_query(const std::vector<Location>& locations,
-                                 const std::string& extra_args) const {
+ValhallaWrapper::get_route_query(const std::vector<Location>& locations) const {
   // Building matrix query for Valhalla.
   std::string query = "GET /" + _route_service + "?json={\"locations\":[";
 
@@ -66,9 +66,7 @@ ValhallaWrapper::get_route_query(const std::vector<Location>& locations,
   query.pop_back(); // Remove trailing ','.
 
   query += R"(],"costing":")" + profile + "\"";
-  if (!extra_args.empty()) {
-    query += "," + extra_args;
-  }
+  query += "," + _routing_args;
   query += "}";
 
   query += " HTTP/1.1\r\n";
@@ -80,12 +78,11 @@ ValhallaWrapper::get_route_query(const std::vector<Location>& locations,
 }
 
 std::string ValhallaWrapper::build_query(const std::vector<Location>& locations,
-                                         const std::string& service,
-                                         const std::string& extra_args) const {
+                                         const std::string& service) const {
   assert(service == _matrix_service or service == _route_service);
 
   return (service == _matrix_service) ? get_matrix_query(locations)
-                                      : get_route_query(locations, extra_args);
+                                      : get_route_query(locations);
 }
 
 void ValhallaWrapper::check_response(const rapidjson::Document& json_result,
@@ -127,12 +124,25 @@ bool ValhallaWrapper::duration_value_is_null(
   return matrix_entry["time"].IsNull();
 }
 
+bool ValhallaWrapper::distance_value_is_null(
+  const rapidjson::Value& matrix_entry) const {
+  assert(matrix_entry.HasMember("distance"));
+  return matrix_entry["distance"].IsNull();
+}
+
 UserDuration ValhallaWrapper::get_duration_value(
   const rapidjson::Value& matrix_entry) const {
   assert(matrix_entry["time"].IsUint());
   return matrix_entry["time"].GetUint();
 }
 
+UserDistance ValhallaWrapper::get_distance_value(
+  const rapidjson::Value& matrix_entry) const {
+  assert(matrix_entry["distance"].IsDouble());
+  return round_cost<UserDistance>(km_to_m *
+                                  matrix_entry["distance"].GetDouble());
+}
+
 double
 ValhallaWrapper::get_total_distance(const rapidjson::Value& result) const {
   return km_to_m * result["trip"]["summary"]["length"].GetDouble();
diff --git a/src/routing/valhalla_wrapper.h b/src/routing/valhalla_wrapper.h
index 186347e46..4644b072b 100644
--- a/src/routing/valhalla_wrapper.h
+++ b/src/routing/valhalla_wrapper.h
@@ -18,12 +18,10 @@ class ValhallaWrapper : public HttpWrapper {
 private:
   std::string get_matrix_query(const std::vector<Location>& locations) const;
 
-  std::string get_route_query(const std::vector<Location>& locations,
-                              const std::string& extra_args = "") const;
+  std::string get_route_query(const std::vector<Location>& locations) const;
 
   std::string build_query(const std::vector<Location>& locations,
-                          const std::string& service,
-                          const std::string& extra_args) const override;
+                          const std::string& service) const override;
 
   void check_response(const rapidjson::Document& json_result,
                       const std::string& service) const override;
@@ -31,9 +29,15 @@ class ValhallaWrapper : public HttpWrapper {
   bool
   duration_value_is_null(const rapidjson::Value& matrix_entry) const override;
 
+  bool
+  distance_value_is_null(const rapidjson::Value& matrix_entry) const override;
+
   UserDuration
   get_duration_value(const rapidjson::Value& matrix_entry) const override;
 
+  UserDistance
+  get_distance_value(const rapidjson::Value& matrix_entry) const override;
+
   double get_total_distance(const rapidjson::Value& result) const override;
 
   unsigned get_legs_number(const rapidjson::Value& result) const override;
diff --git a/src/routing/wrapper.h b/src/routing/wrapper.h
index 754d973c2..1024a9b81 100644
--- a/src/routing/wrapper.h
+++ b/src/routing/wrapper.h
@@ -19,13 +19,19 @@ All rights reserved (see LICENSE).
 
 namespace vroom::routing {
 
+struct Matrices {
+  Matrix<UserDuration> durations;
+  Matrix<UserDistance> distances;
+
+  Matrices(std::size_t n) : durations(n), distances(n){};
+};
+
 class Wrapper {
 
 public:
   std::string profile;
 
-  virtual Matrix<UserCost>
-  get_matrix(const std::vector<Location>& locs) const = 0;
+  virtual Matrices get_matrices(const std::vector<Location>& locs) const = 0;
 
   virtual void add_route_info(Route& route) const = 0;
 
@@ -35,9 +41,9 @@ class Wrapper {
   Wrapper(std::string profile) : profile(std::move(profile)) {
   }
 
-  static UserCost round_cost(double value) {
+  template <typename T> static T round_cost(double value) {
     constexpr double round_increment = 0.5;
-    return static_cast<UserCost>(value + round_increment);
+    return static_cast<T>(value + round_increment);
   }
 
   static inline void
diff --git a/src/structures/typedefs.h b/src/structures/typedefs.h
index 002b2d6f8..15e0e3e5d 100644
--- a/src/structures/typedefs.h
+++ b/src/structures/typedefs.h
@@ -32,9 +32,9 @@ using Id = uint64_t;
 using Index = uint16_t;
 using UserCost = uint32_t;
 using Cost = int64_t;
-using Distance = uint32_t;
 using UserDuration = uint32_t;
 using Duration = int64_t;
+using UserDistance = uint32_t;
 using Coordinate = double;
 using Capacity = int64_t;
 using Skill = uint32_t;
diff --git a/src/structures/vroom/input/input.cpp b/src/structures/vroom/input/input.cpp
index 483647995..bcfda1b25 100644
--- a/src/structures/vroom/input/input.cpp
+++ b/src/structures/vroom/input/input.cpp
@@ -386,6 +386,14 @@ void Input::set_durations_matrix(const std::string& profile,
   _durations_matrices.insert_or_assign(profile, m);
 }
 
+void Input::set_distances_matrix(const std::string& profile,
+                                 Matrix<UserDistance>&& m) {
+  if (m.size() == 0) {
+    throw InputException("Empty distances matrix for " + profile + " profile.");
+  }
+  _distances_matrices.insert_or_assign(profile, m);
+}
+
 void Input::set_costs_matrix(const std::string& profile, Matrix<UserCost>&& m) {
   if (m.size() == 0) {
     throw InputException("Empty costs matrix for " + profile + " profile.");
@@ -808,19 +816,27 @@ void Input::set_matrices(unsigned nb_thread) {
   for (const auto& profile : _profiles) {
     thread_profiles[t_rank % nb_buckets].push_back(profile);
     ++t_rank;
+
+    // Even with custom matrices, we still need routing after
+    // optimization if geometry is requested.
+    bool create_routing_wrapper = _geometry;
+
     if (_durations_matrices.find(profile) == _durations_matrices.end()) {
-      // Durations matrix has not been manually set, create routing
-      // wrapper and empty matrix to allow for concurrent modification
-      // later on.
-      add_routing_wrapper(profile);
+      // Durations matrix has not been manually set, create empty
+      // matrix to allow for concurrent modification later on.
+      create_routing_wrapper = true;
       _durations_matrices.emplace(profile, Matrix<UserDuration>());
-    } else {
-      if (_geometry) {
-        // Even with a custom matrix, we still want routing after
-        // optimization.
-        add_routing_wrapper(profile);
+
+      if (_distances_matrices.find(profile) == _distances_matrices.end()) {
+        // Distances matrix has not been manually set, create empty
+        // matrix to allow for concurrent modification later on.
+        _distances_matrices.emplace(profile, Matrix<UserDistance>());
       }
     }
+
+    if (create_routing_wrapper) {
+      add_routing_wrapper(profile);
+    }
   }
 
   std::exception_ptr ep = nullptr;
@@ -830,14 +846,21 @@ void Input::set_matrices(unsigned nb_thread) {
   auto run_on_profiles = [&](const std::vector<std::string>& profiles) {
     try {
       for (const auto& profile : profiles) {
-        auto d_m = _durations_matrices.find(profile);
-        assert(d_m != _durations_matrices.end());
-
-        if (d_m->second.size() == 0) {
-          // Durations matrix not manually set so defined as empty
-          // above.
+        auto durations_m = _durations_matrices.find(profile);
+        auto distances_m = _distances_matrices.find(profile);
+
+        // Required matrices not manually set have been defined as
+        // empty above.
+        assert(durations_m != _durations_matrices.end());
+        const bool define_durations = (durations_m->second.size() == 0);
+        const bool has_distance_matrix =
+          (distances_m != _distances_matrices.end());
+        const bool define_distances =
+          has_distance_matrix and (distances_m->second.size() == 0);
+        if (define_durations or define_distances) {
           if (_locations.size() == 1) {
-            d_m->second = Matrix<UserCost>(1);
+            durations_m->second = Matrix<UserDuration>(1);
+            distances_m->second = Matrix<UserDistance>(1);
           } else {
             auto rw = std::find_if(_routing_wrappers.begin(),
                                    _routing_wrappers.end(),
@@ -848,36 +871,66 @@ void Input::set_matrices(unsigned nb_thread) {
 
             if (!_has_custom_location_index) {
               // Location indices are set based on order in _locations.
-              d_m->second = (*rw)->get_matrix(_locations);
+              auto matrices = (*rw)->get_matrices(_locations);
+              if (define_durations) {
+                durations_m->second = std::move(matrices.durations);
+              }
+              if (define_distances) {
+                distances_m->second = std::move(matrices.distances);
+              }
             } else {
               // Location indices are provided in input so we need an
               // indirection based on order in _locations.
-              auto m = (*rw)->get_matrix(_locations);
-
-              Matrix<UserDuration> full_m(_max_matrices_used_index + 1);
-              for (Index i = 0; i < _locations.size(); ++i) {
-                const auto& loc_i = _locations[i];
-                for (Index j = 0; j < _locations.size(); ++j) {
-                  full_m[loc_i.index()][_locations[j].index()] = m[i][j];
+              auto matrices = (*rw)->get_matrices(_locations);
+
+              if (define_durations) {
+                Matrix<UserDuration> full_m(_max_matrices_used_index + 1);
+                for (Index i = 0; i < _locations.size(); ++i) {
+                  const auto& loc_i = _locations[i];
+                  for (Index j = 0; j < _locations.size(); ++j) {
+                    full_m[loc_i.index()][_locations[j].index()] =
+                      matrices.durations[i][j];
+                  }
                 }
+
+                durations_m->second = std::move(full_m);
               }
+              if (define_distances) {
+                Matrix<UserDistance> full_m(_max_matrices_used_index + 1);
+                for (Index i = 0; i < _locations.size(); ++i) {
+                  const auto& loc_i = _locations[i];
+                  for (Index j = 0; j < _locations.size(); ++j) {
+                    full_m[loc_i.index()][_locations[j].index()] =
+                      matrices.distances[i][j];
+                  }
+                }
 
-              d_m->second = std::move(full_m);
+                distances_m->second = std::move(full_m);
+              }
             }
           }
         }
 
-        if (d_m->second.size() <= _max_matrices_used_index) {
-          throw InputException("location_index exceeding matrix size for " +
-                               profile + " profile.");
+        if (durations_m->second.size() <= _max_matrices_used_index) {
+          throw InputException(
+            "location_index exceeding durations matrix size for " + profile +
+            " profile.");
+        }
+
+        if (has_distance_matrix and
+            distances_m->second.size() <= _max_matrices_used_index) {
+          throw InputException(
+            "location_index exceeding distances matrix size for " + profile +
+            " profile.");
         }
 
         const auto c_m = _costs_matrices.find(profile);
 
         if (c_m != _costs_matrices.end()) {
           if (c_m->second.size() <= _max_matrices_used_index) {
-            throw InputException("location_index exceeding matrix size for " +
-                                 profile + " profile.");
+            throw InputException(
+              "location_index exceeding costs matrix size for " + profile +
+              " profile.");
           }
 
           // Check for potential overflow in solution cost.
@@ -889,7 +942,7 @@ void Input::set_matrices(unsigned nb_thread) {
           cost_bound_m.unlock();
         } else {
           // Durations matrix will be used for costs.
-          const UserCost current_bound = check_cost_bound(d_m->second);
+          const UserCost current_bound = check_cost_bound(durations_m->second);
           cost_bound_m.lock();
           _cost_upper_bound =
             std::max(_cost_upper_bound,
diff --git a/src/structures/vroom/input/input.h b/src/structures/vroom/input/input.h
index 249143301..705dbebd0 100644
--- a/src/structures/vroom/input/input.h
+++ b/src/structures/vroom/input/input.h
@@ -50,6 +50,7 @@ class Input {
   bool _has_jobs{false};
   bool _has_shipments{false};
   std::unordered_map<std::string, Matrix<UserDuration>> _durations_matrices;
+  std::unordered_map<std::string, Matrix<UserDistance>> _distances_matrices;
   std::unordered_map<std::string, Matrix<UserCost>> _costs_matrices;
   Cost _cost_upper_bound{0};
   std::vector<Location> _locations;
@@ -110,6 +111,10 @@ class Input {
 
   void set_durations_matrix(const std::string& profile,
                             Matrix<UserDuration>&& m);
+
+  void set_distances_matrix(const std::string& profile,
+                            Matrix<UserDistance>&& m);
+
   void set_costs_matrix(const std::string& profile, Matrix<UserCost>&& m);
 
   const Amount& zero_amount() const {
diff --git a/src/structures/vroom/solution/route.h b/src/structures/vroom/solution/route.h
index 3c25f4553..3c93cc420 100644
--- a/src/structures/vroom/solution/route.h
+++ b/src/structures/vroom/solution/route.h
@@ -33,7 +33,7 @@ struct Route {
   Violations violations;
 
   std::string geometry;
-  Distance distance;
+  UserDistance distance;
 
   Route();
 
diff --git a/src/structures/vroom/solution/step.h b/src/structures/vroom/solution/step.h
index 048332707..47c8010b9 100644
--- a/src/structures/vroom/solution/step.h
+++ b/src/structures/vroom/solution/step.h
@@ -32,7 +32,7 @@ struct Step {
   UserDuration arrival;
   UserDuration duration;
   UserDuration waiting_time;
-  Distance distance;
+  UserDistance distance;
 
   Violations violations;
 
diff --git a/src/structures/vroom/solution/summary.h b/src/structures/vroom/solution/summary.h
index 8a8ef363d..6466ec37e 100644
--- a/src/structures/vroom/solution/summary.h
+++ b/src/structures/vroom/solution/summary.h
@@ -29,7 +29,7 @@ struct Summary {
 
   UserDuration duration;
   UserDuration waiting_time;
-  Distance distance;
+  UserDistance distance;
   ComputingTimes computing_times;
 
   Violations violations;
diff --git a/src/utils/input_parser.cpp b/src/utils/input_parser.cpp
index db216e6a5..8d2b752bc 100644
--- a/src/utils/input_parser.cpp
+++ b/src/utils/input_parser.cpp
@@ -603,6 +603,11 @@ void parse(Input& input, const std::string& input_str, bool geometry) {
                                      get_matrix<UserDuration>(
                                        profile_entry.value["durations"]));
         }
+        if (profile_entry.value.HasMember("distances")) {
+          input.set_distances_matrix(profile_entry.name.GetString(),
+                                     get_matrix<UserDistance>(
+                                       profile_entry.value["distances"]));
+        }
         if (profile_entry.value.HasMember("costs")) {
           input.set_costs_matrix(profile_entry.name.GetString(),
                                  get_matrix<UserCost>(