diff --git a/timflow/steady/aquifer.py b/timflow/steady/aquifer.py
index 5f884ca1..8f5a8781 100644
--- a/timflow/steady/aquifer.py
+++ b/timflow/steady/aquifer.py
@@ -12,12 +12,13 @@
 import numpy as np
 import pandas as pd
 
+from timflow.steady.base_io import BaseIO
 from timflow.steady.constant import ConstantStar
 
 __all__ = ["Aquifer", "SimpleAquifer"]
 
 
-class AquiferData:
+class AquiferData(BaseIO):
     def __init__(self, model, kaq, c, z, npor, ltype, model3d=False):
         """Initialize aquifer data.
 
diff --git a/timflow/steady/base_io.py b/timflow/steady/base_io.py
new file mode 100644
index 00000000..644f0717
--- /dev/null
+++ b/timflow/steady/base_io.py
@@ -0,0 +1,136 @@
+from __future__ import annotations
+
+import inspect
+from functools import wraps
+from importlib import import_module
+from typing import TYPE_CHECKING, TypeVar
+
+from numpy import array, ndarray
+
+if TYPE_CHECKING:
+    from timflow.steady import Model
+
+T = TypeVar("T")
+
+
+def store_input(cls: type[T]) -> type[T]:
+
+    original_init = cls.__init__
+
+    @wraps(original_init)
+    def new_init(self, *args, **kwargs) -> None:
+        original_init(self, *args, **kwargs)
+
+        model_instance: Model | None
+        if "Model" in self.__class__.__name__:
+            model_instance = self
+        else:
+            if args != ():
+                model_instance = args[0]
+            else:
+                model_instance = kwargs.pop("model", None)  # remove model ref
+                if model_instance is None:
+                    model_instance = kwargs.pop("ml", None)  # remove model ref
+        if model_instance is not None:
+            # Prevent the reference to the model object from being stored
+            # this is unused and might complicate pickling.
+            if len(args) != 0:
+                args = args[1:] # model ref always first posarg
+                
+            model_instance._obj_registry.append(
+                {
+                    "class": f"{cls.__module__}.{cls.__qualname__}",
+                    "args": args,
+                    "kwargs": kwargs,
+                }
+            )
+
+    cls.__init__ = new_init
+
+    return cls
+
+
+class BaseIO:
+    @classmethod
+    def to_dict(cls, args: tuple, kwargs: dict):
+        """
+        Collect the constructor arguments into a dict.
+
+        :return: Dict with the arguments.
+        """
+        sig = inspect.signature(cls.__init__)
+        if "Model" not in cls.__name__:
+            if "model" not in kwargs or "ml" not in kwargs:
+                args = args + ("model dummy",)  # add dummy for sig.bind
+        bound = sig.bind(cls, *args, **kwargs)
+        # Reference to class for recreation
+        data = {"_type": f"{cls.__module__}.{cls.__qualname__}"}
+        data.update(
+            {
+                k: cls._serialize(v)
+                for k, v in bound.arguments.items()
+                if k not in ("model", "ml", "self")
+            }
+        )
+        return data
+
+    @classmethod
+    def _serialize(cls, value):
+        """Convert python objects to exportable types.
+
+        :param value: Object for export.
+        :return: Object in exportable form.
+        """
+        if isinstance(value, list):
+            return [cls._serialize(v) for v in value]
+        if isinstance(value, dict):
+            return {k: cls._serialize(v) for k, v in value.items()}
+        if isinstance(value, tuple):
+            return {"tuple": [cls._serialize(v) for v in value]}
+        if isinstance(value, ndarray):
+            return {"ndarray": value.tolist()}
+        return value
+
+    @classmethod
+    def from_dict(cls, data: dict):
+        """Factory method to create an instance of this (sub)class.
+
+        :param data: Dict with parameters
+        :return: Instance of this (sub)class.
+        """
+        type_name: str = data["_type"]
+        module_name = ".".join(type_name.split(".")[:-1])
+        class_name = type_name.split(".")[-1]
+        module = import_module(module_name)
+        subclass = getattr(module, class_name)
+        sig = inspect.signature(subclass.__init__)
+        constructor_args = {}
+
+        for name in sig.parameters:
+            if name in ("model", "ml"):
+                constructor_args[name] = cls._setup_model
+            if name != "self" and name in data:
+                constructor_args[name] = cls._deserialize(data.pop(name))
+        obj = subclass(**constructor_args)
+        if cls._setup_model is None:
+            cls._setup_model = obj
+        return obj
+
+    @classmethod
+    def _deserialize(cls, value):
+        """Convert a dict of values to the right python objects.
+
+        :param value: Imported object
+        :return: Object as correct python-type.
+        """
+        if isinstance(value, dict) and "_type" in value:
+            return cls.from_dict(value)
+        if isinstance(value, dict) and "ndarray" in value:
+            return array(value["ndarray"])
+        if isinstance(value, dict) and "tuple" in value:
+            return tuple(cls._deserialize(v) for v in value["tuple"])
+        if isinstance(value, list):
+            return [cls._deserialize(v) for v in value]
+        if isinstance(value, dict):
+            return {k: cls._deserialize(v) for k, v in value.items()}
+        return value
diff --git a/timflow/steady/circareasink.py b/timflow/steady/circareasink.py
index 4b92e85e..07679cec 100644
--- a/timflow/steady/circareasink.py
+++ b/timflow/steady/circareasink.py
@@ -10,11 +10,13 @@
 import numpy as np
 from scipy.special import i0, i1, k0, k1
 
+from timflow.steady.base_io import store_input
 from timflow.steady.element import Element
 
 __all__ = ["CircAreaSink"]
 
 
+@store_input
 class CircAreaSink(Element):
     """Class to create a circular area-sink.
 
