diff --git a/docs/source/io_formats/index.rst b/docs/source/io_formats/index.rst
index 7ae22c0a33..8c89bd2a1d 100644
--- a/docs/source/io_formats/index.rst
+++ b/docs/source/io_formats/index.rst
@@ -45,6 +45,7 @@ Output Files
    statepoint
    source
    summary
+   properties
    depletion_results
    particle_restart
    track
diff --git a/docs/source/io_formats/properties.rst b/docs/source/io_formats/properties.rst
new file mode 100644
index 0000000000..5030e78f35
--- /dev/null
+++ b/docs/source/io_formats/properties.rst
@@ -0,0 +1,36 @@
+.. _io_properties:
+
+======================
+Properties File Format
+======================
+
+The current version of the properties file format is 1.0.
+
+**/**
+
+:Attributes: - **filetype** (*char[]*) -- String indicating the type of file.
+             - **version** (*int[2]*) -- Major and minor version of the
+               statepoint file format.
+             - **openmc_version** (*int[3]*) -- Major, minor, and release
+               version number for OpenMC.
+             - **git_sha1** (*char[40]*) -- Git commit SHA-1 hash.
+             - **date_and_time** (*char[]*) -- Date and time the summary was
+               written.
+             - **path** (*char[]*) -- Path to directory containing input files.
+
+**/geometry/**
+
+:Attributes: - **n_cells** (*int*) -- Number of cells in the problem.
+
+**/geometry/cells/cell <uid>/**
+
+:Datasets: - **temperature** (*double[]*) -- Temperature of the cell in [K].
+
+**/materials/**
+
+:Attributes: - **n_materials** (*int*) -- Number of materials in the problem.
+
+**/materials/material <uid>/**
+
+:Attributes: - **atom_density** (*double*) -- Total density in [atom/b-cm].
+             - **mass_density** (*double*) -- Total density in [g/cm^3].
diff --git a/include/openmc/capi.h b/include/openmc/capi.h
index 2ca9fc430f..f08ec5c02f 100644
--- a/include/openmc/capi.h
+++ b/include/openmc/capi.h
@@ -15,6 +15,7 @@ extern "C" {
   int openmc_cell_get_id(int32_t index, int32_t* id);
   int openmc_cell_get_temperature(int32_t index, const int32_t* instance, double* T);
   int openmc_cell_get_name(int32_t index, const char** name);
+  int openmc_cell_get_num_instances(int32_t index, int32_t* num_instances);
   int openmc_cell_set_name(int32_t index, const char* name);
   int openmc_cell_set_fill(int32_t index, int type, int32_t n, const int32_t* indices);
   int openmc_cell_set_id(int32_t index, int32_t id);
@@ -144,13 +145,13 @@ extern "C" {
   //! \param[in] meshtyally_id id of CMFD Mesh Tally
   //! \param[in] cmfd_indices indices storing spatial and energy dimensions of CMFD problem
   //! \param[in] norm CMFD normalization factor
-  extern "C" void openmc_initialize_mesh_egrid(const int meshtally_id, const int* cmfd_indices,
-                                               const double norm);
+  void openmc_initialize_mesh_egrid(const int meshtally_id, const int* cmfd_indices,
+                                    const double norm);
 
   //! Sets the mesh and energy grid for CMFD reweight
   //! \param[in] feedback whether or not to run CMFD feedback
   //! \param[in] cmfd_src computed CMFD source
-  extern "C" void openmc_cmfd_reweight(const bool feedback, const double* cmfd_src);
+  void openmc_cmfd_reweight(const bool feedback, const double* cmfd_src);
 
   //! Sets the fixed variables that are used for CMFD linear solver
   //! \param[in] indptr CSR format index pointer array of loss matrix
@@ -161,10 +162,10 @@ extern "C" {
   //! \param[in] spectral spectral radius of CMFD matrices and tolerances
   //! \param[in] map coremap for problem, storing accelerated regions
   //! \param[in] use_all_threads whether to use all threads when running CMFD solver
-  extern "C" void openmc_initialize_linsolver(const int* indptr, int len_indptr,
-                                              const int* indices, int n_elements,
-                                              int dim, double spectral,
-                                              const int* map, bool use_all_threads);
+  void openmc_initialize_linsolver(const int* indptr, int len_indptr,
+                                   const int* indices, int n_elements,
+                                   int dim, double spectral,
+                                   const int* map, bool use_all_threads);
 
   //! Runs a Gauss Seidel linear solver to solve CMFD matrix equations
   //! linear solver
@@ -173,9 +174,19 @@ extern "C" {
   //! \param[out] x unknown vector
   //! \param[in] tol tolerance on final error
   //! \return number of inner iterations required to reach convergence
-  extern "C" int openmc_run_linsolver(const double* A_data, const double* b,
+  int openmc_run_linsolver(const double* A_data, const double* b,
                                       double* x, double tol);
 
+  //! Export physical properties for model
+  //! \param[in] filename Filename to write to
+  //! \return Error code
+  int openmc_properties_export(const char* filename);
+
+  //! Import physical properties for model
+  //! \param[in] filename Filename to read from
+  // \return Error code
+  int openmc_properties_import(const char* filename);
+
   // Error codes
   extern int OPENMC_E_UNASSIGNED;
   extern int OPENMC_E_ALLOCATE;
diff --git a/include/openmc/cell.h b/include/openmc/cell.h
index 9e977bc937..634b489de2 100644
--- a/include/openmc/cell.h
+++ b/include/openmc/cell.h
@@ -130,6 +130,14 @@ public:
 
   virtual void to_hdf5_inner(hid_t group_id) const = 0;
 
+  //! Export physical properties to HDF5
+  //! \param[in] group  HDF5 group to read from
+  void export_properties_hdf5(hid_t group) const;
+
+  //! Import physical properties from HDF5
+  //! \param[in] group  HDF5 group to write to
+  void import_properties_hdf5(hid_t group);
+
   //! Get the BoundingBox for this cell.
   virtual BoundingBox bounding_box() const = 0;
 
diff --git a/include/openmc/constants.h b/include/openmc/constants.h
index 89afadeb43..e400c22bca 100644
--- a/include/openmc/constants.h
+++ b/include/openmc/constants.h
@@ -31,6 +31,7 @@ constexpr array<int, 2> VERSION_SUMMARY {6, 0};
 constexpr array<int, 2> VERSION_VOLUME {1, 0};
 constexpr array<int, 2> VERSION_VOXEL {2, 0};
 constexpr array<int, 2> VERSION_MGXS_LIBRARY {1, 0};
+constexpr array<int, 2> VERSION_PROPERTIES {1, 0};
 
 // ============================================================================
 // ADJUSTABLE PARAMETERS
diff --git a/include/openmc/hdf5_interface.h b/include/openmc/hdf5_interface.h
index f601673041..ae91ce2ae9 100644
--- a/include/openmc/hdf5_interface.h
+++ b/include/openmc/hdf5_interface.h
@@ -51,7 +51,7 @@ inline hid_t create_group(hid_t parent_id, const std::stringstream& name)
 
 
 hid_t file_open(const std::string& filename, char mode, bool parallel=false);
-
+hid_t open_group(hid_t group_id, const std::string& name);
 void write_string(hid_t group_id, const char* name, const std::string& buffer,
                   bool indep);
 
diff --git a/include/openmc/material.h b/include/openmc/material.h
index 7e975a1bb6..f0d8103707 100644
--- a/include/openmc/material.h
+++ b/include/openmc/material.h
@@ -69,6 +69,14 @@ public:
   //! Write material data to HDF5
   void to_hdf5(hid_t group) const;
 
+  //! Export physical properties to HDF5
+  //! \param[in] group  HDF5 group to write to
+  void export_properties_hdf5(hid_t group) const;
+
+  //! Import physical properties from HDF5
+  //! \param[in] group  HDF5 group to read from
+  void import_properties_hdf5(hid_t group);
+
   //! Add nuclide to the material
   //
   //! \param[in] nuclide Name of the nuclide
diff --git a/openmc/__init__.py b/openmc/__init__.py
index 50246cdd8a..8985099ba8 100644
--- a/openmc/__init__.py
+++ b/openmc/__init__.py
@@ -31,7 +31,7 @@ from openmc.search import *
 from openmc.polynomial import *
 from . import examples
 
-# Import a few convencience functions that used to be here
-from openmc.model import rectangular_prism, hexagonal_prism
+# Import a few names from the model module
+from openmc.model import rectangular_prism, hexagonal_prism, Model
 
 __version__ = '0.13.0-dev'
diff --git a/openmc/lib/cell.py b/openmc/lib/cell.py
index c69d1e2a40..7fdc4542db 100644
--- a/openmc/lib/cell.py
+++ b/openmc/lib/cell.py
@@ -25,6 +25,9 @@ _dll.openmc_cell_get_fill.argtypes = [
     c_int32, POINTER(c_int), POINTER(POINTER(c_int32)), POINTER(c_int32)]
 _dll.openmc_cell_get_fill.restype = c_int
 _dll.openmc_cell_get_fill.errcheck = _error_handler
+_dll.openmc_cell_get_num_instances.argtypes = [c_int32, POINTER(c_int32)]
+_dll.openmc_cell_get_num_instances.restype = c_int
+_dll.openmc_cell_get_num_instances.errcheck = _error_handler
 _dll.openmc_cell_get_temperature.argtypes = [
     c_int32, POINTER(c_int32), POINTER(c_double)]
 _dll.openmc_cell_get_temperature.restype = c_int
@@ -78,6 +81,14 @@ class Cell(_FortranObjectWithID):
     ----------
     id : int
         ID of the cell
+    fill : openmc.lib.Material or list of openmc.lib.Material
+        Indicates what the region of space is filled with
+    name : str
+        Name of the cell
+    num_instances : int
+        Number of unique cell instances
+    bounding_box : 2-tuple of numpy.ndarray
+        Lower-left and upper-right coordinates of bounding box
 
     """
     __instances = WeakValueDictionary()
@@ -159,6 +170,12 @@ class Cell(_FortranObjectWithID):
             indices = (c_int32*1)(-1)
             _dll.openmc_cell_set_fill(self._index, 0, 1, indices)
 
+    @property
+    def num_instances(self):
+        n = c_int32()
+        _dll.openmc_cell_get_num_instances(self._index, n)
+        return n.value
+
     def get_temperature(self, instance=None):
         """Get the temperature of a cell
 
@@ -211,6 +228,7 @@ class Cell(_FortranObjectWithID):
 
         return llc, urc
 
+
 class _CellMapping(Mapping):
     def __getitem__(self, key):
         index = c_int32()
diff --git a/openmc/lib/core.py b/openmc/lib/core.py
index e740e28699..0da24aab88 100644
--- a/openmc/lib/core.py
+++ b/openmc/lib/core.py
@@ -62,7 +62,13 @@ _dll.openmc_next_batch.argtypes = [POINTER(c_int)]
 _dll.openmc_next_batch.restype = c_int
 _dll.openmc_next_batch.errcheck = _error_handler
 _dll.openmc_plot_geometry.restype = c_int
-_dll.openmc_plot_geometry.restype = _error_handler
+_dll.openmc_plot_geometry.errcheck = _error_handler
+_dll.openmc_properties_export.argtypes = [c_char_p]
+_dll.openmc_properties_export.restype = c_int
+_dll.openmc_properties_export.errcheck = _error_handler
+_dll.openmc_properties_import.argtypes = [c_char_p]
+_dll.openmc_properties_import.restype = c_int
+_dll.openmc_properties_import.errcheck = _error_handler
 _dll.openmc_run.restype = c_int
 _dll.openmc_run.errcheck = _error_handler
 _dll.openmc_reset.restype = c_int
@@ -120,6 +126,24 @@ def current_batch():
     return c_int.in_dll(_dll, 'current_batch').value
 
 
+def export_properties(filename=None):
+    """Export physical properties.
+
+    Parameters
+    ----------
+    filename : str or None
+        Filename to export properties to (defaults to "properties.h5")
+
+    See Also
+    --------
+    openmc.lib.import_properties
+
+    """
+    if filename is not None:
+        filename = c_char_p(filename.encode())
+    _dll.openmc_properties_export(filename)
+
+
 def finalize():
     """Finalize simulation and free memory"""
     _dll.openmc_finalize()
@@ -178,6 +202,22 @@ def hard_reset():
     _dll.openmc_hard_reset()
 
 
+def import_properties(filename):
+    """Import physical properties.
+
+    Parameters
+    ----------
+    filename : str
+        Filename to import properties from
+
+    See Also
+    --------
+    openmc.lib.export_properties
+
+    """
+    _dll.openmc_properties_import(filename.encode())
+
+
 def init(args=None, intracomm=None):
     """Initialize OpenMC
 
diff --git a/openmc/lib/material.py b/openmc/lib/material.py
index 6f80abb5e2..6770721ab6 100644
--- a/openmc/lib/material.py
+++ b/openmc/lib/material.py
@@ -174,15 +174,6 @@ class Material(_FortranObjectWithID):
         return self._get_densities()[0]
         return nuclides
 
-    @property
-    def density(self):
-      density = c_double()
-      try:
-          _dll.openmc_material_get_density(self._index, density)
-      except OpenMCError:
-          return None
-      return density.value
-
     @property
     def densities(self):
         return self._get_densities()[1]
@@ -224,6 +215,29 @@ class Material(_FortranObjectWithID):
         """
         _dll.openmc_material_add_nuclide(self._index, name.encode(), density)
 
+    def get_density(self, units='atom/b-cm'):
+        """Get density of a material.
+
+        Parameters
+        ----------
+        units : {'atom/b-cm', 'g/cm3'}
+            Units for density
+
+        Returns
+        -------
+        float
+            Density in requested units
+
+        """
+        if units == 'atom/b-cm':
+            return self.densities.sum()
+        elif units == 'g/cm3':
+            density = c_double()
+            _dll.openmc_material_get_density(self._index, density)
+            return density.value
+        else:
+            raise ValueError("Units must be 'atom/b-cm' or 'g/cm3'")
+
     def set_density(self, density, units='atom/b-cm'):
         """Set density of a material.
 
diff --git a/openmc/lib/plot.py b/openmc/lib/plot.py
index c51000e0de..da42ea49e0 100644
--- a/openmc/lib/plot.py
+++ b/openmc/lib/plot.py
@@ -228,11 +228,11 @@ def id_map(plot):
     Returns
     -------
     id_map : numpy.ndarray
-        A NumPy array with shape (vertical pixels, horizontal pixels, 2) of
+        A NumPy array with shape (vertical pixels, horizontal pixels, 3) of
         OpenMC property ids with dtype int32
 
     """
-    img_data = np.zeros((plot.v_res, plot.h_res, 2),
+    img_data = np.zeros((plot.v_res, plot.h_res, 3),
                         dtype=np.dtype('int32'))
     _dll.openmc_id_map(plot, img_data.ctypes.data_as(POINTER(c_int32)))
     return img_data
diff --git a/openmc/model/model.py b/openmc/model/model.py
index 0501f03f3a..d6c91e73e1 100644
--- a/openmc/model/model.py
+++ b/openmc/model/model.py
@@ -2,6 +2,8 @@ from collections.abc import Iterable
 from pathlib import Path
 import time
 
+import h5py
+
 import openmc
 from openmc.checkvalue import check_type, check_value
 
@@ -11,11 +13,11 @@ class Model:
 
     This class can be used to store instances of :class:`openmc.Geometry`,
     :class:`openmc.Materials`, :class:`openmc.Settings`,
-    :class:`openmc.Tallies`, :class:`openmc.Plots`, and :class:`openmc.CMFD`,
-    thus making a complete model. The :meth:`Model.export_to_xml` method will
-    export XML files for all attributes that have been set. If the
-    :meth:`Model.materials` attribute is not set, it will attempt to create a
-    ``materials.xml`` file based on all materials appearing in the geometry.
+    :class:`openmc.Tallies`, and :class:`openmc.Plots`, thus making a complete
+    model. The :meth:`Model.export_to_xml` method will export XML files for all
+    attributes that have been set. If the :attr:`Model.materials` attribute is
+    not set, it will attempt to create a ``materials.xml`` file based on all
+    materials appearing in the geometry.
 
     Parameters
     ----------
@@ -124,6 +126,31 @@ class Model:
             for plot in plots:
                 self._plots.append(plot)
 
+    @classmethod
+    def from_xml(cls, geometry='geometry.xml', materials='materials.xml',
+                 settings='settings.xml'):
+        """Create model from existing XML files
+
+        Parameters
+        ----------
+        geometry : str
+            Path to geometry.xml file
+        materials : str
+            Path to materials.xml file
+        settings : str
+            Path to settings.xml file
+
+        Returns
+        -------
+        openmc.model.Model
+            Model created from XML files
+
+        """
+        materials = openmc.Materials.from_xml(materials)
+        geometry = openmc.Geometry.from_xml(geometry, materials)
+        settings = openmc.Settings.from_xml(settings)
+        return cls(geometry, materials, settings)
+
     def deplete(self, timesteps, chain_file=None, method='cecm',
                 fission_q=None, **kwargs):
         """Deplete model using specified timesteps/power
@@ -196,6 +223,51 @@ class Model:
         if self.plots:
             self.plots.export_to_xml(d)
 
+    def import_properties(self, filename):
+        """Import physical properties
+
+        Parameters
+        ----------
+        filename : str
+            Path to properties HDF5 file
+
+        See Also
+        --------
+        openmc.lib.export_properties
+
+        """
+        cells = self.geometry.get_all_cells()
+        materials = self.geometry.get_all_materials()
+
+        with h5py.File(filename, 'r') as fh:
+            cells_group = fh['geometry/cells']
+
+            # Make sure number of cells matches
+            n_cells = fh['geometry'].attrs['n_cells']
+            if n_cells != len(cells):
+                raise ValueError("Number of cells in properties file doesn't "
+                                 "match current model.")
+
+            # Update temperatures for cells filled with materials
+            for name, group in cells_group.items():
+                cell_id = int(name.split()[1])
+                cell = cells[cell_id]
+                if cell.fill_type in ('material', 'distribmat'):
+                    cell.temperature = group['temperature'][()]
+
+            # Make sure number of materials matches
+            mats_group = fh['materials']
+            n_cells = mats_group.attrs['n_materials']
+            if n_cells != len(materials):
+                raise ValueError("Number of materials in properties file doesn't "
+                                 "match current model.")
+
+            # Update material densities
+            for name, group in mats_group.items():
+                mat_id = int(name.split()[1])
+                atom_density = group.attrs['atom_density']
+                materials[mat_id].set_density('atom/b-cm', atom_density)
+
     def run(self, **kwargs):
         """Creates the XML files, runs OpenMC, and returns the path to the last
         statepoint file generated.
diff --git a/src/cell.cpp b/src/cell.cpp
index 4ae5177fa0..a9e9cc3a5e 100644
--- a/src/cell.cpp
+++ b/src/cell.cpp
@@ -309,6 +309,46 @@ Cell::set_temperature(double T, int32_t instance, bool set_contained)
   }
 }
 
+void Cell::export_properties_hdf5(hid_t group) const
+{
+  // Create a group for this cell.
+  auto cell_group = create_group(group, fmt::format("cell {}", id_));
+
+  // Write temperature in [K] for one or more cell instances
+  vector<double> temps;
+  for (auto sqrtkT_val : sqrtkT_)
+    temps.push_back(sqrtkT_val * sqrtkT_val / K_BOLTZMANN);
+  write_dataset(cell_group, "temperature", temps);
+
+  close_group(cell_group);
+}
+
+void Cell::import_properties_hdf5(hid_t group)
+{
+  auto cell_group = open_group(group, fmt::format("cell {}", id_));
+
+  // Read temperatures from file
+  vector<double> temps;
+  read_dataset(cell_group, "temperature", temps);
+
+  // Ensure number of temperatures makes sense
+  auto n_temps = temps.size();
+  if (n_temps > 1 && n_temps != n_instances_) {
+    throw std::runtime_error(fmt::format(
+      "Number of temperatures for cell {} doesn't match number of instances", id_));
+  }
+
+  // Modify temperatures for the cell
+  sqrtkT_.clear();
+  sqrtkT_.resize(temps.size());
+  for (gsl::index i = 0; i < temps.size(); ++i) {
+    this->set_temperature(temps[i], i);
+  }
+
+  close_group(cell_group);
+}
+
+
 void
 Cell::to_hdf5(hid_t cell_group) const {
 
@@ -1259,6 +1299,17 @@ openmc_cell_set_id(int32_t index, int32_t id)
   }
 }
 
+extern "C" int
+openmc_cell_get_num_instances(int32_t index, int32_t* num_instances)
+{
+  if (index < 0 || index >= model::cells.size()) {
+    set_errmsg("Index in cells array is out of bounds.");
+    return OPENMC_E_OUT_OF_BOUNDS;
+  }
+  *num_instances = model::cells[index]->n_instances_;
+  return 0;
+}
+
 //! Extend the cells array by n elements
 extern "C" int
 openmc_extend_cells(int32_t n, int32_t* index_start, int32_t* index_end)
diff --git a/src/hdf5_interface.cpp b/src/hdf5_interface.cpp
index 0f8db2dc5f..522c9c2625 100644
--- a/src/hdf5_interface.cpp
+++ b/src/hdf5_interface.cpp
@@ -225,6 +225,11 @@ file_open(const std::string& filename, char mode, bool parallel)
   return file_open(filename.c_str(), mode, parallel);
 }
 
+hid_t open_group(hid_t group_id, const std::string& name)
+{
+  return open_group(group_id, name.c_str());
+}
+
 void file_close(hid_t file_id)
 {
   H5Fclose(file_id);
diff --git a/src/material.cpp b/src/material.cpp
index 4e6dfbdf1b..c942b10a61 100644
--- a/src/material.cpp
+++ b/src/material.cpp
@@ -1052,6 +1052,23 @@ void Material::to_hdf5(hid_t group) const
   close_group(material_group);
 }
 
+void Material::export_properties_hdf5(hid_t group) const
+{
+  hid_t material_group = create_group(group, "material " + std::to_string(id_));
+  write_attribute(material_group, "atom_density", density_);
+  write_attribute(material_group, "mass_density", density_gpcc_);
+  close_group(material_group);
+}
+
+void Material::import_properties_hdf5(hid_t group)
+{
+  hid_t material_group = open_group(group, "material " + std::to_string(id_));
+  double density;
+  read_attribute(material_group, "atom_density", density);
+  this->set_density(density, "atom/b-cm");
+  close_group(material_group);
+}
+
 void Material::add_nuclide(const std::string& name, double density)
 {
   // Check if nuclide is already in material
diff --git a/src/plot.cpp b/src/plot.cpp
index 615701341e..999cf3d294 100644
--- a/src/plot.cpp
+++ b/src/plot.cpp
@@ -36,7 +36,7 @@ constexpr int32_t NOT_FOUND {-2};
 constexpr int32_t OVERLAP {-3};
 
 IdData::IdData(size_t h_res, size_t v_res)
-  : data_({v_res, h_res, 2}, NOT_FOUND)
+  : data_({v_res, h_res, 3}, NOT_FOUND)
 { }
 
 void
@@ -44,18 +44,20 @@ IdData::set_value(size_t y, size_t x, const Particle& p, int level) {
   // set cell data
   if (p.n_coord() <= level) {
     data_(y, x, 0) = NOT_FOUND;
+    data_(y, x, 1) = NOT_FOUND;
   } else {
     data_(y, x, 0) = model::cells.at(p.coord(level).cell)->id_;
+    data_(y, x, 1) = p.cell_instance();
   }
 
   // set material data
   Cell* c = model::cells.at(p.coord(p.n_coord() - 1).cell).get();
   if (p.material() == MATERIAL_VOID) {
-    data_(y, x, 1) = MATERIAL_VOID;
+    data_(y, x, 2) = MATERIAL_VOID;
     return;
   } else if (c->type_ == Fill::MATERIAL) {
     Material* m = model::materials.at(p.material()).get();
-    data_(y, x, 1) = m->id_;
+    data_(y, x, 2) = m->id_;
   }
 }
 
@@ -155,7 +157,8 @@ void create_ppm(Plot const& pl)
   // assign colors
   for (size_t y = 0; y < height; y++) {
     for (size_t x = 0; x < width; x++) {
-      auto id = ids.data_(y, x, pl.color_by_);
+      int idx = pl.color_by_ == PlotColorBy::cells ? 0 : 2;
+      auto id = ids.data_(y, x, idx);
       // no setting needed if not found
       if (id == NOT_FOUND) { continue; }
       if (id == OVERLAP) {
@@ -838,7 +841,7 @@ void create_voxel(Plot const& pl)
     IdData ids = pltbase.get_map<IdData>();
 
     // select only cell/material ID data and flip the y-axis
-    int idx = pl.color_by_ == PlotColorBy::cells ? 0 : 1;
+    int idx = pl.color_by_ == PlotColorBy::cells ? 0 : 2;
     xt::xtensor<int32_t, 2> data_slice = xt::view(ids.data_, xt::all(), xt::all(), idx);
     xt::xtensor<int32_t, 2> data_flipped = xt::flip(data_slice, 0);
 
diff --git a/src/summary.cpp b/src/summary.cpp
index e43982a6a3..9076f34dbd 100644
--- a/src/summary.cpp
+++ b/src/summary.cpp
@@ -1,6 +1,10 @@
 #include "openmc/summary.h"
 
+#include <fmt/core.h>
+
+#include "openmc/capi.h"
 #include "openmc/cell.h"
+#include "openmc/file_utils.h"
 #include "openmc/hdf5_interface.h"
 #include "openmc/lattice.h"
 #include "openmc/material.h"
@@ -125,4 +129,120 @@ void write_materials(hid_t file)
   close_group(materials_group);
 }
 
+//==============================================================================
+// C API
+//==============================================================================
+
+extern "C" int openmc_properties_export(const char* filename)
+{
+  // Set a default filename if none was passed
+  std::string name = filename ? filename : "properties.h5";
+
+  // Display output message
+  auto msg = fmt::format("Exporting properties to {}...", name);
+  write_message(msg, 5);
+
+  // Create a new file using default properties.
+  hid_t file = file_open(name, 'w');
+
+  // Write metadata
+  write_attribute(file, "filetype", "properties");
+  write_attribute(file, "version", VERSION_STATEPOINT);
+  write_attribute(file, "openmc_version", VERSION);
+#ifdef GIT_SHA1
+  write_attribute(file, "git_sha1", GIT_SHA1);
+#endif
+  write_attribute(file, "date_and_time", time_stamp());
+  write_attribute(file, "path", settings::path_input);
+
+  // Write cell properties
+  auto geom_group = create_group(file, "geometry");
+  write_attribute(geom_group, "n_cells", model::cells.size());
+  auto cells_group = create_group(geom_group, "cells");
+  for (const auto& c : model::cells) {
+    c->export_properties_hdf5(cells_group);
+  }
+  close_group(cells_group);
+  close_group(geom_group);
+
+  // Write material properties
+  hid_t materials_group = create_group(file, "materials");
+  write_attribute(materials_group, "n_materials", model::materials.size());
+  for (const auto& mat : model::materials) {
+    mat->export_properties_hdf5(materials_group);
+  }
+  close_group(materials_group);
+
+  // Terminate access to the file.
+  file_close(file);
+  return 0;
+}
+
+extern "C" int openmc_properties_import(const char* filename)
+{
+  // Display output message
+  auto msg = fmt::format("Importing properties from {}...", filename);
+  write_message(msg, 5);
+
+  // Create a new file using default properties.
+  if (!file_exists(filename)) {
+    set_errmsg(fmt::format("File '{}' does not exist.", filename));
+    return OPENMC_E_INVALID_ARGUMENT;
+  }
+  hid_t file = file_open(filename, 'r');
+
+  // Ensure the filetype is correct
+  std::string filetype;
+  read_attribute(file, "filetype", filetype);
+  if (filetype != "properties") {
+    file_close(file);
+    set_errmsg(fmt::format("File '{}' is not a properties file.", filename));
+    return OPENMC_E_INVALID_ARGUMENT;
+  }
+
+  // Make sure number of cells matches
+  auto geom_group = open_group(file, "geometry");
+  int32_t n;
+  read_attribute(geom_group, "n_cells", n);
+  if (n != openmc::model::cells.size()) {
+    close_group(geom_group);
+    file_close(file);
+    set_errmsg(fmt::format("Number of cells in {} doesn't match current model.", filename));
+    return OPENMC_E_GEOMETRY;
+  }
+
+  // Read cell properties
+  auto cells_group = open_group(geom_group, "cells");
+  try {
+    for (const auto& c : model::cells) {
+      c->import_properties_hdf5(cells_group);
+    }
+  } catch (const std::exception& e) {
+    set_errmsg(e.what());
+    return OPENMC_E_UNASSIGNED;
+  }
+  close_group(cells_group);
+  close_group(geom_group);
+
+  // Make sure number of cells matches
+  auto materials_group = open_group(file, "materials");
+  read_attribute(materials_group, "n_materials", n);
+  if (n != openmc::model::materials.size()) {
+    close_group(materials_group);
+    file_close(file);
+    set_errmsg(fmt::format("Number of materials in {} doesn't match current model.", filename));
+    return OPENMC_E_GEOMETRY;
+  }
+
+  // Read material properties
+  for (const auto& mat : model::materials) {
+    mat->import_properties_hdf5(materials_group);
+  }
+  close_group(materials_group);
+
+  // Terminate access to the file.
+  file_close(file);
+  return 0;
+}
+
 } // namespace openmc
diff --git a/tests/unit_tests/test_lib.py b/tests/unit_tests/test_lib.py
index 87d4843c77..9a4fc8b9fd 100644
--- a/tests/unit_tests/test_lib.py
+++ b/tests/unit_tests/test_lib.py
@@ -1,4 +1,5 @@
 from collections.abc import Mapping
+from ctypes import ArgumentError
 import os
 
 import numpy as np
@@ -115,6 +116,8 @@ def test_cell(lib_init):
     assert cell.name == "Fuel"
     cell.name = "Not fuel"
     assert cell.name == "Not fuel"
+    assert cell.num_instances == 1
+
 
 def test_cell_temperature(lib_init):
     cell = openmc.lib.cells[1]
@@ -124,6 +127,21 @@ def test_cell_temperature(lib_init):
     assert cell.get_temperature() == pytest.approx(200.0)
 
 
+def test_properties_temperature(lib_init):
+    # Cell temperature should be 200 from above test
+    cell = openmc.lib.cells[1]
+    assert cell.get_temperature() == pytest.approx(200.0)
+
+    # Export properties and change temperature
+    openmc.lib.export_properties('properties.h5')
+    cell.set_temperature(300.0)
+    assert cell.get_temperature() == pytest.approx(300.0)
+
+    # Import properties and check that temperature is restored
+    openmc.lib.import_properties('properties.h5')
+    assert cell.get_temperature() == pytest.approx(200.0)
+
+
 def test_new_cell(lib_init):
     with pytest.raises(exc.AllocationError):
         openmc.lib.Cell(1)
@@ -132,6 +150,13 @@ def test_new_cell(lib_init):
     assert len(openmc.lib.cells) == 5
 
 
+def test_properties_fail_cell(lib_init):
+    # The number of cells was changed in the previous test, so the properties
+    # file is no longer valid
+    with pytest.raises(exc.GeometryError, match="Number of cells"):
+        openmc.lib.import_properties("properties.h5")
+
+
 def test_material_mapping(lib_init):
     mats = openmc.lib.materials
     assert isinstance(mats, Mapping)
@@ -162,11 +187,30 @@ def test_material(lib_init):
     assert sum(m.densities) == pytest.approx(rho)
 
     m.set_density(0.1, 'g/cm3')
-    assert m.density == pytest.approx(0.1)
+    assert m.get_density('g/cm3') == pytest.approx(0.1)
     assert m.name == "Hot borated water"
     m.name = "Not hot borated water"
     assert m.name == "Not hot borated water"
 
+
+def test_properties_density(lib_init):
+    m = openmc.lib.materials[1]
+    orig_density = m.get_density('atom/b-cm')
+    orig_density_gpcc = m.get_density('g/cm3')
+
+    # Export properties and change density
+    openmc.lib.export_properties('properties.h5')
+    m.set_density(orig_density_gpcc*2, 'g/cm3')
+    assert m.get_density() == pytest.approx(orig_density*2)
+
+    # Import properties and check that density was restored
+    openmc.lib.import_properties('properties.h5')
+    assert m.get_density() == pytest.approx(orig_density)
+
+    with pytest.raises(ValueError):
+        m.get_density('🥏')
+
+
 def test_material_add_nuclide(lib_init):
     m = openmc.lib.materials[3]
     m.add_nuclide('Xe135', 1e-12)
@@ -182,6 +226,13 @@ def test_new_material(lib_init):
     assert len(openmc.lib.materials) == 5
 
 
+def test_properties_fail_material(lib_init):
+    # The number of materials was changed in the previous test, so the properties
+    # file is no longer valid
+    with pytest.raises(exc.GeometryError, match="Number of materials"):
+        openmc.lib.import_properties("properties.h5")
+
+
 def test_nuclide_mapping(lib_init):
     nucs = openmc.lib.nuclides
     assert isinstance(nucs, Mapping)
@@ -558,9 +609,9 @@ def test_load_nuclide(lib_init):
 
 
 def test_id_map(lib_init):
-    expected_ids = np.array([[(3, 3), (2, 2), (3, 3)],
-                             [(2, 2), (1, 1), (2, 2)],
-                             [(3, 3), (2, 2), (3, 3)]], dtype='int32')
+    expected_ids = np.array([[(3, 0, 3), (2, 0, 2), (3, 0, 3)],
+                             [(2, 0, 2), (1, 0, 1), (2, 0, 2)],
+                             [(3, 0, 3), (2, 0, 2), (3, 0, 3)]], dtype='int32')
 
     # create a plot object
     s = openmc.lib.plot._PlotBase()
@@ -575,6 +626,7 @@ def test_id_map(lib_init):
     ids = openmc.lib.plot.id_map(s)
     assert np.array_equal(expected_ids, ids)
 
+
 def test_property_map(lib_init):
     expected_properties = np.array(
         [[(293.6, 0.740582), (293.6, 6.55), (293.6, 0.740582)],
diff --git a/tests/unit_tests/test_model.py b/tests/unit_tests/test_model.py
new file mode 100644
index 0000000000..4b658d5ab0
--- /dev/null
+++ b/tests/unit_tests/test_model.py
@@ -0,0 +1,34 @@
+import pytest
+import openmc
+import openmc.lib
+
+
+def test_import_properties(run_in_tmpdir, mpi_intracomm):
+    """Test importing properties on the Model class """
+
+    # Create PWR pin cell model and write XML files
+    openmc.reset_auto_ids()
+    model = openmc.examples.pwr_pin_cell()
+    model.export_to_xml()
+
+    # Change fuel temperature and density and export properties
+    openmc.lib.init(intracomm=mpi_intracomm)
+    cell = openmc.lib.cells[1]
+    cell.set_temperature(600.0)
+    cell.fill.set_density(5.0, 'g/cm3')
+    openmc.lib.export_properties()
+    openmc.lib.finalize()
+
+    # Import properties to existing model and re-export to new directory
+    model.import_properties("properties.h5")
+    model.export_to_xml("with_properties")
+
+    # Load model with properties and confirm temperature/density has been changed
+    model_with_properties = openmc.Model.from_xml(
+        'with_properties/geometry.xml',
+        'with_properties/materials.xml',
+        'with_properties/settings.xml'
+    )
+    cell = model_with_properties.geometry.get_all_cells()[1]
+    assert cell.temperature == [600.0]
+    assert cell.fill.get_mass_density() == pytest.approx(5.0)