diff --git a/climatecritters/core/forcing.py b/climatecritters/core/forcing.py
index ba98848..c51e856 100644
--- a/climatecritters/core/forcing.py
+++ b/climatecritters/core/forcing.py
@@ -25,11 +25,11 @@
     f   = seq.compile()              # → Forcing
 
     # Register with a model
-    model.register_forcing('S', f, attachment_style='additive', timing='pre')
+    cc.Model.register_forcing('S', f, attachment_style='additive', timing='pre')
 
     # Superpose two indefinite signals
     combined = cc.Forcing(orbital_func) + cc.Forcing(noise_func)
-    model.register_forcing('S0', combined)
+    cc.Model.register_forcing('S0', combined)
 """
 
 from __future__ import annotations
diff --git a/climatecritters/core/model.py b/climatecritters/core/model.py
index 9e7a776..1e97a61 100644
--- a/climatecritters/core/model.py
+++ b/climatecritters/core/model.py
@@ -76,7 +76,7 @@ def _format_value(v):
 class Model:
     """The overarching model structure for ClimateCritters.
 
-    CCModel serves as the archetype/parent class for models within the
+    Model serves as the archetype/parent class for models within the
     ``model_critters`` directory.  It is not meant to be instantiated directly.
 
     Parameter handling
@@ -580,7 +580,7 @@ def sde_noise(self, t, y):
 
     def integrate(self, t_span=None, y0=None, method='RK45', dt=None,
                   output_time=None, run_name=None, kwargs=None):
-        """Integrate the model over a time span and return a :class:`CCOutput`.
+        """Integrate the model over a time span and return an `Output`.
 
         Parameters
         ----------