diff --git a/timflow/steady/constant.py b/timflow/steady/constant.py
index dde0baba..b7d32e21 100644
--- a/timflow/steady/constant.py
+++ b/timflow/steady/constant.py
@@ -10,6 +10,7 @@
 
 import numpy as np
 
+from timflow.steady.base_io import store_input
 from timflow.steady.element import Element
 from timflow.steady.equation import PotentialEquation
 
@@ -33,6 +34,7 @@ def __init__(
         )
         # Defined here and not in Element as other elements can have multiple parameters
         # per layers:
+        self.layer = layer
         self.nparam = 1
         self.nunknowns = 0
         self.xr = xr
@@ -71,6 +73,7 @@ def disvecinf(self, x, y, aq=None):
         return rv
 
 
+@store_input
 class Constant(ConstantBase, PotentialEquation):
     """Specify the head at one point in the model in one layer.
 
@@ -178,6 +181,7 @@ def setparams(self, sol):
 
 # class ConstantStar(Element, PotentialEquation):
 # I don't think we need the equation
+# @store_input
 class ConstantStar(Element):
     """Constant representing the particular solution inside a semi-confined aquifer.
 
diff --git a/timflow/steady/element.py b/timflow/steady/element.py
index 8095f98e..09e58f4e 100644
--- a/timflow/steady/element.py
+++ b/timflow/steady/element.py
@@ -11,10 +11,12 @@ def initialize(self):
 
 import numpy as np
 
+from timflow.steady.base_io import BaseIO
+
 __all__ = ["Element"]
 
 
-class Element:
+class Element(BaseIO):
     """Base class for all timflow.steady elements.
 
     Elements represent physical features in the aquifer system such as wells,
diff --git a/timflow/steady/inhomogeneity.py b/timflow/steady/inhomogeneity.py
index 9f94158e..e1e1fe46 100644
--- a/timflow/steady/inhomogeneity.py
+++ b/timflow/steady/inhomogeneity.py
@@ -19,6 +19,7 @@
 
 from timflow.steady.aquifer import AquiferData
 from timflow.steady.aquifer_parameters import param_3d, param_maq
+from timflow.steady.base_io import store_input
 from timflow.steady.constant import ConstantInside, ConstantStar
 from timflow.steady.element import Element
 from timflow.steady.intlinesink import (
@@ -145,6 +146,7 @@ def create_elements(self):
             c.inhomelement = True
 
 
+@store_input
 class PolygonInhomMaq(PolygonInhom):
     """Create a polygonal inhomogeneity.
 
