diff --git a/CMakeLists.txt b/CMakeLists.txt
index c57e3b7fe6..d7627e942d 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -401,6 +401,7 @@ set(CORE_SOURCE
         src/Record.cpp
         src/RecordComponent.cpp
         src/Series.cpp
+        src/UnitDimension.cpp
         src/version.cpp
         src/WriteIterations.cpp
         src/auxiliary/Date.cpp
diff --git a/include/openPMD/IO/AbstractIOHandler.hpp b/include/openPMD/IO/AbstractIOHandler.hpp
index 1288a87b21..649252a877 100644
--- a/include/openPMD/IO/AbstractIOHandler.hpp
+++ b/include/openPMD/IO/AbstractIOHandler.hpp
@@ -25,6 +25,7 @@
 #include "openPMD/IO/IOTask.hpp"
 #include "openPMD/IterationEncoding.hpp"
 #include "openPMD/config.hpp"
+#include "openPMD/version.hpp"
 
 #if openPMD_HAVE_MPI
 #include <mpi.h>
@@ -75,6 +76,20 @@ enum class FlushLevel
     CreateOrOpenFiles
 };
 
+enum class OpenpmdStandard
+{
+    v_1_0_0,
+    v_1_0_1,
+    v_1_1_0,
+    v_2_0_0
+};
+
+namespace auxiliary
+{
+    auto parseStandard(std::string const &) -> OpenpmdStandard;
+    auto formatStandard(OpenpmdStandard) -> char const *;
+} // namespace auxiliary
+
 namespace internal
 {
     /**
@@ -189,38 +204,7 @@ class AbstractIOHandler
     friend class detail::ADIOS2File;
 
 private:
-    IterationEncoding m_encoding = IterationEncoding::groupBased;
-
-    void setIterationEncoding(IterationEncoding encoding)
-    {
-        /*
-         * In file-based iteration encoding, the APPEND mode is handled entirely
-         * by the frontend, the backend should just treat it as CREATE mode.
-         * Similar for READ_LINEAR which should be treated as READ_RANDOM_ACCESS
-         * in the backend.
-         */
-        if (encoding == IterationEncoding::fileBased)
-        {
-            switch (m_backendAccess)
-            {
-
-            case Access::READ_LINEAR:
-                // do we really want to have those as const members..?
-                *const_cast<Access *>(&m_backendAccess) =
-                    Access::READ_RANDOM_ACCESS;
-                break;
-            case Access::APPEND:
-                *const_cast<Access *>(&m_backendAccess) = Access::CREATE;
-                break;
-            case Access::READ_RANDOM_ACCESS:
-            case Access::READ_WRITE:
-            case Access::CREATE:
-                break;
-            }
-        }
-
-        m_encoding = encoding;
-    }
+    void setIterationEncoding(IterationEncoding encoding);
 
 public:
 #if openPMD_HAVE_MPI
@@ -284,8 +268,14 @@ class AbstractIOHandler
      */
     Access m_backendAccess;
     Access m_frontendAccess;
-    internal::SeriesStatus m_seriesStatus = internal::SeriesStatus::Default;
     std::queue<IOTask> m_work;
+
+    /**************************************************************************
+     * Since the AbstractIOHandler is linked to every object of the frontend, *
+     * it stores a number of members that are needed by methods traversing    *
+     * the object hierarchy. Those members are found below.                   *
+     **************************************************************************/
+
     /**
      * This is to avoid that the destructor tries flushing again if an error
      * happened. Otherwise, this would lead to confusing error messages.
@@ -294,6 +284,9 @@ class AbstractIOHandler
      * The destructor will only attempt flushing again if this is true.
      */
     bool m_lastFlushSuccessful = false;
+    internal::SeriesStatus m_seriesStatus = internal::SeriesStatus::Default;
+    IterationEncoding m_encoding = IterationEncoding::groupBased;
+    OpenpmdStandard m_standard = auxiliary::parseStandard(getStandardDefault());
 }; // AbstractIOHandler
 
 } // namespace openPMD
diff --git a/include/openPMD/Mesh.hpp b/include/openPMD/Mesh.hpp
index 53274ac7d4..77ef8b2886 100644
--- a/include/openPMD/Mesh.hpp
+++ b/include/openPMD/Mesh.hpp
@@ -20,11 +20,11 @@
  */
 #pragma once
 
+#include "openPMD/UnitDimension.hpp"
 #include "openPMD/backend/Attributable.hpp"
 #include "openPMD/backend/BaseRecord.hpp"
 #include "openPMD/backend/MeshRecordComponent.hpp"
 
-#include <array>
 #include <ostream>
 #include <string>
 #include <type_traits>
@@ -155,7 +155,7 @@ class Mesh : public BaseRecord<MeshRecordComponent>
      */
     template <
         typename T,
-        typename = std::enable_if_t<std::is_floating_point<T>::value> >
+        typename = std::enable_if_t<std::is_floating_point<T>::value>>
     Mesh &setGridSpacing(std::vector<T> const &gridSpacing);
 
     /**
@@ -184,6 +184,10 @@ class Mesh : public BaseRecord<MeshRecordComponent>
      * Mesh::gridSpacing and Mesh::gridGlobalOffset, in order to convert from
      * simulation units to SI units.
      *
+     * Valid for openPMD version 1.*.
+     * In order to specify the gridUnitSI per dimension (openPMD 2.*),
+     * use the vector overload or `setGridUnitSIPerDimension()`.
+     *
      * @param   gridUnitSI  unit-conversion factor to multiply each value in
      * Mesh::gridSpacing and Mesh::gridGlobalOffset, in order to convert from
      * simulation units to SI units.
@@ -191,6 +195,48 @@ class Mesh : public BaseRecord<MeshRecordComponent>
      */
     Mesh &setGridUnitSI(double gridUnitSI);
 
+    /** Alias for `setGridUnitSIPerDimension(std::vector<double>)`.
+     *
+     * Set the unit-conversion factors per axis to multiply each value in
+     * Mesh::gridSpacing and Mesh::gridGlobalOffset, in order to convert from
+     * simulation units to SI units.
+     *
+     * Valid for openPMD 2.*.
+     * The legacy behavior (openPMD 1.*, a scalar gridUnitSI) is implemented
+     * by `setGridUnitSI(double)`.
+     *
+     * @param   gridUnitSI  unit-conversion factor to multiply each value in
+     * Mesh::gridSpacing and Mesh::gridGlobalOffset, in order to convert from
+     * simulation units to SI units.
+     *
+     * @return  Reference to modified mesh.
+     */
+    Mesh &setGridUnitSI(std::vector<double> const &gridUnitSI);
+
+    /**
+     * @return  A vector of the gridUnitSI per grid axis in the order of
+     * the axisLabels. If the gridUnitSI is defined as a scalar
+     * (legacy openPMD), the dimensionality is determined and a vector of
+     * `dimensionality` times the scalar vector is returned.
+     */
+    std::vector<double> gridUnitSIPerDimension() const;
+
+    /* Set the unit-conversion factors per axis to multiply each value in
+     * Mesh::gridSpacing and Mesh::gridGlobalOffset, in order to convert from
+     * simulation units to SI units.
+     *
+     * Valid for openPMD 2.*.
+     * The legacy behavior (openPMD 1.*, a scalar gridUnitSI) is implemented
+     * by `setGridUnitSI(double)`.
+     *
+     * @param   gridUnitSI  unit-conversion factor to multiply each value in
+     * Mesh::gridSpacing and Mesh::gridGlobalOffset, in order to convert from
+     * simulation units to SI units.
+     *
+     * @return  Reference to modified mesh.
+     */
+    Mesh &setGridUnitSIPerDimension(std::vector<double> const &gridUnitSI);
+
     /** Set the powers of the 7 base measures characterizing the record's unit
      * in SI.
      *
@@ -198,8 +244,59 @@ class Mesh : public BaseRecord<MeshRecordComponent>
      * that represent the power of the particular base.
      * @return  Reference to modified mesh.
      */
