diff --git a/NEWS.md b/NEWS.md
index b870d35f0a..004b157664 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -4,6 +4,22 @@ ClimaCore.jl Release Notes
 main
 -------
 
+v0.14.16
+-------
+
+### ![][badge-🐛bugfix] Fix restarting simulations from `Space`s with `enable_bubble = true`
+
+Prior to this change, the `ClimaCore.InputOutput` module did not save whether a
+`Space` was constructed with `enable_bubble = true`. This meant that restarting
+a simulation from a HDF5 file led to inconsistent and incorrect spaces and
+`Field`s. This affected only 2D spectral spaces (and extruded ones that have
+this type of horizontal space).
+
+We now expect `Space`s read from a file to be bitwise identical to the original
+one.
+
+PR [#1999](https://github.com/CliMA/ClimaCore.jl/pull/1999).
+
 v0.14.15
 -------
 
diff --git a/Project.toml b/Project.toml
index 106e10a491..95fe874ec7 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "ClimaCore"
 uuid = "d414da3d-4745-48bb-8d80-42e94e092884"
 authors = ["CliMA Contributors <clima-software@caltech.edu>"]
-version = "0.14.15"
+version = "0.14.16"
 
 [deps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
diff --git a/src/Grids/spectralelement.jl b/src/Grids/spectralelement.jl
index c880c99829..c021421a54 100644
--- a/src/Grids/spectralelement.jl
+++ b/src/Grids/spectralelement.jl
@@ -129,6 +129,7 @@ mutable struct SpectralElementGrid2D{
     local_dss_weights::D
     internal_surface_geometry::IS
     boundary_surface_geometries::BS
+    enable_bubble::Bool
 end
 
 local_geometry_type(
@@ -467,6 +468,7 @@ function _SpectralElementGrid2D(
         dss_local_weights,
         internal_surface_geometry,
         boundary_surface_geometries,
+        enable_bubble,
     )
 end
 
diff --git a/src/InputOutput/readers.jl b/src/InputOutput/readers.jl
index fd8b28c197..744699ecf4 100644
--- a/src/InputOutput/readers.jl
+++ b/src/InputOutput/readers.jl
@@ -333,10 +333,15 @@ function read_grid_new(reader, name)
         quadrature_style =
             _scan_quadrature_style(attrs(group)["quadrature_type"], npts)
         topology = read_topology(reader, attrs(group)["topology"])
+        enable_bubble = get(attrs(group), "bubble", "false") == "true"
         if type == "SpectralElementGrid1D"
             return Grids.SpectralElementGrid1D(topology, quadrature_style)
         else
-            return Grids.SpectralElementGrid2D(topology, quadrature_style)
+            return Grids.SpectralElementGrid2D(
+                topology,
+                quadrature_style;
+                enable_bubble,
+            )
         end
     elseif type == "FiniteDifferenceGrid"
         topology = read_topology(reader, attrs(group)["topology"])
diff --git a/src/InputOutput/writers.jl b/src/InputOutput/writers.jl
index fc78dc9bb6..330f869790 100644
--- a/src/InputOutput/writers.jl
+++ b/src/InputOutput/writers.jl
@@ -343,6 +343,7 @@ function write_new!(
         "quadrature_num_points",
         Quadratures.degrees_of_freedom(Spaces.quadrature_style(grid)),
     )
+    write_attribute(group, "bubble", grid.enable_bubble ? "true" : "false")
     write_attribute(group, "topology", write!(writer, Spaces.topology(grid)))
     return name
 end
diff --git a/test/InputOutput/spectralelement2d.jl b/test/InputOutput/spectralelement2d.jl
index 4dfb6048f8..ae0226f51b 100644
--- a/test/InputOutput/spectralelement2d.jl
+++ b/test/InputOutput/spectralelement2d.jl
@@ -71,20 +71,27 @@ end
     @info "Using device" device
     context = ClimaComms.SingletonCommsContext(device)
     grid_topology = Topologies.Topology2D(context, mesh)
-    space = Spaces.SpectralElementSpace2D(grid_topology, quad)
 
-    y0 = init_state.(Fields.local_geometry_field(space), Ref(parameters))
-    Y = Fields.FieldVector(y0 = y0)
+    for enable_bubble in (true, false)
+        space =
+            Spaces.SpectralElementSpace2D(grid_topology, quad; enable_bubble)
 
-    # write field vector to hdf5 file
-    filename = tempname(pwd())
-    writer = InputOutput.HDF5Writer(filename, context)
-    InputOutput.write!(writer, "Y" => Y) # write field vector from hdf5 file
-    close(writer)
-    reader = InputOutput.HDF5Reader(filename, context)
-    restart_Y = InputOutput.read_field(reader, "Y") # read fieldvector from hdf5 file
-    close(reader)
-    ClimaComms.allowscalar(device) do
-        @test restart_Y == Y # test if restart is exact
+        y0 = init_state.(Fields.local_geometry_field(space), Ref(parameters))
+        Y = Fields.FieldVector(y0 = y0)
+
+        # write field vector to hdf5 file
+        filename = tempname(pwd())
+        writer = InputOutput.HDF5Writer(filename, context)
+        InputOutput.write!(writer, "Y" => Y) # write field vector from hdf5 file
+        close(writer)
+        reader = InputOutput.HDF5Reader(filename, context)
+        restart_Y = InputOutput.read_field(reader, "Y") # read fieldvector from hdf5 file
+        close(reader)
+        ClimaComms.allowscalar(device) do
+            @test restart_Y == Y # test if restart is exact
+        end
+        ClimaComms.allowscalar(device) do
+            @test axes(restart_Y) == axes(Y) # test if restart is exact for space
+        end
     end
 end