@@ -608,7 +608,7 @@ def integrate(self, t_span=None, y0=None, method='RK45', dt=None,
             Fixed timestep for ``euler``, ``euler_maruyama``, ``heun_maruyama``,
             ``milstein``, and ``rk4``.  Required for those methods.
         output_time : array-like, optional
-            If provided, the returned ``CCOutput`` is immediately reframed onto
+            If provided, the returned ``Output`` is immediately reframed onto
             this time axis (e.g. to exclude a spin-up period).
             ``output.model_time`` always retains the raw solver grid.
         run_name : str, optional
diff --git a/climatecritters/core/output.py b/climatecritters/core/output.py
index e7c5357..84ea07e 100644
--- a/climatecritters/core/output.py
+++ b/climatecritters/core/output.py
@@ -1,4 +1,4 @@
-"""CCOutput — container for a single CCModel integration run.
+"""Output — container for a single Model integration run.
 
 Separating run output from model configuration lets the same model instance
 be re-run with different parameters or initial conditions while keeping each
@@ -13,9 +13,9 @@
 
 
 class Output:
-    """Container for the results of one call to ``CCModel.integrate()``.
+    """Container for the results of one call to ``Model.integrate()``.
 
-    ``CCOutput`` carries the full trajectory produced by the solver and
+    ``Output`` carries the full trajectory produced by the solver and
     exposes output-focused operations (noise addition, pyleoclim export,
     time resampling).  Keeping these on the output rather than on the model
     means a single model instance can produce multiple independent outputs
diff --git a/climatecritters/model_critters/box_model.py b/climatecritters/model_critters/box_model.py
index a6f59bb..1d3dd91 100644
--- a/climatecritters/model_critters/box_model.py
+++ b/climatecritters/model_critters/box_model.py
@@ -2,13 +2,13 @@
 
 """Generic box-model utilities for ClimateCritters.
 
-This module provides a small declarative layer on top of :class:`CCModel` for
+This module provides a small declarative layer on top of `Model` for
 building simple ODE box models without writing a full bespoke subclass each
 time.
 
 Two styles are supported:
 
-1. Explicit callable tendencies via :class:`BoxModelSpec.register_tendency`
+1. Explicit callable tendencies via `BoxModelSpec.register_tendency`
 2. Automatic box-network assembly via reciprocal exchange and directed
    transport relations
 
@@ -87,7 +87,7 @@ class BoxModelContext:
     The context exposes a compact, box-model-friendly API:
 
     - ``ctx["A"]`` for direct state lookup
-    - ``ctx.param("k")`` for parameters resolved through ``CCModel.get_param``
+    - ``ctx.param("k")`` for parameters resolved through ``Model.get_param``
     - ``ctx.input("R")`` for external prescribed inputs
     - ``ctx.volume("SP")`` and ``ctx.concentration("SP")`` for box-network
       models
@@ -128,7 +128,7 @@ class BoxModelSpec:
     inputs, diagnostics, and either:
 
     - explicit callable tendencies registered via
-      :meth:`register_tendency` / :meth:`register_relations`, or
+      `BoxModelSpec.register_tendency` / `BoxModelSpec.register_relations`, or
     - an automatic box network assembled from volumes, reciprocal exchange,
       and directed transport terms
 
@@ -140,13 +140,13 @@ class BoxModelSpec:
 
     Notes
     -----
-    Call :meth:`validate` (or :meth:`make_boxmodel`) before integrating.
+    Call `BoxModelSpec.validate` (or `make_boxmodel`) before integrating.
     Validation checks that every state variable has a tendency relation
     (explicit mode) or a registered volume (automatic mode).
 
     See also
     --------
-    GenericBoxModel : The ``CCModel`` subclass produced by :meth:`make_boxmodel`.
+    GenericBoxModel : The ``Model`` subclass produced by `make_boxmodel`.
     BoxModelContext : Evaluation context passed to callable tendencies.
 
     Examples
@@ -218,7 +218,7 @@ def register_parameters(self, **parameters):
         """Register default parameter values.
 
         Values may be constants, callables, or ``Forcing``-like objects
-        compatible with ``CCModel.get_param``.
+        compatible with ``Model.get_param``.
         """
         self.parameter_defaults.update(parameters)
         return self
@@ -280,7 +280,7 @@ def register_input(self, name, fallback_param=None):
     def register_tendency(self, state_variable, relation):
         """Register a callable tendency for one state variable.
 
-        ``relation`` must accept a :class:`BoxModelContext` and return the
+        ``relation`` must accept a `BoxModelContext` and return the
         tendency for ``state_variable``.
         """
         self.tendency_relations[str(state_variable)] = relation
@@ -376,7 +376,7 @@ def validate(self):
             raise ValueError(f"Missing tendency relations for state variables: {missing}")
 
     def make_model(self, var_name=None, **parameter_overrides):
-        """Instantiate a :class:`GenericBoxModel` from this spec."""
+        """Instantiate a `GenericBoxModel` from this spec."""
         self.validate()
         return GenericBoxModel(
             self,
@@ -389,9 +389,9 @@ def make_boxmodel(self, var_name=None, **parameter_overrides):
 
 
 class GenericBoxModel(Model):
-    """``CCModel`` subclass produced by :class:`BoxModelSpec`.
+    """``Model`` subclass produced by `BoxModelSpec`.
 
-    Users construct this via :meth:`BoxModelSpec.make_boxmodel` rather than
+    Users construct this via `make_boxmodel` rather than
     instantiating it directly.  The resulting model integrates with
     ``model.integrate(...)``, exports to Pyleoclim via ``output.to_pyleo()``,
     and computes diagnostics from solved history via
@@ -427,7 +427,7 @@ def __init__(self, spec, var_name=None, *args, **kwargs):
             setattr(self, name, value)
 
     def uses_post_history(self):
-        """Route diagnostics through CCModel's post-history hooks."""
+        """Route diagnostics through Model's post-history hooks."""
         return True
 
     def resolve_input(self, name, t, state):
diff --git a/climatecritters/model_critters/damped_spring.py b/climatecritters/model_critters/damped_spring.py
index a2a82db..1ea0b63 100644
--- a/climatecritters/model_critters/damped_spring.py
+++ b/climatecritters/model_critters/damped_spring.py
@@ -111,7 +111,7 @@ def __init__(
         self.params = ()
 
     # ------------------------------------------------------------------
-    # CCModel interface
+    # Model interface
     # ------------------------------------------------------------------
 
     def uses_post_history(self):
diff --git a/climatecritters/model_critters/ebm.py b/climatecritters/model_critters/ebm.py
index 7ad1076..352646f 100644
--- a/climatecritters/model_critters/ebm.py
+++ b/climatecritters/model_critters/ebm.py
@@ -11,7 +11,7 @@
 
 # ---------------------------------------------------------------------------
 # Module-level physics helpers
-# All callables comply with the CCModel contract: (t), (t, state), or
+# All callables comply with the Model contract: (t), (t, state), or
 # (t, state, model) with the first positional arg named 't' or 'time'.
 # ---------------------------------------------------------------------------
 
@@ -49,7 +49,7 @@ def OLR_func(pRad=650, ps=1000):
     Returns
     -------
     func : callable
-        A ``(t, state)`` callable compliant with the CCModel parameter
+        A ``(t, state)`` callable compliant with the Model parameter
         contract (see ``contracts/signal_model_contract.md``).  Returns OLR
         in W m⁻².
     """
@@ -93,7 +93,7 @@ def albedo_func1D(t, state, model, *, a2=0.25, alpha_ice=0.6, alpha_0=0.1, T1=26
     """P2-corrected latitudinal albedo using a Legendre polynomial parameterization.
 
     Uses the global-mean temperature to set the base albedo via the same
-    quadratic ice-line transition as :func:`albedo_func`, then adds a
+    quadratic ice-line transition as `albedo_func`, then adds a
     second-order Legendre polynomial correction to capture the equator-to-pole
     gradient.  Requires the model to expose a ``phi`` attribute (degrees).
 
@@ -150,7 +150,7 @@ class EBMBase(Model):
 
     See also
     --------
-    climatecritters.core.CCModel : Base class for all ClimateCritters models.
+    climatecritters.core.Model : Base class for all ClimateCritters models.
     EBM0D : Zero-dimensional EBM subclass.
     EBM1DLat : Latitudinally-resolved diffusive EBM subclass.
     """
@@ -477,7 +477,7 @@ def __init__(self, var_name='ebm1d_lat', grid_n=50, C=10.0, D=0.55,
     def validate_initial_state(self, y0):
         """Validate and normalize the initial temperature profile.
 
-        Overrides :meth:`CCModel.validate_initial_state` to accept a scalar
+        Overrides `Model.validate_initial_state` to accept a scalar
         initial condition and broadcast it uniformly across the latitude grid.
 
         Parameters
@@ -503,7 +503,7 @@ def validate_initial_state(self, y0):
     def calc_albedo(self, T, t):
         """Compute the latitudinal ice-albedo with a linear transition zone.
 
-        Overrides :meth:`EBMBase.calc_albedo`.  Each grid point is assigned
+        Overrides `EBMBase.calc_albedo`.  Each grid point is assigned
         an albedo based on local temperature: 0.6 below -10 °C, 0.3 above
         0 °C, and a linear blend in between.
 
@@ -533,7 +533,7 @@ def calc_albedo(self, T, t):
     def calc_OLR(self, T, t):
         """Compute the Budyko linear OLR: ``(A - CO2_forcing) + B * T``.
 
-        Overrides :meth:`EBMBase.calc_OLR`.  Parameters ``A``, ``B``, and
+        Overrides `EBMBase.calc_OLR`.  Parameters ``A``, ``B``, and
         ``CO2_forcing`` are resolved through ``get_param_value``, so they can
         be time-varying or Forcing objects.
 
@@ -699,7 +699,7 @@ def dydt(self, t, state):
 
         This method has **no side effects**: because ``uses_post_history = True``,
         all output is derived from the full solved trajectory in
-        :meth:`populate_diagnostics_from_history` rather than accumulated here.
+        `EBM1DLat.populate_diagnostics_from_history` rather than accumulated here.
 
         Parameters
         ----------
@@ -725,7 +725,7 @@ def dydt(self, t, state):
     def populate_diagnostics_from_history(self, time, history):
         """Compute diagnostic variables from the full solved trajectory.
 
-        Called automatically by :meth:`CCModel.post_integrate` after the
+        Called automatically by `Model.post_integrate` after the
         solver completes.  Populates ``self.diagnostic_variables`` with the
         global-mean temperature and ice-line latitude at every timestep.
 
diff --git a/climatecritters/model_critters/enso_recharge.py b/climatecritters/model_critters/enso_recharge.py
index ea13ac6..3141d6d 100644
--- a/climatecritters/model_critters/enso_recharge.py
+++ b/climatecritters/model_critters/enso_recharge.py
@@ -9,8 +9,8 @@ def seasonal_forcing(A=0.5, period=6.0):
     """Return a sinusoidal seasonal forcing callable for the ENSO recharge oscillator.
 
     The returned function computes ``A * sin(2π t / period)`` and can be
-    wrapped in a :class:`~climatecritters.core.Forcing` or passed directly to
-    :meth:`~climatecritters.core.CCModel.register_forcing`::
+    wrapped in a `Forcing` or passed directly to
+    `Model.register_forcing`::
 
         import climatecritters as cc
         from climatecritters.model_critters.enso_recharge import (
@@ -64,7 +64,7 @@ class ENSORechargeOscillator(Model):
     where ``b = b0*mu`` and ``R = gamma*b - c``.
 
     Seasonal or any other external forcing is added through the standard
-    :meth:`~climatecritters.core.CCModel.register_forcing` interface::
+    `Model.register_forcing` interface::
 
         from climatecritters.utils.forcing import create_sinusoid_forcing
 
diff --git a/climatecritters/model_critters/g24.py b/climatecritters/model_critters/g24.py
index 1ee4af1..ce0eac7 100644
--- a/climatecritters/model_critters/g24.py
+++ b/climatecritters/model_critters/g24.py
@@ -48,7 +48,7 @@ class Model3(Model):
         model.register_forcing('insolation', cc.Forcing(calc_f))
 
     The internal derivative of forcing ``dfdt`` is computed via
-    :func:`calc_df` by default; supply a custom callable or
+    `calc_df` by default; supply a custom callable or
     ``cc.Forcing`` to override.
 
     Parameter defaults are taken from Ganopolski (2024).  Time-varying
diff --git a/climatecritters/model_critters/melcher25.py b/climatecritters/model_critters/melcher25.py
index 292ff9f..0e23567 100644
--- a/climatecritters/model_critters/melcher25.py
+++ b/climatecritters/model_critters/melcher25.py
@@ -110,7 +110,7 @@ class Melcher25(Model):
 
     See also
     --------
-    cc.core.CCModel : Base class.
+    cc.core.Model : Base class.
     classify_bistable_states : Reclassify a Δb signal post-hoc without
         re-running the SDE.
 
@@ -267,7 +267,7 @@ def sde_noise(self, t, y):
     def populate_diagnostics_from_history(self, time, history):
         """Classify states and compute thresholds from the full solved trajectory.
 
-        Called automatically by ``CCModel.post_integrate`` after every
+        Called automatically by ``Model.post_integrate`` after every
         ``integrate()`` call. Populates ``diagnostic_variables['states']``
         with the hysteresis classification and sets
         ``self.stadial_threshold`` / ``self.interstadial_threshold``.
diff --git a/climatecritters/model_critters/pendulum.py b/climatecritters/model_critters/pendulum.py
index c825a9e..052c2fa 100644
--- a/climatecritters/model_critters/pendulum.py
+++ b/climatecritters/model_critters/pendulum.py
@@ -487,7 +487,7 @@ def cartesian_positions(self):
 
 @dataclass
 class PendulumRodBeta:
-    """Properties of one rod in a :class:`MultiPendulumBeta` chain.
+    """Properties of one rod in a `MultiPendulumBeta` chain.
 
     Parameters
     ----------
@@ -504,7 +504,7 @@ class PendulumRodBeta:
     Notes
     -----
     ``forcing`` is resolved at runtime using the same contract as
-    :class:`CCModel` parameters, so constants, ``lambda t: ...`` callables,
+    `Model` parameters, so constants, ``lambda t: ...`` callables,
     ``lambda t, state: ...`` callables, and forcing objects are all accepted.
     """
 
@@ -545,7 +545,7 @@ class MultiPendulumBeta(Model):
     Parameters
     ----------
     rods:
-        Sequence of :class:`PendulumRodBeta` in pivot-to-tip order. Must
+        Sequence of `PendulumRodBeta` in pivot-to-tip order. Must
         contain at least two rods.
     g:
         Gravitational acceleration in m/s^2. Must be > 0.
diff --git a/climatecritters/utils/forcing.py b/climatecritters/utils/forcing.py
index a326ae1..6c6a185 100644
--- a/climatecritters/utils/forcing.py
+++ b/climatecritters/utils/forcing.py
@@ -1,19 +1,19 @@
-"""Convenience factories for common :class:`~climatecritters.core.Forcing` patterns.
+"""Convenience factories for common `Forcing` patterns.
 
-All public functions return either a :class:`~climatecritters.core.Forcing`
-(when no ``duration`` is given) or a :class:`~climatecritters.core.ForcingElement`
-(when ``duration`` is given).  The unified entry point is :func:`create_forcing`;
+All public functions return either a `Forcing`
+(when no ``duration`` is given) or a `ForcingElement`
+(when ``duration`` is given).  The unified entry point is `create_forcing`;
 the named factories are ergonomic aliases that build their callable internally.
 
 .. rubric:: Duration gate
 
 Every factory accepts an optional ``duration`` keyword:
 
-* **No duration** → returns an indefinite :class:`~climatecritters.core.Forcing`
+* **No duration** → returns an indefinite `Forcing`
   backed by a lambda.  Suitable for perpetual signals (orbital forcing, seasonal
   cycle, noise) registered directly with a model.
-* **With duration** → returns a bounded :class:`~climatecritters.core.ForcingElement`
-  that can be composed into a :class:`~climatecritters.core.ForcingSequence` and
+* **With duration** → returns a bounded `ForcingElement`
+  that can be composed into a `ForcingSequence` and
   then compiled::
 
       elem = create_sinusoid_forcing(A=5.0, period=1.0, duration=10.0)
@@ -89,15 +89,15 @@ def create_forcing(func, duration=None):
     func : callable
         Function with signature ``f(t) -> float | ndarray``.
     duration : float, optional
-        If given, returns a bounded :class:`~climatecritters.core.ForcingElement`
+        If given, returns a bounded `ForcingElement`
         lasting ``duration`` time units.  If omitted, returns an indefinite
-        :class:`~climatecritters.core.Forcing`.
+        `Forcing`.
 
     Returns
     -------
     Forcing or ForcingElement
-        * No ``duration`` → :class:`~climatecritters.core.Forcing` (indefinite)
-        * With ``duration`` → :class:`~climatecritters.core.ForcingElement` (bounded)
+        * No ``duration`` → `Forcing` (indefinite)
+        * With ``duration`` → `ForcingElement` (bounded)
 
     Examples
     --------
@@ -134,8 +134,8 @@ def create_constant_forcing(value, duration=None):
     value : float
         The constant value returned for all ``t``.
     duration : float, optional
-        If given, returns a :class:`~climatecritters.core.ForcingElement`.
-        If omitted, returns an indefinite :class:`~climatecritters.core.Forcing`.
+        If given, returns a `ForcingElement`.
+        If omitted, returns an indefinite `Forcing`.
 
     Returns
     -------
@@ -173,8 +173,8 @@ def create_sinusoid_forcing(A, period, y0=0.0, duration=None):
     y0 : float
         Constant offset.  Default 0.0.
     duration : float, optional
-        If given, returns a :class:`~climatecritters.core.ForcingElement`.
-        If omitted, returns an indefinite :class:`~climatecritters.core.Forcing`.
+        If given, returns a `ForcingElement`.
+        If omitted, returns an indefinite `Forcing`.
 
     Returns
     -------
@@ -222,8 +222,8 @@ def create_periodic_forcing(periods_powers, desired_amplitude=1, y0=0, duration=
     y0 : float
         Constant offset.  Default 0.
     duration : float, optional
-        If given, returns a :class:`~climatecritters.core.ForcingElement`.
-        If omitted, returns an indefinite :class:`~climatecritters.core.Forcing`.
+        If given, returns a `ForcingElement`.
+        If omitted, returns an indefinite `Forcing`.
 
     Returns
     -------
@@ -250,17 +250,17 @@ def create_periodic_forcing(periods_powers, desired_amplitude=1, y0=0, duration=
 
 
 def create_piecewise_forcing(elements, label="forcing"):
-    """Build a piecewise forcing from a sequence of :class:`~climatecritters.core.ForcingElement` parts.
+    """Build a piecewise forcing from a sequence of `ForcingElement` parts.
 
     Compiles the sequence immediately and returns a callable
-    :class:`~climatecritters.core.Forcing`.
+    `Forcing`.
 
     Parameters
     ----------
     elements : sequence of ForcingElement
-        Ordered :class:`~climatecritters.core.Hold`, :class:`~climatecritters.core.Ramp`,
-        :class:`~climatecritters.core.Harmonic`, or general
-        :class:`~climatecritters.core.ForcingElement` instances.
+        Ordered `Hold`, `Ramp`,
+        `Harmonic`, or general
+        `ForcingElement` instances.
     label : str
         Human-readable label.  Default ``'forcing'``.
 
@@ -294,12 +294,12 @@ def create_piecewise_forcing(elements, label="forcing"):
 # ---------------------------------------------------------------------------
 
 def make_forcing_element(forcing, duration=None):
-    """Convert a :class:`~climatecritters.core.Forcing` into a bounded
-    :class:`~climatecritters.core.ForcingElement`.
+    """Convert a `Forcing` into a bounded
+    `ForcingElement`.
 
-    This is the reverse of :meth:`~climatecritters.core.ForcingSequence.compile` —
+    This is the reverse of `ForcingSequence.compile` —
     it lets you embed an existing ``Forcing`` as a timed segment inside a
-    :class:`~climatecritters.core.ForcingSequence`.
+    `ForcingSequence`.
 
     Parameters
     ----------
diff --git a/climatecritters/utils/func.py b/climatecritters/utils/func.py
index b3e282c..dd26d88 100644
--- a/climatecritters/utils/func.py
+++ b/climatecritters/utils/func.py
@@ -17,9 +17,9 @@ def make_derivative_func(method='numpy', derivative=None, data=None, time=None):
     1. **Pass-through** — if ``derivative`` is already a callable, return it
        unchanged.
     2. **Numpy mode** — compute a finite-difference derivative via
-       ``np.gradient``, then fit a :class:`scipy.interpolate.CubicSpline` to
+       ``np.gradient``, then fit a ``scipy.interpolate.CubicSpline`` to
        the result.  Suitable for non-uniformly spaced data.
-    3. **Scipy mode** — fit a :class:`scipy.interpolate.CubicSpline` directly
+    3. **Scipy mode** — fit a ``scipy.interpolate.CubicSpline`` directly
        to ``data``, then return its analytical first derivative.  Slightly
        smoother than numpy mode for well-resolved data.
 
@@ -43,7 +43,7 @@ def make_derivative_func(method='numpy', derivative=None, data=None, time=None):
     deriv_func : callable
         Callable with signature ``f(t) -> float`` that returns the derivative
         of the forcing at time ``t``.  In numpy and scipy modes this is a
-        :class:`~scipy.interpolate.CubicSpline` (or its derivative object).
+        ``scipy.interpolate.CubicSpline`` (or its derivative object).
 
     Raises
     ------
diff --git a/climatecritters/utils/noise.py b/climatecritters/utils/noise.py
index a62c4a5..106164a 100644
--- a/climatecritters/utils/noise.py
+++ b/climatecritters/utils/noise.py
@@ -1,6 +1,6 @@
 """Noise generation and surrogate time series utilities.
 
-Thin wrappers around :class:`pyleoclim.SurrogateSeries` that expose surrogate
+Thin wrappers around ``pyleoclim.SurrogateSeries`` that expose surrogate
 and parametric noise generation through a consistent ClimateCritters interface.
 """
 
diff --git a/climatecritters/utils/solver.py b/climatecritters/utils/solver.py
index b322d8c..9ce3a67 100644
--- a/climatecritters/utils/solver.py
+++ b/climatecritters/utils/solver.py
@@ -1,8 +1,8 @@
 """Numerical integrators and solver support utilities.
 
 Provides fixed-step (RK4, Euler, Euler-Maruyama, Heun-Maruyama, Milstein)
-integrators that return a :class:`Solution` object, plus internal helpers
-used by :class:`~climatecritters.core.CCModel` for state validation and history
+integrators that return a `Solution` object, plus internal helpers
+used by `Model` for state validation and history
 reconstruction.
 """
 
@@ -112,7 +112,7 @@ def flux_divergence(face_fluxes, dz):
 
 
 def define_t_eval(t_span, delta_t=None, num_points=None):
-    """Build a ``t_eval`` array for use with :func:`scipy.integrate.solve_ivp`.
+    """Build a ``t_eval`` array for use with ``scipy.integrate.solve_ivp``.
 
     Parameters
     ----------
@@ -313,7 +313,7 @@ def euler_method(f, t_span, y0, dt, args=(), post_step=None):
     Notes
     -----
     Forward Euler is first-order accurate and can be unstable for stiff
-    systems or large ``dt``.  Prefer :func:`rk4_method` for most applications.
+    systems or large ``dt``.  Prefer `rk4_method` for most applications.
     """
     n_steps = int((t_span[1] - t_span[0]) / dt) + 1
     t = np.linspace(t_span[0], t_span[1], n_steps)
@@ -442,7 +442,7 @@ def heun_maruyama_method(f, t_span, y0, dt, si=None, noise_func=None, rng=None,
 
     A predictor-corrector scheme that achieves strong order 1.0 for SDEs with
     additive noise (diffusion independent of state) and weak order 2.0.  This
-    is a meaningful improvement over :func:`euler_maruyama_method` (strong
+    is a meaningful improvement over `euler_maruyama_method` (strong
     order 0.5) when transition timing is the quantity of interest, as in
     bistable climate models.
 
@@ -586,7 +586,7 @@ def milstein_method(f, t_span, y0, dt, si=None, noise_func=None, rng=None, args=
     Achieves strong order 1.0 for both additive *and* multiplicative noise
     (diagonal diffusion), making it the right choice when ``sde_noise``
     depends on the state.  For additive noise it reduces to Euler-Maruyama
-    plus a zero correction; prefer :func:`heun_maruyama_method` in that case
+    plus a zero correction; prefer `heun_maruyama_method` in that case
     as it also improves the drift approximation.
 
     Solves:
diff --git a/docs/_quarto.yml b/docs/_quarto.yml
index 8ef3b8f..f60ed16 100644
--- a/docs/_quarto.yml
+++ b/docs/_quarto.yml
@@ -33,8 +33,12 @@ interlinks:
       inv: objects.txt
   aliases:
     climatecritters.core.forcing: null
-    climatecritters.core.ccmodel: null
-    climatecritters.core.ccoutput: null
+    climatecritters.core.model: null
+    climatecritters.core.output: null
+    climatecritters.utils.forcing: null
+    climatecritters.utils.solver: null
+    climatecritters.utils.func: null
+    climatecritters.utils.noise: null
     pandas: pd
 
 execute:
diff --git a/docs/get-started/figures/quickstart_integrate.png b/docs/get-started/figures/quickstart_integrate.png
new file mode 100644
index 0000000..5b1cbe3
Binary files /dev/null and b/docs/get-started/figures/quickstart_integrate.png differ
diff --git a/docs/get-started/figures/quickstart_pyleo_plot.png b/docs/get-started/figures/quickstart_pyleo_plot.png
new file mode 100644
index 0000000..76b6fcd
Binary files /dev/null and b/docs/get-started/figures/quickstart_pyleo_plot.png differ
diff --git a/docs/get-started/figures/quickstart_pyleo_spectral.png b/docs/get-started/figures/quickstart_pyleo_spectral.png
new file mode 100644
index 0000000..b8f5081
Binary files /dev/null and b/docs/get-started/figures/quickstart_pyleo_spectral.png differ
diff --git a/docs/get-started/figures/quickstart_time_varying.png b/docs/get-started/figures/quickstart_time_varying.png
new file mode 100644
index 0000000..9fe78fc
Binary files /dev/null and b/docs/get-started/figures/quickstart_time_varying.png differ
diff --git a/docs/get-started/models.qmd b/docs/get-started/models.qmd
index e31a5cd..c711717 100644
--- a/docs/get-started/models.qmd
+++ b/docs/get-started/models.qmd
@@ -59,8 +59,6 @@ For method-evaluation purposes, the most relevant initial question is usually st
 
 *Standalone models drawn from the paleoclimate literature.*
 
-**[Daisyworld](../api/Daisyworld.qmd)** — Watson & Lovelock (1983). Two daisy species with different albedos colonize a planet and modulate its temperature. The system self-regulates temperature across a fairly wide range of solar luminosities, illustrating how biological feedbacks can stabilize climate. Multiple equilibria exist at the edges of the habitable range.
-
 **[ENSORechargeOscillator](../api/ENSORechargeOscillator.qmd)** — Jin (1997) recharge oscillator for ENSO. Two variables — sea surface temperature and thermocline heat content — produce quasi-periodic El Niño / La Niña cycles. Seasonal forcing can push the system into a more irregular or chaotic regime.
 
 **[Model3](../api/Model3.qmd)** — Ganopolski (2024) glacial cycle model. Tracks ice volume and a discrete glacial/deglaciation regime under orbital forcing. A regime-switch structure with hysteresis produces the asymmetric sawtooth of glacial cycles; a time-varying critical ice volume (`vc_func`) can reproduce the Mid-Pleistocene Transition from 41 kyr to 100 kyr dominant periodicity.
diff --git a/docs/get-started/quickstart.qmd b/docs/get-started/quickstart.qmd
index 7eadf52..79ace09 100644
--- a/docs/get-started/quickstart.qmd
+++ b/docs/get-started/quickstart.qmd
@@ -7,14 +7,15 @@ This page walks through a complete end-to-end example in about 15 lines of code.
 
 ## 1. Create a model
 
-Every ClimateCritters model takes a `forcing` argument. For the Lorenz 63 system the forcing adds an optional perturbation to the x-equation; here we set it to zero.
+Any model parameter or state variable can be driven by an external `Forcing`. For the Lorenz 63 system we register a (zero, for now) forcing on the x-equation via `register_forcing`.
 
 ```python
 import climatecritters as cc
 from climatecritters.model_critters import Lorenz63
 
+model = Lorenz63(sigma=10.0, rho=28.0, beta=8/3)
 forcing = cc.Forcing(lambda t: 0.0)
-model = Lorenz63(forcing=forcing, sigma=10.0, rho=28.0, beta=8/3)
+model.register_forcing('x', forcing, attachment_style='additive', timing='pre')
 ```
 
 ## 2. Integrate
@@ -29,6 +30,22 @@ output = model.integrate(
 )
 ```
 
+A quick look at the raw trajectory:
+
+```python
+import matplotlib.pyplot as plt
+
+fig, axes = plt.subplots(1, 2, figsize=(9, 3.5))
+axes[0].plot(output.state_variables['x'], output.state_variables['z'], lw=0.3, alpha=0.8)
+axes[0].set_xlabel('x'); axes[0].set_ylabel('z'); axes[0].set_title('phase portrait')
+axes[1].plot(output.time, output.state_variables['x'], lw=0.8)
+axes[1].set_xlabel('time'); axes[1].set_ylabel('x'); axes[1].set_title('x(t)')
+fig.tight_layout()
+plt.savefig('figures/quickstart_integrate.png', dpi=150, bbox_inches='tight')
+```
+
+![Lorenz 63 phase portrait and x(t) time series](figures/quickstart_integrate.png)
+
 ## 3. Inspect the output
 
 ```python
@@ -42,15 +59,20 @@ print(output.time[:5])                       # first five time points
 
 ```python
 ts_x = output.to_pyleo(var_names='x')
-ts_x.plot()
+ts_x.plot(savefig_settings={'path': 'figures/quickstart_pyleo_plot.png', 'dpi': 150})
 ```
 
-For a dashboard with spectral analysis:
+![Pyleoclim Series plot of the x variable](figures/quickstart_pyleo_plot.png)
+
+For a quick spectral analysis:
 
 ```python
-ts_x.dashboard()
+psd = ts_x.spectral()
+psd.plot(savefig_settings={'path': 'figures/quickstart_pyleo_spectral.png', 'dpi': 150})
 ```
 
+![Power spectral density of the x variable](figures/quickstart_pyleo_spectral.png)
+
 ## 5. Try a time-varying parameter
 
 Any parameter can be replaced with a callable at construction time. The callable must have signature `(t)`, `(t, state)`, or `(t, state, model)`.
@@ -58,14 +80,24 @@ Any parameter can be replaced with a callable at construction time. The callable
 ```python
 # rho increases linearly from 20 to 40 over the integration window
 model_tv = Lorenz63(
-    forcing=forcing,
     sigma=10.0,
     rho=lambda t: 20.0 + (20.0 / 50.0) * t,
     beta=8/3,
 )
+model_tv.register_forcing('x', forcing, attachment_style='additive', timing='pre')
 output_tv = model_tv.integrate(t_span=(0, 50), y0=[1, 1, 1], method='RK45')
 ```
 
+```python
+fig, ax = plt.subplots(figsize=(6, 3.5))
+ax.plot(output_tv.time, output_tv.state_variables['x'], lw=0.8)
+ax.set_xlabel('time'); ax.set_ylabel('x'); ax.set_title('x(t) with time-varying rho')
+fig.tight_layout()
+plt.savefig('figures/quickstart_time_varying.png', dpi=150, bbox_inches='tight')
+```
+
+![x(t) trajectory under a linearly increasing rho](figures/quickstart_time_varying.png)
+
 ## Next steps
 
 - [Core concepts](concepts.qmd) — understand `Model`, `Forcing`, and `Output` in depth
diff --git a/docs/scripts/run_quickstart.py b/docs/scripts/run_quickstart.py
new file mode 100644
index 0000000..e674807
--- /dev/null
+++ b/docs/scripts/run_quickstart.py
@@ -0,0 +1,74 @@
+"""Run every Python code block on the Quick Start page as one script.
+
+Extracts each fenced ```python code block from
+``docs/get-started/quickstart.qmd`` in document order and executes them
+sequentially in a single shared namespace, exactly as a reader following
+the page top-to-bottom would. This both smoke-tests that the snippets stay
+runnable and regenerates the figures referenced by the page (each block's
+``plt.savefig(...)`` writes into ``docs/get-started/figures/``).
+
+Usage
+-----
+```
+python scripts/run_quickstart.py
+```
+"""
+# Agg backend must be set before any other matplotlib import.
+import matplotlib
+matplotlib.use('Agg')
+
+import os
+import re
+import sys
+from pathlib import Path
+
+import matplotlib.pyplot as plt
+
+SCRIPT_DIR = Path(__file__).resolve().parent
+DOCS_DIR = SCRIPT_DIR.parent            # docs/
+PROJECT_DIR = DOCS_DIR.parent           # project root — contains climatecritters/
+PAGE = DOCS_DIR / 'get-started' / 'quickstart.qmd'
+
+FENCE_RE = re.compile(r'```python\s*\n(.*?)```', re.DOTALL)
+
+if str(PROJECT_DIR) not in sys.path:
+    sys.path.insert(0, str(PROJECT_DIR))
+
+
+def extract_blocks(qmd_text: str) -> list[str]:
+    """Return the body of every fenced ```python block, in document order."""
+    return FENCE_RE.findall(qmd_text)
+
+
+def main() -> None:
+    if not PAGE.exists():
+        print(f'Page not found: {PAGE}')
+        sys.exit(1)
+
+    blocks = extract_blocks(PAGE.read_text(encoding='utf-8'))
+    if not blocks:
+        print('No python code blocks found.')
+        sys.exit(0)
+
+    print(f'Running {len(blocks)} code block(s) from {PAGE.relative_to(PROJECT_DIR)}...')
+
+    # Figure paths in the page are relative to the page's own directory
+    # (e.g. 'figures/quickstart_integrate.png'), so run from there.
+    namespace = {'__name__': '__main__', '__builtins__': __builtins__}
+    os.chdir(PAGE.parent)
+
+    for i, block in enumerate(blocks, start=1):
+        print(f'  [{i}/{len(blocks)}]')
+        try:
+            exec(compile(block, f'<quickstart.qmd:block {i}>', 'exec'), namespace)  # noqa: S102
+        except Exception:
+            print(f'Block {i} failed:\n{block}')
+            raise
+        finally:
+            plt.close('all')
+
+    print('Done. All blocks ran successfully.')
+
+
+if __name__ == '__main__':
+    main()
diff --git a/docs/setup.md b/docs/setup.md
index 9df6106..66124eb 100644
--- a/docs/setup.md
+++ b/docs/setup.md
@@ -17,6 +17,7 @@ To rebuild the pages:
 cd docs
 python scripts/sync_notebooks.py
 python scripts/generate_tutorials_index.py
+python scripts/run_quickstart.py
 quartodoc build
 quarto render
 ```