+    Mesh &setUnitDimension(unit_representations::AsMap const &unitDimension);
+
+    /** Set the powers of the 7 base measures characterizing the record's unit
+     * in SI.
+     *
+     * @param   unitDimension   array containing seven doubles, each
+     * representing the power of the particular base in order.
+     * @return  Reference to modified mesh.
+     */
+    Mesh &setUnitDimension(unit_representations::AsArray const &unitDimension);
+
+    /**
+     * @brief Set the unitDimension for each axis of the current grid.
+     *
+     * @param gridUnitDimension A vector of the unitDimensions for each
+     * axis of the grid in the order of the axisLabels, in dict representation.
+
+     * Behavior note: This is an updating method, meaning that an SI unit that
+     * has been defined before and is in the next call not explicitly set
+     * in the `std::map<UnitDimension, double>` will keep its previous value.
+     *
+     * @return Reference to modified mesh.
+     */
     Mesh &
-    setUnitDimension(std::map<UnitDimension, double> const &unitDimension);
+    setGridUnitDimension(unit_representations::AsMaps const &gridUnitDimension);
+
+    /**
+     * @brief Set the unitDimension for each axis of the current grid.
+     *
+     * @param gridUnitDimension A vector of the unitDimensions for each
+     * axis of the grid in the order of the axisLabels, in array representation.
+
+     * Behavior note: This is an updating method, meaning that an SI unit that
+     * has been defined before and is in the next call not explicitly set
+     * in the `std::map<UnitDimension, double>` will keep its previous value.
+     *
+     * @return Reference to modified mesh.
+     */
+    Mesh &setGridUnitDimension(
+        unit_representations::AsArrays const &gridUnitDimension);
+
+    /**
+     * @brief Return the physical dimensions of the mesh axes.
+
+     * If the attribute is not defined, the axes are assumed to be spatial
+     * and the return value will be according to this assumption.
+     * If the attribute is defined, the dimensionality of the return value is
+     * not checked against the dimensionality of the mesh.
+     *
+     * @return A vector of arrays, each array representing the SI unit of one
+     * mesh axis.
+     */
+    unit_representations::AsArrays gridUnitDimension() const;
 
     /**
      * @tparam  T   Floating point type of user-selected precision (e.g. float,
@@ -222,7 +319,7 @@ class Mesh : public BaseRecord<MeshRecordComponent>
      */
     template <
         typename T,
-        typename = std::enable_if_t<std::is_floating_point<T>::value> >
+        typename = std::enable_if_t<std::is_floating_point<T>::value>>
     Mesh &setTimeOffset(T timeOffset);
 
 private:
diff --git a/include/openPMD/Record.hpp b/include/openPMD/Record.hpp
index c875389db5..791c4c15f8 100644
--- a/include/openPMD/Record.hpp
+++ b/include/openPMD/Record.hpp
@@ -21,9 +21,9 @@
 #pragma once
 
 #include "openPMD/RecordComponent.hpp"
+#include "openPMD/UnitDimension.hpp"
 #include "openPMD/backend/BaseRecord.hpp"
 
-#include <map>
 #include <string>
 #include <type_traits>
 
@@ -40,7 +40,8 @@ class Record : public BaseRecord<RecordComponent>
     Record &operator=(Record const &) = default;
     ~Record() override = default;
 
-    Record &setUnitDimension(std::map<UnitDimension, double> const &);
+    Record &setUnitDimension(unit_representations::AsMap const &);
+    Record &setUnitDimension(unit_representations::AsArray const &);
 
     template <typename T>
     T timeOffset() const;
diff --git a/include/openPMD/UnitDimension.hpp b/include/openPMD/UnitDimension.hpp
index 232a6ee1f7..480901c0ca 100644
--- a/include/openPMD/UnitDimension.hpp
+++ b/include/openPMD/UnitDimension.hpp
@@ -20,10 +20,16 @@
  */
 #pragma once
 
+#include <array>
 #include <cstdint>
+#include <map>
+#include <vector>
 
 namespace openPMD
 {
+
+using UnitDimensionExponent = double;
+
 /** Physical dimension of a record
  *
  * Dimensional base quantities of the international system of quantities
@@ -38,4 +44,25 @@ enum class UnitDimension : uint8_t
     N, //!< amount of substance
     J //!< luminous intensity
 };
+
+namespace unit_representations
+{
+    using AsMap = std::map<UnitDimension, UnitDimensionExponent>;
+    using AsArray = std::array<UnitDimensionExponent, 7>;
+
+    using AsMaps = std::vector<AsMap>;
+    using AsArrays = std::vector<AsArray>;
+
+    auto asArray(AsMap const &) -> AsArray;
+    auto asMap(AsArray const &, bool skip_zeros = true) -> AsMap;
+
+    auto asArrays(AsMaps const &) -> AsArrays;
+    auto asMaps(AsArrays const &, bool skip_zeros = true) -> AsMaps;
+
+    namespace auxiliary
+    {
+        void fromMapOfUnitDimension(
+            double *cursor, std::map<UnitDimension, double> const &udim);
+    } // namespace auxiliary
+} // namespace unit_representations
 } // namespace openPMD
diff --git a/include/openPMD/backend/Attributable.hpp b/include/openPMD/backend/Attributable.hpp
index 850d3bf35c..6ff5b85ed1 100644
--- a/include/openPMD/backend/Attributable.hpp
+++ b/include/openPMD/backend/Attributable.hpp
@@ -328,6 +328,8 @@ class Attributable
      */
     void touch();
 
+    [[nodiscard]] OpenpmdStandard openPMDStandard() const;
+
     // clang-format off
 OPENPMD_protected
     // clang-format on
diff --git a/include/openPMD/backend/BaseRecord.hpp b/include/openPMD/backend/BaseRecord.hpp
index ba137b10db..87364b46f4 100644
--- a/include/openPMD/backend/BaseRecord.hpp
+++ b/include/openPMD/backend/BaseRecord.hpp
@@ -474,7 +474,7 @@ class BaseRecord
      *
      * @return powers of the 7 base measures in the order specified above
      */
-    std::array<double, 7> unitDimension() const;
+    unit_representations::AsArray unitDimension() const;
 
     void setDatasetDefined(BaseRecordComponent::Data_t &data) override
     {
@@ -928,7 +928,7 @@ auto BaseRecord<T_elem>::emplace(Args &&...args) -> std::pair<iterator, bool>
 }
 
 template <typename T_elem>
-inline std::array<double, 7> BaseRecord<T_elem>::unitDimension() const
+inline unit_representations::AsArray BaseRecord<T_elem>::unitDimension() const
 {
     return this->getAttribute("unitDimension")
         .template get<std::array<double, 7>>();
diff --git a/include/openPMD/binding/python/UnitDimension.hpp b/include/openPMD/binding/python/UnitDimension.hpp
index ab4cbde650..b7a402f0f6 100644
--- a/include/openPMD/binding/python/UnitDimension.hpp
+++ b/include/openPMD/binding/python/UnitDimension.hpp
@@ -24,7 +24,7 @@ namespace openPMD
 {
 namespace python
 {
-    constexpr auto doc_unit_dimension = R"docstr(
+    constexpr auto doc_unit_dimension = &R"docstr(
 Return the physical dimension (quantity) of a record
 
 Annotating the physical dimension of a record allows us to read data
@@ -40,7 +40,29 @@ See https://en.wikipedia.org/wiki/International_System_of_Quantities#Base_quanti
 See https://github.com/openPMD/openPMD-standard/blob/1.1.0/STANDARD.md#required-for-each-record
 
 Returns the powers of the 7 base measures in the order specified above.
-)docstr";
+)docstr"[1];
+
+    constexpr auto doc_mesh_unit_dimension = &R"docstr(
+Return the physical dimension (quantity) of the record axes
+
+Annotating the physical dimension of the record axes allows us to read data
+sets with arbitrary names and understand their purpose simply by
+dimensional analysis. The dimensional base quantities in openPMD are
+in order: length (L), mass (M), time (T), electric current (I),
+thermodynamic temperature (theta), amount of substance (N),
+luminous intensity (J) after the international system of quantities
+(ISQ).
+This attribute may be left out, the axes will then be interpreted as spatial.
+
+See https://en.wikipedia.org/wiki/Dimensional_analysis
+See https://en.wikipedia.org/wiki/International_System_of_Quantities#Base_quantities
+See https://github.com/openPMD/openPMD-standard/blob/1.1.0/STANDARD.md#required-for-each-record
+
+Returns the powers of the 7 base measures in the order specified above, listed
+for each axis in the order of the axisLabels.
+This attribute has been introduced as part of openPMD 2.0.0 in:
+Ref.: https://github.com/openPMD/openPMD-standard/pull/193
+)docstr"[1];
 
 } // namespace python
 } // namespace openPMD
