diff --git a/include/openPMD/Mesh.hpp b/include/openPMD/Mesh.hpp
index edfb86680f..dd830a5613 100644
--- a/include/openPMD/Mesh.hpp
+++ b/include/openPMD/Mesh.hpp
@@ -184,9 +184,8 @@ class Mesh : public BaseRecord<MeshRecordComponent>
      * Mesh::gridSpacing and Mesh::gridGlobalOffset, in order to convert from
      * simulation units to SI units.
      *
-     * Note that this API follows the legacy (openPMD 1.*) definition
-     * of gridUnitSI.
-     * In order to specify the gridUnitSI per dimension,
+     * 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
@@ -196,12 +195,14 @@ class Mesh : public BaseRecord<MeshRecordComponent>
      */
     Mesh &setGridUnitSI(double gridUnitSI);
 
-    /** Set the unit-conversion factor per dimension to multiply each value in
+    /** Alias for `setGridUnitSI(std::vector<double>)`.
+     *
+     * Set the unit-conversion factor per dimension to multiply each value in
      * Mesh::gridSpacing and Mesh::gridGlobalOffset, in order to convert from
      * simulation units to SI units.
      *
-     * Note that this is a feature of openPMD 2.0.
-     * The legacy behavior (a scalar gridUnitSI) is implemented
+     * 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
@@ -220,14 +221,12 @@ class Mesh : public BaseRecord<MeshRecordComponent>
      */
     std::vector<double> gridUnitSIPerDimension() const;
 
-    /** Alias for `setGridUnitSI(std::vector<double>)`.
-     *
-     * Set the unit-conversion factor per dimension to multiply each value in
+    /* Set the unit-conversion factor per dimension to multiply each value in
      * Mesh::gridSpacing and Mesh::gridGlobalOffset, in order to convert from
      * simulation units to SI units.
      *
-     * Note that this is a feature of openPMD 2.0.
-     * The legacy behavior (a scalar gridUnitSI) is implemented
+     * 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
diff --git a/src/Mesh.cpp b/src/Mesh.cpp
index 97e8fa4387..58c3450c69 100644
--- a/src/Mesh.cpp
+++ b/src/Mesh.cpp
@@ -42,7 +42,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
@@ -405,6 +404,18 @@ void Mesh::flush_impl(
                     comp.second.flush(comp.first, flushParams);
             }
         }
+        if (!containsAttribute("gridUnitSI"))
+        {
+            if (auto series = retrieveSeries(); series.openPMD() < "2.")
+            {
+                setGridUnitSI(1);
+            }
+            else
+            {
+                setGridUnitSIPerDimension(
+                    std::vector<double>(retrieveMeshDimensionality(*this), 1));
+            }
+        }
         flushAttributes(flushParams);
     }
 }
diff --git a/test/CoreTest.cpp b/test/CoreTest.cpp
index 62a2012d41..fd9429b02b 100644
--- a/test/CoreTest.cpp
+++ b/test/CoreTest.cpp
@@ -528,11 +528,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 +556,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);