diff --git a/chainladder/adjustments/bootstrap.py b/chainladder/adjustments/bootstrap.py
index a25929712..ae4921dda 100644
--- a/chainladder/adjustments/bootstrap.py
+++ b/chainladder/adjustments/bootstrap.py
@@ -271,20 +271,13 @@ def _get_simulation(self, X, exp_incr_triangle):
         resampled_triangles = (resampled_residual * xp.sqrt(abs(b)) + b).cumsum(2)
         resampled_triangles = resampled_triangles[None, ...].swapaxes(0, 1)
         obj = X.copy()
+        obj.values = resampled_triangles
         if X.key_labels == ["Total"]:
-            obj.kdims = np.arange(self.n_sims)
-            obj.key_labels = ["Simulation_#"]
+            obj.index = pd.DataFrame({"Simulation_#": np.arange(self.n_sims)})
         else:
-            obj.kdims = np.concat(
-                [
-                    np.tile(X.kdims, (self.n_sims, 1)),
-                    np.arange(self.n_sims).reshape(-1, 1),
-                ],
-                axis=1,
-            )
-            obj.key_labels = X.key_labels + ["Simulation_#"]
-        obj.values = resampled_triangles
-        obj._set_slicers()
+            obj_idx = pd.concat([X.index] * self.n_sims, ignore_index=True)
+            obj_idx["Simulation_#"] = np.arange(self.n_sims)
+            obj.index = obj_idx
         return obj, scale_phi
 
     def _get_design_matrix(self, X):
diff --git a/chainladder/adjustments/tests/test_bootstrap.py b/chainladder/adjustments/tests/test_bootstrap.py
index a0d683e59..75ba06ccb 100644
--- a/chainladder/adjustments/tests/test_bootstrap.py
+++ b/chainladder/adjustments/tests/test_bootstrap.py
@@ -24,13 +24,6 @@ def test_bs_multiple_cols():
 def test_multi_index(clrd):
     tri = clrd["CumPaidLoss"].sum()
     resampled_triangles = cl.BootstrapODPSample().fit(tri).resampled_triangles_