@@ -240,6 +242,7 @@ def __init__(
         )
 
 
+@store_input
 class PolygonInhom3D(PolygonInhom):
     """Create a multi-layer model object consisting of many aquifer layers.
 
@@ -545,6 +548,7 @@ def create_elements(self):
             c.inhomelement = True
 
 
+@store_input
 class BuildingPitMaq(BuildingPit):
     """Element to simulate a building pit with an impermeable wall in ModelMaq.
 
@@ -627,6 +631,7 @@ def __init__(
         )
 
 
+@store_input
 class BuildingPit3D(BuildingPit):
     """Element to simulate a building pit with an impermeable wall in Model3D.
 
@@ -917,6 +922,7 @@ def create_elements(self):
             c.inhomelement = True
 
 
+@store_input
 class LeakyBuildingPitMaq(LeakyBuildingPit):
     """Element to simulate a building pit with a leaky wall in ModelMaq.
 
@@ -1005,6 +1011,7 @@ def __init__(
         )
 
 
+@store_input
 class LeakyBuildingPit3D(LeakyBuildingPit):
     """Element to simulate a building pit with a leaky wall in Model3D.
 
diff --git a/timflow/steady/inhomogeneity1d.py b/timflow/steady/inhomogeneity1d.py
index 6baa2e0c..58e96b15 100644
--- a/timflow/steady/inhomogeneity1d.py
+++ b/timflow/steady/inhomogeneity1d.py
@@ -17,6 +17,7 @@
 
 from timflow.steady.aquifer import AquiferData
 from timflow.steady.aquifer_parameters import param_3d, param_maq
+from timflow.steady.base_io import store_input
 from timflow.steady.constant import ConstantStar
 from timflow.steady.linesink1d import FluxDiffLineSink1D, HeadDiffLineSink1D
 from timflow.steady.stripareasink import XsectionAreaSinkInhom
@@ -326,6 +327,7 @@ def plot(
         return ax
 
 
+@store_input
 class XsectionMaq(Xsection):
     """Cross-section inhomogeneity for a multi-aquifer sequence.
 
@@ -379,6 +381,7 @@ def __init__(
         N=None,
         name=None,
     ):
+        self.topboundary = topboundary
         if c is None:
             c = []
         if z is None:
@@ -395,6 +398,7 @@ def __init__(
         )
 
 
+@store_input
 class Xsection3D(Xsection):
     """Cross-section inhomogeneity consisting of stacked aquifer layers.
 
@@ -459,6 +463,7 @@ def __init__(
         N=None,
         name=None,
     ):
+        self.topboundary = topboundary
         if z is None:
             z = [1, 0]
         (
@@ -487,6 +492,7 @@ def __init__(self, model, x1, x2, kaq, c, z, npor, ltype, hstar, N, name=None):
         super().__init__(model, x1, x2, kaq, c, z, npor, ltype, hstar, N, name=name)
 
 
+@store_input
 class StripInhomMaq(XsectionMaq):
     def __init__(
         self,
@@ -510,6 +516,7 @@ def __init__(
         super().__init__(model, x1, x2, kaq, z, c, npor, topboundary, hstar, N, name)
 
 
+@store_input
 class StripInhom3D(Xsection3D):
     def __init__(
         self,
diff --git a/timflow/steady/linedoublet.py b/timflow/steady/linedoublet.py
index 29ad827b..274eeafb 100644
--- a/timflow/steady/linedoublet.py
+++ b/timflow/steady/linedoublet.py
@@ -13,6 +13,7 @@
 import numpy as np
 
 from timflow.bessel.besselnumba import disbesldv, potbesldv
+from timflow.steady.base_io import store_input
 from timflow.steady.controlpoints import controlpoints
 from timflow.steady.element import Element
 from timflow.steady.equation import DisvecEquation, LeakyWallEquation
@@ -185,6 +186,7 @@ def plot(self, ax=None, layer=None):
             ax.plot([self.x1, self.x2], [self.y1, self.y2], "k")
 
 
+@store_input
 class ImpermeableWall(LineDoubletHoBase, DisvecEquation):
     """Create a segment of an impermeable wall, which is simulated with a line-doublet.
 