diff --git a/src/IO/AbstractIOHandler.cpp b/src/IO/AbstractIOHandler.cpp
index 440b663286..3d83978d1d 100644
--- a/src/IO/AbstractIOHandler.cpp
+++ b/src/IO/AbstractIOHandler.cpp
@@ -21,10 +21,79 @@
 
 #include "openPMD/IO/AbstractIOHandler.hpp"
 
+#include "openPMD/Error.hpp"
 #include "openPMD/IO/FlushParametersInternal.hpp"
+#include <utility>
+
+namespace openPMD::auxiliary
+{
+using pair_t = std::pair<OpenpmdStandard, char const *>;
+constexpr pair_t STANDARD_VERSIONS[] = {
+    pair_t{OpenpmdStandard::v_1_0_0, "1.0.0"},
+    pair_t{OpenpmdStandard::v_1_0_1, "1.0.1"},
+    pair_t{OpenpmdStandard::v_1_1_0, "1.1.0"},
+    pair_t{OpenpmdStandard::v_2_0_0, "2.0.0"}};
+
+auto parseStandard(const std::string &str) -> OpenpmdStandard
+{
+    for (auto const &[res, compare] : STANDARD_VERSIONS)
+    {
+        if (str == compare)
+        {
+            return res;
+        }
+    }
+    throw error::IllegalInOpenPMDStandard(
+        "Standard version is not supported: '" + str + "'.");
+}
+
+auto formatStandard(OpenpmdStandard std) -> char const *
+{
+    for (auto const &[compare, res] : STANDARD_VERSIONS)
+    {
+        if (std == compare)
+        {
+            return res;
+        }
+    }
+    throw error::Internal(
+        "[auxiliary::formatStandard] Match should be exhaustive.");
+}
+} // namespace openPMD::auxiliary
 
 namespace openPMD
 {
+void AbstractIOHandler::setIterationEncoding(IterationEncoding encoding)
+{
+    /*
+     * In file-based iteration encoding, the APPEND mode is handled entirely
+     * by the frontend, the backend should just treat it as CREATE mode.
+     * Similar for READ_LINEAR which should be treated as READ_RANDOM_ACCESS
+     * in the backend.
+     */
+    if (encoding == IterationEncoding::fileBased)
+    {
+        switch (m_backendAccess)
+        {
+
+        case Access::READ_LINEAR:
+            // do we really want to have those as const members..?
+            *const_cast<Access *>(&m_backendAccess) =
+                Access::READ_RANDOM_ACCESS;
+            break;
+        case Access::APPEND:
+            *const_cast<Access *>(&m_backendAccess) = Access::CREATE;
+            break;
+        case Access::READ_RANDOM_ACCESS:
+        case Access::READ_WRITE:
+        case Access::CREATE:
+            break;
+        }
+    }
+
+    m_encoding = encoding;
+}
+
 std::future<void> AbstractIOHandler::flush(internal::FlushParams const &params)
 {
     internal::ParsedFlushParams parsedParams{params};
diff --git a/src/Iteration.cpp b/src/Iteration.cpp
index 366fea0de1..28839fea3b 100644
--- a/src/Iteration.cpp
+++ b/src/Iteration.cpp
@@ -493,10 +493,10 @@ void Iteration::read_impl(std::string const &groupPath)
     Series s = retrieveSeries();
 
     Parameter<Operation::LIST_PATHS> pList;
-    std::string version = s.openPMD();
+    auto version = IOHandler()->m_standard;
     bool hasMeshes = false;
     bool hasParticles = false;
-    if (version == "1.0.0" || version == "1.0.1")
+    if (version <= OpenpmdStandard::v_1_0_1)
     {
         IOHandler()->enqueue(IOTask(this, pList));
         IOHandler()->flush(internal::defaultFlushParams);
diff --git a/src/Mesh.cpp b/src/Mesh.cpp
index f977bbe905..c5cdefe483 100644
--- a/src/Mesh.cpp
+++ b/src/Mesh.cpp
@@ -20,9 +20,13 @@
  */
 #include "openPMD/Mesh.hpp"
 #include "openPMD/Error.hpp"
+#include "openPMD/IO/AbstractIOHandler.hpp"
 #include "openPMD/Series.hpp"
+#include "openPMD/ThrowError.hpp"
+#include "openPMD/UnitDimension.hpp"
 #include "openPMD/auxiliary/DerefDynamicCast.hpp"
 #include "openPMD/auxiliary/StringManip.hpp"
+#include "openPMD/backend/Attribute.hpp"
 #include "openPMD/backend/Writable.hpp"
 
 #include <algorithm>
@@ -40,7 +44,6 @@ Mesh::Mesh()
     setAxisLabels({"x"}); // empty strings are not allowed in HDF5
     setGridSpacing(std::vector<double>{1});
     setGridGlobalOffset({0});
-    setGridUnitSI(1);
 }
 
 Mesh::Geometry Mesh::geometry() const
@@ -181,22 +184,168 @@ double Mesh::gridUnitSI() const
 
 Mesh &Mesh::setGridUnitSI(double gusi)
 {
+    if (auto standard = IOHandler()->m_standard;
+        standard >= OpenpmdStandard::v_2_0_0)
+    {
+        std::cerr << "[Mesh::setGridUnitSI] Warning: Setting a scalar "
+                     "`gridUnitSI` in a file with openPMD version '" +
+                std::string(auxiliary::formatStandard(standard)) +
+                "'. Consider specifying a vector instead in order to "
+                "specify the gridUnitSI per axis (ref.: "
+                "https://github.com/openPMD/openPMD-standard/pull/193).\n";
+    }
     setAttribute("gridUnitSI", gusi);
     return *this;
 }
 
-Mesh &Mesh::setUnitDimension(std::map<UnitDimension, double> const &udim)
+Mesh &Mesh::setGridUnitSI(std::vector<double> const &gusi)
+{
+    return setGridUnitSIPerDimension(gusi);
+}
+
+namespace
+{
+    uint64_t retrieveMeshDimensionality(Mesh const &m)
+    {
+        if (m.containsAttribute("axisLabels"))
+        {
+            return m.axisLabels().size();
+        }
+
+        // maybe we have record components and can ask them
+        if (auto it = m.begin(); it != m.end())
+        {
+            return it->second.getDimensionality();
+        }
+        /*
+         * Since some backends cannot distinguish between vector and
+         * scalar values, the most likely answer here is 1.
+         */
+        return 1;
+    }
+} // namespace
+
+std::vector<double> Mesh::gridUnitSIPerDimension() const
+{
+    if (containsAttribute("gridUnitSI"))
+    {
+        if (IOHandler()->m_standard < OpenpmdStandard::v_2_0_0)
+        {
+            // If the openPMD version is lower than 2.0, the gridUnitSI is a
+            // scalar interpreted for all axes. Copy it d times.
+            return std::vector<double>(
+                retrieveMeshDimensionality(*this),
+                getAttribute("gridUnitSI").get<double>());
+        }
+        return getAttribute("gridUnitSI").get<std::vector<double>>();
+    }
+    else
+    {
+        // gridUnitSI is an optional attribute
+        // if it is missing, the mesh is interpreted as unscaled
+        return std::vector<double>(retrieveMeshDimensionality(*this), 1.);
+    }
+}
+
+Mesh &Mesh::setGridUnitSIPerDimension(std::vector<double> const &gridUnitSI)
+{
+    if (auto standard = IOHandler()->m_standard;
+        standard < OpenpmdStandard::v_2_0_0)
+    {
+        throw error::IllegalInOpenPMDStandard(
+            "[Mesh::setGridUnitSI] Setting `gridUnitSI` as a vector in a "
+            "file with openPMD version '" +
+            std::string(auxiliary::formatStandard(standard)) +
+            "', but per-axis specification is only supported as of "
+            "openPMD 2.0. Either upgrade the file to openPMD >= 2.0 "
+            "or specify a scalar that applies to all axes.");
+    }
+    setAttribute("gridUnitSI", gridUnitSI);
+    return *this;
+}
+
+Mesh &Mesh::setUnitDimension(unit_representations::AsMap const &udim)
 {
     if (!udim.empty())
     {
         std::array<double, 7> tmpUnitDimension = this->unitDimension();
-        for (auto const &entry : udim)
-            tmpUnitDimension[static_cast<uint8_t>(entry.first)] = entry.second;
+        unit_representations::auxiliary::fromMapOfUnitDimension(
+            tmpUnitDimension.data(), udim);
         setAttribute("unitDimension", tmpUnitDimension);
     }
     return *this;
 }
 
+Mesh &Mesh::setUnitDimension(unit_representations::AsArray const &udim)
+{
+    return setUnitDimension(
+        unit_representations::asMap(udim, /* skip_zeros = */ false));
+}
+
+Mesh &Mesh::setGridUnitDimension(unit_representations::AsMaps const &udims)
+{
+    auto rawGridUnitDimension = [this]() {
+        try
+        {
+            return this->getAttribute("gridUnitDimension")
+                .get<std::vector<double>>();
+        }
+        catch (no_such_attribute_error const &)
+        {
+            return std::vector<double>();
+        }
+    }();
+    rawGridUnitDimension.resize(7 * udims.size());
+    auto cursor = rawGridUnitDimension.begin();
+    for (auto const &udim : udims)
+    {
+        unit_representations::auxiliary::fromMapOfUnitDimension(&*cursor, udim);
+        cursor += 7;
+    }
+    setAttribute("gridUnitDimension", rawGridUnitDimension);
+    return *this;
+}
+
+Mesh &Mesh::setGridUnitDimension(unit_representations::AsArrays const &udims)
+{
+    return setGridUnitDimension(
+        unit_representations::asMaps(udims, /* skip_zeros = */ false));
+}
+
+unit_representations::AsArrays Mesh::gridUnitDimension() const
+{
+    if (containsAttribute("gridUnitDimension"))
+    {
+        std::vector<double> rawRes =
+            getAttribute("gridUnitDimension").get<std::vector<double>>();
+        if (rawRes.size() % 7 != 0)
+        {
+            throw error::ReadError(
+                error::AffectedObject::Attribute,
+                error::Reason::UnexpectedContent,
+                std::nullopt,
+                "[Mesh::gridUnitDimension()] `gridUnitDimension` attribute "
+                "must have a length equal to a multiple of 7.");
+        }
+        unit_representations::AsArrays res(rawRes.size() / 7);
+        for (size_t dim = 0; dim < res.size(); ++dim)
+        {
+            std::copy_n(rawRes.begin() + dim * 7, 7, res.at(dim).begin());
+        }
+        return res;
+    }
+    else
+    {
+        // gridUnitDimension is an optional attribute
+        // if it is missing, the mesh is interpreted as spatial
+        auto spatialMesh =
+            unit_representations::asArray({{UnitDimension::L, 1}});
+        auto dim = retrieveMeshDimensionality(*this);
+        unit_representations::AsArrays res(dim, spatialMesh);
+        return res;
+    }
+}
+
 template <typename T, typename>
 Mesh &Mesh::setTimeOffset(T to)
 {
@@ -263,6 +412,18 @@ void Mesh::flush_impl(
                     comp.second.flush(comp.first, flushParams);
             }
         }
+        if (!containsAttribute("gridUnitSI"))
+        {
+            if (IOHandler()->m_standard < OpenpmdStandard::v_2_0_0)
+            {
+                setGridUnitSI(1);
+            }
+            else
+            {
+                setGridUnitSIPerDimension(
+                    std::vector<double>(retrieveMeshDimensionality(*this), 1));
+            }
+        }
         flushAttributes(flushParams);
     }
 }
