diff --git a/include/cantera/oneD/Domain1D.h b/include/cantera/oneD/Domain1D.h
index 555d6cd7d59..4425141815c 100644
--- a/include/cantera/oneD/Domain1D.h
+++ b/include/cantera/oneD/Domain1D.h
@@ -182,7 +182,7 @@ class Domain1D
     }
 
     //! index of component with name \a name.
-    size_t componentIndex(const std::string& name) const;
+    virtual size_t componentIndex(const std::string& name) const;
 
     void setBounds(size_t n, doublereal lower, doublereal upper) {
         m_min[n] = lower;
diff --git a/include/cantera/oneD/Sim1D.h b/include/cantera/oneD/Sim1D.h
index 967206e3d9a..8ab374e26a8 100644
--- a/include/cantera/oneD/Sim1D.h
+++ b/include/cantera/oneD/Sim1D.h
@@ -129,15 +129,22 @@ class Sim1D : public OneDim
     int refine(int loglevel=0);
 
     //! Add node for fixed temperature point of freely propagating flame
-    int setFixedTemperature(doublereal t);
+    int setFixedTemperature(double t);
+
+    //! Return temperature at the point used to fix the flame location
+    double fixedTemperature();
+
+    //! Return location of the point where temperature is fixed
+    double fixedTemperatureLocation();
 
     /**
      * Set grid refinement criteria. If dom >= 0, then the settings
      * apply only to the specified domain.  If dom < 0, the settings
      * are applied to each domain.  @see Refiner::setCriteria.
      */
-    void setRefineCriteria(int dom = -1, doublereal ratio = 10.0,
-                           doublereal slope = 0.8, doublereal curve = 0.8, doublereal prune = -0.1);
+    void setRefineCriteria(int dom = -1, double ratio = 10.0,
+                           double slope = 0.8, double curve = 0.8,
+                           double prune = -0.1);
 
     /**
      * Get the grid refinement criteria. dom must be greater than
diff --git a/include/cantera/oneD/StFlow.h b/include/cantera/oneD/StFlow.h
index aeaf0b668f9..fef7d81a33e 100644
--- a/include/cantera/oneD/StFlow.h
+++ b/include/cantera/oneD/StFlow.h
@@ -198,13 +198,24 @@ class StFlow : public Domain1D
         return m_do_radiation;
     }
 
+    //! Return radiative heat loss at grid point j
+    double radiativeHeatLoss(size_t j) const {
+        return m_qdotRadiation[j];
+    }
+
     //! Set the emissivities for the boundary values
     /*!
      * Reads the emissivities for the left and right boundary values in the
      * radiative term and writes them into the variables, which are used for the
      * calculation.
      */
-    void setBoundaryEmissivities(doublereal e_left, doublereal e_right);
+    void setBoundaryEmissivities(double e_left, double e_right);
+
+    //! Return emissivitiy at left boundary
+    double leftEmissivity() const { return m_epsilon_left; }
+
+    //! Return emissivitiy at right boundary
+    double rightEmissivity() const { return m_epsilon_right; }
 
     void fixTemperature(size_t j=npos);
 
@@ -215,7 +226,14 @@ class StFlow : public Domain1D
     //! Change the grid size. Called after grid refinement.
     virtual void resize(size_t components, size_t points);
 
-    virtual void setFixedPoint(int j0, doublereal t0) {}
+    /*!
+     * @deprecated To be removed after Cantera 2.5.
+     */
+    virtual void setFixedPoint(int j0, doublereal t0) {
+        // this does nothing and does not appear to be overloaded
+        warn_deprecated("StFlow::setFixedPoint",
+                        "To be removed after Cantera 2.5.");
+    }
 
     //! Set the gas object state to be consistent with the solution at point j.
     void setGas(const doublereal* x, size_t j);
diff --git a/interfaces/cython/cantera/_cantera.pxd b/interfaces/cython/cantera/_cantera.pxd
index 20746d9d2e1..dd12c9de196 100644
--- a/interfaces/cython/cantera/_cantera.pxd
+++ b/interfaces/cython/cantera/_cantera.pxd
@@ -750,9 +750,12 @@ cdef extern from "cantera/oneD/StFlow.h":
         void setPressure(double)
         void enableRadiation(cbool)
         cbool radiationEnabled()
+        double radiativeHeatLoss(size_t)
         double pressure()
         void setFixedTempProfile(vector[double]&, vector[double]&)
         void setBoundaryEmissivities(double, double)
+        double leftEmissivity()
+        double rightEmissivity()
         void solveEnergyEqn()
         void fixTemperature()
         cbool doEnergy(size_t)
@@ -760,6 +763,7 @@ cdef extern from "cantera/oneD/StFlow.h":
         cbool withSoret()
         void setFreeFlow()
         void setAxisymmetricFlow()
+        string flowType()
 
 
 cdef extern from "cantera/oneD/IonFlow.h":
@@ -815,6 +819,8 @@ cdef extern from "cantera/oneD/Sim1D.h":
         size_t maxGridPoints(size_t) except +translate_exception
         void setGridMin(int, double) except +translate_exception
         void setFixedTemperature(double) except +translate_exception
+        double fixedTemperature()
+        double fixedTemperatureLocation()
         void setInterrupt(CxxFunc1*) except +translate_exception
         void setTimeStepCallback(CxxFunc1*)
         void setSteadyCallback(CxxFunc1*)
diff --git a/interfaces/cython/cantera/composite.py b/interfaces/cython/cantera/composite.py
index 1248bf090c6..2c126e08a08 100644
--- a/interfaces/cython/cantera/composite.py
+++ b/interfaces/cython/cantera/composite.py
@@ -731,12 +731,16 @@ def restore_data(self, data, labels):
                 if s in valid_species:
                     state_data[-1][:, i] = data[:, valid_species[s]]
 
-        # labels may include calculated properties that must not be restored
-        calculated = self._scalar + self._n_species + self._n_reactions
-        exclude = [l for l in labels
-                   if any([v in l for v in calculated])]
-        extra = {l: list(data[:, i]) for i, l in enumerate(labels)
-                 if l not in exclude}
+        # labels may include calculated properties that must not be restored:
+        # compare column labels to names that are reserved for SolutionArray
+        # attributes (see `SolutionArray.collect_data`), i.e. scalar values,
+        # arrays with number of species, and arrays with number of reactions.
+        exclude = [lab for lab in labels
+                   if any([lab in self._scalar,
+                           '_'.join(lab.split('_')[:-1]) in self._n_species,
+                           lab.split(' ')[0] in self._n_reactions])]
+        extra = {lab: list(data[:, i]) for i, lab in enumerate(labels)
+                 if lab not in exclude}
         if len(self._extra):
             extra_lists = {k: extra[k] for k in self._extra}
         else:
@@ -861,7 +865,7 @@ def write_csv(self, filename, cols=None, *args, **kwargs):
         only with 1D `SolutionArray` objects.
         """
         data, labels = self.collect_data(cols=cols, *args, **kwargs)
-        with open(filename, 'w') as outfile:
+        with open(filename, 'w', newline='') as outfile:
             writer = _csv.writer(outfile)
             writer.writerow(labels)
             for row in data:
diff --git a/interfaces/cython/cantera/examples/onedim/adiabatic_flame.py b/interfaces/cython/cantera/examples/onedim/adiabatic_flame.py
index 5c3b3f9f4d4..75890d8135d 100644
--- a/interfaces/cython/cantera/examples/onedim/adiabatic_flame.py
+++ b/interfaces/cython/cantera/examples/onedim/adiabatic_flame.py
@@ -6,6 +6,11 @@
 """
 
 import cantera as ct
+try:
+    import pandas as pd
+    import tables
+except:
+    pd = None
 
 # Simulation parameters
 p = ct.one_atm  # pressure [Pa]
@@ -29,16 +34,25 @@
 f.solve(loglevel=loglevel, auto=True)
 
 # Solve with the energy equation enabled
-f.save('h2_adiabatic.xml', 'mix', 'solution with mixture-averaged transport')
+if pd:
+    f.write_hdf('h2_adiabatic.h5', key='mix', mode='w')
+else:
+    f.save('h2_adiabatic.xml', 'mix',
+           'solution with mixture-averaged transport')
+
 f.show_solution()
-print('mixture-averaged flamespeed = {0:7f} m/s'.format(f.u[0]))
+print('mixture-averaged flamespeed = {0:7f} m/s'.format(f.velocity[0]))
 
 # Solve with multi-component transport properties
 f.transport_model = 'Multi'
 f.solve(loglevel)  # don't use 'auto' on subsequent solves
 f.show_solution()
-print('multicomponent flamespeed = {0:7f} m/s'.format(f.u[0]))
-f.save('h2_adiabatic.xml', 'multi', 'solution with multicomponent transport')
+print('multicomponent flamespeed = {0:7f} m/s'.format(f.velocity[0]))
+if pd:
+    f.write_hdf('h2_adiabatic.h5', key='multi')
+else:
+    f.save('h2_adiabatic.xml', 'multi',
+           'solution with multicomponent transport')
 
 # write the velocity, temperature, density, and mole fractions to a CSV file
 f.write_csv('h2_adiabatic.csv', quiet=False)
diff --git a/interfaces/cython/cantera/examples/onedim/diffusion_flame_batch.py b/interfaces/cython/cantera/examples/onedim/diffusion_flame_batch.py
index 311b7dfb91e..488b0e45a34 100644
--- a/interfaces/cython/cantera/examples/onedim/diffusion_flame_batch.py
+++ b/interfaces/cython/cantera/examples/onedim/diffusion_flame_batch.py
@@ -111,10 +111,10 @@ def interrupt_extinction(t):
     f.fuel_inlet.mdot *= rel_pressure_increase ** exp_mdot_p
     f.oxidizer_inlet.mdot *= rel_pressure_increase ** exp_mdot_p
     # Update velocities
-    f.set_profile('u', normalized_grid,
-                  f.u * rel_pressure_increase ** exp_u_p)
-    f.set_profile('V', normalized_grid,
-                  f.V * rel_pressure_increase ** exp_V_p)
+    f.set_profile('velocity', normalized_grid,
+                  f.velocity * rel_pressure_increase ** exp_u_p)
+    f.set_profile('spread_rate', normalized_grid,
+                  f.spread_rate * rel_pressure_increase ** exp_V_p)
     # Update pressure curvature
     f.set_profile('lambda', normalized_grid,
                   f.L * rel_pressure_increase ** exp_lam_p)
@@ -167,8 +167,10 @@ def interrupt_extinction(t):
     f.fuel_inlet.mdot *= strain_factor ** exp_mdot_a
     f.oxidizer_inlet.mdot *= strain_factor ** exp_mdot_a
     # Update velocities
-    f.set_profile('u', normalized_grid, f.u * strain_factor ** exp_u_a)
-    f.set_profile('V', normalized_grid, f.V * strain_factor ** exp_V_a)
+    f.set_profile('velocity', normalized_grid,
+                  f.velocity * strain_factor ** exp_u_a)
+    f.set_profile('spread_rate', normalized_grid,
+                  f.spread_rate * strain_factor ** exp_V_a)
     # Update pressure curvature
     f.set_profile('lambda', normalized_grid, f.L * strain_factor ** exp_lam_a)
     try:
@@ -203,7 +205,7 @@ def interrupt_extinction(t):
 
     # Plot the axial velocity profiles (normalized by the fuel inlet velocity)
     # for selected pressures
-    ax2.plot(f.grid / f.grid[-1], f.u / f.u[0],
+    ax2.plot(f.grid / f.grid[-1], f.velocity / f.velocity[0],
              label='{0:05.1f} bar'.format(p))
 
 ax1.legend(loc=0)
@@ -231,7 +233,7 @@ def interrupt_extinction(t):
 
     # Plot the axial velocity profiles (normalized by the fuel inlet velocity)
     # for the strain rate loop (selected)
-    ax4.plot(f.grid / f.grid[-1], f.u / f.u[0],
+    ax4.plot(f.grid / f.grid[-1], f.velocity / f.velocity[0],
              label=format(a_max, '.2e') + ' 1/s')
 
 ax3.legend(loc=0)
diff --git a/interfaces/cython/cantera/examples/onedim/diffusion_flame_extinction.py b/interfaces/cython/cantera/examples/onedim/diffusion_flame_extinction.py
index 1183f2a8c70..df9d1e5075d 100644
--- a/interfaces/cython/cantera/examples/onedim/diffusion_flame_extinction.py
+++ b/interfaces/cython/cantera/examples/onedim/diffusion_flame_extinction.py
@@ -88,7 +88,7 @@
 # List of peak temperatures
 T_max = [np.max(f.T)]
 # List of maximum axial velocity gradients
-a_max = [np.max(np.abs(np.gradient(f.u) / np.gradient(f.grid)))]
+a_max = [np.max(np.abs(np.gradient(f.velocity) / np.gradient(f.grid)))]
 
 # Simulate counterflow flames at increasing strain rates until the flame is
 # extinguished. To achieve a fast simulation, an initial coarse strain rate
@@ -109,8 +109,10 @@
     f.fuel_inlet.mdot *= strain_factor ** exp_mdot_a
     f.oxidizer_inlet.mdot *= strain_factor ** exp_mdot_a
     # Update velocities
-    f.set_profile('u', normalized_grid, f.u * strain_factor ** exp_u_a)
-    f.set_profile('V', normalized_grid, f.V * strain_factor ** exp_V_a)
+    f.set_profile('velocity', normalized_grid,
+                  f.velocity * strain_factor ** exp_u_a)
+    f.set_profile('spread_rate', normalized_grid,
+                  f.spread_rate * strain_factor ** exp_V_a)
     # Update pressure curvature
     f.set_profile('lambda', normalized_grid, f.L * strain_factor ** exp_lam_a)
     try:
@@ -125,7 +127,7 @@
                description='Cantera version ' + ct.__version__ +
                ', reaction mechanism ' + reaction_mechanism)
         T_max.append(np.max(f.T))
-        a_max.append(np.max(np.abs(np.gradient(f.u) / np.gradient(f.grid))))
+        a_max.append(np.max(np.abs(np.gradient(f.velocity) / np.gradient(f.grid))))
         # If the temperature difference is too small and the minimum relative
         # strain rate increase is reached, abort
         if ((T_max[-2] - T_max[-1] < delta_T_min) &
diff --git a/interfaces/cython/cantera/examples/onedim/flamespeed_sensitivity.py b/interfaces/cython/cantera/examples/onedim/flamespeed_sensitivity.py
index a4ee73601e9..b477f3e55b4 100644
--- a/interfaces/cython/cantera/examples/onedim/flamespeed_sensitivity.py
+++ b/interfaces/cython/cantera/examples/onedim/flamespeed_sensitivity.py
@@ -24,7 +24,7 @@
 f.set_refine_criteria(ratio=3, slope=0.07, curve=0.14)
 
 f.solve(loglevel=1, auto=True)
-print('\nmixture-averaged flamespeed = {:7f} m/s\n'.format(f.u[0]))
+print('\nmixture-averaged flamespeed = {:7f} m/s\n'.format(f.velocity[0]))
 
 # Use the adjoint method to calculate sensitivities
 sens = f.get_flame_speed_reaction_sensitivities()
diff --git a/interfaces/cython/cantera/examples/onedim/ion_free_flame.py b/interfaces/cython/cantera/examples/onedim/ion_free_flame.py
index 0ab8fc12df9..6dcf715736d 100644
--- a/interfaces/cython/cantera/examples/onedim/ion_free_flame.py
+++ b/interfaces/cython/cantera/examples/onedim/ion_free_flame.py
@@ -31,7 +31,7 @@
 
 f.save('CH4_adiabatic.xml', 'ion', 'solution with ionized gas transport')
 f.show_solution()
-print('mixture-averaged flamespeed = {0:7f} m/s'.format(f.u[0]))
+print('mixture-averaged flamespeed = {0:7f} m/s'.format(f.velocity[0]))
 
 # write the velocity, temperature, density, and mole fractions to a CSV file
 f.write_csv('CH4_adiabatic.csv', quiet=False)
diff --git a/interfaces/cython/cantera/examples/onedim/premixed_counterflow_twin_flame.py b/interfaces/cython/cantera/examples/onedim/premixed_counterflow_twin_flame.py
index 3bf8d05bf50..72039209ce9 100644
--- a/interfaces/cython/cantera/examples/onedim/premixed_counterflow_twin_flame.py
+++ b/interfaces/cython/cantera/examples/onedim/premixed_counterflow_twin_flame.py
@@ -29,12 +29,12 @@ def derivative(x, y):
 
 def computeStrainRates(oppFlame):
     # Compute the derivative of axial velocity to obtain normal strain rate
-    strainRates = derivative(oppFlame.grid, oppFlame.u)
+    strainRates = derivative(oppFlame.grid, oppFlame.velocity)
 
     # Obtain the location of the max. strain rate upstream of the pre-heat zone.
     # This is the characteristic strain rate
     maxStrLocation = abs(strainRates).argmax()
-    minVelocityPoint = oppFlame.u[:maxStrLocation].argmin()
+    minVelocityPoint = oppFlame.velocity[:maxStrLocation].argmin()
 
     # Characteristic Strain Rate = K
     strainRatePoint = abs(strainRates[:minVelocityPoint]).argmax()
@@ -130,15 +130,17 @@ def solveOpposedFlame(oppFlame, massFlux=0.12, loglevel=1,
 
     # Axial Velocity Plot
     plt.subplot(1, 2, 1)
-    plt.plot(oppFlame.grid, oppFlame.u, 'r', lw=2)
+    plt.plot(oppFlame.grid, oppFlame.velocity, 'r', lw=2)
     plt.xlim(oppFlame.grid[0], oppFlame.grid[-1])
     plt.xlabel('Distance (m)')
     plt.ylabel('Axial Velocity (m/s)')
 
     # Identify the point where the strain rate is calculated
-    plt.plot(oppFlame.grid[strainRatePoint], oppFlame.u[strainRatePoint], 'gs')
+    plt.plot(oppFlame.grid[strainRatePoint],
+             oppFlame.velocity[strainRatePoint], 'gs')
     plt.annotate('Strain-Rate point',
-                 xy=(oppFlame.grid[strainRatePoint], oppFlame.u[strainRatePoint]),
+                 xy=(oppFlame.grid[strainRatePoint],
+                     oppFlame.velocity[strainRatePoint]),
                  xytext=(0.001, 0.1),
                  arrowprops={'arrowstyle': '->'})
 
diff --git a/interfaces/cython/cantera/onedim.py b/interfaces/cython/cantera/onedim.py
index 46d75d47c7f..f28e966b87d 100644
--- a/interfaces/cython/cantera/onedim.py
+++ b/interfaces/cython/cantera/onedim.py
@@ -3,14 +3,22 @@
 
 import numpy as np
 from ._cantera import *
-from .composite import Solution
-import csv as _csv
+from .composite import Solution, SolutionArray
 from math import erf
+from os import path
+from email.utils import formatdate
+
+# avoid explicit dependence of cantera on pandas
+try:
+    import pandas as _pandas
+except ImportError as err:
+    _pandas = err
 
 
 class FlameBase(Sim1D):
     """ Base class for flames with a single flow domain """
     __slots__ = ('gas',)
+    _extra = () # extra columns used for saving/restoring of simulation data
 
     def __init__(self, domains, gas, grid=None):
         """
@@ -63,6 +71,87 @@ def get_refine_criteria(self):
         """
         return super().get_refine_criteria(self.flame)
 
+    def set_initial_guess(self, *args, data=None, key=None, **kwargs):
+        """
+        Set the initial guess for the solution, and load restart data if
+        provided. Derived classes extend this function to set approximations
+        for the temperature and composition profiles.
+
+        :param data:
+            Restart data, which are typically based on an earlier simulation
+            result. Restart data may be specified using a `SolutionArray`,
+            pandas' DataFrame, or previously saved CSV or HDF container files. 
+            Note that restart data do not overwrite boundary conditions. 
+            DataFrame and HDF input require working installations of pandas and 
+            PyTables. These packages can be installed using pip (`pandas` and 
+            `tables`) or conda (`pandas` and `pytables`).
+        :param key:
+            Group identifier within a HDF container file (only used in
+            combination with HDF restart data).
+        """
+        super().set_initial_guess(*args, data=data, key=key, **kwargs)
+        if not data:
+            return
+
+        # load restart data into SolutionArray
+        if isinstance(data, SolutionArray):
+            # already a solution array
+            arr = data
+        elif isinstance(data, str):
+            if data.endswith('.hdf5') or data.endswith('.h5'):
+                # data source identifies a HDF file
+                arr = SolutionArray(self.gas, extra=self._extra)
+                arr.read_hdf(data, key=key)
+            elif data.endswith('.csv'):
+                # data source identifies a CSV file
+                arr = SolutionArray(self.gas, extra=self._extra)
+                arr.read_csv(data)
+            else:
+                raise ValueError(
+                    "'{}' does not identify CSV or HDF file.".format(data)
+                )
+        else:
+            # data source is a pandas DataFrame
+            arr = SolutionArray(self.gas, extra=self._extra)
+            arr.from_pandas(data)
+
+        # get left and right boundaries
+        left = self.domains[0]
+        right = self.domains[2]
+
+        if isinstance(left, Inlet1D) and isinstance(right, Inlet1D):
+            # find stagnation plane
+            i = np.flatnonzero(self.velocity > 0)[-1]
+
+            # adjust temperatures
+            T = arr.T
+            xi = arr.grid[1:-1]
+            T[1:-1] += (left.T - T[0]) * (1 - xi) + (right.T - T[-1]) * xi
+            arr.TP = T, self.P
+
+            # adjust velocities
+            u = arr.velocity
+
+            self.gas.TPY = left.T, self.P, left.Y
+            u[:i] = u[:i] * left.mdot / self.gas.density / u[0]
+
+            self.gas.TPY = right.T, self.P, right.Y
+            u[i:] = - u[i:] * right.mdot / self.gas.density / u[-1]
+
+            arr.velocity = u
+
+        elif isinstance(left, Inlet1D):
+            # adjust temperatures
+            arr.TP = arr.T + left.T - arr.T[0], self.P
+
+            # adjust velocities
+            if self.flame.flow_type != "Free Flame":
+                self.gas.TPY = left.T, self.P, left.Y
+                u0 = left.mdot/self.gas.density
+                arr.velocity *= u0 / arr.velocity[0]
+
+        self.from_solution_array(arr, restore_boundaries=False)
+
     def set_profile(self, component, locations, values):
         """
         Set an initial estimate for a profile of one component.
@@ -162,15 +251,46 @@ def T(self):
     def u(self):
         """
         Array containing the velocity [m/s] normal to the flame at each point.
+
+        .. deprecated:: 2.5
+
+             To be deprecated with version 2.5, and removed thereafter.
+             Replaced by property `velocity`.
+        """
+        warnings.warn("Property 'u' to be removed after Cantera 2.5. "
+                      "Replaced by property 'velocity'",
+                      DeprecationWarning)
+        return self.profile(self.flame, 'velocity')
+
+    @property
+    def velocity(self):
+        """
+        Array containing the velocity [m/s] normal to the flame at each point.
         """
-        return self.profile(self.flame, 'u')
+        return self.profile(self.flame, 'velocity')
 
     @property
     def V(self):
         """
         Array containing the tangential velocity gradient [1/s] at each point.
+
+        .. deprecated:: 2.5
+
+             To be deprecated with version 2.5, and removed thereafter.
+             Replaced by property `spread_rate`.
+        """
+        warnings.warn("Property 'V' to be removed after Cantera 2.5. "
+                      "Replaced by property 'spread_rate'",
+                      DeprecationWarning)
+        return self.profile(self.flame, 'spread_rate')
+
+    @property
+    def spread_rate(self):
+        """
+        Array containing the tangential velocity gradient [1/s] (e.g. radial
+        velocity divided by radius) at each point.
         """
-        return self.profile(self.flame, 'V')
+        return self.profile(self.flame, 'spread_rate')
 
     @property
     def L(self):
@@ -262,23 +382,354 @@ def write_csv(self, filename, species='X', quiet=True):
             mole fractions or ``Y`` for mass fractions.
         """
 
-        z = self.grid
-        T = self.T
-        u = self.u
-        V = self.V
-
-        with open(filename, 'w', newline='') as csvfile:
-            writer = _csv.writer(csvfile)
-            writer.writerow(['z (m)', 'u (m/s)', 'V (1/s)',
-                            'T (K)', 'rho (kg/m3)'] + self.gas.species_names)
-            for n in range(self.flame.n_points):
-                self.set_gas_state(n)
-                writer.writerow([z[n], u[n], V[n], T[n], self.gas.density] +
-                                list(getattr(self.gas, species)))
+        # save data
+        cols = ('extra', 'T', 'D', species)
+        self.to_solution_array().write_csv(filename, cols=cols)
+
+        if not quiet:
+            print("Solution saved to '{0}'.".format(filename))
+
+    def to_solution_array(self):
+        """
+        Return the solution vector as a Cantera `SolutionArray` object.
+
+        The `SolutionArray` has the following ``extra`` entries:
+         * ``grid``: grid point positions along the flame [m]
+         * ``velocity``: normal velocity [m/s]
+         * ``spread_rate``: tangential velocity gradient [1/s] (if applicable)
+         * ``lambda``: radial pressure gradient [N/m^4] (if applicable)
+         * ``eField``: electric field strength (if applicable)
+        """
+        # create extra columns
+        extra = {}
+        for e in self._extra:
+            if e == 'grid':
+                val = self.grid
+            else:
+                val = self.profile(self.flame, e)
+            extra[e] = np.hstack([np.nan, val, np.nan])
+        if self.radiation_enabled:
+            qdot = self.flame.radiative_heat_loss
+            extra['qdot'] = np.hstack([np.nan, qdot, np.nan])
+
+        # consider inlet boundaries
+        left = self.domains[0]  # left boundary is always an inlet
+        if isinstance(self.domains[2], Inlet1D):
+            right = self.domains[2]
+            n_arr = self.flame.n_points + 2
+        else:
+            right = None
+            n_arr = self.flame.n_points + 1
+            for e in extra:
+                extra[e] = extra[e][:-1]
+
+        # create solution array object
+        arr = SolutionArray(self.gas, n_arr, extra=extra)
+
+        # add left boundary
+        self.gas.TPY = left.T, self.P, left.Y
+        arr[0].TPY = self.gas.TPY
+        arr._extra['grid'][0] = -np.inf
+        arr._extra['velocity'][0] = left.mdot/self.gas.density
+
+        # retrieve species concentrations and set states
+        for n in range(self.flame.n_points):
+            self.set_gas_state(n)
+            arr[n + 1].TPY = self.gas.T, self.P, self.gas.Y
+
+        # add right boundary
+        if right:
+            self.gas.TPY = right.T, self.P, right.Y
+            arr[-1].TPY = self.gas.TPY
+            arr._extra['grid'][-1] = np.inf
+            arr._extra['velocity'][-1] = -right.mdot/self.gas.density
+
+        return arr
+
+    def from_solution_array(self, arr, restore_boundaries=True,
+                            settings=None):
+        """
+        Restore the solution vector from a Cantera `SolutionArray` object.
+
+        :param arr:
+            SolutionArray to be restored
+        :param restore_boundaries:
+            Boolean flag to indicate whether boundaries should be restored
+            (default is ``True``)
+        :param settings:
+            dictionary containing simulation settings
+            (see `FlameBase.settings`)
+
+        The `SolutionArray` requires the following ``extra`` entries:
+         * ``grid``: grid point positions along the flame [m]
+         * ``velocity``: normal velocity [m/s]
+         * ``spread_rate``: tangential velocity gradient [1/s] (if applicable)
+         * ``lambda``: radial pressure gradient [N/m^4] (if applicable)
+         * ``eField``: electric field strength (if applicable)
+        """
+        # extent (indices) of flame domain
+        idx = np.isfinite(arr.grid)
+        grid = arr.grid[idx]
+
+        # restore grid
+        self.flame.grid = grid
+        self._get_initial_solution()
+        xi = (grid - grid[0]) / (grid[-1] - grid[0])
+
+        # restore temperature and 'extra' profiles
+        self.set_profile('T', xi, arr.T[idx])
+        for e in self._extra:
+            val = getattr(arr, e)[idx]
+            if e in ['grid', 'qdot']:
+                pass
+            else:
+                self.set_profile(e, xi, val)
+
+        # restore species profiles
+        X = arr.X[idx, :]
+        for i, spc in enumerate(self.gas.species_names):
+            self.set_profile(spc, xi, X[:, i])
+
+        # restore pressure
+        self.P = arr.P[0]
+
+        # restore boundaries
+        if restore_boundaries:
+
+            # left boundary
+            left = self.domains[0]
+            left.T = arr[0].T
+            left.Y = arr[0].Y
+            left.mdot = arr.velocity[0] * arr[0].density
+
+            # right boundary
+            if np.isinf(arr.grid[-1]):
+                right = self.domains[2]
+                right.T = arr[-1].T
+                right.Y = arr[-1].Y
+                right.mdot = -arr.velocity[-1] * arr[-1].density
+
+        if settings:
+            self.settings = settings
+
+    def to_pandas(self, species='X'):
+        """
+        Return the solution vector as a `pandas.DataFrame`.
+
+        :param species:
+            Attribute to use obtaining species profiles, e.g. ``X`` for
+            mole fractions or ``Y`` for mass fractions.
+
+        This method uses `to_solution_array` and requires a working pandas
+        installation. Use pip or conda to install `pandas` to enable this
+        method.
+        """
+        cols = ('extra', 'T', 'D', species)
+        return self.to_solution_array().to_pandas(cols=cols)
+
+    def from_pandas(self, df, restore_boundaries=True, settings=None):
+        """
+        Restore the solution vector from a `pandas.DataFrame`.
+
+        :param df:
+            `pandas.DataFrame` containing data to be restored
+        :param restore_boundaries:
+            Boolean flag to indicate whether boundaries should be restored
+            (default is ``True``)
+        :param settings:
+            dictionary containing simulation settings
+            (see `FlameBase.settings`)
+
+        This method is intendend for loading of data that were previously
+        exported by `to_pandas`. The method uses `from_solution_array` and
+        requires a working pandas installation. The package 'pandas' can be
+        installed using pip or conda.
+        """
+        arr = SolutionArray(self.gas, extra=self._extra)
+        arr.from_pandas(df)
+        self.from_solution_array(arr, restore_boundaries=restore_boundaries,
+                                 settings=settings)
+
+    def write_hdf(self, filename, key=None, species='X',
+                  mode=None, complevel=None, quiet=True):
+        """
+        Write the solution vector to a HDF container file. Note that it is
+        possible to write multiple data entries to a single HDF container file.
+        Simulation settings are stored in tabular form as a separate HDF group
+        named ``settings``.
+
+        :param filename:
+            HDF container file containing data to be restored
+        :param key:
+            String identifying the HDF group containing the data. The default
+            is the name of the simulated configuration, e.g. ``FreeFlame``.
+        :param species:
+            Attribute to use obtaining species profiles, e.g. ``X`` for
+            mole fractions or ``Y`` for mass fractions.
+        :param mode:
+            Mode to open file (see `pandas.DataFrame.to_hdf`)
+        :param complevel:
+            Compression level (see `pandas.DataFrame.to_hdf`)
+
+        The method exports data using `SolutionArray.write_hdf` via
+        `to_solution_array` and requires working installations of pandas and
+        PyTables. These packages can be installed using pip (`pandas` and
+        `tables`) or conda (`pandas` and `pytables`).
+        """
+        # ensure key identifying HDF group within container is unique
+        if not key:
+            key = 'data'
+        if path.exists(filename) and mode != 'w':
+            with _pandas.HDFStore(filename) as hdf:
+                 keys = hdf.keys()
+            count = sum([k.startswith(key, 1) for k in keys])
+            if count:
+                key = '{}_{}'.format(key, count)
+            row = len(keys) - 1
+        else:
+            row = 0
+
+        # save data
+        cols = ('extra', 'T', 'D', species)
+        self.to_solution_array().write_hdf(filename, cols=cols, key=key,
+                                           mode=mode, complevel=complevel)
+
+        # convert simulation settings to tabular format
+        df = _pandas.DataFrame()
+        df['key'] = [key]
+        df['date'] = formatdate(localtime=True)
+        for key, val in self.settings.items():
+            df[key] = [val]
+        df.index = _pandas.RangeIndex(start=row, stop=row+1, step=1)
+
+        # store settings to HDF container file as a separate group
+        df.to_hdf(filename, key='settings', format='table',
+                  append=True)
 
         if not quiet:
             print("Solution saved to '{0}'.".format(filename))
 
+    def read_hdf(self, filename, key=None,
+                 restore_boundaries=True, restore_settings=True):
+        """
+        Restore the solution vector from a HDF container file.
+
+        :param filename:
+            HDF container file containing data to be restored
+        :param key:
+            String identifying the HDF group containing the data
+        :param restore_boundaries:
+            Boolean flag to indicate whether boundaries should be restored
+            (default is ``True``)
+        :param restore_settings:
+            Boolean flag to indicate whether simulation settings should be
+            restored (default is ``True``)
+
+        The method imports data using `SolutionArray.read_hdf` via
+        `from_solution_array` and requires working installations of pandas and
+        PyTables. These packages can be installed using pip (`pandas` and
+        `tables`) or conda (`pandas` and `pytables`).
+        """
+        if key is None:
+            with _pandas.HDFStore(filename) as hdf:
+                 keys = hdf.keys()
+            key = [k.lstrip('/') for k in keys
+                   if k != '/settings'][0]
+
+        # retrieve data
+        arr = SolutionArray(self.gas, extra=self._extra)
+        arr.read_hdf(filename, key=key)
+
+        # retrieve settings
+        if restore_settings:
+            df = _pandas.read_hdf(filename, key='settings')
+            df = df.set_index('key')
+            series = df.loc[key]
+            settings = {k: v for k, v in series.items()
+                        if k not in ['key', 'date']}
+        else:
+            settings = None
+
+        self.from_solution_array(arr, restore_boundaries=restore_boundaries,
+                                 settings=settings)
+
+    @property
+    def settings(self):
+        """ Return a dictionary listing simulation settings """
+        out = {'configuration': type(self).__name__}
+        out['transport_model'] = self.transport_model
+        out['energy_enabled'] = self.energy_enabled
+        out['soret_enabled'] = self.soret_enabled
+        out['radiation_enabled'] = self.radiation_enabled
+        epsilon = self.flame.boundary_emissivities
+        out['emissivity_left'] = epsilon[0]
+        out['emissivity_right'] = epsilon[1]
+        out['fixed_temperature'] = self.fixed_temperature
+        out.update(self.get_refine_criteria())
+        out['max_time_step_count'] = self.max_time_step_count
+        out['max_grid_points'] = self.get_max_grid_points(self.flame)
+
+        # add tolerance settings
+        tols = {'steady_reltol': self.flame.steady_reltol(),
+                'steady_abstol': self.flame.steady_abstol(),
+                'transient_reltol': self.flame.transient_reltol(),
+                'transient_abstol': self.flame.transient_abstol()}
+        comp = np.array(self.flame.component_names)
+        for tname, tol in tols.items():
+            # add mode (most frequent tolerance setting)
+            values, counts = np.unique(tol, return_counts=True)
+            ix = np.argmax(counts)
+            out.update({tname: values[ix]})
+
+            # add values deviating from mode
+            ix = np.logical_not(np.isclose(tol, values[ix]))
+            out.update({'{}_{}'.format(tname, c): t
+                        for c, t in zip(comp[ix], tol[ix])})
+
+        return out
+
+    @settings.setter
+    def settings(self, s):
+        # simple setters
+        attr = {'transport_model',
+                'energy_enabled', 'soret_enabled', 'radiation_enabled',
+                'fixed_temperature',
+                'max_time_step_count', 'max_grid_points'}
+        attr = attr & set(s.keys())
+        for key in attr:
+            setattr(self, key, s[key])
+
+        # boundary emissivities
+        if 'emissivity_left' in s or 'emissivity_right' in s:
+            epsilon = self.flame.boundary_emissivities
+            epsilon = (s.get('emissivity_left', epsilon[0]),
+                       s.get('emissivity_right', epsilon[1]))
+            self.flame.boundary_emissivities = epsilon
+
+        # refine criteria
+        refine = {k: v for k, v in s.items()
+                  if k in ['ratio', 'slope', 'curve', 'prune']}
+        if refine:
+            self.set_refine_criteria(**refine)
+
+        # tolerances
+        tols = ['steady_reltol', 'steady_abstol',
+                'transient_reltol', 'transient_abstol']
+        tols = [t for t in tols if t in s]
+        comp = np.array(self.flame.component_names)
+        for tname in tols:
+            mode = tname.split('_')
+            tol = s[tname] * np.ones(len(comp))
+            for i, c in enumerate(comp):
+                key = '{}_{}'.format(tname, c)
+                if key in s:
+                    tol[i] = s[key]
+            tol = {mode[1][:3]: tol}
+            if mode[0] == 'steady':
+                self.flame.set_steady_tolerances(**tol)
+            else:
+                self.flame.set_transient_tolerances(**tol)
+
 
 def _trim(docstring):
     """Remove block indentation from a docstring."""
@@ -364,6 +815,7 @@ def getter(self):
 class FreeFlame(FlameBase):
     """A freely-propagating flat flame."""
     __slots__ = ('inlet', 'outlet', 'flame')
+    _extra = ('grid', 'velocity')
 
     def __init__(self, gas, grid=None, width=None):
         """
@@ -395,18 +847,28 @@ def __init__(self, gas, grid=None, width=None):
         self.inlet.T = gas.T
         self.inlet.X = gas.X
 
-    def set_initial_guess(self, locs=[0.0, 0.3, 0.5, 1.0]):
+    def set_initial_guess(self, locs=[0.0, 0.3, 0.5, 1.0], data=None, key=None):
         """
-        Set the initial guess for the solution. The adiabatic flame
+        Set the initial guess for the solution. By default, the adiabatic flame
         temperature and equilibrium composition are computed for the inlet gas
-        composition.
+        composition. Alternatively, a previously calculated result can be
+        supplied as an initial guess via 'data' and 'key' inputs (see
+        `FlameBase.set_initial_guess`).
 
         :param locs:
-            A list of four locations to define the temperature and mass fraction profiles.
-            Profiles rise linearly between the second and third location.
-            Locations are given as a fraction of the entire domain
+            A list of four locations to define the temperature and mass fraction
+            profiles. Profiles rise linearly between the second and third
+            location. Locations are given as a fraction of the entire domain
         """
-        super().set_initial_guess()
+        super().set_initial_guess(data=data, key=key)
+        if data:
+            # set fixed temperature
+            Tmid = .75 * self.T[0] + .25 * self.T[-1]
+            i = np.flatnonzero(data.T < Tmid)[-1]
+            self.fixed_temperature = data.T[i]
+
+            return
+
         self.gas.TPY = self.inlet.T, self.P, self.inlet.Y
 
         if not self.inlet.mdot:
@@ -423,7 +885,7 @@ def set_initial_guess(self, locs=[0.0, 0.3, 0.5, 1.0]):
         Yeq = self.gas.Y
         u1 = self.inlet.mdot/self.gas.density
 
-        self.set_profile('u', locs, [u0, u0, u1, u1])
+        self.set_profile('velocity', locs, [u0, u0, u1, u1])
         self.set_profile('T', locs, [T0, T0, Teq, Teq])
 
         # Pick the location of the fixed temperature point, using an existing
@@ -432,11 +894,11 @@ def set_initial_guess(self, locs=[0.0, 0.3, 0.5, 1.0]):
         Tmid = 0.75 * T0 + 0.25 * Teq
         i = np.flatnonzero(T < Tmid)[-1] # last point less than Tmid
         if Tmid - T[i] < 0.2 * (Tmid - T0):
-            self.set_fixed_temperature(T[i])
+            self.fixed_temperature = T[i]
         elif T[i+1] - Tmid < 0.2 * (Teq - Tmid):
-            self.set_fixed_temperature(T[i+1])
+            self.fixed_temperature = T[i+1]
         else:
-            self.set_fixed_temperature(Tmid)
+            self.fixed_temperature = Tmid
 
         for n in range(self.gas.n_species):
             self.set_profile(self.gas.species_name(n),
@@ -517,12 +979,12 @@ def get_flame_speed_reaction_sensitivities(self):
         """
 
         def g(sim):
-            return sim.u[0]
+            return sim.velocity[0]
 
         Nvars = sum(D.n_components * D.n_points for D in self.domains)
 
         # Index of u[0] in the global solution vector
-        i_Su = self.inlet.n_components + self.flame.component_index('u')
+        i_Su = self.inlet.n_components + self.flame.component_index('velocity')
 
         dgdx = np.zeros(Nvars)
         dgdx[i_Su] = 1
@@ -537,34 +999,6 @@ def perturb(sim, i, dp):
 
 class IonFlameBase(FlameBase):
 
-    def write_csv(self, filename, species='X', quiet=True):
-        """
-        Write the velocity, temperature, density, electric potential,
-        electric field strength, and species profiles to a CSV file.
-        :param filename:
-            Output file name
-        :param species:
-            Attribute to use obtaining species profiles, e.g. ``X`` for
-            mole fractions or ``Y`` for mass fractions.
-        """
-        z = self.grid
-        T = self.T
-        u = self.u
-        V = self.V
-        E = self.E
-
-        with open(filename, 'w', newline='') as csvfile:
-            writer = _csv.writer(csvfile)
-            writer.writerow(['z (m)', 'u (m/s)', 'V (1/s)', 'T (K)',
-                            'E (V/m)', 'rho (kg/m3)'] + self.gas.species_names)
-            for n in range(self.flame.n_points):
-                self.set_gas_state(n)
-                writer.writerow([z[n], u[n], V[n], T[n], E[n], self.gas.density] +
-                                list(getattr(self.gas, species)))
-
-        if not quiet:
-            print("Solution saved to '{0}'.".format(filename))
-
     @property
     def electric_field_enabled(self):
         """ Get/Set whether or not to solve the Poisson's equation."""
@@ -593,6 +1027,7 @@ def solve(self, loglevel=1, refine_grid=True, auto=False, stage=1, enable_energy
 class IonFreeFlame(IonFlameBase, FreeFlame):
     """A freely-propagating flame with ionized gas."""
     __slots__ = ('inlet', 'outlet', 'flame')
+    _extra = ('grid', 'velocity', 'eField')
 
     def __init__(self, gas, grid=None, width=None):
         if not hasattr(self, 'flame'):
@@ -606,6 +1041,7 @@ def __init__(self, gas, grid=None, width=None):
 class BurnerFlame(FlameBase):
     """A burner-stabilized flat flame."""
     __slots__ = ('burner', 'flame', 'outlet')
+    _extra = ('grid', 'velocity')
 
     def __init__(self, gas, grid=None, width=None):
         """
@@ -641,15 +1077,19 @@ def __init__(self, gas, grid=None, width=None):
         self.burner.T = gas.T
         self.burner.X = gas.X
 
-    def set_initial_guess(self):
+    def set_initial_guess(self, data=None, key=None):
         """
-        Set the initial guess for the solution. The adiabatic flame
+        Set the initial guess for the solution. By default, the adiabatic flame
         temperature and equilibrium composition are computed for the burner
         gas composition. The temperature profile rises linearly in the first
         20% of the flame to Tad, then is flat. The mass fraction profiles are
-        set similarly.
+        set similarly. Alternatively, a previously calculated result can be
+        supplied as an initial guess  via 'data' and 'key' inputs (see
+        `FlameBase.set_initial_guess`).
         """
-        super().set_initial_guess()
+        super().set_initial_guess(data=data, key=key)
+        if data:
+            return
 
         self.gas.TPY = self.burner.T, self.P, self.burner.Y
         Y0 = self.burner.Y
@@ -663,7 +1103,7 @@ def set_initial_guess(self):
         u1 = self.burner.mdot/self.gas.density
 
         locs = [0.0, 0.2, 1.0]
-        self.set_profile('u', locs, [u0, u1, u1])
+        self.set_profile('velocity', locs, [u0, u1, u1])
         self.set_profile('T', locs, [T0, Teq, Teq])
         for n in range(self.gas.n_species):
             self.set_profile(self.gas.species_name(n),
@@ -731,6 +1171,7 @@ def check_blowoff(t):
 class IonBurnerFlame(IonFlameBase, BurnerFlame):
     """A burner-stabilized flat flame with ionized gas."""
     __slots__ = ('burner', 'flame', 'outlet')
+    _extra = ('grid', 'velocity', 'eField')
 
     def __init__(self, gas, grid=None, width=None):
         if not hasattr(self, 'flame'):
@@ -744,6 +1185,7 @@ def __init__(self, gas, grid=None, width=None):
 class CounterflowDiffusionFlame(FlameBase):
     """ A counterflow diffusion flame """
     __slots__ = ('fuel_inlet', 'flame', 'oxidizer_inlet')
+    _extra = ('grid', 'velocity', 'spread_rate', 'lambda')
 
     def __init__(self, gas, grid=None, width=None):
         """
@@ -777,13 +1219,16 @@ def __init__(self, gas, grid=None, width=None):
 
         super().__init__((self.fuel_inlet, self.flame, self.oxidizer_inlet), gas, grid)
 
-    def set_initial_guess(self):
+    def set_initial_guess(self, data=None, key=None):
         """
-        Set the initial guess for the solution. The initial guess is generated
-        by assuming infinitely-fast chemistry.
+        Set the initial guess for the solution. By default, the initial guess
+        is generated by assuming infinitely-fast chemistry. Alternatively, a
+        previously calculated result can be supplied as an initial guess via
+        'data' and 'key' inputs (see `FlameBase.set_initial_guess`).
         """
-
-        super().set_initial_guess()
+        super().set_initial_guess(data=data, key=key)
+        if data:
+            return
 
         moles = lambda el: (self.gas.elemental_mass_fraction(el) /
                             self.gas.atomic_weight(el))
@@ -851,8 +1296,8 @@ def set_initial_guess(self):
         T[-1] = T0o
         zrel = (zz - zz[0])/dz
 
-        self.set_profile('u', [0.0, 1.0], [u0f, -u0o])
-        self.set_profile('V', [0.0, x0/dz, 1.0], [0.0, a, 0.0])
+        self.set_profile('velocity', [0.0, 1.0], [u0f, -u0o])
+        self.set_profile('spread_rate', [0.0, x0/dz, 1.0], [0.0, a, 0.0])
         self.set_profile('T', zrel, T)
         for k,spec in enumerate(self.gas.species_names):
             self.set_profile(spec, zrel, Y[:,k])
@@ -962,10 +1407,10 @@ def strain_rate(self, definition, fuel=None, oxidizer='O2', stoich=None):
             .. math:: a_{o} = \sqrt{-\frac{\Lambda}{\rho_{o}}}
         """
         if definition == 'mean':
-            return - (self.u[-1] - self.u[0]) / self.grid[-1]
+            return - (self.velocity[-1] - self.velocity[0]) / self.grid[-1]
 
         elif definition == 'max':
-            return np.max(np.abs(np.gradient(self.u) / np.gradient(self.grid)))
+            return np.max(np.abs(np.gradient(self.velocity) / np.gradient(self.grid)))
 
         elif definition == 'stoichiometric':
             if fuel is None:
@@ -981,7 +1426,7 @@ def strain_rate(self, definition, fuel=None, oxidizer='O2', stoich=None):
                 if 'C' in self.gas.element_names:
                     stoich += self.gas.n_atoms(fuel, 'C')
 
-            d_u_d_z = np.gradient(self.u) / np.gradient(self.grid)
+            d_u_d_z = np.gradient(self.velocity) / np.gradient(self.grid)
             phi = (self.X[self.gas.species_index(fuel)] * stoich /
                    np.maximum(self.X[self.gas.species_index(oxidizer)], 1e-20))
             z_stoich = np.interp(-1., -phi, self.grid)
@@ -1022,6 +1467,7 @@ def mixture_fraction(self, m):
 class ImpingingJet(FlameBase):
     """An axisymmetric flow impinging on a surface at normal incidence."""
     __slots__ = ('inlet', 'flame', 'surface')
+    _extra = ('grid', 'velocity', 'spread_rate', 'lambda')
 
     def __init__(self, gas, grid=None, width=None, surface=None):
         """
@@ -1064,14 +1510,18 @@ def __init__(self, gas, grid=None, width=None, surface=None):
         self.inlet.T = gas.T
         self.inlet.X = gas.X
 
-    def set_initial_guess(self, products='inlet'):
+    def set_initial_guess(self, products='inlet', data=None, key=None):
         """
         Set the initial guess for the solution. If products = 'equil', then
         the equilibrium composition at the adiabatic flame temperature will be
         used to form the initial guess. Otherwise the inlet composition will
-        be used.
+        be used. Alternatively, a previously calculated result can be supplied
+        as an initial guess via 'data' and 'key' inputs (see
+        `FlameBase.set_initial_guess`).
         """
-        super().set_initial_guess(products=products)
+        super().set_initial_guess(data=data, key=key, products=products)
+        if data:
+            return
 
         Y0 = self.inlet.Y
         T0 = self.inlet.T
@@ -1095,13 +1545,14 @@ def set_initial_guess(self, products='inlet'):
                                  [Y0[k], Y0[k]])
 
         locs = np.array([0.0, 1.0])
-        self.set_profile('u', locs, [u0, 0.0])
-        self.set_profile('V', locs, [0.0, 0.0])
+        self.set_profile('velocity', locs, [u0, 0.0])
+        self.set_profile('spread_rate', locs, [0.0, 0.0])
 
 
 class CounterflowPremixedFlame(FlameBase):
     """ A premixed counterflow flame """
     __slots__ = ('reactants', 'flame', 'products')
+    _extra = ('grid', 'velocity', 'spread_rate', 'lambda')
 
     def __init__(self, gas, grid=None, width=None):
         """
@@ -1138,15 +1589,19 @@ def __init__(self, gas, grid=None, width=None):
         # Setting X needs to be deferred until linked to the flow domain
         self.reactants.X = gas.X
 
-    def set_initial_guess(self, equilibrate=True):
+    def set_initial_guess(self, equilibrate=True, data=None, key=None):
         """
         Set the initial guess for the solution.
 
         If `equilibrate` is True, then the products composition and temperature
         will be set to the equilibrium state of the reactants mixture.
+        Alternatively, a previously calculated result can be supplied as an
+        initial guess via 'data' and 'key' inputs (see
+        `FlameBase.set_initial_guess`).
         """
-
-        super().set_initial_guess(equilibrate=equilibrate)
+        super().set_initial_guess(data=data, key=key, equilibrate=equilibrate)
+        if data:
+            return
 
         Yu = self.reactants.Y
         Tu = self.reactants.T
@@ -1185,8 +1640,8 @@ def set_initial_guess(self, equilibrate=True):
         # estimate stagnation point
         x0 = rhou*uu * dz / (rhou*uu + rhob*ub)
 
-        self.set_profile('u', [0.0, 1.0], [uu, -ub])
-        self.set_profile('V', [0.0, x0/dz, 1.0], [0.0, a, 0.0])
+        self.set_profile('velocity', [0.0, 1.0], [uu, -ub])
+        self.set_profile('spread_rate', [0.0, x0/dz, 1.0], [0.0, a, 0.0])
 
 
 class CounterflowTwinPremixedFlame(FlameBase):
@@ -1195,6 +1650,7 @@ class CounterflowTwinPremixedFlame(FlameBase):
     shooting into each other.
     """
     __slots__ = ('reactants', 'flame', 'products')
+    _extra = ('grid', 'velocity', 'spread_rate', 'lambda')
 
     def __init__(self, gas, grid=None, width=None):
         """
@@ -1231,11 +1687,14 @@ def __init__(self, gas, grid=None, width=None):
         # Setting X needs to be deferred until linked to the flow domain
         self.reactants.X = gas.X
 
-    def set_initial_guess(self):
+    def set_initial_guess(self, data=None, key=None):
         """
-        Set the initial guess for the solution.
+        Set the initial guess for the solution based on an equiibrium solution.
+        Alternatively, a previously calculated result can be supplied as an
+        initial guess via 'data' and 'key' inputs (see
+        `FlameBase.set_initial_guess`).
         """
-        super().set_initial_guess()
+        super().set_initial_guess(data=data, key=key)
 
         Yu = self.reactants.Y
         Tu = self.reactants.T
@@ -1257,5 +1716,5 @@ def set_initial_guess(self):
         dz = zz[-1] - zz[0]
         a = 2 * uu / dz
 
-        self.set_profile('u', [0.0, 1.0], [uu, 0])
-        self.set_profile('V', [0.0, 1.0], [0.0, a])
+        self.set_profile('velocity', [0.0, 1.0], [uu, 0])
+        self.set_profile('spread_rate', [0.0, 1.0], [0.0, a])
diff --git a/interfaces/cython/cantera/onedim.pyx b/interfaces/cython/cantera/onedim.pyx
index 13216763170..78334db6435 100644
--- a/interfaces/cython/cantera/onedim.pyx
+++ b/interfaces/cython/cantera/onedim.pyx
@@ -467,7 +467,26 @@ cdef class _FlowBase(Domain1D):
         del self.flow
 
     def set_boundary_emissivities(self, e_left, e_right):
-        self.flow.setBoundaryEmissivities(e_left, e_right)
+        """
+        .. deprecated:: 2.5
+
+             To be deprecated with version 2.5, and removed thereafter.
+             Replaced by property `boundary_emissivities`.
+        """
+        warnings.warn("Method 'set_boundary_emissivities' to be removed after "
+                      "Cantera 2.5. Replaced by property "
+                      "'boundary_emissivities'", DeprecationWarning)
+        self.boundary_emissivities = e_left, e_right
+
+    property boundary_emissivities:
+        """ Set/get boundary emissivities. """
+        def __get__(self):
+            return self.flow.leftEmissivity(), self.flow.rightEmissivity()
+        def __set__(self, tuple epsilon):
+            if len(epsilon) == 2:
+                self.flow.setBoundaryEmissivities(epsilon[0], epsilon[1])
+            else:
+                raise ValueError("Setter requires tuple of length 2.")
 
     property radiation_enabled:
         """ Determines whether or not to include radiative heat transfer """
@@ -476,6 +495,17 @@ cdef class _FlowBase(Domain1D):
         def __set__(self, do_radiation):
             self.flow.enableRadiation(<cbool>do_radiation)
 
+    property radiative_heat_loss:
+        """
+        Return radiative heat loss (only non-zero if radiation is enabled).
+        """
+        def __get__(self):
+            cdef int j
+            cdef np.ndarray[np.double_t, ndim=1] data = np.empty(self.n_points)
+            for j in range(self.n_points):
+                data[j] = self.flow.radiativeHeatLoss(j)
+            return data
+
     def set_free_flow(self):
         """
         Set flow configuration for freely-propagating flames, using an internal
@@ -491,6 +521,14 @@ cdef class _FlowBase(Domain1D):
         """
         self.flow.setAxisymmetricFlow()
 
+    property flow_type:
+        """
+        Return the type of flow domain being represented, either "Free Flame" or
+        "Axisymmetric Stagnation".
+        """
+        def __get__(self):
+            return pystr(self.flow.flowType())
+
 
 cdef CxxIdealGasPhase* getIdealGasPhase(ThermoPhase phase) except *:
     if pystr(phase.thermo.type()) != "IdealGas":
@@ -507,9 +545,9 @@ cdef class IdealGasFlow(_FlowBase):
     equations for the flow in a finite-height gap of infinite radial extent.
     The solution variables are:
 
-    *u*
+    *velocity*
         axial velocity
-    *V*
+    *spread_rate*
         radial velocity divided by radius
     *T*
         temperature
@@ -832,9 +870,7 @@ cdef class Sim1D:
 
     def set_initial_guess(self, *args, **kwargs):
         """
-        Set the initial guess for the solution. Derived classes extend this
-        function to set approximations for the temperature and composition
-        profiles.
+        Store arguments for initial guess and prepare storage for solution.
         """
         self._initial_guess_args = args
         self._initial_guess_kwargs = kwargs
@@ -920,7 +956,13 @@ cdef class Sim1D:
         flow_domains = [D for D in self.domains if isinstance(D, _FlowBase)]
         zmin = [D.grid[0] for D in flow_domains]
         zmax = [D.grid[-1] for D in flow_domains]
-        nPoints = [len(flow_domains[0].grid), 12, 24, 48]
+
+        # 'data' entry is used for restart
+        data = self._initial_guess_kwargs.get('data')
+        if data:
+           nPoints = [len(flow_domains[0].grid)]
+        else:
+           nPoints = [len(flow_domains[0].grid), 12, 24, 48]
 
         for N in nPoints:
             for i,D in enumerate(flow_domains):
@@ -930,8 +972,9 @@ cdef class Sim1D:
                 if N != len(D.grid):
                     D.grid = np.linspace(zmin[i], zmax[i], N)
 
-            self.set_initial_guess(*self._initial_guess_args,
-                                   **self._initial_guess_kwargs)
+            if not data:
+                self.set_initial_guess(*self._initial_guess_args,
+                                       **self._initial_guess_kwargs)
 
             # Try solving with energy enabled, which usually works
             log('Solving on {} point grid with energy equation enabled', N)
@@ -1124,8 +1167,34 @@ cdef class Sim1D:
         """
         Set the temperature used to fix the spatial location of a freely
         propagating flame.
+
+        .. deprecated:: 2.5
+
+             To be deprecated with version 2.5, and removed thereafter.
+             Replaced by property `fixed_temperature`.
+        """
+        warnings.warn("Method 'set_fixed_temperature' to be removed after "
+                      "Cantera 2.5. Replaced by property 'fixed_temperature'",
+                      DeprecationWarning)
+        self.fixed_temperature = T
+
+    property fixed_temperature:
+        """
+        Set the temperature used to fix the spatial location of a freely
+        propagating flame.
         """
-        self.sim.setFixedTemperature(T)
+        def __get__(self):
+            return self.sim.fixedTemperature()
+        def __set__(self, T):
+            self.sim.setFixedTemperature(T)
+
+    property fixed_temperature_location:
+        """
+        Return the location of the point where temperature is fixed for a freely
+        propagating flame.
+        """
+        def __get__(self):
+            return self.sim.fixedTemperatureLocation()
 
     def save(self, filename='soln.xml', name='solution', description='none',
              loglevel=1):
diff --git a/interfaces/cython/cantera/test/test_onedim.py b/interfaces/cython/cantera/test/test_onedim.py
index 4f6345c38f3..c2934847c6b 100644
--- a/interfaces/cython/cantera/test/test_onedim.py
+++ b/interfaces/cython/cantera/test/test_onedim.py
@@ -157,6 +157,44 @@ def solve_multi(self):
 
         self.assertEqual(self.sim.transport_model, 'Multi')
 
+    def test_flow_type(self):
+        Tin = 300
+        p = ct.one_atm
+        reactants = 'H2:0.65, O2:0.5, AR:2'
+        self.create_sim(p, Tin, reactants, width=0.0001)
+        self.assertEqual(self.sim.flame.flow_type, 'Free Flame')
+
+    def test_fixed_temperature(self):
+        # test setting of fixed temperature
+        Tin = 300
+        p = ct.one_atm
+        reactants = 'H2:0.65, O2:0.5, AR:2'
+        self.create_sim(p, Tin, reactants, width=0.0001)
+        self.sim.set_initial_guess()
+        zvec = self.sim.grid
+        tvec = self.sim.T
+        tfixed = 900.
+        self.sim.fixed_temperature = tfixed
+        zfixed = np.interp(tfixed, tvec, zvec)
+        self.assertNear(self.sim.fixed_temperature, tfixed)
+        self.assertNear(self.sim.fixed_temperature_location, zfixed)
+
+    def test_deprecated(self):
+        Tin = 300
+        p = ct.one_atm
+        reactants = 'H2:0.65, O2:0.5, AR:2'
+        self.create_sim(p, Tin, reactants, width=0.0001)
+        with self.assertWarnsRegex(DeprecationWarning, "Replaced by property"):
+            self.sim.flame.set_boundary_emissivities(0.5, 0.5)
+        with self.assertWarnsRegex(DeprecationWarning, "property 'velocity"):
+            self.sim.u
+        with self.assertWarnsRegex(DeprecationWarning, "property 'spread"):
+            self.sim.V
+        with self.assertRaisesRegex(ct.CanteraError, "renamed to 'velocity"):
+            self.sim.flame.component_index('u')
+        with self.assertRaisesRegex(ct.CanteraError, "renamed to 'spread_rate"):
+            self.sim.flame.component_index('V')
+
     def test_auto_width(self):
         Tin = 300
         p = ct.one_atm
@@ -183,7 +221,7 @@ def test_auto_width2(self):
 
         self.sim.set_refine_criteria(ratio=4, slope=0.8, curve=0.8)
         self.sim.solve(refine_grid=True, auto=True, loglevel=0)
-        self.assertNear(self.sim.u[0], 17.02, 1e-1)
+        self.assertNear(self.sim.velocity[0], 17.02, 1e-1)
         self.assertLess(self.sim.grid[-1], 2.0) # grid should not be too wide
 
 
@@ -230,9 +268,67 @@ def run_mix(self, phi, T, width, p, refine):
         rhou = self.sim.inlet.mdot
 
         # Check continuity
-        for rhou_j in self.sim.density * self.sim.u:
+        for rhou_j in self.sim.density * self.sim.velocity:
+            self.assertNear(rhou_j, rhou, 1e-4)
+
+    def test_solution_array_output(self):
+        self.run_mix(phi=1.0, T=300, width=2.0, p=1.0, refine=False)
+        arr = self.sim.to_solution_array()
+        ix = np.isfinite(arr.grid)
+        self.assertArrayNear(self.sim.grid, arr.grid[ix])
+        self.assertArrayNear(self.sim.T, arr.T[ix])
+        for k in arr._extra.keys():
+            self.assertIn(k, self.sim._extra)
+
+        f2 = ct.FreeFlame(self.gas)
+        f2.from_solution_array(arr)
+        self.assertArrayNear(self.sim.grid, f2.grid)
+        self.assertArrayNear(self.sim.T, f2.T)
+        self.assertNear(self.sim.inlet.T, f2.inlet.T)
+        self.assertNear(self.sim.inlet.mdot, f2.inlet.mdot)
+        self.assertArrayNear(self.sim.inlet.Y, f2.inlet.Y)
+
+    def test_restart(self):
+        self.run_mix(phi=1.0, T=300, width=2.0, p=1.0, refine=False)
+        arr = self.sim.to_solution_array()
+
+        reactants = {'H2': 0.9, 'O2': 0.5, 'AR': 2}
+        self.create_sim(1.1 * ct.one_atm, 500, reactants, 2.0)
+        self.sim.set_initial_guess(data=arr)
+        self.solve_mix(refine=False)
+
+        # Check continuity
+        rhou = self.sim.inlet.mdot
+        for rhou_j in self.sim.density * self.sim.velocity:
             self.assertNear(rhou_j, rhou, 1e-4)
 
+    def test_settings(self):
+        self.run_mix(phi=1.0, T=300, width=2.0, p=1.0, refine=False)
+        settings = self.sim.settings
+
+        keys = ['configuration', 'transport_model',
+                'energy_enabled', 'soret_enabled', 'radiation_enabled',
+                'emissivity_left', 'emissivity_right',
+                'fixed_temperature',
+                'ratio', 'slope', 'curve', 'prune',
+                'max_time_step_count',
+                'max_grid_points',
+                'steady_abstol', 'steady_reltol',
+                'transient_abstol', 'transient_reltol']
+        for k in keys:
+            self.assertIn(k, settings)
+
+        changed = {'fixed_temperature': 900,
+                   'max_time_step_count': 100,
+                   'energy_enabled': False,
+                   'radiation_enabled': True,
+                   'transport_model': 'Multi'}
+        settings.update(changed)
+
+        self.sim.settings = settings
+        for k, v in changed.items():
+            self.assertEqual(getattr(self.sim, k), v)
+
     def test_mixture_averaged_case1(self):
         self.run_mix(phi=0.65, T=300, width=0.03, p=1.0, refine=True)
 
@@ -267,15 +363,15 @@ def test_adjoint_sensitivities(self):
 
         # Forward sensitivities
         dk = 1e-4
-        Su0 = self.sim.u[0]
+        Su0 = self.sim.velocity[0]
         for m in range(self.gas.n_reactions):
             self.gas.set_multiplier(1.0) # reset all multipliers
             self.gas.set_multiplier(1+dk, m) # perturb reaction m
             self.sim.solve(loglevel=0, refine_grid=False)
-            Suplus = self.sim.u[0]
+            Suplus = self.sim.velocity[0]
             self.gas.set_multiplier(1-dk, m) # perturb reaction m
             self.sim.solve(loglevel=0, refine_grid=False)
-            Suminus = self.sim.u[0]
+            Suminus = self.sim.velocity[0]
             fwd = (Suplus-Suminus)/(2*Su0*dk)
             self.assertNear(fwd, dSdk_adj[m], 5e-3)
 
@@ -288,12 +384,12 @@ def test_multicomponent(self):
         self.create_sim(p, Tin, reactants)
         self.solve_fixed_T()
         self.solve_mix(ratio=5, slope=0.5, curve=0.3)
-        Su_mix = self.sim.u[0]
+        Su_mix = self.sim.velocity[0]
 
         # Multicomponent flame speed should be similar but not identical to
         # the mixture-averaged flame speed.
         self.solve_multi()
-        Su_multi = self.sim.u[0]
+        Su_multi = self.sim.velocity[0]
         self.assertFalse(self.sim.soret_enabled)
 
         self.assertNear(Su_mix, Su_multi, 5e-2)
@@ -304,7 +400,7 @@ def test_multicomponent(self):
         self.sim.soret_enabled = True
         self.sim.solve(loglevel=0, refine_grid=True)
         self.assertTrue(self.sim.soret_enabled)
-        Su_soret = self.sim.u[0]
+        Su_soret = self.sim.velocity[0]
 
         self.assertNear(Su_multi, Su_soret, 2e-1)
         self.assertNotEqual(Su_multi, Su_soret)
@@ -387,8 +483,8 @@ def test_save_restore(self):
             os.remove(filename)
 
         Y1 = self.sim.Y
-        u1 = self.sim.u
-        V1 = self.sim.V
+        u1 = self.sim.velocity
+        V1 = self.sim.spread_rate
         P1 = self.sim.P
 
         self.sim.save(filename, 'test', loglevel=0)
@@ -413,8 +509,8 @@ def test_save_restore(self):
         self.assertNear(P1, P2a)
 
         Y2 = self.sim.Y
-        u2 = self.sim.u
-        V2 = self.sim.V
+        u2 = self.sim.velocity
+        V2 = self.sim.spread_rate
 
         self.assertArrayNear(Y1, Y2)
         self.assertArrayNear(u1, u2)
@@ -422,8 +518,8 @@ def test_save_restore(self):
 
         self.solve_fixed_T()
         Y3 = self.sim.Y
-        u3 = self.sim.u
-        V3 = self.sim.V
+        u3 = self.sim.velocity
+        V3 = self.sim.spread_rate
 
         # TODO: These tolereances seem too loose, but the tests fail on some
         # systems with tighter tolerances.
@@ -436,7 +532,12 @@ def test_array_properties(self):
 
         for attr in ct.FlameBase.__dict__:
             if isinstance(ct.FlameBase.__dict__[attr], property):
-                getattr(self.sim, attr)
+                if attr in ['u', 'V']:
+                    msg = "Replaced by property"
+                    with self.assertWarnsRegex(DeprecationWarning, msg):
+                        getattr(self.sim, attr)
+                else:
+                    getattr(self.sim, attr)
 
     def test_save_restore_add_species(self):
         reactants = 'H2:1.1, O2:1, AR:5'
@@ -501,11 +602,12 @@ def test_write_csv(self):
 
         self.create_sim(2e5, 350, 'H2:1.0, O2:2.0', mech='h2o2.xml')
         self.sim.write_csv(filename)
-        data = np.genfromtxt(filename, delimiter=',', skip_header=1)
-        self.assertArrayNear(data[:,0], self.sim.grid)
-        self.assertArrayNear(data[:,3], self.sim.T)
+        data = ct.SolutionArray(self.gas)
+        data.read_csv(filename)
+        self.assertArrayNear(data.grid[1:], self.sim.grid)
+        self.assertArrayNear(data.T[1:], self.sim.T)
         k = self.gas.species_index('H2')
-        self.assertArrayNear(data[:,5+k], self.sim.X[k,:])
+        self.assertArrayNear(data.X[1:,k], self.sim.X[k,:])
 
     def test_refine_criteria_boundscheck(self):
         self.create_sim(ct.one_atm, 300.0, 'H2:1.1, O2:1, AR:5')
@@ -529,7 +631,7 @@ def test_refine_criteria(self):
     def test_replace_grid(self):
         self.create_sim(ct.one_atm, 300.0, 'H2:1.1, O2:1, AR:5')
         self.solve_fixed_T()
-        ub = self.sim.u[-1]
+        ub = self.sim.velocity[-1]
 
         # add some points to the grid
         grid = list(self.sim.grid)
@@ -539,7 +641,7 @@ def test_replace_grid(self):
         self.sim.set_initial_guess()
 
         self.solve_fixed_T()
-        self.assertNear(self.sim.u[-1], ub, 1e-2)
+        self.assertNear(self.sim.velocity[-1], ub, 1e-2)
 
     def test_exceed_max_grid_points(self):
         self.create_sim(ct.one_atm, 400.0, 'H2:1.1, O2:1, AR:5')
@@ -615,8 +717,8 @@ def test_mixture_averaged(self, saveReference=False):
         self.solve_mix()
         data = np.empty((self.sim.flame.n_points, self.gas.n_species + 4))
         data[:,0] = self.sim.grid
-        data[:,1] = self.sim.u
-        data[:,2] = self.sim.V
+        data[:,1] = self.sim.velocity
+        data[:,2] = self.sim.spread_rate
         data[:,3] = self.sim.T
         data[:,4:] = self.sim.Y.T
 
@@ -654,8 +756,8 @@ def time_step_func(dt):
 
         data = np.empty((self.sim.flame.n_points, self.gas.n_species + 4))
         data[:,0] = self.sim.grid
-        data[:,1] = self.sim.u
-        data[:,2] = self.sim.V
+        data[:,1] = self.sim.velocity
+        data[:,2] = self.sim.spread_rate
         data[:,3] = self.sim.T
         data[:,4:] = self.sim.Y.T
 
@@ -693,6 +795,22 @@ def test_extinction_case6(self):
     def test_extinction_case7(self):
         self.run_extinction(mdot_fuel=0.2, mdot_ox=2.0, T_ox=600, width=0.2, P=0.05)
 
+    def test_restart(self):
+        self.run_extinction(mdot_fuel=0.5, mdot_ox=3.0, T_ox=300, width=0.018, P=1.0)
+
+        arr = self.sim.to_solution_array()
+
+        self.create_sim(mdot_fuel=5.5, mdot_ox=3.3, T_ox=400, width=0.018,
+                        p=ct.one_atm*1.1)
+        self.sim.set_initial_guess(data=arr)
+        self.sim.solve(loglevel=0, auto=True)
+
+        # Check inlet
+        mdot = self.sim.density * self.sim.velocity
+        self.assertNear(mdot[0], self.sim.fuel_inlet.mdot, 1e-4)
+        self.assertNear(self.sim.T[0], self.sim.fuel_inlet.T, 1e-4)
+        self.assertNear(mdot[-1], -self.sim.oxidizer_inlet.mdot, 1e-4)
+
     def test_mixture_averaged_rad(self, saveReference=False):
         referenceFile = pjoin(self.test_data_dir, 'DiffusionFlameTest-h2-mix-rad.csv')
         self.create_sim(p=ct.one_atm)
@@ -705,16 +823,19 @@ def test_mixture_averaged_rad(self, saveReference=False):
         self.assertFalse(self.sim.radiation_enabled)
         self.sim.radiation_enabled = True
         self.assertTrue(self.sim.radiation_enabled)
-        self.sim.set_boundary_emissivities(0.25,0.15)
+        self.sim.flame.boundary_emissivities = 0.25, 0.15
 
         self.solve_mix()
         data = np.empty((self.sim.flame.n_points, self.gas.n_species + 4))
         data[:,0] = self.sim.grid
-        data[:,1] = self.sim.u
-        data[:,2] = self.sim.V
+        data[:,1] = self.sim.velocity
+        data[:,2] = self.sim.spread_rate
         data[:,3] = self.sim.T
         data[:,4:] = self.sim.Y.T
 
+        qdot = self.sim.flame.radiative_heat_loss
+        self.assertEqual(len(qdot), self.sim.flame.n_points)
+
         if saveReference:
             np.savetxt(referenceFile, data, '%11.6e', ', ')
         else:
@@ -796,8 +917,8 @@ def test_mixture_averaged(self, saveReference=False):
 
         data = np.empty((sim.flame.n_points, gas.n_species + 4))
         data[:,0] = sim.grid
-        data[:,1] = sim.u
-        data[:,2] = sim.V
+        data[:,1] = sim.velocity
+        data[:,2] = sim.spread_rate
         data[:,3] = sim.T
         data[:,4:] = sim.Y.T
 
@@ -818,7 +939,7 @@ def run_case(self, phi, T, width, P):
         sim.set_refine_criteria(ratio=6, slope=0.7, curve=0.8, prune=0.4)
         sim.solve(loglevel=0, auto=True)
         self.assertTrue(all(sim.T >= T - 1e-3))
-        self.assertTrue(all(sim.V >= -1e-9))
+        self.assertTrue(all(sim.spread_rate >= -1e-9))
         return sim
 
     def test_solve_case1(self):
@@ -836,6 +957,20 @@ def test_solve_case4(self):
     def test_solve_case5(self):
         self.run_case(phi=2.0, T=300, width=0.2, P=0.2)
 
+    def test_restart(self):
+        sim = self.run_case(phi=2.0, T=300, width=0.2, P=0.2)
+
+        arr = sim.to_solution_array()
+        sim.reactants.mdot *= 1.1
+        sim.products.mdot *= 1.1
+        sim.set_initial_guess(data=arr)
+        sim.solve(loglevel=0, auto=True)
+
+        # Check inlet / outlet
+        mdot = sim.density * sim.velocity
+        self.assertNear(mdot[0], sim.reactants.mdot, 1e-4)
+        self.assertNear(mdot[-1], -sim.products.mdot, 1e-4)
+
 
 class TestBurnerFlame(utilities.CanteraTest):
     def solve(self, phi, T, width, P):
@@ -885,9 +1020,27 @@ def test_blowoff(self):
         sim.solve(loglevel=0, auto=True)
         # nonreacting solution
         self.assertNear(sim.T[-1], sim.T[0], 1e-6)
-        self.assertNear(sim.u[-1], sim.u[0], 1e-6)
+        self.assertNear(sim.velocity[-1], sim.velocity[0], 1e-6)
         self.assertArrayNear(sim.Y[:,0], sim.Y[:,-1], 1e-6, atol=1e-6)
 
+    def test_restart(self):
+        gas = ct.Solution('h2o2.cti')
+        gas.set_equivalence_ratio(0.4, 'H2', 'O2:1.0, AR:5')
+        gas.TP = 300, ct.one_atm
+        sim = ct.BurnerFlame(gas=gas, width=0.1)
+        sim.burner.mdot = 1.2
+        sim.set_refine_criteria(ratio=3, slope=0.3, curve=0.5, prune=0)
+        sim.solve(loglevel=0, auto=True)
+
+        arr = sim.to_solution_array()
+        sim.burner.mdot = 1.1
+        sim.set_initial_guess(data=arr)
+        sim.solve(loglevel=0, auto=True)
+
+        # Check continuity
+        rhou = sim.burner.mdot
+        for rhou_j in sim.density * sim.velocity:
+            self.assertNear(rhou_j, rhou, 1e-4)
 
 class TestImpingingJet(utilities.CanteraTest):
     def run_reacting_surface(self, xch4, tsurf, mdot, width):
@@ -940,10 +1093,26 @@ def solve(self, phi, T, width, P):
         sim.reactants.mdot = gas.density * axial_velocity
         sim.solve(loglevel=0, auto=True)
         self.assertGreater(sim.T[-1], T + 100)
+        return sim
 
     def test_case1(self):
         self.solve(phi=0.4, T=300, width=0.05, P=0.1)
 
+    def test_restart(self):
+        sim = self.solve(phi=0.4, T=300, width=0.05, P=0.1)
+
+        arr = sim.to_solution_array()
+        axial_velocity = 2.2
+        sim.reactants.mdot *= 1.1
+        sim.reactants.T = sim.reactants.T + 100
+        sim.set_initial_guess(data=arr)
+        sim.solve(loglevel=0, auto=True)
+
+        # Check inlet
+        mdot = sim.density * sim.velocity
+        self.assertNear(mdot[0], sim.reactants.mdot, 1e-4)
+        self.assertNear(sim.T[0], sim.reactants.T, 1e-4)
+
 
 class TestIonFreeFlame(utilities.CanteraTest):
     def test_ion_profile(self):
diff --git a/samples/cxx/flamespeed/flamespeed.cpp b/samples/cxx/flamespeed/flamespeed.cpp
index d687671d7c6..c6a219a1011 100644
--- a/samples/cxx/flamespeed/flamespeed.cpp
+++ b/samples/cxx/flamespeed/flamespeed.cpp
@@ -101,7 +101,7 @@ int flamespeed(double phi)
 
         double uout = inlet.mdot()/rho_out;
         value = {uin, uin, uout, uout};
-        flame.setInitialGuess("u",locs,value);
+        flame.setInitialGuess("velocity",locs,value);
         value = {temp, temp, Tad, Tad};
         flame.setInitialGuess("T",locs,value);
 
@@ -135,21 +135,24 @@ int flamespeed(double phi)
         bool refine_grid = true;
 
         flame.solve(loglevel,refine_grid);
-        double flameSpeed_mix = flame.value(flowdomain,flow.componentIndex("u"),0);
+        double flameSpeed_mix = flame.value(flowdomain,
+                                            flow.componentIndex("velocity"),0);
         print("Flame speed with mixture-averaged transport: {} m/s\n",
               flameSpeed_mix);
 
         // now switch to multicomponent transport
         flow.setTransport(*trmulti);
         flame.solve(loglevel, refine_grid);
-        double flameSpeed_multi = flame.value(flowdomain,flow.componentIndex("u"),0);
+        double flameSpeed_multi = flame.value(flowdomain,
+                                              flow.componentIndex("velocity"),0);
         print("Flame speed with multicomponent transport: {} m/s\n",
               flameSpeed_multi);
 
         // now enable Soret diffusion
         flow.enableSoret(true);
         flame.solve(loglevel, refine_grid);
-        double flameSpeed_full = flame.value(flowdomain,flow.componentIndex("u"),0);
+        double flameSpeed_full = flame.value(flowdomain,
+                                             flow.componentIndex("velocity"),0);
         print("Flame speed with multicomponent transport + Soret: {} m/s\n",
               flameSpeed_full);
 
@@ -159,9 +162,12 @@ int flamespeed(double phi)
               "z (m)", "T (K)", "U (m/s)", "Y(CO)");
         for (size_t n = 0; n < flow.nPoints(); n++) {
             Tvec.push_back(flame.value(flowdomain,flow.componentIndex("T"),n));
-            COvec.push_back(flame.value(flowdomain,flow.componentIndex("CO"),n));
-            CO2vec.push_back(flame.value(flowdomain,flow.componentIndex("CO2"),n));
-            Uvec.push_back(flame.value(flowdomain,flow.componentIndex("u"),n));
+            COvec.push_back(flame.value(flowdomain,
+                                        flow.componentIndex("CO"),n));
+            CO2vec.push_back(flame.value(flowdomain,
+                                         flow.componentIndex("CO2"),n));
+            Uvec.push_back(flame.value(flowdomain,
+                                       flow.componentIndex("velocity"),n));
             zvec.push_back(flow.grid(n));
             print("{:9.6f}\t{:8.3f}\t{:5.3f}\t{:7.5f}\n",
                   flow.grid(n), Tvec[n], Uvec[n], COvec[n]);
diff --git a/src/oneD/Sim1D.cpp b/src/oneD/Sim1D.cpp
index 5231bdf35ed..70972062b54 100644
--- a/src/oneD/Sim1D.cpp
+++ b/src/oneD/Sim1D.cpp
@@ -12,6 +12,7 @@
 #include "cantera/numerics/funcs.h"
 #include "cantera/base/xml.h"
 #include "cantera/numerics/Func1.h"
+#include <limits>
 
 using namespace std;
 
@@ -424,7 +425,7 @@ int Sim1D::refine(int loglevel)
     return np;
 }
 
-int Sim1D::setFixedTemperature(doublereal t)
+int Sim1D::setFixedTemperature(double t)
 {
     int np = 0;
     vector_fp znew, xnew;
@@ -511,8 +512,34 @@ int Sim1D::setFixedTemperature(doublereal t)
     return np;
 }
 
-void Sim1D::setRefineCriteria(int dom, doublereal ratio,
-                              doublereal slope, doublereal curve, doublereal prune)
+double Sim1D::fixedTemperature()
+{
+    double t_fixed = std::numeric_limits<double>::quiet_NaN();
+    for (size_t n = 0; n < nDomains(); n++) {
+        StFlow* d = dynamic_cast<StFlow*>(&domain(n));
+        if (d && d->domainType() == cFreeFlow && d->m_tfixed > 0) {
+            t_fixed = d->m_tfixed;
+            break;
+        }
+    }
+    return t_fixed;
+}
+
+double Sim1D::fixedTemperatureLocation()
+{
+    double z_fixed = std::numeric_limits<double>::quiet_NaN();
+    for (size_t n = 0; n < nDomains(); n++) {
+        StFlow* d = dynamic_cast<StFlow*>(&domain(n));
+        if (d && d->domainType() == cFreeFlow && d->m_tfixed > 0) {
+            z_fixed = d->m_zfixed;
+            break;
+        }
+    }
+    return z_fixed;
+}
+
+void Sim1D::setRefineCriteria(int dom, double ratio,
+                              double slope, double curve, double prune)
 {
     if (dom >= 0) {
         Refiner& r = domain(dom).refiner();
diff --git a/src/oneD/StFlow.cpp b/src/oneD/StFlow.cpp
index 1185cc8d827..df9b1e7d65f 100644
--- a/src/oneD/StFlow.cpp
+++ b/src/oneD/StFlow.cpp
@@ -28,7 +28,7 @@ StFlow::StFlow(ThermoPhase* ph, size_t nsp, size_t points) :
     m_kExcessLeft(0),
     m_kExcessRight(0),
     m_zfixed(Undef),
-    m_tfixed(Undef)
+    m_tfixed(-1.)
 {
     if (ph->type() == "IdealGas") {
         m_thermo = static_cast<IdealGasPhase*>(ph);
@@ -574,9 +574,9 @@ string StFlow::componentName(size_t n) const
 {
     switch (n) {
     case 0:
-        return "u";
+        return "velocity";
     case 1:
-        return "V";
+        return "spread_rate";
     case 2:
         return "T";
     case 3:
@@ -595,8 +595,18 @@ string StFlow::componentName(size_t n) const
 size_t StFlow::componentIndex(const std::string& name) const
 {
     if (name=="u") {
+        warn_deprecated("StFlow::componentIndex",
+                        "To be changed after Cantera 2.5. "
+                        "Solution component 'u' renamed to 'velocity'");
+        return 0;
+    } else if (name=="velocity") {
         return 0;
     } else if (name=="V") {
+        warn_deprecated("StFlow::componentIndex",
+                        "To be changed after Cantera 2.5. "
+                        "Solution component 'V' renamed to 'spread_rate'");
+        return 1;
+    } else if (name=="spread_rate") {
         return 1;
     } else if (name=="T") {
         return 2;
@@ -610,8 +620,9 @@ size_t StFlow::componentIndex(const std::string& name) const
                 return n;
             }
         }
+        throw CanteraError("StFlow1D::componentIndex",
+                           "no component named " + name);
     }
-    return npos;
 }
 
 void StFlow::restore(const XML_Node& dom, doublereal* soln, int loglevel)