@@ -251,6 +253,7 @@ def setparams(self, sol):
         self.parameters[:, 0] = sol
 
 
+@store_input
 class LeakyWall(LineDoubletHoBase, LeakyWallEquation):
     """Create a segment of a leaky wall, which is simulated with a line-doublet.
 
@@ -423,6 +426,7 @@ def plot(self, ax=None, layer=None):
             ax.plot(self.x, self.y, "k")
 
 
+@store_input
 class ImpermeableWallString(LineDoubletStringBase, DisvecEquation):
     """Create a string of impermeable wall segments consisting of line-doublets.
 
@@ -473,6 +477,7 @@ def setparams(self, sol):
         self.parameters[:, 0] = sol
 
 
+@store_input
 class LeakyWallString(LineDoubletStringBase, LeakyWallEquation):
     """Create a string of leaky wall segments consisting of line-doublets.
 
@@ -525,6 +530,7 @@ def setparams(self, sol):
         self.parameters[:, 0] = sol
 
 
+@store_input
 class ImpLineDoublet(ImpermeableWall):
     """Deprecated alias for :class:`.ImpermeableWall`.
 
@@ -542,6 +548,7 @@ def __init__(self, *args, **kwargs):
         super().__init__(*args, **kwargs)
 
 
+@store_input
 class ImpLineDoubletString(ImpermeableWallString):
     """Deprecated alias for :class:`.ImpermeableWallString`.
 
@@ -559,6 +566,7 @@ def __init__(self, *args, **kwargs):
         super().__init__(*args, **kwargs)
 
 
+@store_input
 class LeakyLineDoublet(LeakyWall):
     """Deprecated alias for :class:`.LeakyWall`.
 
@@ -576,6 +584,7 @@ def __init__(self, *args, **kwargs):
         super().__init__(*args, **kwargs)
 
 
+@store_input
 class LeakyLineDoubletString(LeakyWallString):
     """Deprecated alias for :class:`.LeakyWallString`.
 
diff --git a/timflow/steady/linedoublet1d.py b/timflow/steady/linedoublet1d.py
index 3fcd5e62..b001e114 100644
--- a/timflow/steady/linedoublet1d.py
+++ b/timflow/steady/linedoublet1d.py
@@ -13,6 +13,7 @@
 import matplotlib.pyplot as plt
 import numpy as np
 
+from timflow.steady.base_io import store_input
 from timflow.steady.element import Element
 from timflow.steady.equation import DisvecEquation, LeakyWallEquation
 
@@ -122,6 +123,7 @@ def plot(self, ax=None):
         )
 
 
+@store_input
 class ImpermeableWall1D(LineDoublet1D, DisvecEquation):
     """Create 1D impermeable wall."""
 
@@ -147,6 +149,7 @@ def setparams(self, sol):
         self.parameters[:, 0] = sol
 
 
+@store_input
 class LeakyWall1D(LineDoublet1D, LeakyWallEquation):
     """Create an infinitely long leaky or impermeable wall.
 
@@ -196,6 +199,7 @@ def setparams(self, sol):
         self.parameters[:, 0] = sol
 
 
+@store_input
 class ImpLineDoublet1D(ImpermeableWall1D):
     """Deprecated alias for :class:`.ImpermeableWall1D`.
 
@@ -213,6 +217,7 @@ def __init__(self, *args, **kwargs):
         super().__init__(*args, **kwargs)
 
 
+@store_input
 class LeakyLineDoublet1D(LeakyWall1D):
     """Deprecated alias for :class:`.LeakyWall1D`.
 