@@ -385,17 +546,35 @@ void Mesh::read()
     aRead.name = "gridUnitSI";
     IOHandler()->enqueue(IOTask(this, aRead));
     IOHandler()->flush(internal::defaultFlushParams);
-    if (auto val = Attribute(*aRead.resource).getOptional<double>();
-        val.has_value())
-        setGridUnitSI(val.value());
+    if (IOHandler()->m_standard >= OpenpmdStandard::v_2_0_0)
+    {
+        if (auto val =
+                Attribute(*aRead.resource).getOptional<std::vector<double>>();
+            val.has_value())
+            setGridUnitSIPerDimension(val.value());
+        else
+            throw error::ReadError(
+                error::AffectedObject::Attribute,
+                error::Reason::UnexpectedContent,
+                {},
+                "Unexpected Attribute datatype for 'gridUnitSI' "
+                "(expected vector of double, found " +
+                    datatypeToString(Attribute(*aRead.resource).dtype) + ")");
+    }
     else
-        throw error::ReadError(
-            error::AffectedObject::Attribute,
-            error::Reason::UnexpectedContent,
-            {},
-            "Unexpected Attribute datatype for 'gridUnitSI' (expected double, "
-            "found " +
-                datatypeToString(Attribute(*aRead.resource).dtype) + ")");
+    {
+        if (auto val = Attribute(*aRead.resource).getOptional<double>();
+            val.has_value())
+            setGridUnitSI(val.value());
+        else
+            throw error::ReadError(
+                error::AffectedObject::Attribute,
+                error::Reason::UnexpectedContent,
+                {},
+                "Unexpected Attribute datatype for 'gridUnitSI' "
+                "(expected double, found " +
+                    datatypeToString(Attribute(*aRead.resource).dtype) + ")");
+    }
 
     if (scalar())
     {
diff --git a/src/Record.cpp b/src/Record.cpp
index 3bcac4d7e1..7d41fce5c2 100644
--- a/src/Record.cpp
+++ b/src/Record.cpp
@@ -20,6 +20,7 @@
  */
 #include "openPMD/Record.hpp"
 #include "openPMD/RecordComponent.hpp"
+#include "openPMD/UnitDimension.hpp"
 #include "openPMD/backend/BaseRecord.hpp"
 
 #include <iostream>
@@ -31,7 +32,7 @@ Record::Record()
     setTimeOffset(0.f);
 }
 
-Record &Record::setUnitDimension(std::map<UnitDimension, double> const &udim)
+Record &Record::setUnitDimension(unit_representations::AsMap const &udim)
 {
     if (!udim.empty())
     {
@@ -42,6 +43,11 @@ Record &Record::setUnitDimension(std::map<UnitDimension, double> const &udim)
     }
     return *this;
 }
