diff --git a/src/include/serif/geometry/surface/automatic.hpp b/src/include/serif/geometry/surface/automatic.hpp
new file mode 100644
index 00000000..bdef120b
--- /dev/null
+++ b/src/include/serif/geometry/surface/automatic.hpp
@@ -0,0 +1,83 @@
+#pragma once
+
+#include "serif/geometry/surface/detail/automatic_relations.hpp"
+
+namespace serif::geometry::surface {
+    enum class PressureSurfaceSelection {
+        supported,
+        unavailable,
+        ambiguous
+    };
+
+    namespace detail {
+        /** @brief Strong boundary enforcement requires an actual solved, scalar H1 equation row. */
+        template <typename Field, typename Rows, typename SolvedTerms, typename PreferredField>
+        concept EligibleSurfaceRow = discretization::blocks::SpatialTerm<Field> && discretization::blocks::HasResidualBlock<Field> &&
+            utils::dependencies::ContainsKey<Field, Rows> && utils::dependencies::ContainsKey<Field, SolvedTerms> &&
+            (std::same_as<PreferredField, void> || std::same_as<PreferredField, Field>);
+
+        /** @brief Search one precedence level, without asking the EOS for a relation when pressure is already an unknown. */
+        template <int Priority, typename EOS, typename Field, typename Fields, typename Rows, typename SolvedTerms, typename PreferredField>
+        [[nodiscard]] consteval auto field_candidates() {
+            if constexpr (!EligibleSurfaceRow<Field, Rows, SolvedTerms, PreferredField>) {
+                return std::tuple<>{};
+            }
+
+            else if constexpr (Priority == 0) {
+                if constexpr (std::same_as<typename Field::Quantity, dimensions::Pressure>) {
+                    using Bindings = utils::dependencies::Bindings<utils::dependencies::Binding<dimensions::Pressure, Field>>;
+                    return std::tuple<std::type_identity<PressureFieldSurfaceFormulation<EOS, Bindings>>>{};
+                }
+                else {
+                    return std::tuple<>{};
+                }
+            }
+            else if constexpr (eos::models::EOSModel<EOS>) {
+                return relation_candidates<Priority == 1, EOS, Field, Fields>(typename EOS::Relations{});
+            }
+            else {
+                return std::tuple<>{};
+            }
+        }
+
+        /** @brief Prefer a pressure unknown, then an explicit inverse, then a direct pressure residual. */
+        template <int Priority, typename EOS, typename Rows, typename SolvedTerms, typename PreferredField, typename... Fields>
+        [[nodiscard]] consteval auto select_candidates(utils::dependencies::KeyList<Fields...> fields) {
+            using FieldsList = utils::dependencies::KeyList<Fields...>;
+            using Candidates = decltype(std::tuple_cat(field_candidates<Priority, EOS, Fields, FieldsList, Rows, SolvedTerms, PreferredField>()...));
+
+            if constexpr (std::tuple_size_v<Candidates> != 0 || Priority == 2) {
+                return Candidates{};
+            }
+            else {
+                return select_candidates<Priority + 1, EOS, Rows, SolvedTerms, PreferredField>(fields);
+            }
+        }
+
+        template <typename Candidates>
+        struct SelectedPressureSurface {
+            static constexpr PressureSurfaceSelection status = std::tuple_size_v<Candidates> == 0 ? PressureSurfaceSelection::unavailable : PressureSurfaceSelection::ambiguous;
+        };
+
+        template <typename Candidate>
+        struct SelectedPressureSurface<std::tuple<Candidate>> {
+            static constexpr PressureSurfaceSelection status = PressureSurfaceSelection::supported;
+            using Formulation = typename Candidate::type;
+        };
+    }
+
+    /**
+     * @brief Compile an isobaric condition from available physical fields, actual equation rows, and the EOS catalog.
+     * @note This searches declared algebraic relations, never performs inversion, and rejects equal-priority choices.
+     * Global thermodynamic inputs may opt in with thermodynamic_state = true; units alone do not establish that role.
+     * An optional preferred field narrows the row search without choosing an EOS relation for the model user.
+     */
+    template <typename EOS, utils::dependencies::IsKeyList AvailableTerms, utils::dependencies::IsKeyList ActualRows,
+        utils::dependencies::IsKeyList SolvedTerms = AvailableTerms, typename PreferredField = void>
+    struct AutomaticPressureSurface : detail::SelectedPressureSurface<decltype(detail::select_candidates<0, std::remove_cvref_t<EOS>, ActualRows, SolvedTerms, PreferredField>(detail::thermodynamic_fields(AvailableTerms{})))> {
+    };
+
+    /** @brief Whether the available equations and numerical backend select one implementable isobaric residual. */
+    template <typename EOS, typename AvailableTerms, typename ActualRows, typename SolvedTerms = AvailableTerms, typename PreferredField = void>
+    concept AutomaticPressureSurfaceCompilable = requires { typename AutomaticPressureSurface<EOS, AvailableTerms, ActualRows, SolvedTerms, PreferredField>::Formulation; };
+}
diff --git a/src/include/serif/geometry/surface/compilation.hpp b/src/include/serif/geometry/surface/compilation.hpp
new file mode 100644
index 00000000..2e5d3f25
--- /dev/null
+++ b/src/include/serif/geometry/surface/compilation.hpp
@@ -0,0 +1,108 @@
+#pragma once
+
+#include "serif/geometry/surface/detail/relations.hpp"
+#include "serif/geometry/surface/isobaric.hpp"
+
+namespace serif::geometry::surface
+{
+    /**
+     * @brief Whether the available state selects exactly one pressure-to-carrier EOS relation.
+     * @tparam Carrier Optional physical quantity whose equation enforces the boundary condition.
+     * @note An ambiguous catalog never selects its first entry. Specifying a carrier narrows the same search.
+     */
+    template <typename EOS, typename Bindings, typename Carrier = void>
+    concept PressureSurfaceCompilable = eos::models::EOSModel<EOS> && SurfaceStateBindings<Bindings> &&
+        (detail::MatchingPressureRelations<std::remove_cvref_t<EOS>, std::remove_cvref_t<Bindings>, Carrier, typename std::remove_cvref_t<EOS>::Relations>::count == 1);
+
+    template <typename EOS, typename Bindings, typename Carrier = void>
+    requires PressureSurfaceCompilable<EOS, Bindings, Carrier>
+    using pressure_relation_t = std::tuple_element_t<0, typename detail::MatchingPressureRelations<std::remove_cvref_t<EOS>, std::remove_cvref_t<Bindings>, Carrier, typename std::remove_cvref_t<EOS>::Relations>::Matches>;
+
+    namespace detail
+    {
+        /** @brief Common physical metadata; numerical values and preparation belong to the dependency context and cache. */
+        template <typename EOS, typename Bindings, typename RowQuantity, typename RelationT, bool ResolveCarrier>
+        struct PressureSurfaceDescription
+        {
+            using EquationOfState = std::remove_cvref_t<EOS>;
+            using StateBindings = std::remove_cvref_t<Bindings>;
+            using Relation = RelationT;
+            using PhysicalQuantity = dimensions::Pressure;
+            using ResidualQuantity = std::conditional_t<ResolveCarrier, RowQuantity, dimensions::Pressure>;
+            using CarrierQuantity = RowQuantity;
+            using CarrierTerm = utils::dependencies::bound_type_t<StateBindings, RowQuantity>;
+            using Test = CarrierTerm;
+            using Support = discretization::domain::StellarSurfaceBoundary;
+            using EOSStateQuantities = std::conditional_t<ResolveCarrier, typename RelationQuantities<Relation>::NonPressureInputs, typename RelationQuantities<Relation>::AllInputs>;
+            using StateQuantities = utils::dependencies::CombinedKeyList<utils::dependencies::KeyList<RowQuantity>, EOSStateQuantities>;
+            using StateTerms = BoundTerms<StateBindings, StateQuantities>;
+            using EOSStateTerms = BoundTerms<StateBindings, EOSStateQuantities>;
+            using Operands = typename SurfaceOperands<StateTerms>::Type;
+
+            static constexpr bool resolves_carrier = ResolveCarrier;
+            static constexpr bool uses_eos = true;
+        };
+    }
+
+    /** @brief The boundary residual is carrier minus the EOS value at prescribed pressure; no carrier quantity is privileged. */
+    template <typename EOS, typename Bindings, typename Carrier = void>
+    requires PressureSurfaceCompilable<EOS, Bindings, Carrier>
+    struct PressureSurfaceFormulation final : detail::PressureSurfaceDescription<EOS, Bindings, typename pressure_relation_t<EOS, Bindings, Carrier>::OutputQuantity, pressure_relation_t<EOS, Bindings, Carrier>, true>
+    {
+    };
+
+    /** @brief A selected explicit inverse relation; choosing it requires no numerical EOS inversion. */
+    template <eos::models::EOSModel EOS, SurfaceStateBindings Bindings, typename Relation>
+    requires eos::relations::CatalogContainingRelation<typename EOS::Relations, Relation> && (detail::matches_pressure_relation<EOS, Bindings, void, Relation>())
+    struct InversePressureSurfaceFormulation final : detail::PressureSurfaceDescription<EOS, Bindings, typename Relation::OutputQuantity, Relation, true>
+    {
+    };
+
+    /** @brief A direct pressure residual requires an explicitly chosen pressure relation and an input quantity whose equation row it replaces. */
+    template <typename EOS, typename Bindings, typename RowQuantity, typename Relation>
+    concept DirectPressureSurfaceCompilable = eos::models::EOSModel<EOS> && SurfaceStateBindings<Bindings> && eos::relations::ValidThermodynamicRelation<Relation> &&
+        eos::relations::CatalogContainingRelation<typename std::remove_cvref_t<EOS>::Relations, Relation> &&
+        std::same_as<typename Relation::OutputQuantity, dimensions::Pressure> && eos::relations::RelationContainingInput<Relation, RowQuantity> &&
+        SurfaceResidualRow<Bindings, RowQuantity> &&
+        utils::dependencies::SatisfiesKeyList<Bindings, typename detail::RelationQuantities<Relation>::AllInputs> &&
+        detail::supports_partials<EOS, Relation>(typename detail::RelationQuantities<Relation>::AllInputs{});
+
+    /**
+     * @brief Enforce P_EOS(state) - P_surface = 0 without inventing an inverse EOS relation.
+     * @note The explicit row and relation make the choice of residual, units, and Newton scaling visible.
+     */
+    template <typename EOS, typename Bindings, typename RowQuantity, typename Relation>
+    requires DirectPressureSurfaceCompilable<EOS, Bindings, RowQuantity, Relation>
+    struct DirectPressureSurfaceFormulation final : detail::PressureSurfaceDescription<EOS, Bindings, RowQuantity, Relation, false>
+    {
+    };
+
+    /** @brief An existing pressure unknown enforces P - P_surface directly, without evaluating the EOS. */
+    template <typename EOS, SurfaceStateBindings Bindings>
+    requires SurfaceResidualRow<Bindings, dimensions::Pressure>
+    struct PressureFieldSurfaceFormulation final : detail::PressureSurfaceDescription<EOS, Bindings, dimensions::Pressure, eos::relations::Relation<dimensions::Pressure, dimensions::Pressure>, true>
+    {
+        static constexpr bool uses_eos = false;
+    };
+
+    namespace detail
+    {
+        template <typename Candidate>
+        inline constexpr bool is_pressure_surface_formulation = false;
+
+        template <typename EOS, typename Bindings, typename Carrier>
+        inline constexpr bool is_pressure_surface_formulation<PressureSurfaceFormulation<EOS, Bindings, Carrier>> = true;
+
+        template <typename EOS, typename Bindings, typename Relation>
+        inline constexpr bool is_pressure_surface_formulation<InversePressureSurfaceFormulation<EOS, Bindings, Relation>> = true;
+
+        template <typename EOS, typename Bindings, typename RowQuantity, typename Relation>
+        inline constexpr bool is_pressure_surface_formulation<DirectPressureSurfaceFormulation<EOS, Bindings, RowQuantity, Relation>> = true;
+
+        template <typename EOS, typename Bindings>
+        inline constexpr bool is_pressure_surface_formulation<PressureFieldSurfaceFormulation<EOS, Bindings>> = true;
+    }
+
+    template <typename Candidate>
+    concept PressureSurfaceFormulationType = detail::is_pressure_surface_formulation<Candidate>;
+}
diff --git a/src/include/serif/geometry/surface/computations.hpp b/src/include/serif/geometry/surface/computations.hpp
new file mode 100644
index 00000000..4d432a2f
--- /dev/null
+++ b/src/include/serif/geometry/surface/computations.hpp
@@ -0,0 +1,108 @@
+#pragma once
+
+#include <tuple>
+#include <type_traits>
+#include <utility>
+
+#include "serif/geometry/surface/detail/batches.hpp"
+#include "serif/utils/dependencies/computation_graph.hpp"
+
+namespace serif::geometry::surface
+{
+    namespace detail
+    {
+        /** @brief Declare the exact EOS inputs from the formulation, so derivative metadata and cache dependencies cannot diverge. */
+        template <PressureSurfaceFormulationType Formulation, template <typename> typename SampleKey, typename... Terms>
+        [[nodiscard]] auto make_surface_eos_computation(utils::dependencies::KeyList<Terms...>)
+        {
+            using namespace utils::dependencies;
+            using EOS = typename Formulation::EquationOfState;
+
+            if constexpr (!Formulation::uses_eos)
+            {
+                // Reuse the carrier residual stage with one constant target; pressure unknowns need no EOS dependency.
+                return make_computation<keys::EOSResponse, InputList<keys::PressureCondition>>(
+                    [](const IsobaricSurface& condition)
+                    {
+                        return std::vector<SurfaceEOSResponse<Formulation>>{{.value = condition.target_pressure()}};
+                    });
+            }
+            else if constexpr (Formulation::resolves_carrier)
+            {
+                return make_computation<keys::EOSResponse, InputList<keys::EquationOfState, keys::PressureCondition, SampleKey<Terms>...>>(
+                    [](const EOS& equation_of_state, const IsobaricSurface& condition, const BoundaryValues<Terms>&... values)
+                    {
+                        return evaluate_response_samples<Formulation>(equation_of_state, condition.target_pressure(), values...);
+                    });
+            }
+            else
+            {
+                return make_computation<keys::EOSResponse, InputList<keys::EquationOfState, SampleKey<Terms>...>>(
+                    [](const EOS& equation_of_state, const BoundaryValues<Terms>&... values)
+                    {
+                        return evaluate_response_samples<Formulation>(equation_of_state, dimensions::PressureValue{0.0}, values...);
+                    });
+            }
+        }
+
+        /** @brief Keep the inexpensive boundary residual separate from the reusable EOS response and its partials. */
+        template <PressureSurfaceFormulationType Formulation, template <typename> typename SampleKey>
+        [[nodiscard]] auto make_surface_linearization_computation()
+        {
+            using namespace utils::dependencies;
+            using Responses = std::vector<SurfaceEOSResponse<Formulation>>;
+
+            if constexpr (Formulation::resolves_carrier)
+            {
+                using Carrier = typename Formulation::CarrierTerm;
+                return make_computation<keys::Linearization, InputList<keys::EOSResponse, SampleKey<Carrier>>>(
+                    [](const Responses& responses, const BoundaryValues<Carrier>& carrier)
+                    {
+                        return linearize_carrier_samples<Formulation>(responses, carrier);
+                    });
+            }
+            else
+            {
+                return make_computation<keys::Linearization, InputList<keys::EOSResponse, keys::PressureCondition>>(
+                    [](const Responses& responses, const IsobaricSurface& condition)
+                    {
+                        return linearize_pressure_samples<Formulation>(condition, responses);
+                    });
+            }
+        }
+    }
+
+    /** @brief Expose the existing EOS and residual stages for composition with backend sampling or other physics. */
+    template <PressureSurfaceFormulationType Formulation, template <typename> typename SampleKey = std::type_identity_t>
+    [[nodiscard]] auto make_pressure_surface_computations()
+    {
+        return std::tuple{
+            detail::make_surface_eos_computation<Formulation, SampleKey>(typename Formulation::EOSStateTerms{}),
+            detail::make_surface_linearization_computation<Formulation, SampleKey>()};
+    }
+
+    /**
+     * @brief Describe the isobaric calculation as an ordinary dependency graph, without evaluating the EOS.
+     * @tparam Formulation A carrier-resolved or explicitly selected direct-pressure surface formulation.
+     * @return Two computations: the EOS response and the boundary residual with analytical derivatives.
+     * @note Supply pressure_condition and boundary_values, and bind keys::EquationOfState with a tracking policy when the formulation uses the EOS.
+     * A pressure-field formulation needs no EOS source. Use make_computation_cache<SurfaceFailure> and its normal try_prepare/get interface;
+     * the cache owns the complete preparation lifecycle.
+     */
+    template <PressureSurfaceFormulationType Formulation>
+    [[nodiscard]] auto make_pressure_surface_graph()
+    {
+        return std::apply([](auto... computations)
+        {
+            return utils::dependencies::make_computation_graph(std::move(computations)...);
+        }, make_pressure_surface_computations<Formulation>());
+    }
+
+    /** @brief Infer a unique carrier-resolved formulation directly from an EOS and the available physical state. */
+    template <typename EOS, typename Bindings, typename Carrier = void>
+    requires PressureSurfaceCompilable<EOS, Bindings, Carrier>
+    [[nodiscard]] auto make_pressure_surface_graph()
+    {
+        return make_pressure_surface_graph<PressureSurfaceFormulation<EOS, Bindings, Carrier>>();
+    }
+}
diff --git a/src/include/serif/geometry/surface/detail/automatic_relations.hpp b/src/include/serif/geometry/surface/detail/automatic_relations.hpp
new file mode 100644
index 00000000..a50d7e5c
--- /dev/null
+++ b/src/include/serif/geometry/surface/detail/automatic_relations.hpp
@@ -0,0 +1,87 @@
+#pragma once
+
+#include "serif/geometry/surface/compilation.hpp"
+#include "serif/discretization/blocks/traits.hpp"
+
+namespace serif::geometry::surface::detail {
+    /** @brief Globals require an explicit physical role: matching units alone would mistake C or a multiplier for enthalpy. */
+    template <typename Term>
+    concept ThermodynamicBoundaryField = discretization::blocks::ScalarTraceTerm<Term> &&
+        discretization::blocks::HasQuantity<Term> && dimensions::ThermodynamicQuantityType<typename Term::Quantity> &&
+        (discretization::blocks::SpatialTerm<Term> || requires { requires Term::thermodynamic_state; });
+
+    /** @brief Collect physical fields, retaining distinct terms with equal quantities so ambiguity remains visible. */
+    template <typename... Terms>
+    [[nodiscard]] consteval auto thermodynamic_fields(utils::dependencies::KeyList<Terms...>) {
+        return utils::dependencies::CombinedKeyList<std::conditional_t<ThermodynamicBoundaryField<Terms>, utils::dependencies::KeyList<Terms>, utils::dependencies::KeyList<>>...>{};
+    }
+
+    /** @brief Return a unique field for a quantity; the chosen residual carrier supplies its own quantity unambiguously. */
+    template <typename Quantity, typename Carrier, typename... Fields>
+    [[nodiscard]] consteval auto quantity_field(utils::dependencies::KeyList<Fields...>) {
+        using Matches = decltype(std::tuple_cat(std::declval<std::conditional_t<std::same_as<Quantity, typename Fields::Quantity>, std::tuple<Fields>, std::tuple<>>>()...));
+        if constexpr (std::same_as<Quantity, typename Carrier::Quantity>) {
+            return std::type_identity<Carrier>{};
+        }
+        else if constexpr (std::tuple_size_v<Matches> == 1) {
+            return std::type_identity<std::tuple_element_t<0, Matches>>{};
+        }
+        else {
+            return std::type_identity<void>{};
+        }
+    }
+
+    /** @brief Bind only the quantities read by this candidate, ignoring unrelated fields and quantities. */
+    template <typename Carrier, typename Fields, typename... Quantities>
+    [[nodiscard]] consteval auto relation_bindings(utils::dependencies::KeyList<Quantities...>) {
+        using namespace utils::dependencies;
+
+        if constexpr ((!std::same_as<void, typename decltype(quantity_field<Quantities, Carrier>(Fields{}))::type> && ...)) {
+            return std::type_identity<Bindings<Binding<Quantities, typename decltype(quantity_field<Quantities, Carrier>(Fields{}))::type>...>>{};
+        }
+        else {
+            return std::type_identity<void>{};
+        }
+    }
+
+    /** @brief Test a declared inverse and every required analytic partial before publishing its formulation. */
+    template <typename EOS, typename Carrier, typename Fields, typename Relation>
+    [[nodiscard]] consteval auto inverse_candidate() {
+        using namespace utils::dependencies;
+        using Quantities = CombinedKeyList<KeyList<typename Carrier::Quantity>, typename surface::detail::RelationQuantities<Relation>::NonPressureInputs>;
+        using Bindings = typename decltype(relation_bindings<Carrier, Fields>(Quantities{}))::type;
+
+        if constexpr (SurfaceStateBindings<Bindings> && surface::detail::matches_pressure_relation<EOS, Bindings, typename Carrier::Quantity, Relation>()) {
+            return std::tuple<std::type_identity<InversePressureSurfaceFormulation<EOS, Bindings, Relation>>>{};
+        }
+        else {
+            return std::tuple<>{};
+        }
+    }
+
+    /** @brief A direct pressure relation may replace only a row belonging to one of its solved inputs. */
+    template <typename EOS, typename Carrier, typename Fields, typename Relation>
+    [[nodiscard]] consteval auto direct_candidate() {
+        using namespace utils::dependencies;
+        using Quantities = CombinedKeyList<KeyList<typename Carrier::Quantity>, typename surface::detail::RelationQuantities<Relation>::AllInputs>;
+        using Bindings = typename decltype(relation_bindings<Carrier, Fields>(Quantities{}))::type;
+
+        if constexpr (DirectPressureSurfaceCompilable<EOS, Bindings, typename Carrier::Quantity, Relation>) {
+            return std::tuple<std::type_identity<DirectPressureSurfaceFormulation<EOS, Bindings, typename Carrier::Quantity, Relation>>>{};
+        }
+        else {
+            return std::tuple<>{};
+        }
+    }
+
+    /** @brief Keep every admissible relation; catalog order never selects a physical implementation. */
+    template <bool Inverse, typename EOS, typename Carrier, typename Fields, typename... Relations>
+    [[nodiscard]] consteval auto relation_candidates(eos::relations::RelationCatalog<Relations...>) {
+        if constexpr (Inverse) {
+            return std::tuple_cat(inverse_candidate<EOS, Carrier, Fields, Relations>()...);
+        }
+        else {
+            return std::tuple_cat(direct_candidate<EOS, Carrier, Fields, Relations>()...);
+        }
+    }
+}
diff --git a/src/include/serif/geometry/surface/detail/batches.hpp b/src/include/serif/geometry/surface/detail/batches.hpp
new file mode 100644
index 00000000..44f78813
--- /dev/null
+++ b/src/include/serif/geometry/surface/detail/batches.hpp
@@ -0,0 +1,114 @@
+#pragma once
+
+#include <array>
+
+#include "serif/geometry/surface/detail/evaluation.hpp"
+#include "serif/geometry/surface/samples.hpp"
+
+namespace serif::geometry::surface::detail
+{
+    /** @brief All EOS columns must describe the same ordered, nonempty set of boundary points. */
+    template <typename... Terms>
+    [[nodiscard]] std::expected<std::size_t, SurfaceFailure> surface_sample_count(const BoundaryValues<Terms>&... columns)
+    {
+        if constexpr (sizeof...(Terms) == 0)
+        {
+            // A pressure-only relation has one constant response, broadcast over the carrier samples later.
+            return 1;
+        }
+        else
+        {
+            const std::array<std::size_t, sizeof...(Terms)> counts{columns.values.size()...};
+            if (counts.front() == 0)
+            {
+                return std::unexpected(SurfaceFailure{.reason = SurfaceFailureReason::empty_surface, .message = "A surface calculation requires at least one boundary sample."});
+            }
+            for (const std::size_t count : counts)
+            {
+                if (count != counts.front())
+                {
+                    return std::unexpected(SurfaceFailure{.reason = SurfaceFailureReason::inconsistent_sample_count, .message = "Surface terms must be sampled at the same ordered boundary points; their sample counts differ."});
+                }
+            }
+            return counts.front();
+        }
+    }
+
+    /** @brief Assemble the conventional equation point from typed columns; no mesh assumptions enter the physical calculation. */
+    template <PressureSurfaceFormulationType Formulation, typename... Terms>
+    [[nodiscard]] SurfacePoint<Formulation> surface_point(const std::size_t sample, const BoundaryValues<Terms>&... columns)
+    {
+        SurfacePoint<Formulation> point{};
+        ((point.fields.template get<Terms>() = columns.values[sample]), ...);
+        return point;
+    }
+
+    /** @brief Collect one ordered batch, attaching the first failure's sample index before any result is published. */
+    template <typename Value, typename Evaluate>
+    [[nodiscard]] std::expected<std::vector<Value>, SurfaceFailure> collect_surface_samples(std::size_t count, Evaluate evaluate)
+    {
+        std::vector<Value> result;
+        result.reserve(count);
+        for (std::size_t sample = 0; sample < count; ++sample)
+        {
+            auto value = evaluate(sample);
+            if (!value)
+            {
+                SurfaceFailure failure = value.error();
+                failure.sample = sample;
+                return std::unexpected(failure);
+            }
+            result.push_back(std::move(*value));
+        }
+        return result;
+    }
+
+    /** @brief Evaluate the EOS once per state sample, or once for a genuinely pressure-only relation. */
+    template <PressureSurfaceFormulationType Formulation, typename... Terms>
+    [[nodiscard]] std::expected<std::vector<SurfaceEOSResponse<Formulation>>, SurfaceFailure> evaluate_response_samples(const typename Formulation::EquationOfState& equation_of_state, const dimensions::PressureValue pressure, const BoundaryValues<Terms>&... columns)
+    {
+        const auto count = surface_sample_count(columns...);
+        if (!count)
+        {
+            return std::unexpected(count.error());
+        }
+        return collect_surface_samples<SurfaceEOSResponse<Formulation>>(*count, [&](std::size_t sample)
+        {
+            return evaluate_surface_response<Formulation>(equation_of_state, pressure, surface_point<Formulation>(sample, columns...));
+        });
+    }
+
+    /** @brief Broadcast only a pressure-only relation; every state-dependent carrier sample needs its own response. */
+    template <PressureSurfaceFormulationType Formulation>
+    [[nodiscard]] std::expected<std::vector<SurfaceLinearization<Formulation>>, SurfaceFailure> linearize_carrier_samples(const std::vector<SurfaceEOSResponse<Formulation>>& responses, const BoundaryValues<typename Formulation::CarrierTerm>& carrier)
+    {
+        constexpr bool constant_response = std::same_as<typename Formulation::EOSStateTerms, utils::dependencies::KeyList<>>;
+        const auto count = surface_sample_count(carrier);
+        if (!count)
+        {
+            return std::unexpected(count.error());
+        }
+        if (responses.size() != (constant_response ? 1 : *count))
+        {
+            return std::unexpected(SurfaceFailure{.reason = SurfaceFailureReason::inconsistent_sample_count, .message = "The carrier and EOS inputs must refer to the same boundary points."});
+        }
+        return collect_surface_samples<SurfaceLinearization<Formulation>>(*count, [&](std::size_t sample) -> std::expected<SurfaceLinearization<Formulation>, SurfaceFailure>
+        {
+            if (!std::isfinite(carrier.values[sample]))
+            {
+                return std::unexpected(SurfaceFailure{.reason = SurfaceFailureReason::non_finite_state, .message = "A boundary carrier value is non-finite."});
+            }
+            return checked_linearization<Formulation>(dimensions::PressureValue{0.0}, surface_point<Formulation>(sample, carrier), responses[constant_response ? 0 : sample]);
+        });
+    }
+
+    /** @brief Changing prescribed pressure reuses all values and partials of a direct-pressure EOS relation. */
+    template <PressureSurfaceFormulationType Formulation>
+    [[nodiscard]] std::expected<std::vector<SurfaceLinearization<Formulation>>, SurfaceFailure> linearize_pressure_samples(const IsobaricSurface& condition, const std::vector<SurfaceEOSResponse<Formulation>>& responses)
+    {
+        return collect_surface_samples<SurfaceLinearization<Formulation>>(responses.size(), [&](std::size_t sample)
+        {
+            return checked_linearization<Formulation>(condition.target_pressure(), SurfacePoint<Formulation>{}, responses[sample]);
+        });
+    }
+} // namespace serif::geometry::surface::detail
diff --git a/src/include/serif/geometry/surface/detail/evaluation.hpp b/src/include/serif/geometry/surface/detail/evaluation.hpp
new file mode 100644
index 00000000..5bb2eaca
--- /dev/null
+++ b/src/include/serif/geometry/surface/detail/evaluation.hpp
@@ -0,0 +1,135 @@
+#pragma once
+
+#include <cmath>
+#include <expected>
+
+#include "serif/eos/evaluation.hpp"
+#include "serif/geometry/surface/linearization.hpp"
+
+namespace serif::geometry::surface::detail
+{
+    /** @brief Read one thermodynamic role from its bound term, preserving the EOS argument's physical type. */
+    template <typename Quantity, PressureSurfaceFormulationType Formulation>
+    [[nodiscard]] dimensions::QuantityValue<Quantity> state_value(const SurfacePoint<Formulation>& point)
+    {
+        using Term = utils::dependencies::bound_type_t<typename Formulation::StateBindings, Quantity>;
+        return dimensions::QuantityValue<Quantity>{point.fields.template get<Term>()};
+    }
+
+    /** @brief Only the carrier-resolved formulation substitutes prescribed pressure for an EOS input. */
+    template <typename Quantity, PressureSurfaceFormulationType Formulation>
+    [[nodiscard]] dimensions::QuantityValue<Quantity> relation_input(const dimensions::PressureValue pressure, const SurfacePoint<Formulation>& point)
+    {
+        if constexpr (Formulation::resolves_carrier && std::same_as<Quantity, dimensions::Pressure>)
+        {
+            return pressure;
+        }
+        else
+        {
+            return state_value<Quantity, Formulation>(point);
+        }
+    }
+
+    /** @brief Validate precisely the state read by this calculation, leaving unrelated state outside its dependencies. */
+    template <PressureSurfaceFormulationType Formulation, typename... Quantities>
+    [[nodiscard]] bool finite_state(const SurfacePoint<Formulation>& point, utils::dependencies::KeyList<Quantities...>)
+    {
+        return (std::isfinite(state_value<Quantities, Formulation>(point).value()) && ...);
+    }
+
+    /** @brief Add each quantity's chain-rule contribution, including distinct roles bound to the same unknown. */
+    template <typename Quantity, PressureSurfaceFormulationType Formulation, typename... Inputs>
+    void accumulate_partial(const typename Formulation::EquationOfState& equation_of_state, const dimensions::PressureValue pressure, const SurfacePoint<Formulation>& point, SurfaceEOSResponse<Formulation>& response, dimensions::QuantityList<Inputs...>)
+    {
+        using Output = typename Formulation::Relation::OutputQuantity;
+        using Term = utils::dependencies::bound_type_t<typename Formulation::StateBindings, Quantity>;
+
+        const dimensions::PartialDerivative<Output, Quantity> derivative = eos::partial_derivative<Output, Quantity>(equation_of_state, relation_input<Inputs, Formulation>(pressure, point)...);
+        response.field_derivatives.template get<Term>() += derivative.value();
+    }
+
+    /** @brief Evaluate one selected relation and all required partials without any hidden EOS inversion. */
+    template <PressureSurfaceFormulationType Formulation, typename... Inputs, typename... StateQuantities>
+    [[nodiscard]] SurfaceEOSResponse<Formulation> evaluate_relation(const typename Formulation::EquationOfState& equation_of_state, const dimensions::PressureValue pressure, const SurfacePoint<Formulation>& point, dimensions::QuantityList<Inputs...> inputs, utils::dependencies::KeyList<StateQuantities...>)
+    {
+        using Output = typename Formulation::Relation::OutputQuantity;
+
+        SurfaceEOSResponse<Formulation> response{};
+        response.value = eos::evaluate<Output>(equation_of_state, relation_input<Inputs, Formulation>(pressure, point)...);
+        (accumulate_partial<StateQuantities, Formulation>(equation_of_state, pressure, point, response, inputs), ...);
+        return response;
+    }
+
+    /**
+     * @brief Reject non-finite thermodynamic inputs, EOS outputs, and partials as values.
+     * @note Existing EOS implementations can still throw their documented domain errors; these propagate to the caller.
+     */
+    template <PressureSurfaceFormulationType Formulation>
+    [[nodiscard]] std::expected<SurfaceEOSResponse<Formulation>, SurfaceFailure> evaluate_surface_response(const typename Formulation::EquationOfState& equation_of_state, const dimensions::PressureValue pressure, const SurfacePoint<Formulation>& point)
+    {
+        if constexpr (!Formulation::uses_eos)
+        {
+            // A pressure unknown uses the prescribed pressure itself as its target, without an EOS relation.
+            return SurfaceEOSResponse<Formulation>{.value = pressure};
+        }
+        else
+        {
+            if (!finite_state<Formulation>(point, typename Formulation::EOSStateQuantities{}))
+            {
+                return std::unexpected(SurfaceFailure{.reason = SurfaceFailureReason::non_finite_state, .message = "The surface EOS inputs must be finite; supply an admissible boundary state."});
+            }
+
+            const SurfaceEOSResponse<Formulation> response = evaluate_relation<Formulation>(equation_of_state, pressure, point, typename Formulation::Relation::InputQuantities{}, typename Formulation::EOSStateQuantities{});
+            bool finite = std::isfinite(response.value.value());
+            for (const double derivative : response.field_derivatives.values)
+            {
+                finite = finite && std::isfinite(derivative);
+            }
+            if (!finite)
+            {
+                return std::unexpected(SurfaceFailure{.reason = SurfaceFailureReason::non_finite_eos_response, .message = "The surface EOS returned a non-finite value or partial derivative; check its thermodynamic domain."});
+            }
+            return response;
+        }
+    }
+
+    /** @brief Form the declared residual from an already evaluated EOS response, keeping the carrier identity explicit. */
+    template <PressureSurfaceFormulationType Formulation>
+    [[nodiscard]] SurfaceLinearization<Formulation> linearize_surface(const dimensions::PressureValue pressure, const SurfacePoint<Formulation>& point, const SurfaceEOSResponse<Formulation>& response)
+    {
+        SurfaceLinearization<Formulation> result{};
+        result.field_derivatives = response.field_derivatives;
+
+        if constexpr (Formulation::resolves_carrier)
+        {
+            result.residual = state_value<typename Formulation::CarrierQuantity, Formulation>(point).value() - response.value.value();
+            for (double& derivative : result.field_derivatives.values)
+            {
+                derivative = -derivative;
+            }
+            result.field_derivatives.template get<typename Formulation::CarrierTerm>() += 1.0;
+        }
+        else
+        {
+            result.residual = response.value.value() - pressure.value();
+        }
+        return result;
+    }
+
+    /** @brief Catch overflow in residual subtraction or derivative accumulation before publishing a linearization. */
+    template <PressureSurfaceFormulationType Formulation>
+    [[nodiscard]] std::expected<SurfaceLinearization<Formulation>, SurfaceFailure> checked_linearization(const dimensions::PressureValue pressure, const SurfacePoint<Formulation>& point, const SurfaceEOSResponse<Formulation>& response)
+    {
+        const SurfaceLinearization<Formulation> result = linearize_surface<Formulation>(pressure, point, response);
+        bool finite = std::isfinite(result.residual);
+        for (const double derivative : result.field_derivatives.values)
+        {
+            finite = finite && std::isfinite(derivative);
+        }
+        if (!finite)
+        {
+            return std::unexpected(SurfaceFailure{.reason = SurfaceFailureReason::non_finite_linearization, .message = "The surface residual or derivative is non-finite; reduce the trial state before using this linearization."});
+        }
+        return result;
+    }
+}
diff --git a/src/include/serif/geometry/surface/detail/relations.hpp b/src/include/serif/geometry/surface/detail/relations.hpp
new file mode 100644
index 00000000..f378d7d1
--- /dev/null
+++ b/src/include/serif/geometry/surface/detail/relations.hpp
@@ -0,0 +1,97 @@
+#pragma once
+
+#include <tuple>
+#include <type_traits>
+
+#include "serif/discretization/forms/value.hpp"
+#include "serif/eos/models/concepts.hpp"
+#include "serif/utils/dependencies/key_list_operations.hpp"
+
+namespace serif::geometry::surface {
+    /** @brief Scalar thermodynamic data may come from any declared scalar space; boundary enforcement is a separate backend capability. */
+    template <typename Term>
+    concept SurfaceScalarTerm = discretization::blocks::Term<Term> && requires {
+        typename discretization::forms::OperandValue<discretization::forms::Operand<Term>>;
+        requires std::same_as<discretization::forms::OperandValue<discretization::forms::Operand<Term>>, double>;
+    };
+
+    namespace detail {
+        template <typename Bindings>
+        struct AreSurfaceBindings : std::false_type {
+        };
+
+        template <typename... Entries>
+        struct AreSurfaceBindings<utils::dependencies::Bindings<Entries...>> : std::bool_constant<
+            sizeof...(Entries) != 0 && (dimensions::ThermodynamicQuantityType<typename Entries::Key> && ...) &&
+            (SurfaceScalarTerm<typename Entries::Target> && ...)> {
+        };
+    }
+
+    /** @brief Reuse dependency bindings, retaining distinct quantity roles even when they refer to one term. */
+    template <typename Bindings>
+    concept SurfaceStateBindings = utils::dependencies::IsBindings<Bindings> && detail::AreSurfaceBindings<std::remove_cvref_t<Bindings>>::value;
+
+    /** @brief The quantity enforcing the boundary must have an actual equation row, not just sampled values. */
+    template <typename Bindings, typename Quantity>
+    concept SurfaceResidualRow = utils::dependencies::BindingsWithBinding<Bindings, Quantity> &&
+        discretization::blocks::HasResidualBlock<utils::dependencies::bound_type_t<Bindings, Quantity>>;
+
+    namespace detail {
+        template <typename Relation>
+        struct RelationQuantities;
+
+        template <typename Output, typename... Inputs>
+        struct RelationQuantities<eos::relations::Relation<Output, Inputs...>> {
+            using AllInputs = utils::dependencies::KeyList<Inputs...>;
+            using NonPressureInputs = utils::dependencies::CombinedKeyList<
+                std::conditional_t<std::same_as<Inputs, dimensions::Pressure>, utils::dependencies::KeyList<>, utils::dependencies::KeyList<Inputs>>...>;
+        };
+
+        /** @brief Check only the derivatives of unknown state; prescribed pressure has no Newton column. */
+        template <typename EOS, typename Relation, typename... Quantities>
+        [[nodiscard]] consteval bool supports_partials(utils::dependencies::KeyList<Quantities...>) {
+            return (eos::models::EOSSupportingPartialDerivative<EOS, Relation, Quantities> && ...);
+        }
+
+        /** @brief A resolved carrier must be supplied by one pressure-input relation, with all remaining inputs and partials available. */
+        template <typename EOS, typename Bindings, typename Carrier, typename Relation>
+        [[nodiscard]] consteval bool matches_pressure_relation() {
+            using Output = typename Relation::OutputQuantity;
+            using Inputs = typename RelationQuantities<Relation>::NonPressureInputs;
+
+            return !std::same_as<Output, dimensions::Pressure> &&
+                (std::same_as<Carrier, void> || std::same_as<Carrier, Output>) &&
+                eos::relations::RelationContainingInput<Relation, dimensions::Pressure> &&
+                SurfaceResidualRow<Bindings, Output> &&
+                utils::dependencies::SatisfiesKeyList<Bindings, Inputs> && supports_partials<EOS, Relation>(Inputs{});
+        }
+
+        template <typename EOS, typename Bindings, typename Carrier, typename Catalog>
+        struct MatchingPressureRelations;
+
+        template <typename EOS, typename Bindings, typename Carrier, typename... Relations>
+        struct MatchingPressureRelations<EOS, Bindings, Carrier, eos::relations::RelationCatalog<Relations...>> {
+            using Matches = decltype(std::tuple_cat(std::declval<std::conditional_t<matches_pressure_relation<EOS, Bindings, Carrier, Relations>(), std::tuple<Relations>, std::tuple<>>>()...));
+            static constexpr std::size_t count = std::tuple_size_v<Matches>;
+        };
+
+        template <typename Terms>
+        struct SurfaceOperands;
+
+        template <typename... Terms>
+        struct SurfaceOperands<utils::dependencies::KeyList<Terms...>> {
+            using Type = utils::types::TypeList<discretization::forms::Operand<Terms>...>;
+        };
+
+        template <typename Terms>
+        struct UniqueTerms;
+
+        template <typename... Terms>
+        struct UniqueTerms<utils::types::TypeList<Terms...>> {
+            using Type = utils::dependencies::CombinedKeyList<utils::dependencies::KeyList<Terms>...>;
+        };
+
+        template <typename Bindings, typename Quantities>
+        using BoundTerms = typename UniqueTerms<utils::dependencies::ResolvedTargets<Bindings, Quantities>>::Type;
+    }
+}
diff --git a/src/include/serif/geometry/surface/evaluation.hpp b/src/include/serif/geometry/surface/evaluation.hpp
new file mode 100644
index 00000000..805b6058
--- /dev/null
+++ b/src/include/serif/geometry/surface/evaluation.hpp
@@ -0,0 +1,26 @@
+#pragma once
+
+#include "serif/geometry/surface/detail/evaluation.hpp"
+
+namespace serif::geometry::surface {
+    /**
+     * @brief Evaluate the chosen isobaric residual and its analytical derivatives at one boundary point.
+     * @param equation_of_state The EOS declared by Formulation.
+     * @param condition The prescribed surface pressure, independent of the choice of carrier.
+     * @param point Values of the formulation's bound scalar terms, already sampled at a common boundary point.
+     * @return A conventional equation linearization, or a numerical rejection. No preparation or borrowed state is retained.
+     * @note Domain exceptions from an existing throwing EOS propagate; no exception is used for surface validation.
+     */
+    template <PressureSurfaceFormulationType Formulation>
+    [[nodiscard]] std::expected<SurfaceLinearization<Formulation>, SurfaceFailure> evaluate_pressure_surface(const typename Formulation::EquationOfState& equation_of_state, const IsobaricSurface& condition, const SurfacePoint<Formulation>& point) {
+        if (!detail::finite_state<Formulation>(point, typename Formulation::StateQuantities{})) {
+            return std::unexpected(SurfaceFailure{.reason = SurfaceFailureReason::non_finite_state, .message = "The surface state must be finite; supply finite values for every required thermodynamic quantity."});
+        }
+
+        const auto response = detail::evaluate_surface_response<Formulation>(equation_of_state, condition.target_pressure(), point);
+        if (!response) {
+            return std::unexpected(response.error());
+        }
+        return detail::checked_linearization<Formulation>(condition.target_pressure(), point, *response);
+    }
+}
diff --git a/src/include/serif/geometry/surface/isobaric.hpp b/src/include/serif/geometry/surface/isobaric.hpp
index e69de29b..8b902a64 100644
--- a/src/include/serif/geometry/surface/isobaric.hpp
+++ b/src/include/serif/geometry/surface/isobaric.hpp
@@ -0,0 +1,40 @@
+#pragma once
+
+#include <array>
+#include <type_traits>
+
+#include "serif/dimensions/dimensions.hpp"
+#include "serif/dimensions/type_alias.hpp"
+
+namespace serif::geometry::surface {
+    struct PressureSurfaceDescriptor final {
+        double target_pressure;
+    };
+
+    class IsobaricSurface final {
+    public:
+        using PhysicalQuantity = dimensions::Pressure;
+        using TargetValue = dimensions::PressureValue;
+
+    public:
+        struct Parameters final {
+            dimensions::PressureValue PSurf;
+        };
+
+    public:
+        explicit IsobaricSurface(const Parameters parameters);
+        explicit IsobaricSurface (const TargetValue targetPressure);
+
+        [[nodiscard]] TargetValue target_pressure() const noexcept;
+        [[nodiscard]] std::array<double, 1> continuation_identity() const noexcept;
+        [[nodiscard]] PressureSurfaceDescriptor descriptor() const noexcept;
+
+    private:
+        TargetValue m_targetPressure;
+    };
+
+    template <typename Candidate>
+    concept IsobaricSurfaceType = std::same_as<std::remove_cvref_t<Candidate>, IsobaricSurface>;
+
+    using Isobaric = IsobaricSurface;
+}
\ No newline at end of file
diff --git a/src/include/serif/geometry/surface/linearization.hpp b/src/include/serif/geometry/surface/linearization.hpp
new file mode 100644
index 00000000..db6d6c8c
--- /dev/null
+++ b/src/include/serif/geometry/surface/linearization.hpp
@@ -0,0 +1,41 @@
+#pragma once
+
+#include <cstddef>
+#include <limits>
+#include <string_view>
+
+#include "serif/equations/point.hpp"
+#include "serif/geometry/surface/compilation.hpp"
+
+namespace serif::geometry::surface {
+    enum class SurfaceFailureReason {
+        non_finite_state,
+        non_finite_eos_response,
+        non_finite_linearization,
+        inconsistent_sample_count,
+        empty_surface
+    };
+
+    /** @brief A rejected thermodynamic candidate, with a sample index when evaluating a boundary batch. */
+    struct SurfaceFailure {
+        SurfaceFailureReason reason;
+        std::string_view message;
+        std::size_t sample = std::numeric_limits<std::size_t>::max();
+    };
+
+    template <PressureSurfaceFormulationType Formulation>
+    using SurfacePoint = equations::EquationPoint<typename Formulation::Operands>;
+
+    template <PressureSurfaceFormulationType Formulation>
+    using SurfaceLinearization = equations::PointLinearization<typename Formulation::Operands>;
+
+    /**
+     * @brief An EOS value and all of its local derivatives, evaluated together at one state.
+     * @note Derivatives sharing a term are accumulated; the later residual calculation supplies the sign and carrier identity.
+     */
+    template <PressureSurfaceFormulationType Formulation>
+    struct SurfaceEOSResponse {
+        dimensions::QuantityValue<typename Formulation::Relation::OutputQuantity> value{0.0};
+        equations::OperandValues<typename Formulation::Operands, double> field_derivatives{};
+    };
+}
diff --git a/src/include/serif/geometry/surface/samples.hpp b/src/include/serif/geometry/surface/samples.hpp
new file mode 100644
index 00000000..fc1ab3ee
--- /dev/null
+++ b/src/include/serif/geometry/surface/samples.hpp
@@ -0,0 +1,70 @@
+#pragma once
+
+#include <utility>
+#include <vector>
+
+#include "serif/geometry/surface/linearization.hpp"
+#include "serif/utils/dependencies/resource.hpp"
+
+namespace serif::geometry::surface {
+    namespace keys {
+        struct EquationOfState {
+        };
+        struct PressureCondition {
+        };
+        struct EOSResponse {
+        };
+        struct Linearization {
+        };
+    }
+
+    /**
+     * @brief Values of one scalar term at an ordered collection of boundary points.
+     * @note Every term in a preparation must use the same point ordering. Sampling and finite-element trace choices belong to the backend.
+     */
+    template <SurfaceScalarTerm Term>
+    struct BoundaryValues {
+        std::vector<double> values;
+
+        /// Compare the actual sampled values, including their count, without a lossy checksum.
+        [[nodiscard]] bool operator==(const BoundaryValues&) const = default;
+    };
+
+    namespace detail {
+        struct BoundaryValueTracking : utils::dependencies::ValueSnapshot {
+            /// Numerical samples have value semantics; rebuilding an equal context need not reevaluate the physics.
+            template <SurfaceScalarTerm Term>
+            [[nodiscard]] const void* identity(const BoundaryValues<Term>&) const {
+                return nullptr;
+            }
+        };
+
+        struct PressureConditionTracking {
+            /// Only prescribed pressure can affect the isobaric calculation.
+            [[nodiscard]] dimensions::PressureValue capture(const IsobaricSurface& condition) const {
+                return condition.target_pressure();
+            }
+
+            /// Compare pressure exactly, retaining the dependency engine's normal change semantics.
+            [[nodiscard]] bool matches(const dimensions::PressureValue captured, const IsobaricSurface& condition) const {
+                return captured == condition.target_pressure();
+            }
+
+            /// An equivalent pressure condition may be supplied by a newly constructed context.
+            [[nodiscard]] const void* identity(const IsobaricSurface&) const {
+                return nullptr;
+            }
+        };
+    }
+
+    /** @brief Supply owned boundary samples to an ordinary dependency context, keyed by their physical term. */
+    template <SurfaceScalarTerm Term>
+    [[nodiscard]] auto boundary_values(std::vector<double> values) {
+        return utils::dependencies::own<Term>(BoundaryValues<Term>{.values = std::move(values)}, detail::BoundaryValueTracking{});
+    }
+
+    /** @brief Supply the isobaric condition as an owned dependency source, with pressure tracked by value. */
+    [[nodiscard]] inline auto pressure_condition(IsobaricSurface condition) {
+        return utils::dependencies::own<keys::PressureCondition>(std::move(condition), detail::PressureConditionTracking{});
+    }
+}
diff --git a/src/include/serif/geometry/surface/surface.hpp b/src/include/serif/geometry/surface/surface.hpp
new file mode 100644
index 00000000..ebb1b5e5
--- /dev/null
+++ b/src/include/serif/geometry/surface/surface.hpp
@@ -0,0 +1,17 @@
+#pragma once
+
+/**
+ * @file surface.hpp
+ * @brief Isobaric conditions, EOS-selected formulations, and reusable boundary calculations.
+ *
+ * IsobaricSurface describes the prescribed pressure. PressureSurfaceFormulation selects one compatible EOS relation and its carrier;
+ * DirectPressureSurfaceFormulation explicitly chooses a forward pressure residual and its equation row. Neither selects a finite-element trace.
+ *
+ * Use evaluate_pressure_surface for one sampled point, or make_pressure_surface_graph for an ordered batch. The latter returns an ordinary
+ * dependency graph: the caller supplies the EOS with its tracking policy, the pressure condition, and typed boundary values. Results use the
+ * same PointLinearization representation as other equations. Finite-element row enforcement remains a separate numerical capability.
+ */
+
+#include "serif/geometry/surface/compilation.hpp"
+#include "serif/geometry/surface/computations.hpp"
+#include "serif/geometry/surface/evaluation.hpp"
diff --git a/src/lib/serif/geometry/surface/isobaric.cpp b/src/lib/serif/geometry/surface/isobaric.cpp
index e69de29b..fed8ac4a 100644
--- a/src/lib/serif/geometry/surface/isobaric.cpp
+++ b/src/lib/serif/geometry/surface/isobaric.cpp
@@ -0,0 +1,31 @@
+#include <cmath>
+#include <stdexcept>
+#include <array>
+
+#include "serif/geometry/surface/isobaric.hpp"
+#include "serif/utils/error/errors.hpp"
+
+namespace serif::geometry::surface {
+    using namespace utils::errors;
+
+    IsobaricSurface::IsobaricSurface(const Parameters parameters) : IsobaricSurface(parameters.PSurf) {}
+
+    IsobaricSurface::IsobaricSurface(const TargetValue targetPressure) : m_targetPressure(targetPressure) {
+        serif_ensure<SERiFErrorCode::NON_FINITE_VALUE_ERROR, std::invalid_argument>(std::isfinite(targetPressure.value()), "The target surface pressure must be finite, got {}", targetPressure.value());
+        serif_ensure<SERiFErrorCode::NEGATIVE_VALUE_ERROR, std::invalid_argument>(targetPressure.value() >= 0, "The target surface pressure must be non-negative, got {}", targetPressure.value());
+    }
+
+    IsobaricSurface::TargetValue IsobaricSurface::target_pressure() const noexcept {
+        return m_targetPressure;
+    }
+
+    std::array<double, 1> IsobaricSurface::continuation_identity() const noexcept {
+        return {m_targetPressure.value()};
+    }
+
+    PressureSurfaceDescriptor IsobaricSurface::descriptor() const noexcept {
+        return PressureSurfaceDescriptor{
+            .target_pressure = m_targetPressure.value()
+        };
+    }
+}