diff --git a/timflow/steady/linesink.py b/timflow/steady/linesink.py
index 145accd9..b5073f83 100644
--- a/timflow/steady/linesink.py
+++ b/timflow/steady/linesink.py
@@ -13,6 +13,7 @@
 import numpy as np
 
 from timflow.bessel.besselnumba import disbeslsv, potbeslsv
+from timflow.steady.base_io import store_input
 from timflow.steady.controlpoints import controlpoints, strengthinf_controlpoints
 from timflow.steady.element import Element
 from timflow.steady.equation import HeadEquation
@@ -500,6 +501,7 @@ def headinside(self, icp=0):
         return hinside
 
 
+@store_input
 class River(LineSinkHoBase, HeadEquation):
     """Head-specified line-sink which may optionally have a width and resistance.
 
@@ -609,6 +611,7 @@ def setparams(self, sol):
         self.parameters[:, 0] = sol
 
 
+@store_input
 class Ditch(River):
     """Line-sink with specified total discharge, and uniform but unknown head.
 
@@ -872,6 +875,7 @@ def plot(self, ax=None, layer=None, **kwargs):
             ls.plot(layer=layer, ax=ax, **kwargs)
 
 
+@store_input
 class RiverString(LineSinkStringBase2):
     """String of head-specified line-sinks with optional width and resistance.
 
@@ -1029,6 +1033,7 @@ def equation(self):
         return mat, rhs
 
 
+@store_input
 class DitchString(RiverString):
     """String of Ditches with specified discharge and uniform unknown head.
 
@@ -1117,6 +1122,7 @@ def equation(self):
         return mat, rhs
 
 
+@store_input
 class LineSinkDitch(Ditch):
     """Deprecated alias for :class:`.Ditch`.
 
@@ -1133,6 +1139,7 @@ def __init__(self, *args, **kwargs):
         super().__init__(*args, **kwargs)
 
 
+@store_input
 class LineSinkDitchString(DitchString):
     """Deprecated alias for :class:`.DitchString`.
 
@@ -1150,6 +1157,7 @@ def __init__(self, *args, **kwargs):
         super().__init__(*args, **kwargs)
 
 
+@store_input
 class HeadLineSink(River):
     """Deprecated alias for :class:`.River`.
 
@@ -1167,6 +1175,7 @@ def __init__(self, *args, **kwargs):
         super().__init__(*args, **kwargs)
 
 
+@store_input
 class HeadLineSinkString(RiverString):
     """Deprecated alias for :class:`.RiverString`.
 
@@ -1184,6 +1193,7 @@ def __init__(self, *args, **kwargs):
         super().__init__(*args, **kwargs)
 
 
+@store_input
 class CollectorWell(DitchString):
     """Collector well: collection of line sinks with a specified total discharge.
 
@@ -1246,6 +1256,7 @@ def __init__(
         self.name = "CollectorWell"
 
 
+@store_input
 class RadialCollectorWell(CollectorWell):
     """Radial collector well.
 
diff --git a/timflow/steady/linesink1d.py b/timflow/steady/linesink1d.py
index bfff0c15..6101f524 100644
--- a/timflow/steady/linesink1d.py
+++ b/timflow/steady/linesink1d.py
@@ -12,6 +12,7 @@
 import matplotlib.pyplot as plt
 import numpy as np
 
+from timflow.steady.base_io import store_input
 from timflow.steady.element import Element
 from timflow.steady.equation import (
     DisvecDiffEquation,
@@ -136,6 +137,7 @@ def plot(self, ax=None):
             )
 
 
+@store_input
 class LineSink1D(LineSink1DBase, MscreenWellEquation):
     """Create an infinitely long line-sink with a given discharge per unit length.
 
@@ -195,7 +197,7 @@ def initialize(self):
     def setparams(self, sol):
         self.parameters[:, 0] = sol
 
-
+@store_input
 class River1D(LineSink1DBase, HeadEquation):
     """Create an infinitely long line-sink with a given head.
 