+Record &Record::setUnitDimension(unit_representations::AsArray const &udim)
+{
+    return setUnitDimension(
+        unit_representations::asMap(udim, /* skip_zeros = */ false));
+}
 
 void Record::flush_impl(
     std::string const &name, internal::FlushParams const &flushParams)
diff --git a/src/Series.cpp b/src/Series.cpp
index a3749c5e73..2d10575c17 100644
--- a/src/Series.cpp
+++ b/src/Series.cpp
@@ -139,12 +139,13 @@ std::string Series::openPMD() const
 
 Series &Series::setOpenPMD(std::string const &o)
 {
-    if (o >= "2.0")
+    setAttribute("openPMD", o);
+    auto standard = auxiliary::parseStandard(o);
+    IOHandler()->m_standard = standard;
+    if (standard >= OpenpmdStandard::v_2_0_0)
     {
-        std::cerr << "[Warning] openPMD 2.0 is still under development."
-                  << std::endl;
+        std::cerr << "[Warning] openPMD 2.0 is still under development.\n";
     }
-    setAttribute("openPMD", o);
     return *this;
 }
 
@@ -166,11 +167,15 @@ std::string Series::basePath() const
 
 Series &Series::setBasePath(std::string const &bp)
 {
-    std::string version = openPMD();
-    if (version == "1.0.0" || version == "1.0.1" || version == "1.1.0" ||
-        version == "2.0.0")
+    switch (IOHandler()->m_standard)
+    {
+    case OpenpmdStandard::v_1_0_0:
+    case OpenpmdStandard::v_1_0_1:
+    case OpenpmdStandard::v_1_1_0:
+    case OpenpmdStandard::v_2_0_0:
         throw std::runtime_error(
             "Custom basePath not allowed in openPMD <=2.0");
+    }
 
     setAttribute("basePath", bp);
     return *this;
@@ -630,7 +635,8 @@ Series &Series::setIterationFormat(std::string const &i)
             setBasePath(i);
         }
         else if (
-            basePath() != i && (openPMD() == "1.0.1" || openPMD() == "1.0.0"))
+            basePath() != i &&
+            IOHandler()->m_standard <= OpenpmdStandard::v_1_0_1)
             throw std::invalid_argument(
                 "iterationFormat must not differ from basePath " + basePath() +
                 " for group- or variableBased data");
diff --git a/src/UnitDimension.cpp b/src/UnitDimension.cpp
new file mode 100644
index 0000000000..3c8c7755a3
--- /dev/null
+++ b/src/UnitDimension.cpp
@@ -0,0 +1,59 @@
+#include "openPMD/UnitDimension.hpp"
+#include <algorithm>
+#include <iterator>
+
+namespace openPMD::unit_representations
+{
+auto asArray(AsMap const &udim) -> AsArray
+{
+    AsArray res{};
+    auxiliary::fromMapOfUnitDimension(res.data(), udim);
+    return res;
+}
+auto asMap(AsArray const &array, bool skip_zeros) -> AsMap
+{
+    AsMap udim;
+    for (size_t i = 0; i < array.size(); ++i)
+    {
+        if (!skip_zeros || array[i] != 0)
+        {
+            udim[static_cast<UnitDimension>(i)] = array[i];
+        }
+    }
+    return udim;
+}
+
+auto asArrays(AsMaps const &vec) -> AsArrays
+{
+    AsArrays res;
+    res.reserve(vec.size());
+    std::transform(
+        vec.begin(), vec.end(), std::back_inserter(res), [](auto const &map) {
+            return asArray(map);
+        });
+    return res;
+}
+auto asMaps(AsArrays const &vec, bool skip_zeros) -> AsMaps
+{
+    AsMaps res;
+    res.reserve(vec.size());
+    std::transform(
+        vec.begin(),
+        vec.end(),
+        std::back_inserter(res),
+        [&](auto const &array) { return asMap(array, skip_zeros); });
+    return res;
+}
+
+namespace auxiliary
+{
+    void fromMapOfUnitDimension(
+        double *cursor, std::map<UnitDimension, double> const &udim)
+    {
+        for (auto [unit, exponent] : udim)
+        {
+            cursor[static_cast<uint8_t>(unit)] = exponent;
+        }
+    }
+} // namespace auxiliary
+} // namespace openPMD::unit_representations
diff --git a/src/backend/Attributable.cpp b/src/backend/Attributable.cpp
index 9ffbe25ee3..e7b48efddf 100644
--- a/src/backend/Attributable.cpp
+++ b/src/backend/Attributable.cpp
@@ -19,6 +19,7 @@
  * If not, see <http://www.gnu.org/licenses/>.
  */
 #include "openPMD/backend/Attributable.hpp"
+#include "openPMD/IO/AbstractIOHandler.hpp"
 #include "openPMD/Iteration.hpp"
 #include "openPMD/ParticleSpecies.hpp"
 #include "openPMD/RecordComponent.hpp"
@@ -250,6 +251,11 @@ void Attributable::touch()
     setDirtyRecursive(true);
 }
 
