diff --git a/.gitignore b/.gitignore
index e98536dfb0..aeb98c555c 100644
--- a/.gitignore
+++ b/.gitignore
@@ -47,3 +47,4 @@ Testing
 *.pdf
 *.gv
 *.json
+*.gz
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 325dfaabc7..aefe007b14 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -273,6 +273,15 @@ if(NEML2_DOC)
       FetchContent_MakeAvailable(doxygen-awesome-css)
 endif()
 
+# ----------------------------------------------------------------------------
+# zlib
+# ----------------------------------------------------------------------------
+find_package(ZLIB)
+
+if(NOT ZLIB_FOUND)
+      message(WARNING "ZLIB not found, model packaging features will be disabled")
+endif()
+
 # ----------------------------------------------------------------------------
 # base neml2 library
 # ----------------------------------------------------------------------------
diff --git a/include/neml2/base/Factory.h b/include/neml2/base/Factory.h
index 622439075e..4fcb208202 100644
--- a/include/neml2/base/Factory.h
+++ b/include/neml2/base/Factory.h
@@ -26,6 +26,7 @@
 
 #include <filesystem>
 #include <iostream>
+#include <memory>
 
 #include "neml2/misc/errors.h"
 #include "neml2/base/InputFile.h"
@@ -41,6 +42,7 @@ class Data;
 class Model;
 class Driver;
 class WorkScheduler;
+class BundledModel;
 
 /**
  * @brief A convenient function to parse all options from an input file
@@ -73,7 +75,7 @@ class Factory
   const InputFile & input_file() const { return _input_file; }
 
   /// Global settings
-  const std::shared_ptr<Settings> & settings() const { return _input_file.settings(); }
+  const std::shared_ptr<Settings> & settings() const { return _settings; }
 
   /// Check if an object with the given name exists under the given section.
   bool has_object(const std::string & section, const std::string & name);
@@ -116,6 +118,17 @@ class Factory
   template <class T = WorkScheduler>
   std::shared_ptr<T> get_scheduler(const std::string & name);
 
+  /**
+   * @brief Serialize an object to an input file. The returned input file contains the exact
+   * information needed to reconstruct the object.
+   *
+   * @note Behind the scenes, this method calls the get_object method with \p force_create set to
+   * true, which has the side effect of creating the object if it does not already exist.
+   */
+  std::unique_ptr<InputFile> serialize_object(const std::string & section,
+                                              const std::string & name,
+                                              const OptionSet & additional_options = OptionSet());
+
   /// @brief Delete all factories and destruct all the objects.
   void clear();
 
@@ -139,14 +152,26 @@ class Factory
   /// Check if the options are compatible with the object
   bool options_compatible(const std::shared_ptr<NEML2Object> & obj, const OptionSet & opts) const;
 
+  /// BundledModel will need to squeeze the unpacked model into the factory
+  friend class BundledModel;
+
   /// The input file
   InputFile _input_file;
 
+  /// Global settings of the input file
+  const std::shared_ptr<Settings> _settings;
+
   /**
    * Manufactured objects. The key of the outer map is the section name, and the key of the inner
    * map is the object name.
    */
   std::map<std::string, std::map<std::string, std::vector<std::shared_ptr<NEML2Object>>>> _objects;
+
+  /// Whether the factory is currently serializing an object
+  bool _serializing = false;
+
+  /// The output serialized input file (used by the serialize_object method)
+  std::unique_ptr<InputFile> _serialized_file;
 };
 
 template <class T>