@@ -250,6 +252,7 @@ def setparams(self, sol):
         self.parameters[:, 0] = sol
 
 
+@store_input
 class HeadLineSink1D(River1D):
     """Deprecated alias for :class:`.River1D`.
 
diff --git a/timflow/steady/model.py b/timflow/steady/model.py
index 30caf3f4..a0342be3 100644
--- a/timflow/steady/model.py
+++ b/timflow/steady/model.py
@@ -10,8 +10,11 @@
     ml.solve()
 """
 
+import json
 import multiprocessing as mp
 import warnings
+from importlib import import_module
+from typing import Any
 
 import numpy as np
 import pandas as pd
@@ -19,6 +22,7 @@
 
 from timflow.steady.aquifer import Aquifer, SimpleAquifer
 from timflow.steady.aquifer_parameters import param_3d, param_maq
+from timflow.steady.base_io import BaseIO, store_input
 from timflow.steady.constant import ConstantStar
 from timflow.steady.plots import PlotSteady
 from timflow.version import check_tqdm_parallel
@@ -42,7 +46,7 @@ def _compute_velocity_mp(args):
     return i, vv
 
 
-class Model:
+class Model(BaseIO):
     """Create a model consisting of an arbitrary sequence of aquifers and leaky layers.
 
     Notes
@@ -80,8 +84,51 @@ def __init__(self, kaq, z, c, npor, ltype, model3d=False):
 
         self.plots = PlotSteady(self)
 
+        self._obj_registry: list[dict[str, Any]] = []
+
         self.initialized = False
 
+    def to_json(self, filepath) -> None:
+        """
+        Write the constructor arguments to a JSON-file.
+
+        :param filepath: Filepath for the to be created JSON-file.
+        """
+        data = {}
+        i = 0
+        for item in self._obj_registry:
+            type_name:str = item["class"]
+            args: list[Any] = item.get("args", [])
+            kwargs: dict[str, Any] = item.get("kwargs", {})
+            module_name: str = ".".join(type_name.split(".")[:-1])
+            class_name: str = type_name.split(".")[-1]
+            module = import_module(module_name)
+            subclass = getattr(module, class_name)
+            data.update({f"object{i}": subclass.to_dict(args, kwargs)})
+            i += 1
+        
+        with open(filepath, "w") as f:
+            f.write(json.dumps(data, indent=4))
+
+    @classmethod
+    def from_json(cls, filepath):
+        """
+        Read the constructor arguments and potential addition attributes from a JSON-file.
+
+        :param filepath: Filepath to the to be created JSON-file.
+        """
+        cls._setup_model = None
+        with open(filepath, "r") as f:
+            data: dict = json.load(f)
+        for k, v in data.items():
+            if k == "object0":  # Model object is always first created.
+                obj = cls.from_dict(v)
+                continue
+            if "obj" not in locals():  # No model in json
+                raise ImportError("No main model found in the JSON-file.")
+            cls.from_dict(v)
+        return obj
+
     def initialize(self):
         # remove inhomogeneity elements (they are added again)
         self.elementlist = [e for e in self.elementlist if not e.inhomelement]
@@ -943,6 +990,7 @@ def vcontoursf1D(self, *args, **kwargs):
         return self.plots.vcontour_stream_function(*args, **kwargs)
 
 
+@store_input
 class ModelMaq(Model):
     """Create a model by specifying a multi-aquifer sequence of aquifer-leaky layer.
 
@@ -982,6 +1030,7 @@ class ModelMaq(Model):
     """
 
     def __init__(self, kaq=1, z=None, c=None, npor=0.3, topboundary="conf", hstar=None):
+        self.topboundary = topboundary
         if c is None:
             c = []
         if z is None:
@@ -993,6 +1042,7 @@ def __init__(self, kaq=1, z=None, c=None, npor=0.3, topboundary="conf", hstar=No
             ConstantStar(self, hstar, aq=self.aq)
 
 
+@store_input
 class Model3D(Model):
     """Create a multi-layer model object consisting of stacked aquifer layers.
 