-    assert np.all(
-        resampled_triangles.index
-        == pd.DataFrame(
-            np.concat(
-                (np.array([["(All)", "(All)"]] * 1000), np.arange(1000).reshape(-1, 1)),
-                axis=1,
-            ),
-            columns=["GRNAME", "LOB", "Simulation_#"],
-        )
-    )
+    expected = pd.concat([tri.index] * 1000, ignore_index=True)
+    expected["Simulation_#"] = np.arange(1000)
+    pd.testing.assert_frame_equal(resampled_triangles.index, expected)
diff --git a/chainladder/core/axis.py b/chainladder/core/axis.py
new file mode 100644
index 000000000..86ac3c434
--- /dev/null
+++ b/chainladder/core/axis.py
@@ -0,0 +1,164 @@
+"""Axis descriptor for Triangle dimensions, analogous to Pandas AxisProperty."""
+
+from __future__ import annotations
+
+from typing import TYPE_CHECKING, Any, Callable
+import pandas as pd
+
+from chainladder.utils.utility_functions import _get_axis_name, _get_axis_number
+
+if TYPE_CHECKING:
+    from chainladder.core.triangle import Triangle
+
+
+class TriangleAxis:
+    """
+    Generalized descriptor for representing a Triangle dimension,
+    analogous to Pandas AxisProperty.
+
+    Parameters
+    ----------
+    axis : int or str
+        Axis number (0..3) or canonical name ('index', 'columns', 'origin',
+        'development') used to store and access raw axis data in ``obj._axes``.
+    fget : callable, optional
+        Transformation function taking ``(obj, raw)`` and returning the
+        public axis representation. If None, the raw value is returned.
+    fset : callable, optional
+        Transformation function taking ``(obj, value)`` and returning the
+        raw axis representation to be stored in ``obj._axes``. If None,
+        the value is stored as-is.
+    doc : str, optional
+        Docstring for the descriptor attribute on the class.
+    """
+
+    def __init__(
+        self,
+        axis: int | str,
+        *,
+        fget: Callable[[Triangle, Any], Any] | None = None,
+        fset: Callable[[Triangle, Any], Any] | None = None,
+        doc: str | None = None,
+    ) -> None:
+        """
+        Initialize a TriangleAxis descriptor.
+
+        Parameters
+        ----------
+        axis : int or str
+            Axis number or name mapped to an integer in ``obj._axes``.
+        fget : callable, optional
+            Getter callable mapping ``(obj, raw)`` to public representation.
+        fset : callable, optional
+            Setter callable mapping ``(obj, value)`` to internal representation.
+        doc : str, optional
+            Docstring for the attribute.
+        """
+        self.axis: int = _get_axis_number(axis)
+        self.key: int = self.axis
+        self.fget = fget  # raw -> public; identity if None
+        self.fset = fset  # (obj, public) -> raw; identity if None
+        self.__doc__ = doc
+
+    def __get__(self, obj: Triangle | None, objtype: type | None = None) -> Any:
+        """
+        Access the axis value on a Triangle instance or return descriptor on class.
+
+        Parameters
+        ----------
+        obj : Triangle or None
+            Instance of Triangle or None if accessed via class.
+        objtype : type, optional
+            Type of the owning class.
+
+        Returns
+        -------
+        Any
+            The descriptor itself if accessed on class, or the transformed
+            axis value if accessed on an instance.
+
+        Raises
+        ------
+        AttributeError
+            If the axis key has not been initialized in ``obj._axes``.
+        """
+        if obj is None:
+            return self
+        if not hasattr(obj, "_axes") or self.axis not in obj._axes:
+            axis_name = _get_axis_name(self.axis)
+            raise AttributeError(
+                f"'{type(obj).__name__}' object has no attribute '{axis_name}'"
+            )
+        raw = obj._axes[self.axis]
+        return self.fget(obj, raw) if self.fget else raw
+
+    def __set__(self, obj: Triangle, value: Any) -> None:
+        """
+        Set the axis value on a Triangle instance.
+
+        Transforms ``value`` using ``fset`` if defined, stores the raw value
+        in ``obj._axes[self.axis]``, and updates slicers if necessary.
+
+        Parameters
+        ----------
+        obj : Triangle
+            Instance of Triangle.
+        value : Any
+            New axis value provided by the user.
+        """
+        raw = self.fset(obj, value) if self.fset else value
+        if not hasattr(obj, "_axes"):
+            obj._axes = {}
+        obj._axes[self.axis] = raw
+        if hasattr(obj, "virtual_columns"):
+            obj._set_slicers()
+
+
+def _set_columns(obj: Triangle, value: Any) -> pd.Index:
+    """
+    Validate and transform assigned columns into a pandas Index.
+
+    Parameters
+    ----------
+    obj : Triangle
+        The Triangle instance being updated.
+    value : Any
+        The new column labels. Can be a string, list, or pandas Index.
+
+    Returns
+    -------
+    pd.Index
+        A pandas Index named ``"columns"``.
+    """
+    if isinstance(value, str):
+        value = [value]
+    if hasattr(obj, "values") and obj.values is not None:
+        obj._len_check(range(obj.values.shape[1]), value)
+    elif hasattr(obj, "_axes") and (1 in obj._axes or "columns" in obj._axes):
+        obj._len_check(obj.columns, value)
+    return pd.Index(value, name="columns")
+
+
+def _set_index(obj: Triangle, value: Any) -> pd.DataFrame:
+    """
+    Validate and transform assigned index into a pandas DataFrame.
+
+    Parameters
+    ----------
+    obj : Triangle
+        The Triangle instance being updated.
+    value : Any
+        The new index labels. Must be a pandas DataFrame.
+
+    Returns
+    -------
+    pd.DataFrame
+        A reset-index pandas DataFrame.
+    """
+    if not isinstance(value, pd.DataFrame):
+        raise TypeError("index must be a pandas DataFrame")
+    if hasattr(obj, "values") and obj.values is not None:
+        obj._len_check(range(obj.values.shape[0]), value)
+    elif hasattr(obj, "_axes") and (0 in obj._axes or "index" in obj._axes):
+        obj._len_check(obj.index, value)
+    return value.copy().reset_index(drop=True)
diff --git a/chainladder/core/base.py b/chainladder/core/base.py
index d15d6cd16..4f30769a1 100644
--- a/chainladder/core/base.py
+++ b/chainladder/core/base.py
@@ -282,7 +282,9 @@ def _aggregate_data(
         return data_agg
 
     @staticmethod
-    def _set_kdims(data_agg: DataFrame, index: list) -> tuple[np.ndarray, np.ndarray]:
+    def _factorize_index(
+        data_agg: DataFrame, index: list
+    ) -> tuple[np.ndarray, np.ndarray]:
         """
         Sets the key dimension of the triangle.
 
@@ -302,15 +304,18 @@ def _set_kdims(data_agg: DataFrame, index: list) -> tuple[np.ndarray, np.ndarray
         """
 
         # Get unique values of index and assign an integer value to each one.
-        kdims: DataFrame = (
+        unique_keys: DataFrame = (
             data_agg[index].drop_duplicates().reset_index(drop=True).reset_index()
         )
 
         # Map these integers back to the agg data to generate a key index.
-        key_idx: np.ndarray = (
-            data_agg[index].merge(kdims, how="left", on=index)["index"].values[None].T
+        index_idx: np.ndarray = (
+            data_agg[index]
+            .merge(unique_keys, how="left", on=index)["index"]
+            .values[None]
+            .T
         )
-        return kdims.drop(labels="index", axis=1).values, key_idx
+        return unique_keys.drop(labels="index", axis=1).values, index_idx
 
     @staticmethod
     def _set_odims(
@@ -353,7 +358,7 @@ def _set_ddims(
     @staticmethod
     def _set_values(
         data_agg: DataFrame,
-        key_idx: np.ndarray,
+        index_idx: np.ndarray,
         columns: list,
         orig_idx: np.ndarray,
         dev_idx: np.ndarray,
@@ -368,7 +373,7 @@ def _set_values(
             tuple([np.concatenate((orig_idx, dev_idx), 1)] * len(columns)), 0
         )
         coords: np.ndarray = np.concatenate(
-            (np.concatenate(tuple([key_idx] * len(columns)), 0), val_idx, coords), 1
+            (np.concatenate(tuple([index_idx] * len(columns)), 0), val_idx, coords), 1
         )
         amts: np.ndarray = np.concatenate([
             data_agg[col].fillna(0).values for col in data_agg[columns]
diff --git a/chainladder/core/common.py b/chainladder/core/common.py
index 9c3d44761..8b73648b8 100644
--- a/chainladder/core/common.py
+++ b/chainladder/core/common.py
@@ -198,7 +198,7 @@ def ibnr_(self):
         else:
             ld = self.latest_diagonal if self.is_cumulative else self.sum(axis=3)
         ibnr = self.ultimate_ - ld
-        ibnr.vdims = self.ultimate_.vdims
+        ibnr.columns = self.ultimate_.columns
         return ibnr
 
     @property
diff --git a/chainladder/core/correlation.py b/chainladder/core/correlation.py
index dfa88981b..ed832c509 100644
--- a/chainladder/core/correlation.py
+++ b/chainladder/core/correlation.py
@@ -173,7 +173,7 @@ def __init__(self, triangle, p_critical: float = 0.5):
 
         # final big T
         self.t_expectation = pd.DataFrame(
-            t_expectation[..., 0, 0], columns=triangle.vdims, index=idx
+            t_expectation[..., 0, 0], columns=triangle.columns, index=idx
         )
 
         # table of Spearman's rank coefficients Tk, can be used to verify consistency with paper
@@ -350,17 +350,17 @@ def p_z_lower(z: int, n: int, p: float = 0.5) -> float:
                     (self.range[0] > VarZ.sum(axis=-1))
                     | (VarZ.sum(axis=-1) > self.range[1])
                 )[..., 0],
-                columns=triangle.vdims,
+                columns=triangle.columns,
                 index=idx,
             )
             self.z = pd.DataFrame(
-                z.sum(axis=-1)[..., 0], columns=triangle.vdims, index=idx
+                z.sum(axis=-1)[..., 0], columns=triangle.columns, index=idx
             )
             self.z_expectation = pd.DataFrame(
-                EZ.sum(axis=-1)[..., 0], columns=triangle.vdims, index=idx
+                EZ.sum(axis=-1)[..., 0], columns=triangle.columns, index=idx
             )
             self.z_variance = pd.DataFrame(
-                VarZ.sum(axis=-1)[..., 0], columns=triangle.vdims, index=idx
+                VarZ.sum(axis=-1)[..., 0], columns=triangle.columns, index=idx
             )
 
 
diff --git a/chainladder/core/display.py b/chainladder/core/display.py
index e2015b5f5..1c6b5bc13 100644
--- a/chainladder/core/display.py
+++ b/chainladder/core/display.py
@@ -52,7 +52,7 @@ def _summary_frame(self) -> DataFrame:
                 "O" + self.origin_grain + "D" + self.development_grain,
                 self.shape,
                 self.key_labels,
-                list(self.vdims),
+                list(self.columns),
             ],
             index=["Valuation:", "Grain:", "Shape:", "Index:", "Columns:"],
             name="Triangle Summary",
diff --git a/chainladder/core/dunders.py b/chainladder/core/dunders.py
index e340366e1..da39ae1a6 100644
--- a/chainladder/core/dunders.py
+++ b/chainladder/core/dunders.py
@@ -72,21 +72,28 @@ def _compatibility_check(self, x, y):
         return x, y
 
     def _prep_index(self, x, y):
-        if x.kdims.shape[0] == 1 and y.kdims.shape[0] > 1:
-            x.kdims = y.kdims
-            x.key_labels = y.key_labels
+        if len(x.index) == 1 and len(y.index) > 1:
+            x._axes[0] = y.index.copy()
+            x.key_labels = list(y.key_labels)
+            x._set_slicers()
             return x, y
-        if x.kdims.shape[0] > 1 and y.kdims.shape[0] == 1:
-            y.kdims = x.kdims
-            y.key_labels = x.key_labels
+        if len(x.index) > 1 and len(y.index) == 1:
+            y._axes[0] = x.index.copy()
+            y.key_labels = list(x.key_labels)
+            y._set_slicers()
             return x, y
-        if x.kdims.shape[0] == y.kdims.shape[0] == 1 and x.key_labels != y.key_labels:
-            kdims = x.kdims if len(x.key_labels) > len(y.key_labels) else y.kdims
-            key_labels = (
-                x.key_labels if len(x.key_labels) > len(y.key_labels) else y.key_labels
+        if len(x.index) == len(y.index) == 1 and x.key_labels != y.key_labels:
+            index = (
+                x.index.copy()
+                if len(x.key_labels) > len(y.key_labels)
+                else y.index.copy()
             )
-            x.kdims = y.kdims = kdims
-            x.key_labels = y.key_labels = key_labels
+            x._axes[0] = index.copy()
+            x.key_labels = list(index.columns)
+            x._set_slicers()
+            y._axes[0] = index.copy()
+            y.key_labels = list(index.columns)
+            y._set_slicers()
             return x, y
 
         # Use sets for faster operations
@@ -95,14 +102,10 @@ def _prep_index(self, x, y):
         common = x_labels.intersection(y_labels)
 
         if common == x_labels or common == y_labels:
-            if x_labels != y_labels or x.kdims.shape[0] != y.kdims.shape[0]:
+            if x_labels != y_labels or len(x.index) != len(y.index):
                 x = x.groupby(list(common))
                 y = y.groupby(list(common))
-            elif (
-                x.kdims.shape[0] > 1
-                and not np.array_equal(x.kdims, y.kdims)
-                and not x.index.equals(y.index)
-            ):
+            elif len(x.index) > 1 and not x.index.equals(y.index):
                 x = x.sort_index()
                 try:
                     y = y.loc[x.index]
@@ -124,12 +127,14 @@ def _prep_index(self, x, y):
 
     def _prep_columns(self, x, y):
         if len(x.columns) == 1 and len(y.columns) > 1:
-            x.vdims = y.vdims
+            x._axes[1] = y.columns
+            x._set_slicers()
         elif len(y.columns) == 1 and len(x.columns) > 1:
-            y.vdims = x.vdims
-        elif len(y.columns) == len(x.columns) == 1 and x.columns != y.columns:
-            y.vdims = x.vdims
-        elif x.shape[1] == y.shape[1] and np.array_equal(x.columns, y.columns):
+            y._axes[1] = x.columns
+            y._set_slicers()
+        elif len(y.columns) == len(x.columns) == 1 and not x.columns.equals(y.columns):
+            y.columns = x.columns
+        elif x.shape[1] == y.shape[1] and x.columns.equals(y.columns):
             return x, y
         else:
             # Find columns to add to each triangle
@@ -260,8 +265,7 @@ def _slice_or_nan(obj, other, k):
                 .reset_index()
             )
             new_idx = new_idx[new_obj.key_labels].iloc[-1:]
-            new_obj.kdims = new_idx.values
-            new_obj.key_labels = list(new_idx.columns)
+            new_obj.index = new_idx
             return new_obj
 
     @staticmethod
diff --git a/chainladder/core/io.py b/chainladder/core/io.py
index c41d7f16b..26a0b0c94 100644
--- a/chainladder/core/io.py
+++ b/chainladder/core/io.py
@@ -104,7 +104,11 @@ def to_json(self):
         json_dict["sub_tris"] = {
             sub_tri: getattr(self, sub_tri).to_json() for sub_tri in sub_tris
         }
-        dfs = [k for k, v in vars(self).items() if isinstance(v, pd.DataFrame)]
+        dfs = [
+            k
+            for k, v in vars(self).items()
+            if isinstance(v, pd.DataFrame) and k != "_index"
+        ]
         json_dict["dfs"] = {df: getattr(self, df).to_json() for df in dfs}
         dfs = [k for k, v in vars(self).items() if isinstance(v, pd.Series)]
         json_dict["dfs"].update({
diff --git a/chainladder/core/pandas.py b/chainladder/core/pandas.py
index 56e300999..061f47390 100644
--- a/chainladder/core/pandas.py
+++ b/chainladder/core/pandas.py
@@ -170,8 +170,8 @@ def to_frame(
             elif len(axes) in [1, 2]:
                 tri: ndarray = np.squeeze(self.set_backend("numpy").values)
                 axes_lookup: dict = {
-                    0: self.kdims,
-                    1: self.vdims,
+                    0: self.index.values,
+                    1: self.columns.values,
                     2: self.origin,
                     3: self.development,
                 }
@@ -1141,10 +1141,17 @@ def sort_index(self, *args, **kwargs) -> Triangle:
         -------
         Triangle
         """
+        inplace = kwargs.pop("inplace", False)
         sorted_index: DataFrame = cast(
-            "DataFrame", self.index.sort_values(self.key_labels, *args, **kwargs)
+            "DataFrame",
+            self.index.copy().sort_values(self.key_labels, *args, **kwargs),
         )
-        return self.iloc[sorted_index.index]
+        result = self.iloc[sorted_index.index]
+        if inplace:
+            self.values = result.values
+            self.index = result.index
+            return self
+        return result
 
     def exp(self) -> Triangle:
         """
@@ -1300,10 +1307,12 @@ def agg_func(
 
         # Aggregation function will collapse a dimension, so
         # adjust the dimensions of the original object to match that of the aggregation.
-        if axis == 0 and obj.values.shape[axis] == 1 and len(obj.kdims) > 1:
-            obj.kdims = np.array([["(All)"] * len(obj.key_labels)])
-        if axis == 1 and obj.values.shape[axis] == 1 and len(obj.vdims) > 1:
-            obj.vdims = np.array([0])
+        if axis == 0 and obj.values.shape[axis] == 1 and len(obj.index) > 1:
+            obj.index = pd.DataFrame(
+                [["(All)"] * len(obj.key_labels)], columns=obj.key_labels
+            )
+        if axis == 1 and obj.values.shape[axis] == 1 and len(obj.columns) > 1:
+            obj.columns = pd.Index([0], name="columns")
         if axis == 2 and obj.values.shape[axis] == 1 and len(obj.odims) > 1:
             obj.odims = obj.odims[0:1]
         # If axis is development, set the ddims to be the valuation date.
@@ -1366,10 +1375,9 @@ def aggregate(i, obj, axis, v):
                 obj.index = index
             else:
                 index = pd.DataFrame(group_index)
-                obj.key_labels = index.columns.tolist()
-                obj.kdims = index.values
+                obj.index = index
         if self.axis == 1:
-            obj.vdims = pd.DataFrame(group_index).values[:, 0]
+            obj.columns = group_index
         if self.axis == 2:
             odims = self.obj._to_datetime(
                 pd.Series(self.groups.indices.keys()).to_frame(), [0]
diff --git a/chainladder/core/slice.py b/chainladder/core/slice.py
index 1b464205b..440870e2c 100644
--- a/chainladder/core/slice.py
+++ b/chainladder/core/slice.py
@@ -9,6 +9,7 @@
 
 import numpy as np
 import pandas as pd
+from pandas.core.indexing import check_bool_indexer
 
 from chainladder.core.typing import _AxisKey, _LabelKey, TriangleProtocol
 
@@ -63,8 +64,8 @@ def get_idx(self, idx: tuple[_AxisKey, _AxisKey, _AxisKey, _AxisKey]) -> Triangl
         else:
             obj.values = obj.values[i_idx, :, o_idx, d_idx][:, c_idx, ...]
         # Set the new dimension values.
-        obj.kdims = obj.kdims[i_idx]
-        obj.vdims = obj.vdims[c_idx]
+        obj.index = obj.index.iloc[i_idx].reset_index(drop=True)
+        obj.columns = obj.columns[c_idx]
         obj.odims, obj.ddims = obj.odims[o_idx], obj.ddims[d_idx]
         obj.valuation_date = cast(
             pd.Timestamp, np.minimum(obj.valuation.max(), obj.valuation_date)
@@ -383,7 +384,7 @@ def index_key(self, key: _LabelKey) -> np.ndarray:
             ).values.flatten()
         # Case scalar, locate position in first level of index.
         else:
-            idx = np.where(self.obj.kdims[:, 0] == key)[0]
+            idx = np.where(self.obj.index.iloc[:, 0].values == key)[0]
         return idx
 
     def other_key(
@@ -797,6 +798,13 @@ def __getitem__(
             return self._slice(key, "odims")
         # Case index.
         if isinstance(key, pd.Series):
+            if pd.api.types.is_bool_dtype(key):
+                if len(key) != len(self.index):
+                    raise IndexError(
+                        f"Boolean index has wrong length: {len(key)} instead of {len(self.index)}"
+                    )
+                mask = check_bool_indexer(self.index.index, key)
+                return self.iloc[np.where(mask)[0]]
             return self.iloc[self.index[key].index]
         elif key in self.key_labels:
             return self.index[key]
@@ -811,7 +819,7 @@ def __getitem__(
                 [key] if isinstance(key, (str, int, float, np.generic)) else key
             )
             # Identify the position of each requested element within the valuation dimension.
-            idx = [list(self.vdims).index(item) for item in keys]
+            idx = [list(self.columns).index(item) for item in keys]
             return self.iloc[:, idx]
 
     def __setitem__(
@@ -825,7 +833,7 @@ def __setitem__(
         Parameters
         ----------
         key: str | int
-            The vdims label of the column to set.
+            The column label of the column to set.
         value: int | float | TriangleSlicer | Callable[[Triangle], TriangleSlicer]
             The value(s) to assign to the column. A callable defines a virtual
             (lazily-computed) column.
@@ -873,23 +881,32 @@ def __setitem__(
         if callable(value):
             self.virtual_columns[key] = value
             if self.array_backend == "sparse":
-                if key not in self.vdims:
+                if key not in self.columns:
                     k, v, o, d = self.values.shape
                     self.values.shape = k, v + 1, o, d
-                    self.vdims = np.append(self.vdims, key)
+                    self.columns = self.columns.append(pd.Index([key]))
                 return
             # Create a placeholder column.
             value = (self.iloc[:, 0].copy() * xp.nan).set_backend(self.array_backend)
         # Otherwise, remove column from virtual columns if previously reserved for lazy evaluation.
         else:
             self.virtual_columns.pop(key)
+        if isinstance(value, (int, float, np.number)):
+            # Broadcast scalar across the Triangle's shape.
+            value = self.iloc[:, 0] * 0 + value
         if isinstance(value, TriangleSlicer):
             value = cast("TriangleProtocol", cast(object, value))
             if value.array_backend != self.array_backend:
                 value = value.set_backend(self.array_backend)
         # Key exists in columns, replace data.
-        if key in self.vdims:
-            i = np.where(self.vdims == key)[0][0]
+        if key in self.columns:
+            loc = self.columns.get_loc(key)
+            if isinstance(loc, slice):
+                i = loc.start or 0
+            elif isinstance(loc, np.ndarray):
+                i = int(np.where(loc)[0][0])
+            else:
+                i = int(loc)
             # Case sparse backend.
             if self.array_backend == "sparse":
                 # Unwrap a Triangle-valued assignment to its raw array. A raw
@@ -927,10 +944,6 @@ def __setitem__(
                 cast(np.ndarray, self.values)[:, i : i + 1] = value
         # Key is new, create a column and update data.
         else:
-            self.vdims = np.append(self.vdims, key)
-            if isinstance(value, (int, float, np.number)):
-                # Broadcast scalar across the Triangle's shape.
-                value = self.iloc[:, 0] * 0 + value
             try:
                 self.values = xp.concatenate((self.values, value.values), axis=1)
             except (ValueError, AttributeError, AssertionError):
@@ -940,6 +953,7 @@ def __setitem__(
                     (self.iloc[:, 0] * 0 + cast("Triangle", value)).values,
                 )
                 self.values = xp.concatenate(conc, axis=1)
+            self.columns = self.columns.append(pd.Index([key]))
 
     def _slice_valuation(self: TriangleProtocol, key: np.ndarray) -> Triangle:
         """
diff --git a/chainladder/core/tests/test_triangle.py b/chainladder/core/tests/test_triangle.py
index cd3515747..542579f39 100644
--- a/chainladder/core/tests/test_triangle.py
+++ b/chainladder/core/tests/test_triangle.py
@@ -10,7 +10,12 @@
 import warnings
 
 from chainladder.core.common import Common
-from chainladder.utils.utility_functions import date_delta_adjustment
+from chainladder.utils.utility_functions import (
+    _axis_orders,
+    _get_axis_name,
+    _get_axis_number,
+    date_delta_adjustment,
+)
 from chainladder.utils.sparse import COO
 
 from io import StringIO
@@ -174,8 +179,7 @@ def test_sum_of_diff_eq_diff_of_sum(clrd):
 
 def test_append(raa):
     raa2 = raa.copy()
-    raa2.kdims = np.array([["P2"]])
-    raa.append(raa2).sum() == raa * 2
+    raa2.index = pd.DataFrame([["P2"]], columns=raa.key_labels)
     assert raa.append(raa2).sum() == 2 * raa
 
 
@@ -328,7 +332,7 @@ def test_origin_and_value_setters(raa):
         np.all(raa2.origin == raa.origin),
         np.all(raa2.development == raa.development),
         np.all(raa2.odims == raa.odims),
-        np.all(raa2.vdims == raa.vdims),
+        np.all(raa2.columns == raa.columns),
     ))
 
 
@@ -353,9 +357,9 @@ def test_index_setter_with_dataframe(clrd: Triangle) -> None:
     tri.index = new_index
 
     assert tri.key_labels == ["Company"]
-    np.testing.assert_array_equal(tri.kdims, new_index.values)
+    assert tri.index.equals(new_index)
     # _set_slicers() must have rebuilt .loc against the new key label.
-    assert tri.loc["A"].kdims.tolist() == [["A"]]
+    assert tri.loc["A"].index.values.tolist() == [["A"]]
     assert tri.loc["A"] == clrd.iloc[:1]
 
 
@@ -401,8 +405,8 @@ def test_index_setter_non_dataframe_raises(clrd: Triangle) -> None:
 
 def test_set_index_inplace(clrd: Triangle) -> None:
     """
-    Triangle.set_index(value, inplace=True) should mutate the calling
-    Triangle's index via the index setter and return that same object.
+    Triangle.set_index(value, inplace=True) should mutate the Triangle and
+    return self.
 
     Parameters
     ----------
@@ -420,14 +424,14 @@ def test_set_index_inplace(clrd: Triangle) -> None:
 
     assert result is tri
     assert tri.key_labels == ["Company"]
-    np.testing.assert_array_equal(tri.kdims, new_index.values)
+    assert tri.index.equals(new_index)
 
 
 def test_set_index_not_inplace(clrd: Triangle) -> None:
     """
     Triangle.set_index(value) with the default inplace=False should operate
     on a copy: it returns a distinct Triangle with the new index applied,
-    leaving the original Triangle's kdims/key_labels untouched.
+    leaving the original Triangle's index/key_labels untouched.
 
     Parameters
     ----------
@@ -439,7 +443,7 @@ def test_set_index_not_inplace(clrd: Triangle) -> None:
     None
     """
     tri = clrd.iloc[:3]
-    original_kdims = tri.kdims.copy()
+    original_index = tri.index.copy()
     original_key_labels = list(tri.key_labels)
     new_index = pd.DataFrame({"Company": ["A", "B", "C"]})
 
@@ -447,9 +451,370 @@ def test_set_index_not_inplace(clrd: Triangle) -> None:
 
     assert result is not tri
     assert result.key_labels == ["Company"]
-    np.testing.assert_array_equal(result.kdims, new_index.values)
+    assert result.index.equals(new_index)
     assert tri.key_labels == original_key_labels
-    np.testing.assert_array_equal(tri.kdims, original_kdims)
+    assert tri.index.equals(original_index)
+
+
+def test_triangle_columns_setter(raa):
+    """Setting Triangle.columns should update columns and slicers, and validate lengths."""
+    tri = raa.copy()
+    tri.columns = ["NewColumn"]
+    assert list(tri.columns) == ["NewColumn"]
+    assert tri["NewColumn"].shape == tri.shape
+
+    tri.columns = "SingleColumn"
+    assert list(tri.columns) == ["SingleColumn"]
+    assert tri["SingleColumn"].shape == tri.shape
+
+    with pytest.raises(ValueError, match="Length mismatch"):
+        tri.columns = ["Col1", "Col2"]
+
+
+def test_triangle_axes(raa):
+    """Triangle.axes should return a list of [index, columns, origin, development]."""
+    axes = raa.axes
+    assert isinstance(axes, list)
+    assert len(axes) == 4
+    assert axes[0].equals(raa.index)
+    assert axes[1].equals(raa.columns)
+    assert axes[2].equals(raa.origin)
+    assert axes[3].equals(raa.development)
+
+
+def test_axis_mappings():
+    """Verify axis number and name resolutions in chainladder.core.typing."""
+    for idx, name in enumerate(_axis_orders):
+        assert _get_axis_number(idx) == idx
+        assert _get_axis_number(idx - 4) == idx
+        assert _get_axis_number(name) == idx
+        assert _get_axis_name(idx) == name
+        assert _get_axis_name(idx - 4) == name
+        assert _get_axis_name(name) == name
+
+    with pytest.raises(
+        ValueError, match="No axis named invalid for object type Triangle"
+    ):
+        _get_axis_number("invalid")
+
+    with pytest.raises(ValueError, match=r"No axis named 5 for object type Triangle"):
+        _get_axis_number(5)
+
+
+def test_triangle_axis_descriptor():
+    """Verify TriangleAxis descriptor correctly uses integer axis keys for index and columns."""
+    from chainladder.core.axis import TriangleAxis, _set_columns, _set_index
+
+    assert isinstance(cl.Triangle.index, TriangleAxis)
+    assert cl.Triangle.index.__doc__ is not None
+    assert isinstance(cl.Triangle.columns, TriangleAxis)
+    assert cl.Triangle.columns.__doc__ is not None
+
+    class MockTriangle:
+        def __init__(self):
+            self.values = None
+            self._axes = {}
+
+        def _len_check(self, target, obj):
+            if len(target) != len(obj):
+                raise ValueError("length mismatch")
+
+        index = TriangleAxis(0, fset=_set_index)
+        columns = TriangleAxis(1, fset=_set_columns)
+
+    obj = MockTriangle()
+    with pytest.raises(TypeError, match="index must be a pandas DataFrame"):
+        obj.index = "not_a_df"
+
+    df = pd.DataFrame({"Company": ["A", "B"]})
+    obj.index = df
+    assert 0 in obj._axes
+    pd.testing.assert_frame_equal(obj.index, df)
+
+    # When self.values is None and 0 in _axes, _set_index checks against obj.index
+    df_same_len = pd.DataFrame({"Company": ["C", "D"]})
+    obj.index = df_same_len
+    pd.testing.assert_frame_equal(obj.index, df_same_len)
+
+    with pytest.raises(ValueError, match="length mismatch"):
+        obj.index = pd.DataFrame({"Company": ["A", "B", "C"]})
+
+    # When self.values is not None, _set_index checks against range(self.values.shape[0])
+    obj.values = np.zeros((2, 3))
+    obj.index = df
+    with pytest.raises(ValueError, match="length mismatch"):
+        obj.index = pd.DataFrame({"Company": ["A"]})
+
+    obj.columns = ["a", "b", "c"]
+    assert 1 in obj._axes
+    assert list(obj.columns) == ["a", "b", "c"]
+    assert obj.columns.name == "columns"
+
+    # Custom fget mapping
+    class CustomAxisTriangle:
+        col = TriangleAxis(1, fget=lambda o, raw: [x.upper() for x in raw])
+
+    c_obj = CustomAxisTriangle()
+    c_obj.col = pd.Index(["x", "y"], name="columns")
+    assert c_obj.col == ["X", "Y"]
+
+    # _set_columns with single string and len_check without values
+    class MockColumnsTriangle:
+        def __init__(self):
+            self.values = None
+            self._axes = {}
+
+        def _len_check(self, target, val):
+            if len(target) != len(val):
+                raise ValueError("length mismatch")
+
+        columns = TriangleAxis(1, fset=_set_columns)
+
+    m = MockColumnsTriangle()
+    m.columns = "single_col"
+    assert list(m.columns) == ["single_col"]
+
+    m.columns = ["a"]
+    with pytest.raises(ValueError, match="length mismatch"):
+        m.columns = ["a", "b"]
+
+
+def test_triangle_index_property(raa):
+    """Verify Triangle.index descriptor and slicer reset."""
+    tri = raa.copy()
+    tri.index = tri.index.copy()
+    assert tri.index.equals(raa.index)
+
+    tri.index = pd.DataFrame({"Total": ["A"]})
+    assert list(tri.index["Total"]) == ["A"]
+
+
+def test_prep_index_single_element_different_key_labels(raa):
+    """Verify arithmetic index broadcasting when both triangles have 1 row but differing key_labels."""
+    x = raa.copy()
+    x.key_labels = ["LabelA"]
+    y = raa.copy()
+    y.index = pd.DataFrame([["B", "C"]], columns=["LabelB", "LabelC"])
+
+    z1 = x + y
+    assert z1.key_labels == ["LabelB", "LabelC"]
+    assert 0 in z1._axes
+
+    z2 = y + x
+    assert z2.key_labels == ["LabelB", "LabelC"]
+    assert 0 in z2._axes
+
+
+def test_concat_axis0_index(raa):
+    """Verify cl.concat along axis 0 updates index and key_labels."""
+    out_ignore = cl.concat([raa, raa], axis=0, ignore_index=True)
+    assert len(out_ignore.index) == 2
+    assert out_ignore.key_labels == ["Index"]
+    assert 0 in out_ignore._axes
+
+    r1 = raa.copy()
+    r1.index = pd.DataFrame({"Item": ["A"]})
+    r2 = raa.copy()
+    r2.index = pd.DataFrame({"Item": ["B"]})
+    out_keep = cl.concat([r1, r2], axis=0, ignore_index=False)
+    assert len(out_keep.index) == 2
+    assert list(out_keep.index["Item"]) == ["A", "B"]
+    assert 0 in out_keep._axes
+
+
+def test_key_labels_setter(raa):
+    """Setting key_labels should update Triangle.index columns and slicers."""
+    tri = raa.copy()
+    tri.key_labels = ["NewCompany"]
+    assert tri.key_labels == ["NewCompany"]
+    assert list(tri.index.columns) == ["NewCompany"]
+    indexed = tri.index.set_index(tri.key_labels)
+    assert list(indexed.index.names) == ["NewCompany"]
+
+    tri.key_labels = "SingleCompany"
+    assert tri.key_labels == ["SingleCompany"]
+    assert list(tri.index.columns) == ["SingleCompany"]
+
+
+def test_series_indexing(raa):
+    """Indexing Triangle by boolean Series or label Series should work."""
+    s_bool = pd.Series([True], index=raa.index.index)
+    assert raa[s_bool].shape == raa.shape
+
+    s_label = pd.Series(["Total"])
+    assert raa[s_label].shape == raa.shape
+
+    # Length mismatch raises IndexError
+    s_short = pd.Series([], dtype=bool)
+    with pytest.raises(IndexError, match="Boolean index has wrong length"):
+        _ = raa[s_short]
+
+    s_long = pd.Series([True, False])
+    with pytest.raises(IndexError, match="Boolean index has wrong length"):
+        _ = raa[s_long]
+
+    # Index misalignment raises IndexingError
+    s_misaligned = pd.Series([True], index=[99])
+    with pytest.raises(pd.errors.IndexingError):
+        _ = raa[s_misaligned]
+
+    # Permuted index aligns correctly to instance index
+    two_row = cl.concat([raa, raa], axis=0, ignore_index=True)
+    two_row.index = pd.DataFrame({"Company": ["A", "B"]})
+    s_perm = pd.Series([False, True], index=[1, 0])
+    sliced = two_row[s_perm]
+    assert len(sliced.index) == 1
+    assert sliced.index.iloc[0]["Company"] == "A"
+
+
+def test_legacy_pickle_compatibility(raa):
+    """Pickles saved prior to kdims/vdims migration should unpickle cleanly with warning."""
+    import pickle
+
+    # 1. Simulate oldest serialized Triangle containing kdims and vdims
+    state = raa.__dict__.copy()
+    state.pop("_axes", None)
+    state.pop("_index", None)
+    state.pop("_columns", None)
+    state["kdims"] = raa.index.values
+    state["vdims"] = raa.columns.values
+
+    restored = cl.Triangle.__new__(cl.Triangle)
+    with pytest.warns(
+        UserWarning, match="legacy pickled Triangle instance was detected"
+    ):
+        restored.__setstate__(state)
+
+    assert isinstance(restored.columns, pd.Index)
+    assert isinstance(restored.index, pd.DataFrame)
+    assert restored.index.equals(raa.index)
+    assert list(restored.columns) == list(raa.columns)
+    assert restored.loc["Total"].shape == raa.loc["Total"].shape
+    assert restored == raa
+    assert 0 in restored._axes
+    assert restored._axes[0].equals(raa.index)
+    assert 1 in restored._axes
+    assert restored._axes[1].equals(raa.columns)
+
+    # 2. Simulate intermediate serialized Triangle containing _kdims and _vdims
+    state_mid = raa.__dict__.copy()
+    state_mid.pop("_axes", None)
+    state_mid.pop("_index", None)
+    state_mid.pop("_columns", None)
+    state_mid["_kdims"] = raa.index.values
+    state_mid["_vdims"] = raa.columns.values
+
+    restored_mid = cl.Triangle.__new__(cl.Triangle)
+    with pytest.warns(
+        UserWarning, match="legacy pickled Triangle instance was detected"
+    ):
+        restored_mid.__setstate__(state_mid)
+
+    assert isinstance(restored_mid.columns, pd.Index)
+    assert isinstance(restored_mid.index, pd.DataFrame)
+    assert restored_mid.index.equals(raa.index)
+    assert list(restored_mid.columns) == list(raa.columns)
+    assert restored_mid == raa
+    assert 0 in restored_mid._axes
+    assert restored_mid._axes[0].equals(raa.index)
+    assert 1 in restored_mid._axes
+    assert restored_mid._axes[1].equals(raa.columns)
+
+    # 3. Simulate unprefixed kdims/vdims legacy pickle
+    state_unprefixed = raa.__dict__.copy()
+    state_unprefixed.pop("_axes", None)
+    state_unprefixed.pop("_index", None)
+    state_unprefixed.pop("_columns", None)
+    state_unprefixed.pop("_kdims", None)
+    state_unprefixed.pop("_vdims", None)
+    state_unprefixed["kdims"] = raa.index.values
+    state_unprefixed["vdims"] = raa.columns.values
+
+    restored_unprefixed = cl.Triangle.__new__(cl.Triangle)
+    with pytest.warns(
+        UserWarning, match="legacy pickled Triangle instance was detected"
+    ):
+        restored_unprefixed.__setstate__(state_unprefixed)
+    assert 0 in restored_unprefixed._axes
+    assert restored_unprefixed.index.equals(raa.index)
+
+    # 4. Simulate serialized Triangle containing _columns
+    state_col = raa.__dict__.copy()
+    state_col.pop("_axes", None)
+    state_col.pop("_index", None)
+    state_col.pop("_columns", None)
+    state_col["_kdims"] = raa.index.values
+    state_col["_columns"] = list(raa.columns)
+
+    restored_col = cl.Triangle.__new__(cl.Triangle)
+    with pytest.warns(
+        UserWarning, match="legacy pickled Triangle instance was detected"
+    ):
+        restored_col.__setstate__(state_col)
+    assert list(restored_col.columns) == list(raa.columns)
+    assert 1 in restored_col._axes
+
+    # 5. Simulate string-keyed _axes
+    state_str_axes = raa.__dict__.copy()
+    state_str_axes.pop("_index", None)
+    state_str_axes.pop("_columns", None)
+    state_str_axes["_axes"] = {
+        "index": raa.index.copy(),
+        "columns": raa.columns.copy(),
+    }
+    restored_str = cl.Triangle.__new__(cl.Triangle)
+    restored_str.__setstate__(state_str_axes)
+    assert 0 in restored_str._axes
+    assert 1 in restored_str._axes
+    assert restored_str.index.equals(raa.index)
+
+    # 6. Simulate redundant _index when 0 already in _axes
+    state_redundant = raa.__dict__.copy()
+    state_redundant["_axes"] = {0: raa.index.copy(), 1: raa.columns.copy()}
+    state_redundant["_index"] = raa.index.copy()
+    restored_red = cl.Triangle.__new__(cl.Triangle)
+    restored_red.__setstate__(state_redundant)
+    assert 0 in restored_red._axes
+    assert "_index" not in restored_red.__dict__
+
+    # 7. Simulate fallback when no kdims/vdims keys are found (no warning)
+    state_empty = raa.__dict__.copy()
+    state_empty.pop("_axes", None)
+    state_empty.pop("_index", None)
+    state_empty.pop("_columns", None)
+    state_empty.pop("kdims", None)
+    state_empty.pop("_kdims", None)
+    state_empty.pop("vdims", None)
+    state_empty.pop("_vdims", None)
+
+    restored_empty = cl.Triangle.__new__(cl.Triangle)
+    with warnings.catch_warnings():
+        warnings.simplefilter("error", UserWarning)
+        restored_empty.__setstate__(state_empty)
+    assert list(restored_empty.columns) == ["values"]
+    assert list(restored_empty.index.columns) == ["Total"]
+    assert 0 in restored_empty._axes
+    assert 1 in restored_empty._axes
+
+    # 8. Verify re-pickling the migrated instance works without warnings
+    with warnings.catch_warnings():
+        warnings.simplefilter("error", UserWarning)
+        roundtripped = pickle.loads(pickle.dumps(restored))
+    assert roundtripped == raa
+    assert 0 in roundtripped._axes
+    assert 1 in roundtripped._axes
+
+
+def test_sort_axis_columns_out_of_order(raa):
+    """Verify sort_axis(1) reorders values when columns are unsorted."""
+    import copy
+
+    t1 = copy.deepcopy(raa).rename("columns", ["B"])
+    t2 = copy.deepcopy(raa).rename("columns", ["A"])
+    combined = cl.concat([t1, t2], axis=1)
+    assert list(combined.columns) == ["B", "A"]
+    sorted_tri = combined.sort_axis(1)
+    assert list(sorted_tri.columns) == ["A", "B"]
 
 
 def test_valdev1(qtr):
@@ -3450,3 +3815,57 @@ def test_full_fill(raa: Triangle) -> None:
     full_tri = cl.Chainladder().fit(raa).full_triangle_
     fill_full_tri = full_tri.fill(200)
     assert np.all(fill_full_tri.values == np.broadcast_to([200], (1, 1, 10, 12)))
+
+
+def test_json_roundtrip_preserves_dataframe_index(raa, clrd) -> None:
+    """JSON serialization roundtrip must preserve Triangle.index as a pd.DataFrame."""
+    r_single = cl.read_json(raa.to_json())
+    assert isinstance(r_single.index, pd.DataFrame)
+    assert r_single.index.iloc[:, 0].tolist() == ["Total"]
+    assert r_single == raa
+
+    r_multi = cl.read_json(clrd.to_json())
+    assert isinstance(r_multi.index, pd.DataFrame)
+    pd.testing.assert_frame_equal(r_multi.index, clrd.index)
+    assert r_multi == clrd
+
+
+def test_triangle_copy_isolated_index(raa) -> None:
+    """Triangle.copy() must not share the mutable index DataFrame."""
+    copied = raa.copy()
+    assert copied.index is not raa.index
+    copied.index.iloc[0, 0] = "Modified"
+    assert raa.index.iloc[0, 0] == "Total"
+
+
+def test_triangle_index_setter_resets_row_index(raa) -> None:
+    """Assigning a DataFrame with non-default index to Triangle.index resets row index."""
+    custom_df = pd.DataFrame({"Total": ["Total"]}, index=[42])
+    raa1 = raa.copy()
+    raa1.index = custom_df
+    assert list(raa1.index.index) == [0]
+    # Boolean Series slicing should work cleanly
+    mask = raa1.index["Total"] == "Total"
+    sliced = raa1[mask]
+    assert sliced == raa1
+
+
+def test_sort_index_inplace(clrd: Triangle) -> None:
+    """sort_index with inplace=True updates values and index in sync without crashing."""
+    shuffled = clrd.iloc[::-1].copy()
+    result = shuffled.sort_index(inplace=True)
+    assert result is shuffled
+    pd.testing.assert_frame_equal(shuffled.index, clrd.index)
+    assert shuffled == clrd
+
+
+def test_duplicate_columns_assignment(raa) -> None:
+    """Assigning to duplicate column label updates first matching column without error."""
+    r1 = raa.copy()
+    r1.columns = ["a"]
+    r2 = raa.copy()
+    r2.columns = ["b"]
+    two = cl.concat([r1, r2], axis=1)
+    two.columns = ["col", "col"]
+    two["col"] = 42
+    assert two.values[0, 0, 0, 0] == 42
diff --git a/chainladder/core/triangle.py b/chainladder/core/triangle.py
index 38efe3324..850e47185 100644
--- a/chainladder/core/triangle.py
+++ b/chainladder/core/triangle.py
@@ -11,6 +11,9 @@
 import numpy as np
 import warnings
 from chainladder.core.base import TriangleBase
+from chainladder.core.axis import TriangleAxis
+from chainladder.core.axis import _set_columns
+from chainladder.core.axis import _set_index
 from chainladder.utils.sparse import sp
 from chainladder.core.slice import VirtualColumns
 from chainladder.core.correlation import DevelopmentCorrelation, ValuationCorrelation
@@ -28,7 +31,13 @@
 except ImportError:
     db = None
 
-from typing import cast, Optional, TYPE_CHECKING, Self
+from typing import (
+    Any,
+    cast,
+    Optional,
+    TYPE_CHECKING,
+    Self,
+)
 
 if TYPE_CHECKING:
     from pandas import DataFrame, Series
@@ -465,6 +474,17 @@ class Triangle(TriangleBase):
         1982  12000.0
     """
 
+    index = TriangleAxis(
+        0,
+        fset=_set_index,
+        doc="Represents the index axis of the triangle.",
+    )
+    columns = TriangleAxis(
+        1,
+        fset=_set_columns,
+        doc="Represents the column axis of the triangle.",
+    )
+
     @_deprecated_rename_argument("development", "valuation", remove_in_version="v2.0")
     @_deprecated_rename_argument(
         "development_format",
@@ -489,6 +509,7 @@ def __init__(
         *args,
         **kwargs,
     ):
+        self._axes: dict[int, Any] = {}
 
         # If data are present, validate the dimensions.
         if data is None:
@@ -616,23 +637,21 @@ def __init__(
             self.index_label: list = index
             data_agg[index[0]] = "Total"
 
-        self.kdims: np.ndarray
-        key_idx: np.ndarray
-        self.vdims: np.ndarray
+        index_idx: np.ndarray
         self.odims: np.ndarray
         orig_idx: np.ndarray
         self.ddims: ArrayLike
         dev_idx: np.ndarray
 
-        self.kdims, key_idx = self._set_kdims(data_agg, index)
-        self.vdims = np.array(columns)
+        index_values, index_idx = self._factorize_index(data_agg, index)
+        self.index = pd.DataFrame(list(index_values), columns=index)
+        self.columns = columns
         self.odims, orig_idx = self._set_odims(data_agg, date_axes)
         self.ddims, dev_idx = self._set_ddims(data_agg, date_axes)
 
         # Set remaining triangle properties.
         val_date: Timestamp = data_agg["__development__"].max()
         val_date = val_date.compute() if hasattr(val_date, "compute") else val_date
-        self.key_labels: list = index
         self.valuation_date: Timestamp = val_date
 
         if cumulative is None:
@@ -690,7 +709,7 @@ def __init__(
 
         coords, amts = self._set_values(
             data_agg=data_agg,
-            key_idx=key_idx,
+            index_idx=index_idx,
             columns=columns,
             orig_idx=orig_idx,
             dev_idx=dev_idx,
@@ -707,8 +726,8 @@ def __init__(
                     has_duplicates=False,
                     sorted=True,
                     shape=(
-                        len(self.kdims),
-                        len(self.vdims),
+                        len(self.index),
+                        len(self.columns),
                         len(self.odims),
                         len(self.ddims),
                     ),
@@ -779,33 +798,27 @@ def _split_ult(
         return data, ult
 
     @property
-    def index(self) -> DataFrame:
+    def key_labels(self) -> list:
         """
-        Returns a DataFrame of the unique values of the index.
+        Returns a list of the labels corresponding to the levels of the index.
         """
-        return pd.DataFrame(list(self.kdims), columns=self.key_labels)
-
-    @index.setter
-    def index(self, value) -> None:
-        self._len_check(self.index, value)
-        if type(value) is pd.DataFrame:
-            self.kdims = value.values
-            self.key_labels = list(value.columns)
-            self._set_slicers()
+        return list(self.index.columns)
+
+    @key_labels.setter
+    def key_labels(self, value) -> None:
+        if isinstance(value, str):
+            value = [value]
         else:
-            raise TypeError("index must be a pandas DataFrame")
+            value = list(value)
+        self.index.columns = value
+        self._set_slicers()
 
     @property
-    def columns(self):
-        return pd.Index(self.vdims, name="columns")
-
-    @columns.setter
-    def columns(self, value):
-        self._len_check(self.columns, value)
-        self.vdims = [value] if type(value) is str else value
-        if type(self.vdims) is list:
-            self.vdims = np.array(self.vdims)
-        self._set_slicers()
+    def axes(self) -> list:
+        """
+        Return a list representing the axes of the Triangle.
+        """
+        return [self.index, self.columns, self.origin, self.development]
 
     @property
     def columns_label(self) -> list:
@@ -2148,10 +2161,60 @@ def copy(self):
         """
         X = object.__new__(self.__class__)
         X.__dict__.update(vars(self))
+        if hasattr(self, "_axes"):
+            X._axes = {k: v.copy() for k, v in self._axes.items()}
         X._set_slicers()
         X.values = X.values.copy()
         return X
 
+    def __setstate__(self, state: dict) -> None:
+        """Migrate legacy pickled instances with 'kdims'/'_kdims' to '_index' and 'vdims'/'_vdims'/'_columns' to '_axes'."""
+        has_legacy_keys = any(
+            k in state for k in ("kdims", "_kdims", "vdims", "_vdims", "_columns")
+        )
+        if has_legacy_keys:
+            warnings.warn(
+                "A legacy pickled Triangle instance was detected. Please re-save your "
+                "triangle with chainladder 0.11.0+ as legacy pickle support will be "
+                "removed in a future release.",
+                UserWarning,
+                stacklevel=2,
+            )
+        key_labels = state.pop("key_labels", ["Total"])
+        if "_axes" not in state:
+            state["_axes"] = {}
+        if 0 not in state["_axes"]:
+            if "_index" in state:
+                state["_axes"][0] = state.pop("_index")
+            elif "index" in state["_axes"]:
+                state["_axes"][0] = state["_axes"].pop("index")
+            else:
+                raw_kdims = state.pop("_kdims", None)
+                if raw_kdims is None:
+                    raw_kdims = state.pop("kdims", None)
+                if raw_kdims is not None:
+                    state["_axes"][0] = pd.DataFrame(
+                        list(raw_kdims), columns=key_labels
+                    )
+                else:
+                    state["_axes"][0] = pd.DataFrame([["Total"]], columns=["Total"])
+        else:
+            state.pop("_index", None)
+        if 1 not in state["_axes"]:
+            if "columns" in state["_axes"]:
+                state["_axes"][1] = state["_axes"].pop("columns")
+            else:
+                raw_cols = state.pop("_columns", None)
+                if raw_cols is None:
+                    raw_cols = state.pop("_vdims", None)
+                if raw_cols is None:
+                    raw_cols = state.pop("vdims", None)
+                if raw_cols is not None:
+                    state["_axes"][1] = pd.Index(raw_cols, name="columns")
+                else:
+                    state["_axes"][1] = pd.Index(["values"], name="columns")
+        self.__dict__.update(state)
+
     def development_correlation(self, p_critical=0.5):
         """
         Mack (1997) test for correlations between subsequent development
@@ -2426,10 +2489,10 @@ def sort_axis(self, axis):
             return self.sort_index()
         obj = self.copy()
         if axis == 1:
-            sort = pd.Series(self.vdims).sort_values().index
-            if np.any(sort != pd.Series(self.vdims).index):
+            sort = pd.Series(self.columns).sort_values().index
+            if np.any(sort != pd.Series(self.columns).index):
                 obj.values = obj.values[:, list(sort), ...]
-                obj.vdims = obj.vdims[list(sort)]
+                obj.columns = self.columns[list(sort)]
         if axis == 2:
             sort = pd.Series(self.odims).sort_values().index
             if np.any(sort != pd.Series(self.odims).index):
diff --git a/chainladder/development/clark.py b/chainladder/development/clark.py
index adda5e414..f254a6490 100644
--- a/chainladder/development/clark.py
+++ b/chainladder/development/clark.py
@@ -329,10 +329,16 @@ def ldf(age):
         self.ldf_ = obj
         self.ldf_.valuation_date = pd.to_datetime(options.ULT_VAL)
         rows = X.index.set_index(X.key_labels).index
-        self.omega_ = pd.DataFrame(params[..., 0, 0], index=rows, columns=X.vdims)
-        self.theta_ = pd.DataFrame(params[..., 0, 1], index=rows, columns=X.vdims)
+        self.omega_ = pd.DataFrame(
+            params[..., 0, 0], index=rows, columns=X.columns.to_list()
+        )
+        self.theta_ = pd.DataFrame(
+            params[..., 0, 1], index=rows, columns=X.columns.to_list()
+        )
         if sample_weight:
-            self.elr_ = pd.DataFrame(params[..., 0, 2], index=rows, columns=X.vdims)
+            self.elr_ = pd.DataFrame(
+                params[..., 0, 2], index=rows, columns=X.columns.to_list()
+            )
         ultimate_ = (
             self._G(age=(latest_age - age_offset)[::-1]).swapaxes(-1, -2) * ld.values
         )
diff --git a/chainladder/development/munich.py b/chainladder/development/munich.py
index 87f34b3f3..cad2efd34 100644
--- a/chainladder/development/munich.py
+++ b/chainladder/development/munich.py
@@ -364,11 +364,11 @@ def _get_mcl_cdf(self, X, munich_full_triangle_):
         cdf_triangle = cdf_triangle[..., -1:] / cdf_triangle[..., :-1]
         paid = [item[0] for item in p_to_i]
         for n, item in enumerate(paid):
-            idx = np.where(X.cdf_.vdims == item)[0][0]
+            idx = np.where(X.cdf_.columns == item)[0][0]
             obj.values[:, idx : idx + 1, ...] = cdf_triangle[0, :, n : n + 1, ...]
         incurred = [item[1] for item in p_to_i]
         for n, item in enumerate(incurred):
-            idx = np.where(X.cdf_.vdims == item)[0][0]
+            idx = np.where(X.cdf_.columns == item)[0][0]
             obj.values[:, idx : idx + 1, ...] = cdf_triangle[1, :, n : n + 1, ...]
         obj._set_slicers()
         return obj
diff --git a/chainladder/methods/benktander.py b/chainladder/methods/benktander.py
index 7ae68b8a6..add140083 100644
--- a/chainladder/methods/benktander.py
+++ b/chainladder/methods/benktander.py
@@ -271,8 +271,7 @@ def _get_benktander_aprioris(self, X, sample_weight):
             apriori = random_state.lognormal(mu_log, sigma_log, X.shape[0])
             apriori = apriori.reshape(X.shape[0], -1)[..., None, None]
             apriori = sample_weight * apriori
-            apriori.kdims = X.kdims
-            apriori.key_labels = X.key_labels
+            apriori.index = X.index.copy()
         else:
             apriori = sample_weight * self.apriori
         apriori.columns = sample_weight.columns
diff --git a/chainladder/tails/bondy.py b/chainladder/tails/bondy.py
index ae32f1cd8..5c1072d88 100644
--- a/chainladder/tails/bondy.py
+++ b/chainladder/tails/bondy.py
@@ -161,7 +161,7 @@ def fit(self, X, y=None, sample_weight=None):
         obj = Development().fit_transform(X) if "ldf_" not in X else X
         b_optimized = []
         initial = xp.where(obj.ddims == earliest_age)[0][0] if earliest_age else 0
-        for num in range(len(obj.vdims)):
+        for num in range(len(obj.columns)):
             b0 = (xp.ones(obj.shape[0]) * 0.5)[:, None]
             data = xp.log(obj.ldf_.values[:, num, 0, initial:])
             b0 = xp.concatenate((b0, data[..., 0:1]), axis=1)
@@ -192,9 +192,11 @@ def fit(self, X, y=None, sample_weight=None):
         )
         fitted = xp.repeat(fitted, self.ldf_.shape[2], axis=2)
         rows = X.index.set_index(X.key_labels).index
-        self.b_ = pd.DataFrame(self.b_[..., 0, 0], index=rows, columns=X.vdims)
+        self.b_ = pd.DataFrame(
+            self.b_[..., 0, 0], index=rows, columns=X.columns.to_list()
+        )
         self.earliest_ldf_ = pd.DataFrame(
-            self.earliest_ldf_[..., 0, 0], index=rows, columns=X.vdims
+            self.earliest_ldf_[..., 0, 0], index=rows, columns=X.columns.to_list()
         )
         self.ldf_.values = xp.concatenate(
             (
diff --git a/chainladder/tails/curve.py b/chainladder/tails/curve.py
index 7782f460b..bb90db9fa 100644
--- a/chainladder/tails/curve.py
+++ b/chainladder/tails/curve.py
@@ -174,7 +174,6 @@ def fit(self, X, y=None, sample_weight=None):
         self : object
             Returns the instance itself.
         """
-
         X = X.copy()
         xp = X.get_array_module()
         if type(self.fit_period) is slice:
@@ -290,7 +289,7 @@ def slope_(self):
         """Does not work with munich"""
         rows = self.ldf_.index.set_index(self.ldf_.key_labels).index
         return pd.DataFrame(
-            self._slope_[..., 0, 0], index=rows, columns=self.ldf_.vdims
+            self._slope_[..., 0, 0], index=rows, columns=self.ldf_.columns.to_list()
         )
 
     @property
@@ -298,5 +297,7 @@ def intercept_(self):
         """Does not work with munich"""
         rows = self.ldf_.index.set_index(self.ldf_.key_labels).index
         return pd.DataFrame(
-            self._intercept_[..., 0, 0], index=rows, columns=self.ldf_.vdims
+            self._intercept_[..., 0, 0],
+            index=rows,
+            columns=self.ldf_.columns.to_list(),
         )
diff --git a/chainladder/utils/tests/test_utilities.py b/chainladder/utils/tests/test_utilities.py
index 96e5e7243..ae355634e 100644
--- a/chainladder/utils/tests/test_utilities.py
+++ b/chainladder/utils/tests/test_utilities.py
@@ -59,8 +59,8 @@ def from_sequence(seq):
 def test_triangle_json_io(clrd):
     clrd2 = cl.read_json(clrd.to_json(), array_backend=clrd.array_backend)
     assert clrd == clrd2
-    assert np.all(clrd.kdims == clrd2.kdims)
-    assert np.all(clrd.vdims == clrd2.vdims)
+    assert clrd.index.equals(clrd2.index)
+    assert np.all(clrd.columns == clrd2.columns)
     assert np.all(clrd.odims == clrd2.odims)
     assert np.all(clrd.ddims == clrd2.ddims)
     assert np.all(clrd.valuation == clrd2.valuation)
@@ -370,7 +370,7 @@ def test_load_sample_uspp() -> None:
         "friedland_uspp_increasing_claim_case",
     ]:
         tri = cl.load_sample(key)
-        assert set(str(c) for c in tri.vdims) == {
+        assert set(str(c) for c in tri.columns) == {
             "Reported Claims",
             "Paid Claims",
             "Earned Premium",
@@ -396,7 +396,7 @@ def test_load_sample_clrd2025() -> None:
         "EarnedPremCeded",
         "EarnedPremNet",
     }
-    assert set(str(c) for c in tri.vdims) == expected_columns
+    assert set(str(c) for c in tri.columns) == expected_columns
 
     # Accident years span 1998-2007.
     assert str(tri.origin.min()) == "1998"
@@ -612,6 +612,49 @@ def test_concat_axis1_duplicate_columns(raa: Triangle) -> None:
         cl.concat([raa, raa], axis=1)
 
 
+def test_concat_ignore_index_axes(raa: Triangle) -> None:
+    """Test concat with ignore_index=True along axes 0, 1, 2, and 3."""
+    res0 = cl.concat([raa, raa], axis=0, ignore_index=True)
+    assert len(res0.index) == 2
+    assert res0.key_labels == ["Index"]
+    assert list(res0.index["Index"]) == [0, 1]
+
+    t1 = copy.deepcopy(raa).rename("columns", ["A"])
+    t2 = copy.deepcopy(raa).rename("columns", ["B"])
+    res1 = cl.concat([t1, t2], axis=1, ignore_index=True)
+    assert list(res1.columns) == [0, 1]
+
+    o1 = raa.iloc[:, :, :5, :]
+    o2 = raa.iloc[:, :, 5:, :]
+    res2 = cl.concat([o1, o2], axis=2, ignore_index=True)
+    assert len(res2.odims) == len(raa.odims)
+
+    d1 = raa.iloc[:, :, :, :5]
+    d2 = raa.iloc[:, :, :, 5:]
+    res3 = cl.concat([d1, d2], axis=3, ignore_index=True)
+    assert len(res3.ddims) == len(raa.ddims)
+
+
+def test_concat_axis_1_index_dtype_relaxed(raa: Triangle) -> None:
+    """Verify concat on axis 1 succeeds when index values match but row index labels differ."""
+    t1 = copy.deepcopy(raa).rename("columns", ["A"])
+    t2 = copy.deepcopy(raa).rename("columns", ["B"])
+    t2.index = pd.DataFrame(t2.index.values, columns=t2.index.columns, index=[10])
+    res = cl.concat([t1, t2], axis=1)
+    assert list(res.columns) == ["A", "B"]
+
+
+def test_concat_axis_1_index_nan_preserved(raa: Triangle) -> None:
+    """Verify concat on axis 1 succeeds when index contains NaN entries."""
+    t1 = copy.deepcopy(raa).rename("columns", ["A"])
+    t2 = copy.deepcopy(raa).rename("columns", ["B"])
+    t1.index = pd.DataFrame([["Total", np.nan]], columns=["Key1", "Key2"])
+    t2.index = pd.DataFrame([["Total", np.nan]], columns=["Key1", "Key2"])
+    res = cl.concat([t1, t2], axis=1)
+    assert list(res.columns) == ["A", "B"]
+    assert pd.isna(res.index.iloc[0, 1])
+
+
 def test_maximum_2(raa: Triangle) -> None:
     """
     Run cl.maximum(raa, 5000) and check if each element in the resulting triangle is at least 5000.
diff --git a/chainladder/utils/utility_functions.py b/chainladder/utils/utility_functions.py
index a7822e062..0b9401392 100644
--- a/chainladder/utils/utility_functions.py
+++ b/chainladder/utils/utility_functions.py
@@ -428,10 +428,14 @@ def read_json(json_str, array_backend=None):
                 setattr(getattr(tri, k), "development_grain", tri.development_grain)
         if "dfs" in json_dict.keys():
             for k, v in json_dict["dfs"].items():
+                if k == "_index":
+                    continue
                 df = pd.read_json(StringIO(v))
                 if len(df.columns) == 1:
                     df = df.iloc[:, 0]
                 setattr(tri, k, df)
+        if isinstance(tri.index, pd.Series):
+            tri.index = tri.index.to_frame()
         if array_backend:
             return tri.set_backend(array_backend)
         else:
@@ -703,31 +707,67 @@ def concat(
             if list(objs[num].columns) != all_columns:
                 objs[num] = objs[num][all_columns]
     objs = set_common_backend(objs)
-    mapper = {0: "kdims", 1: "vdims", 2: "odims", 3: "ddims"}
-    for k in mapper.keys():
-        if k != axis and k != 1:  # All non-concat axes must be identical
-            a = np.array([getattr(obj, mapper[k]) for obj in objs])
-            assert np.all(a == a[0])
-        else:  # All elements of concat axis must be unique
-            if ignore_index:
-                new_axis = np.arange(
-                    np.sum([len(getattr(obj, mapper[axis])) for obj in objs])
+    if axis != 0:
+        for obj in objs[1:]:
+            assert (
+                list(obj.index.columns) == list(objs[0].index.columns)
+                and obj.index.shape == objs[0].index.shape
+                and bool(
+                    (
+                        (obj.index.values == objs[0].index.values)
+                        | (pd.isna(obj.index.values) & pd.isna(objs[0].index.values))
+                    ).all()
                 )
-                new_axis = new_axis[:, None] if axis == 0 else new_axis
-            else:
-                new_axis = np.concatenate([getattr(obj, mapper[axis]) for obj in objs])
-            if axis == 0:
-                assert len(pd.DataFrame(new_axis).drop_duplicates()) == len(new_axis)
-            else:
-                assert len(new_axis) == len(set(new_axis))
+            )
+    if axis != 2:
+        a = np.array([obj.odims for obj in objs])
+        assert np.all(a == a[0])
+    if axis != 3:
+        a = np.array([obj.ddims for obj in objs])
+        assert np.all(a == a[0])
+
     out = copy.deepcopy(objs[0])
     out.values = xp.concatenate([obj.values for obj in objs], axis=axis)
-    setattr(out, mapper[axis], new_axis)
-    if ignore_index and axis == 0:
-        out.key_labels = ["Index"]
+
+    if axis == 0:
+        if ignore_index:
+            new_axis = np.arange(sum([len(obj.index) for obj in objs]))[:, None]
+            out.index = pd.DataFrame(new_axis, columns=["Index"])
+            out.key_labels = ["Index"]
+        else:
+            new_axis = pd.concat([obj.index for obj in objs], ignore_index=True)
+            assert len(new_axis.drop_duplicates()) == len(new_axis)
+            out.index = new_axis
+            out.key_labels = list(new_axis.columns)
+    elif axis == 1:
+        if ignore_index:
+            new_axis = pd.Index(
+                np.arange(sum([len(obj.columns) for obj in objs])), name="columns"
+            )
+        else:
+            new_axis = pd.Index(
+                np.concatenate([obj.columns for obj in objs]), name="columns"
+            )
+            assert len(new_axis) == len(set(new_axis))
+        out.columns = new_axis
+    elif axis == 2:
+        if ignore_index:
+            new_axis = np.arange(sum([len(obj.odims) for obj in objs]))
+        else:
+            new_axis = np.concatenate([obj.odims for obj in objs])
+            assert len(new_axis) == len(set(new_axis))
+        out.odims = new_axis
+    elif axis == 3:
+        if ignore_index:
+            new_axis = np.arange(sum([len(obj.ddims) for obj in objs]))
+        else:
+            new_axis = np.concatenate([obj.ddims for obj in objs])
+            assert len(new_axis) == len(set(new_axis))
+        out.ddims = new_axis
+        if out.ddims.dtype == __dt64_dtype__ and type(out.ddims) is np.ndarray:
+            out.ddims = pd.DatetimeIndex(out.ddims)
+
     out.valuation_date = pd.Series([obj.valuation_date for obj in objs]).max()
-    if out.ddims.dtype == __dt64_dtype__ and type(out.ddims) is np.ndarray:
-        out.ddims = pd.DatetimeIndex(out.ddims)
     out._set_slicers()
     if sort:
         return out.sort_axis(axis)
@@ -1212,3 +1252,53 @@ def warn_exclusions_ignored(preserve):
             + " link ratio(s) is required for development estimation."
         )
     warnings.warn(warning)
+
+
+# Canonical ordering of Triangle dimensions.
+_axis_orders: tuple[str, ...] = ("index", "columns", "origin", "development")
+
+# Lookup mapping supporting integers (positive and negative) and string names.
+_axis_to_axis_number: dict[int | str, int] = {
+    **{i: i for i in range(4)},
+    **{i - 4: i for i in range(4)},
+    **{name: i for i, name in enumerate(_axis_orders)},
+}
+
+
+def _get_axis_number(axis: int | str) -> int:
+    """
+    Convert an axis name or integer (-4..3) to its canonical integer 0..3.
+
+    Parameters
+    ----------
+    axis : int or str
+        Axis name or integer representation.
+
+    Returns
+    -------
+    int
+        Canonical integer axis in 0..3.
+    """
+    if axis not in _axis_to_axis_number:
+        raise ValueError(
+            f"No axis named {axis} for object type Triangle. "
+            f"Valid axes are {_axis_orders} or integers 0..3."
+        )
+    return _axis_to_axis_number[axis]
+
+
+def _get_axis_name(axis: int | str) -> str:
+    """
+    Convert an axis name or integer to its canonical name.
+
+    Parameters
+    ----------
+    axis : int or str
+        Axis name or integer representation.
+
+    Returns
+    -------
+    str
+        Canonical axis name ('index', 'columns', 'origin', 'development').
+    """
+    return _axis_orders[_get_axis_number(axis)]