+OpenpmdStandard Attributable::openPMDStandard() const
+{
+    return IOHandler()->m_standard;
+}
+
 template <bool flush_entire_series>
 void Attributable::seriesFlush_impl(internal::FlushParams const &flushParams)
 {
diff --git a/src/binding/python/Mesh.cpp b/src/binding/python/Mesh.cpp
index 9d53e15591..bdc6f182f6 100644
--- a/src/binding/python/Mesh.cpp
+++ b/src/binding/python/Mesh.cpp
@@ -19,6 +19,9 @@
  * If not, see <http://www.gnu.org/licenses/>.
  */
 #include "openPMD/Mesh.hpp"
+#include "openPMD/Error.hpp"
+#include "openPMD/IO/AbstractIOHandler.hpp"
+#include "openPMD/UnitDimension.hpp"
 #include "openPMD/backend/Attributable.hpp"
 #include "openPMD/backend/BaseRecord.hpp"
 #include "openPMD/backend/MeshRecordComponent.hpp"
@@ -29,6 +32,7 @@
 #include "openPMD/binding/python/UnitDimension.hpp"
 
 #include <string>
+#include <variant>
 #include <vector>
 
 void init_Mesh(py::module &m)
@@ -36,7 +40,7 @@ void init_Mesh(py::module &m)
     auto py_m_cont =
         declare_container<PyMeshContainer, Attributable>(m, "Mesh_Container");
 
-    py::class_<Mesh, BaseRecord<MeshRecordComponent> > cl(m, "Mesh");
+    py::class_<Mesh, BaseRecord<MeshRecordComponent>> cl(m, "Mesh");
 
     py::enum_<Mesh::Geometry>(m, "Geometry") // TODO: m -> cl
         .value("cartesian", Mesh::Geometry::cartesian)
@@ -58,9 +62,33 @@ void init_Mesh(py::module &m)
         .def_property(
             "unit_dimension",
             &Mesh::unitDimension,
-            &Mesh::setUnitDimension,
+            [](Mesh &self,
+               std::variant<
+                   unit_representations::AsMap,
+                   unit_representations::AsArray> const &arg) -> Mesh & {
+                return std::visit(
+                    [&](auto const &arg_resolved) -> Mesh & {
+                        return self.setUnitDimension(arg_resolved);
+                    },
+                    arg);
+            },
             python::doc_unit_dimension)
 
+        .def_property(
+            "grid_unit_dimension",
+            &Mesh::gridUnitDimension,
+            [](Mesh &self,
+               std::variant<
+                   unit_representations::AsMaps,
+                   unit_representations::AsArrays> const &arg) -> Mesh & {
+                return std::visit(
+                    [&](auto const &arg_resolved) -> Mesh & {
+                        return self.setGridUnitDimension(arg_resolved);
+                    },
+                    arg);
+            },
+            python::doc_mesh_unit_dimension)
+
         .def_property(
             "geometry",
             &Mesh::geometry,
@@ -95,14 +123,65 @@ void init_Mesh(py::module &m)
             "grid_global_offset",
             &Mesh::gridGlobalOffset,
             &Mesh::setGridGlobalOffset)
-        .def_property("grid_unit_SI", &Mesh::gridUnitSI, &Mesh::setGridUnitSI)
+        .def_property(
+            "grid_unit_SI",
+            /*
+             * Using pybind11's support for std::variant in order to implement a
+             * polymorphic type for this property. Will be a scalar double in
+             * openPMD 1.*, a list of double in openPMD 2.*.
+             * Unlike in the C++ API, this means that no new API calls
+             * such as gridUnitSIPerDimension() must be added.
+             */
+            [](Mesh &self) {
+                using return_t = std::variant<double, std::vector<double>>;
+                if (self.openPMDStandard() < OpenpmdStandard::v_2_0_0)
+                {
+                    return return_t(self.gridUnitSI());
+                }
+                else
+                {
+                    return return_t(self.gridUnitSIPerDimension());
+                }
+            },
+            [](Mesh &self,
+               std::variant<double, std::vector<double>> const &arg) -> Mesh & {
+                return std::visit(
+                    [&](auto const &arg_resolved) -> Mesh & {
+                        return self.setGridUnitSI(arg_resolved);
+                    },
+                    arg);
+            },
+            &R"(
+For openPMD versions 1.*:
+
+Set the unit-conversion factor to multiply each value in
+Mesh.grid_spacing and Mesh.grid_global_offset, in order to convert from
+simulation units to SI units.
+The type is a scalar floating point.
+
+For openPMD versions 2.*:
+
+Set the unit-conversion **factors per axis** in the order of the axisLabels
+to multiply each value in Mesh.grid_spacing and Mesh.grid_global_offset,
+in order to convert from simulation units to SI units.
+The type is a list of floating points.
+
+When writing a scalar value to an openPMD 2.* file, a warning will be printed
+(for enabling a more comfortable migration to openPMD 2.*).
+When writing a list value to an openPMD 1.* file, an error will be thrown,
+since most openPMD 1.*-based readers will not be able to interpret this
+properly.
+Ref.: https://github.com/openPMD/openPMD-standard/pull/193)"[1])
         .def_property(
             "time_offset",
             &Mesh::timeOffset<double>,
             &Mesh::setTimeOffset<double>)
 
         // TODO remove in future versions (deprecated)
-        .def("set_unit_dimension", &Mesh::setUnitDimension)
+        .def(
+            "set_unit_dimension",
+            py::overload_cast<unit_representations::AsMap const &>(
+                &Mesh::setUnitDimension))
         .def(
             "set_geometry",
             py::overload_cast<Mesh::Geometry>(&Mesh::setGeometry))
@@ -113,7 +192,9 @@ void init_Mesh(py::module &m)
         .def("set_grid_spacing", &Mesh::setGridSpacing<double>)
         .def("set_grid_spacing", &Mesh::setGridSpacing<long double>)
         .def("set_grid_global_offset", &Mesh::setGridGlobalOffset)
-        .def("set_grid_unit_SI", &Mesh::setGridUnitSI);
+        .def(
+            "set_grid_unit_SI",
+            py::overload_cast<double>(&Mesh::setGridUnitSI));
     add_pickle(
         cl, [](openPMD::Series series, std::vector<std::string> const &group) {
             uint64_t const n_it = std::stoull(group.at(1));
diff --git a/src/binding/python/Record.cpp b/src/binding/python/Record.cpp
index b4f732a83d..0d8a0e94bb 100644
--- a/src/binding/python/Record.cpp
+++ b/src/binding/python/Record.cpp
@@ -20,6 +20,7 @@
  */
 #include "openPMD/Record.hpp"
 #include "openPMD/RecordComponent.hpp"
+#include "openPMD/UnitDimension.hpp"
 #include "openPMD/backend/Attributable.hpp"
 #include "openPMD/backend/BaseRecord.hpp"
 
@@ -50,7 +51,16 @@ void init_Record(py::module &m)
         .def_property(
             "unit_dimension",
             &Record::unitDimension,
-            &Record::setUnitDimension,
+            [](Record &self,
+               std::variant<
+                   unit_representations::AsMap,
+                   unit_representations::AsArray> const &arg) -> Record & {
+                return std::visit(
+                    [&](auto const &arg_resolved) -> Record & {
+                        return self.setUnitDimension(arg_resolved);
+                    },
+                    arg);
+            },
             python::doc_unit_dimension)
 
         .def_property(
@@ -67,7 +77,10 @@ void init_Record(py::module &m)
             &Record::setTimeOffset<long double>)
 
         // TODO remove in future versions (deprecated)
-        .def("set_unit_dimension", &Record::setUnitDimension)
+        .def(
+            "set_unit_dimension",
+            py::overload_cast<unit_representations::AsMap const &>(
+                &Record::setUnitDimension))
         .def("set_time_offset", &Record::setTimeOffset<float>)
         .def("set_time_offset", &Record::setTimeOffset<double>)
         .def("set_time_offset", &Record::setTimeOffset<long double>);
diff --git a/src/binding/python/UnitDimension.cpp b/src/binding/python/UnitDimension.cpp
index 6e46a6cfcf..abd119c5b3 100644
--- a/src/binding/python/UnitDimension.cpp
+++ b/src/binding/python/UnitDimension.cpp
@@ -31,5 +31,25 @@ void init_UnitDimension(py::module &m)
         .value("I", UnitDimension::I)
         .value("theta", UnitDimension::theta)
         .value("N", UnitDimension::N)
-        .value("J", UnitDimension::J);
+        .value("J", UnitDimension::J)
+        .def(
+            "as_index",
+            [](UnitDimension ud) -> uint8_t {
+                return static_cast<uint8_t>(ud);
+            })
+        .def(
+            "from_index",
+            [](uint8_t idx) -> UnitDimension {
+                return static_cast<UnitDimension>(idx);
+            })
+        .def("as_array", &unit_representations::asArray)
+        .def(
+            "as_map",
+            &unit_representations::asMap,
+            py::arg("skip_zeros") = true)
+        .def("as_arrays", &unit_representations::asArrays)
+        .def(
+            "as_maps",
+            &unit_representations::asMaps,
+            py::arg("skip_zeros") = true);
 }
diff --git a/test/CoreTest.cpp b/test/CoreTest.cpp
index d27d68a8c5..2c4a0b1680 100644
--- a/test/CoreTest.cpp
+++ b/test/CoreTest.cpp
@@ -1,4 +1,5 @@
 // expose private and protected members for invasive testing
+#include "openPMD/UnitDimension.hpp"
 #if openPMD_USE_INVASIVE_TESTS
 #define OPENPMD_private public:
 #define OPENPMD_protected public:
@@ -24,6 +25,21 @@
 #include <string>
 #include <vector>
 
+// On Windows, REQUIRE() might not be able to print more complex data structures
+// upon failure:
+// CoreTest.obj : error LNK2001: unresolved external symbol
+// "class std::string const Catch::Detail::unprintableString" (...)
+#ifdef _WIN32
+#define OPENPMD_REQUIRE_GUARD_WINDOWS(...)                                     \
+    do                                                                         \
+    {                                                                          \
+        bool guarded_require_boolean = __VA_ARGS__;                            \
+        REQUIRE(guarded_require_boolean);                                      \
+    } while (0);
+#else
+#define OPENPMD_REQUIRE_GUARD_WINDOWS(...) REQUIRE(__VA_ARGS__)
+#endif
+
 using namespace openPMD;
 
 Dataset globalDataset(Datatype::CHAR, {1});
@@ -528,11 +544,10 @@ TEST_CASE("mesh_constructor_test", "[core]")
     REQUIRE(m.gridSpacing<double>() == gs);
     std::vector<double> ggo{0};
     REQUIRE(m.gridGlobalOffset() == ggo);
-    REQUIRE(m.gridUnitSI() == static_cast<double>(1));
     REQUIRE(
         m.numAttributes() ==
-        8); /* axisLabels, dataOrder, geometry, gridGlobalOffset, gridSpacing,
-               gridUnitSI, timeOffset, unitDimension */
+        7); /* axisLabels, dataOrder, geometry, gridGlobalOffset, gridSpacing,
+               timeOffset, unitDimension */
 
     o.flush();
     REQUIRE(m["x"].unitSI() == 1);
@@ -557,22 +572,22 @@ TEST_CASE("mesh_modification_test", "[core]")
 
     m.setGeometry(Mesh::Geometry::spherical);
     REQUIRE(m.geometry() == Mesh::Geometry::spherical);
-    REQUIRE(m.numAttributes() == 8);
+    REQUIRE(m.numAttributes() == 7);
     m.setDataOrder(Mesh::DataOrder::F);
     REQUIRE(m.dataOrder() == Mesh::DataOrder::F);
-    REQUIRE(m.numAttributes() == 8);
+    REQUIRE(m.numAttributes() == 7);
     std::vector<std::string> al{"z_", "y_", "x_"};
     m.setAxisLabels({"z_", "y_", "x_"});
     REQUIRE(m.axisLabels() == al);
-    REQUIRE(m.numAttributes() == 8);
+    REQUIRE(m.numAttributes() == 7);
     std::vector<double> gs{1e-5, 2e-5, 3e-5};
     m.setGridSpacing(gs);
     REQUIRE(m.gridSpacing<double>() == gs);
-    REQUIRE(m.numAttributes() == 8);
+    REQUIRE(m.numAttributes() == 7);
     std::vector<double> ggo{1e-10, 2e-10, 3e-10};
     m.setGridGlobalOffset({1e-10, 2e-10, 3e-10});
     REQUIRE(m.gridGlobalOffset() == ggo);
-    REQUIRE(m.numAttributes() == 8);
+    REQUIRE(m.numAttributes() == 7);
     m.setGridUnitSI(42.0);
     REQUIRE(m.gridUnitSI() == static_cast<double>(42));
     REQUIRE(m.numAttributes() == 8);
@@ -813,10 +828,10 @@ TEST_CASE("wrapper_test", "[core]")
     REQUIRE(copy.openPMDextension() == 42);
     REQUIRE(copy.iterationEncoding() == IterationEncoding::fileBased);
     REQUIRE(copy.name() == "new_openpmd_output_%T");
-    copy.setOpenPMD("1.2.0");
+    copy.setOpenPMD("1.1.0");
     copy.setIterationEncoding(IterationEncoding::groupBased);
     copy.setName("other_name");
-    REQUIRE(o.openPMD() == "1.2.0");
+    REQUIRE(o.openPMD() == "1.1.0");
     REQUIRE(o.iterationEncoding() == IterationEncoding::groupBased);
     REQUIRE(o.name() == "other_name");
 
@@ -1281,6 +1296,16 @@ TEST_CASE("custom_geometries", "[core]")
         Series write("../samples/custom_geometry.json", Access::CREATE);
         auto E = write.iterations[0].meshes["E"];
         E.setAttribute("geometry", "other:customGeometry");
+        // gridUnitDimension is technically an openPMD 2.0 addition, but since
+        // it's a non-breaking addition, we can also use it in openPMD 1.*
+        // files. However, it only really makes sense to use along with per-axis
+        // gridUnitSI definitions, which are in fact breaking in comparison to
+        // openPMD 1.*.
+        E.setGridUnitDimension(
+            {{{UnitDimension::theta, 1}},
+             {{UnitDimension::M, 1},
+              {UnitDimension::L, 1},
+              {UnitDimension::T, 2}}});
         auto E_x = E["x"];
         E_x.resetDataset({Datatype::INT, {10}});
         E_x.storeChunk(sampleData, {0}, {10});
@@ -1307,6 +1332,13 @@ TEST_CASE("custom_geometries", "[core]")
     {
         Series read("../samples/custom_geometry.json", Access::READ_ONLY);
         auto E = read.iterations[0].meshes["E"];
+        auto compare = unit_representations::AsMaps{
+            {{UnitDimension::theta, 1}},
+            {{UnitDimension::M, 1},
+             {UnitDimension::L, 1},
+             {UnitDimension::T, 2}}};
+        OPENPMD_REQUIRE_GUARD_WINDOWS(
+            unit_representations::asMaps(E.gridUnitDimension()) == compare);
         REQUIRE(
             E.getAttribute("geometry").get<std::string>() ==
             "other:customGeometry");
diff --git a/test/SerialIOTest.cpp b/test/SerialIOTest.cpp
index 3ff3f7572c..7374da7e30 100644
--- a/test/SerialIOTest.cpp
+++ b/test/SerialIOTest.cpp
@@ -2,6 +2,7 @@
 #include "openPMD/ChunkInfo_internal.hpp"
 #include "openPMD/Datatype.hpp"
 #include "openPMD/IO/Access.hpp"
+#include "openPMD/UnitDimension.hpp"
 #if openPMD_USE_INVASIVE_TESTS
 #define OPENPMD_private public:
 #define OPENPMD_protected public:
@@ -47,6 +48,21 @@
 #undef max
 #endif
 
+// On Windows, REQUIRE() might not be able to print more complex data structures
+// upon failure:
+// CoreTest.obj : error LNK2001: unresolved external symbol
+// "class std::string const Catch::Detail::unprintableString" (...)
+#ifdef _WIN32
+#define OPENPMD_REQUIRE_GUARD_WINDOWS(...)                                     \
+    do                                                                         \
+    {                                                                          \
+        bool guarded_require_boolean = __VA_ARGS__;                            \
+        REQUIRE(guarded_require_boolean);                                      \
+    } while (0);
+#else
+#define OPENPMD_REQUIRE_GUARD_WINDOWS(...) REQUIRE(__VA_ARGS__)
+#endif
+
 using namespace openPMD;
 
 struct BackendSelection
@@ -936,15 +952,28 @@ inline void constant_scalar(std::string const &file_ending)
 
         // store a number of predefined attributes in E
         Mesh &E_mesh = s.iterations[1].meshes["E"];
+        // test that these can be defined successively
+        E_mesh.setGridUnitDimension({{{UnitDimension::L, 1}}, {}, {}});
+        E_mesh.setGridUnitDimension(
+            {{}, {{UnitDimension::L, 1}}, {{UnitDimension::L, 1}}});
+        // let's modify the last dimension for the test's purpose now
+        E_mesh.setGridUnitDimension({{}, {}, {{UnitDimension::M, 1}}});
+        // should now all be [1,0,0,0,0,0,0], except the last which should be
+        // [1,1,0,0,0,0,0]
         E_mesh.setGeometry(geometry);
         E_mesh.setGeometryParameters(geometryParameters);
         E_mesh.setDataOrder(dataOrder);
         E_mesh.setGridSpacing(gridSpacing);
         E_mesh.setGridGlobalOffset(gridGlobalOffset);
-        E_mesh.setGridUnitSI(gridUnitSI);
         E_mesh.setAxisLabels(axisLabels);
         E_mesh.setUnitDimension(unitDimensions);
         E_mesh.setTimeOffset(timeOffset);
+        OPENPMD_REQUIRE_GUARD_WINDOWS(
+            E_mesh.gridUnitSIPerDimension() == std::vector<double>{1., 1., 1.});
+        E_mesh.setGridUnitSI(std::vector(3, gridUnitSI));
+        OPENPMD_REQUIRE_GUARD_WINDOWS(
+            E_mesh.gridUnitSIPerDimension() ==
+            std::vector<double>{gridUnitSI, gridUnitSI, gridUnitSI});
 
         // constant scalar
         auto pos =
@@ -1105,12 +1134,18 @@ inline void constant_scalar(std::string const &file_ending)
             Extent{3, 2, 1});
 
         Mesh &E_mesh = s.iterations[1].meshes["E"];
+        OPENPMD_REQUIRE_GUARD_WINDOWS(
+            E_mesh.gridUnitDimension() ==
+            std::vector{
+                std::array<double, 7>{1., 0., 0., 0., 0, .0, 0.},
+                std::array<double, 7>{1., 0., 0., 0., 0, .0, 0.},
+                std::array<double, 7>{1., 1., 0., 0., 0, .0, 0.}});
         REQUIRE(E_mesh.geometry() == geometry);
         REQUIRE(E_mesh.geometryParameters() == geometryParameters);
         REQUIRE(E_mesh.dataOrder() == dataOrder);
         REQUIRE(E_mesh.gridSpacing<double>() == gridSpacing);
         REQUIRE(E_mesh.gridGlobalOffset() == gridGlobalOffset);
-        REQUIRE(E_mesh.gridUnitSI() == gridUnitSI);
+        REQUIRE(E_mesh.gridUnitSIPerDimension() == std::vector(3, gridUnitSI));
         REQUIRE(E_mesh.axisLabels() == axisLabels);
         // REQUIRE( E_mesh.unitDimension() == unitDimensions );
         REQUIRE(E_mesh.timeOffset<double>() == timeOffset);
diff --git a/test/python/unittest/API/APITest.py b/test/python/unittest/API/APITest.py
index 07809760e5..be1609d577 100644
--- a/test/python/unittest/API/APITest.py
+++ b/test/python/unittest/API/APITest.py
@@ -100,6 +100,17 @@ def testFieldData(self):
         self.assertEqual(len(i.meshes), 2)
         for m in i.meshes:
             self.assertTrue(m in ["E", "rho"])
+            self.assertEqual(
+                i.meshes[m].unit_dimension,
+                io.Unit_Dimension.as_array(
+                    io.Unit_Dimension.as_map(i.meshes[m].unit_dimension)))
+            self.assertEqual(
+                io.Unit_Dimension.as_maps(i.meshes[m].grid_unit_dimension),
+                [{io.Unit_Dimension.L: 1}, {io.Unit_Dimension.L: 1}])
+        self.assertEqual(io.Unit_Dimension.from_index(0), io.Unit_Dimension.L)
+        self.assertEqual(io.Unit_Dimension.L.as_index(), 0)
+        for idx in range(7):
+            self.assertEqual(idx, io.Unit_Dimension.from_index(idx).as_index())
 
         # Check entries.
         self.assertEqual(len(i.meshes), 2)
@@ -403,6 +414,97 @@ def testAttributes(self):
         for ext in tested_file_extensions:
             self.attributeRoundTrip(ext)
 
+    def testOpenPMD_2_0(self):
+        write_2_0 = io.Series("../samples/openpmd_2_0.json", io.Access.create)
+        write_2_0.openPMD = "2.0.0"
+        meshes = write_2_0.write_iterations()[100].meshes
+
+        E = meshes["E"]
+        E.reset_dataset(io.Dataset(io.Datatype.DOUBLE, [10, 10, 10]))
+        E.axis_labels = ["x", "y", "z"]
+        E.grid_unit_SI = [1, 2, 3]
+        E.grid_unit_dimension = [
+            {io.Unit_Dimension.L: 1},
+            {io.Unit_Dimension.L: 1},
+            {io.Unit_Dimension.L: 1, io.Unit_Dimension.T: -1}]
+        E.make_constant(17)
+
+        B = meshes["B"]
+        B.reset_dataset(io.Dataset(io.Datatype.DOUBLE, [10, 10, 10]))
+        B.axis_labels = ["x", "y", "z"]
+        # This is deprecated for openPMD 2.0, a warning will be printed.
+        B.grid_unit_SI = 3
+        B.make_constant(18)
+
+        write_2_0.close()
+
+        read_2_0 = io.Series(
+            "../samples/openpmd_2_0.json", io.Access.read_only)
+        meshes = read_2_0.iterations[100].meshes
+
+        E = meshes["E"]
+        self.assertEqual(E.grid_unit_SI, [1, 2, 3])
+        self.assertEqual(E.grid_unit_dimension, io.Unit_Dimension.as_arrays([
+            {io.Unit_Dimension.L: 1},
+            {io.Unit_Dimension.L: 1},
+            {io.Unit_Dimension.L: 1, io.Unit_Dimension.T: -1}]))
+
+        B = meshes["B"]
+        # Will return a list due to openPMD standard being set to 2.0.0
+        self.assertEqual(B.grid_unit_SI, [3])
+        # If the attribute is not defined, the mesh is implicitly spatial
+        self.assertEqual(io.Unit_Dimension.as_maps(B.grid_unit_dimension), [
+                         {io.Unit_Dimension.L: 1} for _ in range(3)])
+        read_2_0.close()
+
+        write_1_1 = io.Series("../samples/openpmd_1_1.json", io.Access.create)
+        write_1_1.openPMD = "1.1.0"
+        meshes = write_1_1.write_iterations()[100].meshes
+
+        E = meshes["E"]
+        E.reset_dataset(io.Dataset(io.Datatype.DOUBLE, [10, 10, 10]))
+
+        def unsupported_in_1_1():
+            E.grid_unit_SI = [1, 2, 3]
+        # self.assertRaises(
+        #     io.ErrorIllegalInOpenPMDStandard, unsupported_in_1_1)
+        E.axis_labels = ["x", "y", "z"]
+        E.grid_unit_dimension = [
+            {io.Unit_Dimension.L: 1},
+            {io.Unit_Dimension.L: 1},
+            {io.Unit_Dimension.L: 1, io.Unit_Dimension.T: -1}]
+        E.make_constant(17)
+
+        B = meshes["B"]
+        B.reset_dataset(io.Dataset(io.Datatype.DOUBLE, [10, 10, 10]))
+        B.axis_labels = ["x", "y", "z"]
+        # This is deprecated for openPMD 2.0, a warning will be printed.
+        B.grid_unit_SI = 3
+        B.make_constant(18)
+
+        write_1_1.close()
+
+        read_1_1 = io.Series(
+            "../samples/openpmd_1_1.json", io.Access.read_only)
+        meshes = read_1_1.iterations[100].meshes
+
+        E = meshes["E"]
+        # Will return a default value due to the failed attempt at setting
+        #  a list at write time
+        self.assertEqual(E.grid_unit_SI, 1)
+        self.assertEqual(E.grid_unit_dimension, io.Unit_Dimension.as_arrays([
+            {io.Unit_Dimension.L: 1},
+            {io.Unit_Dimension.L: 1},
+            {io.Unit_Dimension.L: 1, io.Unit_Dimension.T: -1}]))
+
+        B = meshes["B"]
+        # Will return a scalar due to openPMD standard being set to 2.0.0
+        self.assertEqual(B.grid_unit_SI, 3)
+        # If the attribute is not defined, the mesh is implicitly spatial
+        self.assertEqual(io.Unit_Dimension.as_maps(B.grid_unit_dimension), [
+                         {io.Unit_Dimension.L: 1} for _ in range(3)])
+        read_1_1.close()
+
     def makeConstantRoundTrip(self, file_ending):
         # write
         series = io.Series(