@@ -1077,6 +1127,7 @@ def __init__(
             ConstantStar(self, hstar, aq=self.aq)
 
 
+@store_input
 class ModelXsection(Model):
     r"""Model for cross-section (2D vertical slice) problems.
 
@@ -1097,6 +1148,7 @@ class ModelXsection(Model):
     """
 
     def __init__(self, naq=1):
+        self.naq = naq
         self.elementlist = []
         self.elementdict = {}  # only elements that have a label
         self.aq = SimpleAquifer(self, naq)
diff --git a/timflow/steady/stripareasink.py b/timflow/steady/stripareasink.py
index 0e3d18c5..645d6893 100644
--- a/timflow/steady/stripareasink.py
+++ b/timflow/steady/stripareasink.py
@@ -12,6 +12,7 @@
 import matplotlib.pyplot as plt
 import numpy as np
 
+from timflow.steady.base_io import store_input
 from timflow.steady.element import Element
 
 __all__ = ["XsectionAreaSinkInhom", "XsectionAreaSink"]
@@ -126,6 +127,7 @@ def plot(self, ax=None, n_arrows=10, **kwargs):
         return ax
 
 
+@store_input
 class XsectionAreaSink(Element):
     """Cross-section area-sink for testing purposes only.
 
diff --git a/timflow/steady/uflow.py b/timflow/steady/uflow.py
index 1242426e..b42e5e2a 100644
--- a/timflow/steady/uflow.py
+++ b/timflow/steady/uflow.py
@@ -9,11 +9,13 @@
 
 import numpy as np
 
+from timflow.steady.base_io import store_input
 from timflow.steady.element import Element
 
 __all__ = ["Uflow"]
 
 
+@store_input
 class Uflow(Element):
     """Add uniform flow to the model.
 
diff --git a/timflow/steady/well.py b/timflow/steady/well.py
index a5cd27df..dbf733d4 100644
--- a/timflow/steady/well.py
+++ b/timflow/steady/well.py
@@ -14,6 +14,7 @@
 import numpy as np
 from scipy.special import k0, k1
 
+from timflow.steady.base_io import store_input
 from timflow.steady.element import Element
 from timflow.steady.equation import HeadEquation, MscreenWellNoflowEquation
 from timflow.steady.trace import tracelines
@@ -341,6 +342,7 @@ def plotcapzone(self, *args, **kwargs):
         )
 
 
+@store_input
 class Well(WellBase):
     r"""Well Class to create a well with a specified discharge.
 
@@ -451,6 +453,7 @@ def setparams(self, sol):
         self.parameters[:, 0] = sol
 
 
+@store_input
 class HeadWell(WellBase, HeadEquation):
     r"""HeadWell Class to create a well with a specified head inside the well.
 
@@ -523,6 +526,7 @@ def setparams(self, sol):
         self.parameters[:, 0] = sol
 
 
+@store_input
 class TargetHeadWell(WellBase):
     r"""TargetHeadWell is a well with a specified head at (layer, x, y).
 
@@ -924,6 +928,7 @@ def headinside(self):
         return self.wlist[0].headinside()[0]
 
 
+@store_input
 class WellString(WellStringBase):
     """
     WellString is a string of wells for which the total discharge is specified.
@@ -1012,6 +1017,7 @@ def setparams(self, sol):
             i += w.nparam
 
 
+@store_input
 class HeadWellString(WellStringBase):
     """
     HeadWellString is a string of wells for which the head is specified in the wells.
@@ -1080,6 +1086,7 @@ def setparams(self, sol):
             i += w.nparam
 
 
+@store_input
 class TargetHeadWellString(WellStringBase):
     """
     A string of wells for which the head is specified at a point.