@@ -160,7 +185,7 @@ Factory::get_object(const std::string & section,
     throw FactoryException("The input file is empty.");
 
   // Easy if it already exists
-  if (!force_create)
+  if (!force_create && !_serializing)
     if (_objects.count(section) && _objects.at(section).count(name))
       for (const auto & neml2_obj : _objects[section][name])
       {
@@ -183,8 +208,8 @@ Factory::get_object(const std::string & section,
     if (options.first == name)
     {
       auto new_options = options.second;
-      new_options.set<Factory *>("_factory") = this;
-      new_options.set<std::shared_ptr<Settings>>("_settings") = settings();
+      new_options.set<Factory *>("factory") = this;
+      new_options.set<std::shared_ptr<Settings>>("settings") = settings();
       new_options += additional_options;
       create_object(section, new_options);
       break;
diff --git a/include/neml2/base/HITParser.h b/include/neml2/base/HITParser.h
index e713265b1a..c6eaefad06 100644
--- a/include/neml2/base/HITParser.h
+++ b/include/neml2/base/HITParser.h
@@ -29,6 +29,7 @@
 namespace neml2
 {
 class OptionSet;
+class OptionBase;
 
 /**
  * @copydoc neml2::Parser
@@ -47,11 +48,14 @@ class HITParser : public Parser
   ~HITParser() override = default;
 
   /// Parse a HIT input file from a filename.
-  InputFile parse(const std::filesystem::path & filename,
-                  const std::string & additional_input = "") const override;
+  InputFile parse_from_string(const std::string & input,
+                              const std::string & additional_input = "") const override;
 
   /// Parse a HIT input file from a root node.
-  InputFile parse(hit::Node * root) const;
+  InputFile parse_from_hit_node(hit::Node * root) const;
+
+  /// Serialize an input file to a string.
+  std::string serialize(const InputFile & inp) const override;
 
 private:
   /**
@@ -61,10 +65,12 @@ class HITParser : public Parser
    * @param section The current section node.
    * @return OptionSet The options of the object.
    */
-  virtual OptionSet extract_object_options(hit::Node * object, hit::Node * section) const;
+  OptionSet extract_object_options(hit::Node * object, hit::Node * section) const;
 
   void extract_options(hit::Node * object, OptionSet & options) const;
   void extract_option(hit::Node * node, OptionSet & options) const;
+
+  void serialize_options(hit::Node * node, const OptionSet & options) const;
 };
 
 } // namespace neml2
diff --git a/include/neml2/base/InputFile.h b/include/neml2/base/InputFile.h
index 5b255d661f..41aeaa899c 100644
--- a/include/neml2/base/InputFile.h
+++ b/include/neml2/base/InputFile.h
@@ -28,18 +28,16 @@
 
 namespace neml2
 {
-class Settings;
-
 /**
  * @brief A data structure that holds options of multiple objects.
  */
 class InputFile
 {
 public:
-  InputFile(const OptionSet & settings);
+  InputFile(OptionSet settings);
 
   /// Get global settings
-  const std::shared_ptr<Settings> & settings() const { return _settings; }
+  const OptionSet & settings() const { return _settings; }
 
   /// Get all the object options under a specific section.
   std::map<std::string, OptionSet> & operator[](const std::string & section);
@@ -55,7 +53,7 @@ class InputFile
 
 private:
   /// Global settings specified under the [Settings] section
-  const std::shared_ptr<Settings> _settings;
+  const OptionSet _settings;
 
   /// Collection of options for all manufacturable objects
   std::map<std::string, std::map<std::string, OptionSet>> _data;
diff --git a/include/neml2/base/Option.h b/include/neml2/base/Option.h
index e9ff8821a7..19fba25ba4 100644
--- a/include/neml2/base/Option.h
+++ b/include/neml2/base/Option.h
@@ -24,6 +24,7 @@
 
 #pragma once
 
+#include <cstddef>
 #include <vector>
 
 #include "neml2/base/OptionBase.h"
@@ -49,6 +50,9 @@ template <typename P>
 void _print_helper(std::ostream & os, const std::vector<P> *);
 template <typename P>
 void _print_helper(std::ostream & os, const std::vector<std::vector<P>> *);
+/// bool
+template <>
+void _print_helper(std::ostream & os, const bool *);
 /// The evil vector of bool :/
 template <>
 void _print_helper(std::ostream & os, const std::vector<bool> *);
@@ -58,6 +62,9 @@ void _print_helper(std::ostream & os, const char *);
 /// Specialization so that we don't print out unprintable characters
 template <>
 void _print_helper(std::ostream & os, const unsigned char *);
+/// Specialization for tensor shape
+template <>
+void _print_helper(std::ostream & os, const TensorShape *);
 ///@}
 }
 
@@ -108,16 +115,24 @@ template <typename P>
 void
 _print_helper(std::ostream & os, const std::vector<P> * option)
 {
-  for (const auto & p : *option)
-    os << p << " ";
+  for (std::size_t i = 0; i < option->size(); i++)
+  {
+    if (i > 0)
+      os << " ";
+    _print_helper(os, &(*option)[i]);
+  }
 }
 
 template <typename P>
 void
 _print_helper(std::ostream & os, const std::vector<std::vector<P>> * option)
 {
-  for (const auto & pv : *option)
-    _print_helper(os, &pv);
+  for (std::size_t i = 0; i < option->size(); i++)
+  {
+    if (i > 0)
+      os << "; ";
+    _print_helper(os, &(*option)[i]);
+  }
 }
 } // namespace details
 // LCOV_EXCL_STOP
diff --git a/include/neml2/base/Parser.h b/include/neml2/base/Parser.h
index 6cf52b181f..ee750bd502 100644
--- a/include/neml2/base/Parser.h
+++ b/include/neml2/base/Parser.h
@@ -66,7 +66,13 @@ class Parser
    * @return InputFile The extracted object options.
    */
   virtual InputFile parse(const std::filesystem::path & filename,
-                          const std::string & additional_input = "") const = 0;
+                          const std::string & additional_input = "") const;
+
+  virtual InputFile parse_from_string(const std::string & input,
+                                      const std::string & additional_input = "") const = 0;
+
+  /// @brief Serialize an input file to a string
+  virtual std::string serialize(const InputFile & inp) const = 0;
 };
 
 namespace utils
@@ -117,7 +123,7 @@ parse_vector_(std::vector<T> & vals, const std::string & raw_str)
     vals.resize(tokens.size(), kCPU);
   else
     vals.resize(tokens.size());
-  for (size_t i = 0; i < tokens.size(); i++)
+  for (std::size_t i = 0; i < tokens.size(); i++)
   {
     auto success = parse_<T>(vals[i], tokens[i]);
     if (!success)
@@ -144,7 +150,7 @@ parse_vector_vector_(std::vector<std::vector<T>> & vals, const std::string & raw
 {
   auto token_vecs = split(raw_str, ";");
   vals.resize(token_vecs.size());
-  for (size_t i = 0; i < token_vecs.size(); i++)
+  for (std::size_t i = 0; i < token_vecs.size(); i++)
   {
     auto success = parse_vector_<T>(vals[i], token_vecs[i]);
     if (!success)
diff --git a/include/neml2/base/Registry.h b/include/neml2/base/Registry.h
index b1eb75f8ae..5120a863c8 100644
--- a/include/neml2/base/Registry.h
+++ b/include/neml2/base/Registry.h
@@ -83,6 +83,9 @@ class Registry
   /// Load registry from a dynamic library
   static void load(const std::filesystem::path &);
 
+  /// Check if an object is registered in the registry.
+  static bool is_registered(const std::string & name);
+
   /// Get information of all registered objects.
   static const std::map<std::string, NEML2ObjectInfo> & info();
 
diff --git a/include/neml2/drivers/ModelDriver.h b/include/neml2/drivers/ModelDriver.h
index 8f12d4d4c5..f3a4fdb660 100644
--- a/include/neml2/drivers/ModelDriver.h
+++ b/include/neml2/drivers/ModelDriver.h
@@ -61,12 +61,8 @@ class ModelDriver : public Driver
   /// The device on which to evaluate the model
   const Device _device;
 
-  /// Set to true to list all the model parameters at the beginning
-  const bool _show_params;
-  /// Set to true to show model's input axis at the beginning
-  const bool _show_input;
-  /// Set to true to show model's output axis at the beginning
-  const bool _show_output;
+  /// Set to true to show model summary at the beginning
+  const bool _show_model;
 
 #ifdef NEML2_HAS_DISPATCHER
   /// The work scheduler to use
diff --git a/include/neml2/models/BundledModel.h b/include/neml2/models/BundledModel.h
new file mode 100644
index 0000000000..313c4522b1
--- /dev/null
+++ b/include/neml2/models/BundledModel.h
@@ -0,0 +1,80 @@
+// Copyright 2024, UChicago Argonne, LLC
+// All Rights Reserved
+// Software Name: NEML2 -- the New Engineering material Model Library, version 2
+// By: Argonne National Laboratory
+// OPEN SOURCE LICENSE (MIT)
+//
+// Permission is hereby granted, free of charge, to any person obtaining a copy
+// of this software and associated documentation files (the "Software"), to deal
+// in the Software without restriction, including without limitation the rights
+// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the Software is
+// furnished to do so, subject to the following conditions:
+//
+// The above copyright notice and this permission notice shall be included in
+// all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+// THE SOFTWARE.
+
+#pragma once
+
+#include <filesystem>
+
+#ifdef NEML2_HAS_ZLIB
+#include <zlib.h>
+#endif
+
+#ifdef NEML2_HAS_JSON
+#include "nlohmann/json.hpp"
+#endif
+
+#define NEML2_CAN_BUNDLE_MODEL defined(NEML2_HAS_ZLIB) && defined(NEML2_HAS_JSON)
+
+#include "neml2/models/Model.h"
+
+namespace neml2
+{
+#ifdef NEML2_CAN_BUNDLE_MODEL
+void bundle_model(const std::string & file,
+                  const std::string & name,
+                  const std::string & cliargs = "",
+                  const nlohmann::json & config = nlohmann::json(),
+                  std::filesystem::path output_path = std::filesystem::path());
+
+std::pair<std::shared_ptr<Model>, nlohmann::json> unbundle_model(const std::filesystem::path & pkg,
+                                                                 NEML2Object * host = nullptr);
+#endif // NEML2_CAN_BUNDLE_MODEL
+
+class BundledModel : public Model
+{
+public:
+  static OptionSet expected_options();
+
+  BundledModel(const OptionSet & options);
+
+  const nlohmann::json & config() const { return _config; }
+
+  ///@{
+  /// Methods for retrieving descriptions
+  std::string description() const;
+  std::string input_description(const VariableName & name) const;
+  std::string output_description(const VariableName & name) const;
+  std::string param_description(const std::string & name) const;
+  std::string buffer_description(const std::string & name) const;
+  ///@}
+
+protected:
+  void link_output_variables() override;
+  void set_value(bool, bool, bool) override;
+
+  std::shared_ptr<Model> _bundled_model;
+
+  nlohmann::json _config;
+};
+} // namespace neml2
diff --git a/include/neml2/models/Data.h b/include/neml2/models/Data.h
index c78f49acd6..17b60c1069 100644
--- a/include/neml2/models/Data.h
+++ b/include/neml2/models/Data.h
@@ -59,11 +59,11 @@ class Data : public NEML2Object, public BufferStore
   T & register_data(const std::string & name)
   {
     OptionSet extra_opts;
-    extra_opts.set<NEML2Object *>("_host") = host();
-    if (!host()->factory())
-      throw SetupException("Internal error: Host object '" + host()->name() +
-                           "' does not have a factory set.");
-    auto data = host()->factory()->get_object<T>("Data", name, extra_opts, /*force_create=*/false);
+    extra_opts.set<NEML2Object *>("host") = host();
+    if (!factory())
+      throw SetupException("Internal error: Object '" + this->name() +
+                           "' does not have a factory.");
+    auto data = factory()->get_object<T>("Data", name, extra_opts, /*force_create=*/false);
 
     if (std::find(_registered_data.begin(), _registered_data.end(), data) != _registered_data.end())
       throw SetupException("Data named '" + name + "' has already been registered.");
diff --git a/include/neml2/models/Model.h b/include/neml2/models/Model.h
index 87022ee10d..404ead42c6 100644
--- a/include/neml2/models/Model.h
+++ b/include/neml2/models/Model.h
@@ -192,7 +192,7 @@ class Model : public std::enable_shared_from_this<Model>,
   /// Declaration of nonlinear parameters may require manipulation of input
   friend class ParameterStore;
 
-  /// ComposedModel's set_value need to call submodel's set_value
+  /// ComposedModel::set_value need to call submodel's set_value
   friend class ComposedModel;
 
 protected:
@@ -260,13 +260,13 @@ class Model : public std::enable_shared_from_this<Model>,
                            "' is trying to register itself as a sub-model. This is not allowed.");
 
     OptionSet extra_opts;
-    extra_opts.set<NEML2Object *>("_host") = host();
-    extra_opts.set<bool>("_nonlinear_system") = nonlinear;
+    extra_opts.set<NEML2Object *>("host") = host();
+    extra_opts.set<bool>("nonlinear_system") = nonlinear;
 
-    if (!host()->factory())
-      throw SetupException("Internal error: Host object '" + host()->name() +
-                           "' does not have a factory set.");
-    auto model = host()->factory()->get_object<T>("Models", model_name, extra_opts);
+    if (!factory())
+      throw SetupException("Internal error: Object '" + this->name() +
+                           "' does not have a factory.");
+    auto model = factory()->get_object<T>("Models", model_name, extra_opts);
     if (std::find(_registered_models.begin(), _registered_models.end(), model) !=
         _registered_models.end())
       throw SetupException("Model named '" + model_name + "' has already been registered.");
diff --git a/include/neml2/models/solid_mechanics/elasticity/ElasticityInterface.h b/include/neml2/models/solid_mechanics/elasticity/ElasticityInterface.h
index 52b0d89a68..bbd70f461a 100644
--- a/include/neml2/models/solid_mechanics/elasticity/ElasticityInterface.h
+++ b/include/neml2/models/solid_mechanics/elasticity/ElasticityInterface.h
@@ -26,10 +26,13 @@
 
 #include "neml2/models/solid_mechanics/elasticity/ElasticityConverter.h"
 #include "neml2/base/MultiEnumSelection.h"
+#include "neml2/base/OptionSet.h"
+#include "neml2/base/TensorName.h"
 #include "neml2/tensors/Scalar.h"
 
 namespace neml2
 {
+
 /**
  * @brief Interface for objects defining elasticity tensors in terms of other parameters
  *
diff --git a/runner/src/main.cxx b/runner/src/main.cxx
index 96dfc321be..f1fd4ecb3e 100644
--- a/runner/src/main.cxx
+++ b/runner/src/main.cxx
@@ -26,9 +26,11 @@
 #include "neml2/base/guards.h"
 #include "neml2/drivers/Driver.h"
 #include "neml2/models/Model.h"
+#include "neml2/models/BundledModel.h"
 #include "neml2/misc/errors.h"
 
 #include <argparse/argparse.hpp>
+#include <stdexcept>
 
 std::string get_additional_cliargs(const argparse::ArgumentParser & program);
 
@@ -86,11 +88,34 @@ main(int argc, char * argv[])
       .remaining()
       .help("additional command-line arguments to pass to the input file parser");
 
+// sub-command: bundle
+#ifdef NEML2_CAN_BUNDLE_MODEL
+
+  argparse::ArgumentParser bundle_command("bundle");
+  bundle_command.add_description("Bundle a model into a single archive.");
+  bundle_command.add_argument("input").help("path to the input file");
+  bundle_command.add_argument("model").help("name of the model in the input file to bundle");
+  bundle_command.add_argument("additional_args")
+      .remaining()
+      .help("additional command-line arguments to pass to the input file parser");
+  bundle_command.add_argument("-c", "--config")
+      .default_value("")
+      .help("path to the bundling configuration file");
+  bundle_command.add_argument("-o", "--output")
+      .default_value("")
+      .help(
+          "path to the output file where the bundled model will be saved. If not provided, the "
+          "output file will be named <input>_<model>.gz in the same directory as the input file.");
+#endif
+
   // Add sub-commands to the main program
   program.add_subparser(run_command);
   program.add_subparser(diagnose_command);
   program.add_subparser(inspect_command);
   program.add_subparser(time_command);
+#ifdef NEML2_CAN_BUNDLE_MODEL
+  program.add_subparser(bundle_command);
+#endif
 
   try
   {
@@ -187,6 +212,27 @@ main(int argc, char * argv[])
           std::cout << "    " << object << ": " << time << " ms" << std::endl;
       }
     }
+
+    // sub-command: bundle
+#ifdef NEML2_CAN_BUNDLE_MODEL
+    if (program.is_subcommand_used("bundle"))
+    {
+      const auto input = bundle_command.get<std::string>("input");
+      const auto additional_cliargs = get_additional_cliargs(bundle_command);
+      const auto modelname = bundle_command.get<std::string>("model");
+      const auto outpath = bundle_command.get<std::string>("--output");
+      nlohmann::json config;
+      if (bundle_command.is_used("--config"))
+      {
+        const auto config_file = bundle_command.get<std::string>("--config");
+        std::ifstream f(config_file);
+        if (!f)
+          throw std::runtime_error("Could not open config file: " + config_file);
+        f >> config;
+      }
+      neml2::bundle_model(input, modelname, additional_cliargs, config, outpath);
+    }
+#endif
   }
   catch (const std::exception & err)
   {
diff --git a/src/neml2/CMakeLists.txt b/src/neml2/CMakeLists.txt
index fe878c824d..23f26c7f04 100644
--- a/src/neml2/CMakeLists.txt
+++ b/src/neml2/CMakeLists.txt
@@ -121,6 +121,11 @@ if(NEML2_PCH)
   )
 endif()
 
+if(ZLIB_FOUND)
+  target_compile_definitions(neml2_model PUBLIC NEML2_HAS_ZLIB)
+  target_link_libraries(neml2_model PUBLIC ZLIB::ZLIB)
+endif()
+
 # libneml2_driver
 neml2_add_submodule(neml2_driver SHARED drivers)
 target_link_libraries(neml2_driver PUBLIC neml2_model neml2_tensor neml2_base)
diff --git a/src/neml2/base/Factory.cxx b/src/neml2/base/Factory.cxx
index e8d35ba0fb..3d44ec2906 100644
--- a/src/neml2/base/Factory.cxx
+++ b/src/neml2/base/Factory.cxx
@@ -23,6 +23,8 @@
 // THE SOFTWARE.
 
 #include "neml2/base/Factory.h"
+#include "neml2/base/InputFile.h"
+#include "neml2/base/Settings.h"
 #include "neml2/base/OptionSet.h"
 #include "neml2/base/NEML2Object.h"
 #include "neml2/base/Registry.h"
@@ -46,6 +48,7 @@ load_input(const std::filesystem::path & path, const std::string & additional_in
 
 Factory::Factory(InputFile inp)
   : _input_file(std::move(inp)),
+    _settings(std::make_shared<Settings>(_input_file.settings())),
     _objects()
 {
 }
@@ -75,6 +78,22 @@ Factory::create_object(const std::string & section, const OptionSet & options)
     throw FactoryException("Failed to setup object '" + name + "' of type '" + type +
                            "' in section '" + section + "':\n" + e.what());
   }
+
+  // Record the object options in the serialized file if we are serializing
+  if (_serializing)
+    (*_serialized_file)[section][name] = object->input_options();
+}
+
+std::unique_ptr<InputFile>
+Factory::serialize_object(const std::string & section,
+                          const std::string & name,
+                          const OptionSet & additional_options)
+{
+  _serialized_file = std::make_unique<InputFile>(_input_file.settings());
+  _serializing = true;
+  get_object<NEML2Object>(section, name, additional_options);
+  _serializing = false;
+  return std::move(_serialized_file);
 }
 
 bool
diff --git a/src/neml2/base/HITParser.cxx b/src/neml2/base/HITParser.cxx
index c78628a6c1..4f8dcfe163 100644
--- a/src/neml2/base/HITParser.cxx
+++ b/src/neml2/base/HITParser.cxx
@@ -23,6 +23,7 @@
 // THE SOFTWARE.
 
 #include "hit/braceexpr.h"
+#include "hit/parse.h"
 
 #include "neml2/base/HITParser.h"
 #include "neml2/base/Registry.h"
@@ -32,6 +33,8 @@
 #include "neml2/base/EnumSelection.h"
 #include "neml2/base/MultiEnumSelection.h"
 #include "neml2/base/LabeledAxisAccessor.h"
+#include "neml2/misc/errors.h"
+#include "neml2/misc/string_utils.h"
 #include "neml2/tensors/tensors.h"
 #include "neml2/misc/assertions.h"
 #include "neml2/misc/types.h"
@@ -39,20 +42,11 @@
 namespace neml2
 {
 InputFile
-HITParser::parse(const std::filesystem::path & filename, const std::string & additional_input) const
+HITParser::parse_from_string(const std::string & input, const std::string & additional_input) const
 {
-  // Open and read the file
-  std::ifstream file(filename);
-  neml_assert(file.is_open(), "Unable to open file ", filename);
-
-  // Read the file into a string
-  std::stringstream buffer;
-  buffer << file.rdbuf();
-  std::string input = buffer.str();
-
   // Let HIT lex the string
-  std::unique_ptr<hit::Node> root(hit::parse(filename, input));
-  neml_assert(root.get(), "HIT failed to lex the input file: ", filename);
+  std::unique_ptr<hit::Node> root(hit::parse("neml2 input file", input));
+  neml_assert(root.get(), "HIT failed to lex the input file.");
 
   // Handle additional input (they could be coming from cli args)
   std::unique_ptr<hit::Node> cli_root(hit::parse("cliargs", additional_input));
@@ -66,11 +60,11 @@ HITParser::parse(const std::filesystem::path & filename, const std::string & add
   expander.registerEvaler("raw", raw);
   root->walk(&expander);
 
-  return parse(root.get());
+  return parse_from_hit_node(root.get());
 }
 
 InputFile
-HITParser::parse(hit::Node * root) const
+HITParser::parse_from_hit_node(hit::Node * root) const
 {
   // Extract global settings
   OptionSet settings = Settings::expected_options();
@@ -102,8 +96,16 @@ HITParser::extract_object_options(hit::Node * object, hit::Node * section) const
   // There is a special field reserved for object type
   std::string type = object->param<std::string>("type");
   // Extract the options
-  auto options = Registry::info(type).expected_options;
-  extract_options(object, options);
+  OptionSet options;
+  if (Registry::is_registered(type))
+  {
+    options = Registry::info(type).expected_options;
+    extract_options(object, options);
+  }
+  else
+    std::cerr << "Warning: Object type '" << type
+              << "' is not registered in the NEML2 registry. "
+                 "This may lead to unexpected behavior or errors.\n";
 
   // Also fill in the metadata
   options.name() = object->path();
@@ -178,4 +180,53 @@ HITParser::extract_option(hit::Node * n, OptionSet & options) const
   }
 }
 // NOLINTEND
+
+std::string
+HITParser::serialize(const InputFile & inp) const
+{
+  // Serialized input
+  std::string output;
+
+  // Serialize settings
+  auto node = std::make_unique<hit::Section>("Settings");
+  serialize_options(node.get(), inp.settings());
+  output += node->render();
+
+  // Serialize each section
+  for (const auto & [section_name, section] : inp.data())
+  {
+    output += "\n\n"; // Add a newline before each section for readability
+    auto section_node = std::make_unique<hit::Section>(section_name);
+    for (const auto & [object_name, options] : section)
+    {
+      // Create a section for the object
+      auto object_node = std::make_unique<hit::Section>(object_name);
+      // Add the special field "type". The destructor will delete the children nodes, so it's
+      // generally safe to use raw pointers for children
+      object_node->addChild(new hit::Field("type", hit::Field::Kind::String, options.type()));
+      // Serialize the options
+      serialize_options(object_node.get(), options);
+      // Add the object node to the section
+      section_node->addChild(object_node.release());
+    }
+    output += section_node->render();
+  }
+
+  return output;
+}
+
+void
+HITParser::serialize_options(hit::Node * node, const OptionSet & options) const
+{
+  for (const auto & [key, val] : options)
+  {
+    if (val->suppressed() || !val->user_specified())
+      continue;
+
+    // auto * field = serialize_option(key, val.get());
+    auto val_str = "'" + utils::stringify(*val) + "'"; // Wrap in single quotes
+    auto * field = new hit::Field(key, hit::Field::Kind::String, val_str);
+    node->addChild(field);
+  }
+}
 } // namespace neml2
diff --git a/src/neml2/base/InputFile.cxx b/src/neml2/base/InputFile.cxx
index 2f676271d9..7103dba730 100644
--- a/src/neml2/base/InputFile.cxx
+++ b/src/neml2/base/InputFile.cxx
@@ -25,12 +25,11 @@
 #include <iostream>
 
 #include "neml2/base/InputFile.h"
-#include "neml2/base/Settings.h"
 
 namespace neml2
 {
-InputFile::InputFile(const OptionSet & settings)
-  : _settings(std::make_shared<Settings>(settings)),
+InputFile::InputFile(OptionSet settings)
+  : _settings(std::move(settings)),
     _data()
 {
 }
diff --git a/src/neml2/base/MultiEnumSelection.cxx b/src/neml2/base/MultiEnumSelection.cxx
index c940a473e4..f1c9fa48d5 100644
--- a/src/neml2/base/MultiEnumSelection.cxx
+++ b/src/neml2/base/MultiEnumSelection.cxx
@@ -32,14 +32,12 @@ std::ostream &
 operator<<(std::ostream & os, const MultiEnumSelection & es)
 {
   auto selections = std::vector<std::string>(es);
-  os << '(';
   for (size_t i = 0; i < selections.size(); i++)
   {
     os << selections[i];
     if (i < selections.size() - 1)
-      os << ", ";
+      os << " ";
   }
-  os << ')';
   return os;
 }
 
diff --git a/src/neml2/base/NEML2Object.cxx b/src/neml2/base/NEML2Object.cxx
index ab70db3b7d..a051df65ed 100644
--- a/src/neml2/base/NEML2Object.cxx
+++ b/src/neml2/base/NEML2Object.cxx
@@ -34,23 +34,23 @@ NEML2Object::expected_options()
 {
   auto options = OptionSet();
 
-  options.set<Factory *>("_factory") = nullptr;
-  options.set("_factory").suppressed() = true;
+  options.set<Factory *>("factory") = nullptr;
+  options.set("factory").suppressed() = true;
 
-  options.set<std::shared_ptr<Settings>>("_settings") = nullptr;
-  options.set("_settings").suppressed() = true;
+  options.set<std::shared_ptr<Settings>>("settings") = nullptr;
+  options.set("settings").suppressed() = true;
 
-  options.set<NEML2Object *>("_host") = nullptr;
-  options.set("_host").suppressed() = true;
+  options.set<NEML2Object *>("host") = nullptr;
+  options.set("host").suppressed() = true;
 
   return options;
 }
 
 NEML2Object::NEML2Object(const OptionSet & options)
   : _input_options(options),
-    _factory(options.get<Factory *>("_factory")),
-    _settings(options.get<std::shared_ptr<Settings>>("_settings")),
-    _host(options.get<NEML2Object *>("_host"))
+    _factory(options.get<Factory *>("factory")),
+    _settings(options.get<std::shared_ptr<Settings>>("settings")),
+    _host(options.get<NEML2Object *>("host"))
 {
 }
 
diff --git a/src/neml2/base/Option.cxx b/src/neml2/base/Option.cxx
index 43c23ea28e..8141e8a8a7 100644
--- a/src/neml2/base/Option.cxx
+++ b/src/neml2/base/Option.cxx
@@ -162,12 +162,37 @@ _print_helper(std::ostream & os, const unsigned char * option)
   os << static_cast<int>(*option);
 }
 
+template <>
+void
+_print_helper(std::ostream & os, const bool * option)
+{
+  os << (*option ? "true" : "false");
+}
+
 template <>
 void
 _print_helper(std::ostream & os, const std::vector<bool> * option)
 {
-  for (const auto p : *option)
-    os << static_cast<bool>(p) << " ";
+  for (std::size_t i = 0; i < option->size(); i++)
+  {
+    if (i > 0)
+      os << " ";
+    os << ((*option)[i] ? "true" : "false");
+  }
+}
+
+template <>
+void
+_print_helper(std::ostream & os, const TensorShape * option)
+{
+  os << "(";
+  for (std::size_t i = 0; i < option->size(); i++)
+  {
+    if (i > 0)
+      os << ",";
+    os << (*option)[i];
+  }
+  os << ")";
 }
 } // namespace details
 } // namespace neml2
diff --git a/src/neml2/base/Parser.cxx b/src/neml2/base/Parser.cxx
index 0ca75e9754..0297c039ed 100644
--- a/src/neml2/base/Parser.cxx
+++ b/src/neml2/base/Parser.cxx
@@ -23,15 +23,33 @@
 // THE SOFTWARE.
 
 #include <iostream>
+#include <fstream>
 
 #include "neml2/base/Parser.h"
+#include "neml2/misc/assertions.h"
 #include "neml2/base/LabeledAxisAccessor.h"
+#include "neml2/base/InputFile.h"
 
 namespace neml2
 {
 const std::vector<std::string> Parser::sections = {
     "Tensors", "Solvers", "Data", "Models", "Drivers", "Schedulers"};
 
+InputFile
+Parser::parse(const std::filesystem::path & filename, const std::string & additional_input) const
+{
+  // Open and read the file
+  std::ifstream file(filename);
+  neml_assert(file.is_open(), "Unable to open file ", filename);
+
+  // Read the file into a string
+  std::stringstream buffer;
+  buffer << file.rdbuf();
+  std::string input = buffer.str();
+
+  return parse_from_string(input, additional_input);
+}
+
 namespace utils
 {
 template <>
diff --git a/src/neml2/base/Registry.cxx b/src/neml2/base/Registry.cxx
index a9a0fa2b69..51de0a26a7 100644
--- a/src/neml2/base/Registry.cxx
+++ b/src/neml2/base/Registry.cxx
@@ -66,6 +66,13 @@ Registry::load(const std::filesystem::path & lib)
               dlerror());
 }
 
+bool
+Registry::is_registered(const std::string & name)
+{
+  const auto & reg = get();
+  return reg._info.count(name) > 0;
+}
+
 const std::map<std::string, NEML2ObjectInfo> &
 Registry::info()
 {
diff --git a/src/neml2/drivers/ModelDriver.cxx b/src/neml2/drivers/ModelDriver.cxx
index 05dbb7c32a..414a652e30 100644
--- a/src/neml2/drivers/ModelDriver.cxx
+++ b/src/neml2/drivers/ModelDriver.cxx
@@ -48,12 +48,8 @@ ModelDriver::expected_options()
       "target compute device to be CPU, and device='cuda:1' sets the target compute device to be "
       "CUDA with device ID 1.";
 
-  options.set<bool>("show_parameters") = false;
-  options.set("show_parameters").doc() = "Whether to show model parameters at the beginning";
-  options.set<bool>("show_input_axis") = false;
-  options.set("show_input_axis").doc() = "Whether to show model input axis at the beginning";
-  options.set<bool>("show_output_axis") = false;
-  options.set("show_output_axis").doc() = "Whether to show model output axis at the beginning";
+  options.set<bool>("show_model") = false;
+  options.set("show_model").doc() = "Display a summary of the model being tested.";
 
 #ifdef NEML2_HAS_DISPATCHER
   options.set<std::string>("scheduler");
@@ -69,9 +65,7 @@ ModelDriver::ModelDriver(const OptionSet & options)
   : Driver(options),
     _model(get_model("model")),
     _device(options.get<std::string>("device")),
-    _show_params(options.get<bool>("show_parameters")),
-    _show_input(options.get<bool>("show_input_axis")),
-    _show_output(options.get<bool>("show_output_axis"))
+    _show_model(options.get<bool>("show_model"))
 #ifdef NEML2_HAS_DISPATCHER
     ,
     _scheduler(options.get("scheduler").user_specified() ? get_scheduler("scheduler") : nullptr),
@@ -125,18 +119,8 @@ ModelDriver::setup()
 #endif
 
   // LCOV_EXCL_START
-  if (_show_input)
-    std::cout << _model->name() << "'s input axis:\n" << _model->input_axis() << std::endl;
-
-  if (_show_output)
-    std::cout << _model->name() << "'s output axis:\n" << _model->output_axis() << std::endl;
-
-  if (_show_params)
-  {
-    std::cout << _model->name() << "'s parameters:" << std::endl;
-    for (auto && [pname, pval] : _model->named_parameters())
-      std::cout << "  " << pname << std::endl;
-  }
+  if (_show_model)
+    std::cout << *_model << std::endl;
   // LCOV_EXCL_STOP
 }
 
diff --git a/src/neml2/drivers/TransientDriver.cxx b/src/neml2/drivers/TransientDriver.cxx
index 2c23b386f6..dbace3d42b 100644
--- a/src/neml2/drivers/TransientDriver.cxx
+++ b/src/neml2/drivers/TransientDriver.cxx
@@ -60,7 +60,7 @@ set_ic(ValueMap & storage,
               " and ",
               vals.size(),
               " respectively.");
-  auto * factory = options.get<Factory *>("_factory");
+  auto * factory = options.get<Factory *>("factory");
   neml_assert(factory, "Internal error: factory != nullptr");
   for (std::size_t i = 0; i < names.size(); i++)
   {
@@ -92,7 +92,7 @@ get_force(std::vector<VariableName> & names,
               " and ",
               vals.size(),
               " respectively.");
-  auto * factory = options.get<Factory *>("_factory");
+  auto * factory = options.get<Factory *>("factory");
   neml_assert(factory, "Internal error: factory != nullptr");
   for (std::size_t i = 0; i < force_names.size(); i++)
   {
diff --git a/src/neml2/models/BundledModel.cxx b/src/neml2/models/BundledModel.cxx
new file mode 100644
index 0000000000..71387a1356
--- /dev/null
+++ b/src/neml2/models/BundledModel.cxx
@@ -0,0 +1,276 @@
+// Copyright 2024, UChicago Argonne, LLC
+// All Rights Reserved
+// Software Name: NEML2 -- the New Engineering material Model Library, version 2
+// By: Argonne National Laboratory
+// OPEN SOURCE LICENSE (MIT)
+//
+// Permission is hereby granted, free of charge, to any person obtaining a copy
+// of this software and associated documentation files (the "Software"), to deal
+// in the Software without restriction, including without limitation the rights
+// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the Software is
+// furnished to do so, subject to the following conditions:
+//
+// The above copyright notice and this permission notice shall be included in
+// all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+// THE SOFTWARE.
+
+#include "neml2/models/BundledModel.h"
+#include "neml2/base/HITParser.h"
+
+namespace neml2
+{
+static std::string
+pack_strings(const std::vector<std::string> & strings)
+{
+  std::string packed;
+  for (const auto & s : strings)
+  {
+    std::uint32_t len = static_cast<std::uint32_t>(s.size());
+    packed.append(reinterpret_cast<const char *>(&len), sizeof(len));
+    packed.append(s);
+  }
+  return packed;
+}
+
+static std::vector<std::string>
+unpack_strings(const std::string & packed)
+{
+  std::vector<std::string> strings;
+  std::size_t offset = 0;
+  while (offset + sizeof(std::uint32_t) <= packed.size())
+  {
+    std::uint32_t len = 0;
+    std::memcpy(&len, packed.data() + offset, sizeof(len));
+    offset += sizeof(len);
+    if (offset + len > packed.size())
+      throw NEMLException("Invalid packed string format");
+    strings.emplace_back(packed.data() + offset, len);
+    offset += len;
+  }
+  return strings;
+}
+
+#ifdef NEML2_CAN_BUNDLE_MODEL
+void
+bundle_model(const std::string & file,
+             const std::string & name,
+             const std::string & cliargs,
+             const nlohmann::json & config,
+             std::filesystem::path output_path)
+{
+  namespace fs = std::filesystem;
+
+  // load the input file and serialize the model
+  auto factory = load_input(file, cliargs);
+  auto inp = factory->serialize_object("Models", name);
+  auto parser = HITParser(); // For now we always use HIT for serialization
+  auto inp_str = parser.serialize(*inp);
+  auto model = factory->get_model(name);
+
+  // default output name
+  if (output_path.empty())
+    output_path =
+        fs::path(file).parent_path() / (fs::path(file).stem().string() + "_" + name + ".gz");
+
+  // create the header
+  nlohmann::json header;
+  header["schema"] = 1;
+  header["name"] = name;
+
+  if (config.contains("desc"))
+    header["desc"] = config["desc"];
+
+  for (auto && [vname, var] : model->input_variables())
+    if (config.contains("inputs") && config["inputs"].contains(vname.str()) &&
+        config["inputs"][vname.str()].contains("desc"))
+      header["inputs"][vname.str()]["desc"] = config["inputs"][vname.str()]["desc"];
+
+  for (auto && [vname, var] : model->output_variables())
+    if (config.contains("outputs") && config["outputs"].contains(vname.str()) &&
+        config["outputs"][vname.str()].contains("desc"))
+      header["outputs"][vname.str()]["desc"] = config["outputs"][vname.str()]["desc"];
+
+  for (auto && [pname, param] : model->named_parameters())
+    if (config.contains("parameters") && config["parameters"].contains(pname) &&
+        config["parameters"][pname].contains("desc"))
+      header["parameters"][pname]["desc"] = config["parameters"][pname]["desc"];
+
+  for (auto && [bname, buffer] : model->named_buffers())
+    if (config.contains("buffers") && config["buffers"].contains(bname) &&
+        config["buffers"][bname].contains("desc"))
+      header["buffers"][bname]["desc"] = config["buffers"][bname]["desc"];
+
+  // pack strings for serialization
+  const auto packed_str = pack_strings({header.dump(), inp_str});
+
+  // stream the serialized model to a file
+  gzFile out_file = gzopen(fs::absolute(output_path).c_str(), "wb");
+  if (!out_file)
+    throw NEMLException("Failed to open output file for writing: " + output_path.string());
+  if (gzwrite(out_file, packed_str.c_str(), packed_str.size()) !=
+      static_cast<int>(packed_str.size()))
+  {
+    gzclose(out_file);
+    throw NEMLException("Failed to write to output file: " + output_path.string());
+  }
+  if (gzclose(out_file) != Z_OK)
+    throw NEMLException("Failed to close output file: " + output_path.string());
+}
+
+std::pair<std::shared_ptr<Model>, nlohmann::json>
+unbundle_model(const std::filesystem::path & pkg, NEML2Object * host)
+{
+  namespace fs = std::filesystem;
+
+  // deserialize from gz
+  gzFile in_file = gzopen(fs::absolute(pkg).c_str(), "rb");
+  if (!in_file)
+    throw NEMLException("Failed to open packaged model file for reading: " + pkg.string());
+  std::string packed_str;
+  std::array<char, 4096> buffer{};
+  int bytes_read = 0;
+  while ((bytes_read = gzread(in_file, buffer.data(), buffer.size())) > 0)
+    packed_str.append(buffer.data(), bytes_read);
+  if (gzclose(in_file) != Z_OK)
+    throw NEMLException("Failed to close input file: " + pkg.string());
+
+  // unpack the strings
+  auto unpacked = unpack_strings(packed_str);
+  if (unpacked.size() != 2)
+    throw NEMLException("Invalid packed model format: expected 2 strings, got " +
+                        std::to_string(unpacked.size()));
+  const auto & header_str = unpacked[0];
+  const auto & inp_str = unpacked[1];
+
+  // parse the header
+  auto config = nlohmann::json::parse(header_str);
+
+  // parse the input string
+  auto parser = HITParser(); // For now we always use HIT for serialization
+  auto inp = parser.parse_from_string(inp_str);
+  auto factory = std::make_unique<Factory>(inp);
+  OptionSet additional_options;
+  additional_options.set<NEML2Object *>("host") = host;
+  auto model = factory->get_object<Model>("Models", config["name"], additional_options);
+
+  return {std::move(model), config};
+}
+#endif // NEML2_CAN_BUNDLE_MODEL
+
+register_NEML2_object(BundledModel);
+
+OptionSet
+BundledModel::expected_options()
+{
+  OptionSet options = Model::expected_options();
+  NonlinearSystem::enable_automatic_scaling(options);
+  options.doc() =
+      "Deserialize a model from an bundle and use it as a new model. The deserialized model "
+      "is a 'black box' and can be used in the same way as any other models.";
+
+  options.set<std::string>("bundle");
+  options.set("bundle").doc() =
+      "Path to the bundled model file. The file must be a valid bundled model file generated by "
+      "NEML2. The path can either be relative or absolute. If a relative path is provided, it will "
+      "be resolved against the current working directory.";
+
+  options.set<bool>("jit") = false;
+  options.set("jit").suppressed() = true;
+
+  return options;
+}
+
+BundledModel::BundledModel(const OptionSet & options)
+  : Model(options)
+{
+#ifdef NEML2_CAN_BUNDLE_MODEL
+  // Unpack the model
+  std::tie(_bundled_model, _config) = unbundle_model(options.get<std::string>("bundle"), host());
+
+  // register this model so that it travels with this model
+  _registered_models.push_back(_bundled_model);
+
+  // clone the input and output variables to the current model
+  for (auto && [name, var] : _bundled_model->input_variables())
+    clone_input_variable(*var);
+  for (auto && [name, var] : _bundled_model->output_variables())
+    clone_output_variable(*var);
+#else
+  throw NEMLException("NEML2 was not built with support for bundled models.");
+#endif // NEML2_CAN_BUNDLE_MODEL
+}
+
+std::string
+BundledModel::description() const
+{
+  if (_config.contains("desc"))
+    return _config["desc"].get<std::string>();
+  return "";
+}
+
+std::string
+BundledModel::input_description(const VariableName & name) const
+{
+  if (_config.contains("inputs") && _config["inputs"].contains(name.str()) &&
+      _config["inputs"][name.str()].contains("desc"))
+    return _config["inputs"][name.str()]["desc"].get<std::string>();
+  return "";
+}
+
+std::string
+BundledModel::output_description(const VariableName & name) const
+{
+  if (_config.contains("outputs") && _config["outputs"].contains(name.str()) &&
+      _config["outputs"][name.str()].contains("desc"))
+    return _config["outputs"][name.str()]["desc"].get<std::string>();
+  return "";
+}
+
+std::string
+BundledModel::param_description(const std::string & name) const
+{
+  if (_config.contains("parameters") && _config["parameters"].contains(name) &&
+      _config["parameters"][name].contains("desc"))
+    return _config["parameters"][name]["desc"].get<std::string>();
+  return "";
+}
+
+std::string
+BundledModel::buffer_description(const std::string & name) const
+{
+  if (_config.contains("buffers") && _config["buffers"].contains(name) &&
+      _config["buffers"][name].contains("desc"))
+    return _config["buffers"][name]["desc"].get<std::string>();
+  return "";
+}
+
+void
+BundledModel::link_output_variables()
+{
+  Model::link_output_variables();
+  for (auto && [name, var] : output_variables())
+    var->ref(_bundled_model->output_variable(name));
+}
+
+void
+BundledModel::set_value(bool out, bool dout_din, bool d2out_din2)
+{
+  _bundled_model->forward_maybe_jit(out, dout_din, d2out_din2);
+
+  // copy derivatives
+  if (dout_din)
+    for (auto && [name, var] : output_variables())
+      var->derivatives() = _bundled_model->output_variable(name).derivatives();
+  if (d2out_din2)
+    for (auto && [name, var] : output_variables())
+      var->second_derivatives() = _bundled_model->output_variable(name).second_derivatives();
+}
+} // namespace neml2
diff --git a/src/neml2/models/Model.cxx b/src/neml2/models/Model.cxx
index 9b148388de..936fb9a81a 100644
--- a/src/neml2/models/Model.cxx
+++ b/src/neml2/models/Model.cxx
@@ -35,6 +35,7 @@
 #include "neml2/models/Model.h"
 #include "neml2/models/Assembler.h"
 #include "neml2/models/map_types_fwd.h"
+#include "neml2/models/BundledModel.h"
 
 namespace neml2
 {
@@ -79,8 +80,8 @@ Model::expected_options()
   // Model defaults to _not_ being part of a nonlinear system
   // Model::get_model will set this to true if the model is expected to be part of a nonlinear
   // system, and additional diagnostics will be performed
-  options.set<bool>("_nonlinear_system") = false;
-  options.set("_nonlinear_system").suppressed() = true;
+  options.set<bool>("nonlinear_system") = false;
+  options.set("nonlinear_system").suppressed() = true;
 
   options.set<bool>("jit") = true;
   options.set("jit").doc() = "Use JIT compilation for the forward operator";
@@ -103,7 +104,7 @@ Model::Model(const OptionSet & options)
     _defines_value(options.get<bool>("define_values")),
     _defines_dvalue(options.get<bool>("define_derivatives")),
     _defines_d2value(options.get<bool>("define_second_derivatives")),
-    _nonlinear_system(options.get<bool>("_nonlinear_system")),
+    _nonlinear_system(options.get<bool>("nonlinear_system")),
     _jit(options.get<bool>("jit")),
     _production(options.get<bool>("production"))
 {
@@ -862,57 +863,72 @@ Model::extract_AD_derivatives(bool dout, bool d2out)
 std::ostream &
 operator<<(std::ostream & os, const Model & model)
 {
+  const auto * bundled_model = dynamic_cast<const BundledModel *>(&model);
   bool first = false;
-  const std::string tab = "            ";
+  const std::string tab = "             ";
 
-  os << "Name:       " << model.name() << '\n';
+  os << "Name:        " << model.name() << '\n';
+  if (bundled_model)
+    os << "Description: " << bundled_model->description() << '\n';
 
   if (!model.input_variables().empty())
   {
-    os << "Input:      ";
+    os << "Input:       ";
     first = true;
     for (auto && [name, var] : model.input_variables())
     {
       os << (first ? "" : tab);
-      os << name << " [" << var->type() << "]\n";
+      os << name << " [" << var->type() << "]";
+      if (bundled_model)
+        os << " (" << bundled_model->input_description(name) << ")";
+      os << '\n';
       first = false;
     }
   }
 
-  if (!model.input_variables().empty())
+  if (!model.output_variables().empty())
   {
-    os << "Output:     ";
+    os << "Output:      ";
     first = true;
     for (auto && [name, var] : model.output_variables())
     {
       os << (first ? "" : tab);
-      os << name << " [" << var->type() << "]\n";
+      os << name << " [" << var->type() << "]";
+      if (bundled_model)
+        os << " (" << bundled_model->output_description(name) << ")";
+      os << '\n';
       first = false;
     }
   }
 
   if (!model.named_parameters().empty())
   {
-    os << "Parameters: ";
+    os << "Parameters:  ";
     first = true;
     for (auto && [name, param] : model.named_parameters())
     {
       os << (first ? "" : tab);
       os << name << " [" << param->type() << "][" << Tensor(*param).scalar_type() << "]["
-         << Tensor(*param).device() << "]\n";
+         << Tensor(*param).device() << "]";
+      if (bundled_model)
+        os << " (" << bundled_model->param_description(name) << ")";
+      os << '\n';
       first = false;
     }
   }
 
   if (!model.named_buffers().empty())
   {
-    os << "Buffers:    ";
+    os << "Buffers:     ";
     first = true;
     for (auto && [name, buffer] : model.named_buffers())
     {
       os << (first ? "" : tab);
       os << name << " [" << buffer->type() << "][" << Tensor(*buffer).scalar_type() << "]["
-         << Tensor(*buffer).device() << "]\n";
+         << Tensor(*buffer).device() << "]";
+      if (bundled_model)
+        os << " (" << bundled_model->buffer_description(name) << ")";
+      os << '\n';
       first = false;
     }
   }
diff --git a/src/neml2/models/ParameterStore.cxx b/src/neml2/models/ParameterStore.cxx
index 15f2889d11..5dd9ed6f02 100644
--- a/src/neml2/models/ParameterStore.cxx
+++ b/src/neml2/models/ParameterStore.cxx
@@ -119,7 +119,7 @@ resolve_tensor_name(const TensorName<T> & tn, Model * caller, const std::string
     // When we retrieve a model, we want it to register its own parameters and buffers in the
     // host of the caller.
     OptionSet extra_opts;
-    extra_opts.set<NEML2Object *>("_host") = caller->host();
+    extra_opts.set<NEML2Object *>("host") = caller->host();
 
     // The raw string is interpreted as a _variable specifier_ which takes three possible forms
     // 1. "model_name.variable_name"
diff --git a/src/neml2/models/crystallography/CrystalGeometry.cxx b/src/neml2/models/crystallography/CrystalGeometry.cxx
index 147f0a672b..0be47bfa27 100644
--- a/src/neml2/models/crystallography/CrystalGeometry.cxx
+++ b/src/neml2/models/crystallography/CrystalGeometry.cxx
@@ -65,7 +65,7 @@ CrystalGeometry::expected_options()
 }
 
 CrystalGeometry::CrystalGeometry(const OptionSet & options)
-  : CrystalGeometry(options, options.get<Factory *>("_factory"))
+  : CrystalGeometry(options, options.get<Factory *>("factory"))
 {
 }
 
diff --git a/src/neml2/models/crystallography/CubicCrystal.cxx b/src/neml2/models/crystallography/CubicCrystal.cxx
index 49023c8f1b..1b4cfe2675 100644
--- a/src/neml2/models/crystallography/CubicCrystal.cxx
+++ b/src/neml2/models/crystallography/CubicCrystal.cxx
@@ -50,7 +50,7 @@ CubicCrystal::expected_options()
 }
 
 CubicCrystal::CubicCrystal(const OptionSet & options)
-  : CubicCrystal(options, options.get<Factory *>("_factory"))
+  : CubicCrystal(options, options.get<Factory *>("factory"))
 {
 }
 
diff --git a/tests/include/ModelUnitTest.h b/tests/include/ModelUnitTest.h
index c31bf8fd45..3d5b07cf79 100644
--- a/tests/include/ModelUnitTest.h
+++ b/tests/include/ModelUnitTest.h
@@ -64,8 +64,6 @@ class ModelUnitTest : public Driver
   double _param_rtol;
   double _param_atol;
 
-  const bool _show_params;
-  const bool _show_input;
-  const bool _show_output;
+  const bool _show_model;
 };
 } // namespace neml2
diff --git a/tests/src/ModelUnitTest.cxx b/tests/src/ModelUnitTest.cxx
index 9da26396cc..6c6f065fba 100644
--- a/tests/src/ModelUnitTest.cxx
+++ b/tests/src/ModelUnitTest.cxx
@@ -44,7 +44,7 @@ set_variable(ValueMap & storage,
               " values to ",
               vars.size(),
               " variables.");
-  auto * factory = options.get<Factory *>("_factory");
+  auto * factory = options.get<Factory *>("factory");
   neml_assert(factory, "Internal error: factory != nullptr");
   for (size_t i = 0; i < vars.size(); i++)
     storage[vars[i]] = vals[i].resolve(factory);
@@ -79,12 +79,8 @@ ModelUnitTest::expected_options()
   options.set<std::vector<TensorName<T>>>("output_" #T "_values")
   FOR_ALL_TENSORBASE(OPTION_SET_);
 
-  options.set<bool>("show_parameters") = false;
-  options.set("show_parameters").doc() = "Whether to show model parameters at the beginning";
-  options.set<bool>("show_input_axis") = false;
-  options.set("show_input_axis").doc() = "Whether to show model input axis at the beginning";
-  options.set<bool>("show_output_axis") = false;
-  options.set("show_output_axis").doc() = "Whether to show model output axis at the beginning";
+  options.set<bool>("show_model") = false;
+  options.set("show_model").doc() = "Display a summary of the model being tested.";
 
   return options;
 }
@@ -106,9 +102,7 @@ ModelUnitTest::ModelUnitTest(const OptionSet & options)
     _param_rtol(options.get<double>("parameter_derivative_rel_tol")),
     _param_atol(options.get<double>("parameter_derivative_abs_tol")),
 
-    _show_params(options.get<bool>("show_parameters")),
-    _show_input(options.get<bool>("show_input_axis")),
-    _show_output(options.get<bool>("show_output_axis"))
+    _show_model(options.get<bool>("show_model"))
 {
 #define SET_VARIABLE_(T)                                                                           \
   set_variable<T>(_in, options, "input_" #T "_names", "input_" #T "_values");                      \
@@ -120,18 +114,8 @@ bool
 ModelUnitTest::run()
 {
   // LCOV_EXCL_START
-  if (_show_params)
-  {
-    std::cout << _model->name() << "'s parameters:\n";
-    for (auto && [pname, pval] : _model->named_parameters())
-      std::cout << "  " << pname << std::endl;
-  }
-
-  if (_show_input)
-    std::cout << _model->name() << "'s input axis:\n" << _model->input_axis() << std::endl;
-
-  if (_show_output)
-    std::cout << _model->name() << "'s output axis:\n" << _model->output_axis() << std::endl;
+  if (_show_model)
+    std::cout << *_model << std::endl;
   // LCOV_EXCL_STOP
 
   check_all();
diff --git a/tests/unit/base/test_HITParser.cxx b/tests/unit/base/test_HITParser.cxx
index 9d0c1835d6..24702980dd 100644
--- a/tests/unit/base/test_HITParser.cxx
+++ b/tests/unit/base/test_HITParser.cxx
@@ -76,10 +76,10 @@ TEST_CASE("HITParser", "[base]")
 
       SECTION("global settings")
       {
-        auto & settings = inp.settings();
-        REQUIRE(settings->buffer_name_separator() == "::");
-        REQUIRE(settings->parameter_name_separator() == "::");
-        REQUIRE(!settings->require_double_precision());
+        const auto & settings = inp.settings();
+        REQUIRE(settings.get<std::string>("buffer_name_separator") == "::");
+        REQUIRE(settings.get<std::string>("parameter_name_separator") == "::");
+        REQUIRE(!settings.get<bool>("require_double_precision"));
       }
 
       SECTION("booleans")
@@ -160,6 +160,13 @@ TEST_CASE("HITParser", "[base]")
       }
     }
 
+    SECTION("serialize")
+    {
+      auto inp = parser.parse("base/test_HITParser1.i");
+      std::string out = parser.serialize(inp);
+      std::cout << out << std::endl;
+    }
+
     SECTION("error")
     {
       SECTION("setting a suppressed option")
diff --git a/tests/unit/base/test_Settings.cxx b/tests/unit/base/test_Settings.cxx
index 3edb091743..f26fa613f9 100644
--- a/tests/unit/base/test_Settings.cxx
+++ b/tests/unit/base/test_Settings.cxx
@@ -25,21 +25,19 @@
 #include <catch2/catch_test_macros.hpp>
 #include <catch2/matchers/catch_matchers_all.hpp>
 
+#include "neml2/base/Factory.h"
 #include "neml2/base/Settings.h"
-#include "neml2/base/HITParser.h"
-#include "neml2/base/InputFile.h"
 
 using namespace neml2;
 
 TEST_CASE("Settings", "[Settings]")
 {
   // Parse input file
-  HITParser parser;
-  auto input = parser.parse("base/test_HITParser1.i");
+  auto factory = neml2::load_input("base/test_HITParser1.i");
 
   // After applying the global settings
-  const auto settings = *input.settings();
-  REQUIRE(settings.buffer_name_separator() == "::");
-  REQUIRE(settings.parameter_name_separator() == "::");
-  REQUIRE(!settings.require_double_precision());
+  const auto & settings = factory->settings();
+  REQUIRE(settings->buffer_name_separator() == "::");
+  REQUIRE(settings->parameter_name_separator() == "::");
+  REQUIRE(!settings->require_double_precision());
 }
diff --git a/tests/unit/models/BundledModel.i b/tests/unit/models/BundledModel.i
new file mode 100644
index 0000000000..5a450b5b81
--- /dev/null
+++ b/tests/unit/models/BundledModel.i
@@ -0,0 +1,32 @@
+[Tensors]
+  [T_vals]
+    type = Scalar
+    values = '0 100 200'
+    batch_shape = '(3)'
+  []
+  [c_A_vals]
+    type = Scalar
+    values = '1 2 3'
+    batch_shape = '(3)'
+  []
+[]
+
+[Drivers]
+  [unit]
+    type = ModelUnitTest
+    model = 'model'
+    input_Scalar_names = 'forces/T state/A state/B'
+    input_Scalar_values = '100 0.5 1.5'
+    output_Scalar_names = 'state/C'
+    output_Scalar_values = '9.25'
+    check_AD_parameter_derivatives = false
+    show_model = true
+  []
+[]
+
+[Models]
+  [model]
+    type = BundledModel
+    bundle = 'ComposedModel5_model.gz'
+  []
+[]
diff --git a/tests/unit/models/ComposedModel5.json b/tests/unit/models/ComposedModel5.json
new file mode 100644
index 0000000000..7b366f1006
--- /dev/null
+++ b/tests/unit/models/ComposedModel5.json
@@ -0,0 +1,35 @@
+{
+  "desc": "A composed model for testing",
+  "inputs": {
+    "forces/T": {
+      "desc": "Prescribed temperature used to interpolate the prefactor for the first input variable"
+    },
+    "state/A": {
+      "desc": "The first input variable"
+    },
+    "state/B": {
+      "desc": "The second input variable"
+    }
+  },
+  "outputs": {
+    "state/C": {
+      "desc": "The only output variable"
+    }
+  },
+  "parameters": {
+    "c_A_X": {
+      "desc": "Abscissa values of the linear interpolation"
+    },
+    "c_A_Y": {
+      "desc": "Ordinate values of the linear interpolation"
+    },
+    "model0_c_1": {
+      "desc": "Prefactor of the second input variable"
+    }
+  },
+  "buffers": {
+    "model0_c0": {
+      "desc": "Constant offset"
+    }
+  }
+}
diff --git a/tests/unit/models/ComposedModel5_model.gz b/tests/unit/models/ComposedModel5_model.gz
new file mode 100644
index 0000000000..4d0ad23166
Binary files /dev/null and b/tests/unit/models/ComposedModel5_model.gz differ
diff --git a/tests/unit/models/test_BundledModel.cxx b/tests/unit/models/test_BundledModel.cxx
new file mode 100644
index 0000000000..f4aa928bb3
--- /dev/null
+++ b/tests/unit/models/test_BundledModel.cxx
@@ -0,0 +1,46 @@
+// Copyright 2024, UChicago Argonne, LLC
+// All Rights Reserved
+// Software Name: NEML2 -- the New Engineering material Model Library, version 2
+// By: Argonne National Laboratory
+// OPEN SOURCE LICENSE (MIT)
+//
+// Permission is hereby granted, free of charge, to any person obtaining a copy
+// of this software and associated documentation files (the "Software"), to deal
+// in the Software without restriction, including without limitation the rights
+// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the Software is
+// furnished to do so, subject to the following conditions:
+//
+// The above copyright notice and this permission notice shall be included in
+// all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+// THE SOFTWARE.
+
+#include <catch2/catch_test_macros.hpp>
+#include <filesystem>
+
+#include "neml2/models/BundledModel.h"
+
+using namespace neml2;
+
+TEST_CASE("BundledModel", "[models]")
+{
+#if !defined(NEML2_HAS_ZLIB)
+  SKIP("NEML2 was not compiled with zlib support, cannot package models.");
+#elif !defined(NEML2_HAS_JSON)
+  SKIP("NEML2 was not compiled with json support, cannot package models.");
+#else
+  std::ifstream f("models/ComposedModel5.json");
+  auto config = nlohmann::json::parse(f);
+  bundle_model("models/ComposedModel5.i", "model", "", config, "models/ComposedModel5_model.gz");
+  // assert file exists
+  std::filesystem::path output_path = "models/ComposedModel5_model.gz";
+  REQUIRE(std::filesystem::exists(output_path));
+#endif
+}