diff --git a/h2integrate/converters/natural_gas/natural_gas_cc_ct.py b/h2integrate/converters/natural_gas/natural_gas_cc_ct.py
index 4f3b34324..ce4b401df 100644
--- a/h2integrate/converters/natural_gas/natural_gas_cc_ct.py
+++ b/h2integrate/converters/natural_gas/natural_gas_cc_ct.py
@@ -2,6 +2,7 @@
 from attrs import field, define, validators
 
 from h2integrate.core.utilities import BaseConfig, merge_shared_inputs
+from h2integrate.reliability.models import WeibullReliabilityModel
 from h2integrate.core.model_baseclasses import (
     CostModelBaseClass,
     CostModelBaseConfig,
@@ -66,6 +67,7 @@ def initialize(self):
         self.commodity = "electricity"
         self.commodity_rate_units = "MW"
         self.commodity_amount_units = "MW*h"
+        self.reliability_model = None
 
     def setup(self):
         super().setup()
@@ -74,6 +76,12 @@ def setup(self):
             merge_shared_inputs(self.options["tech_config"]["model_inputs"], "performance"),
             additional_cls_name=self.__class__.__name__,
         )
+        if use_reliability := "reliability" in self.options["tech_config"]["model_inputs"]:
+            self.reliability_model = WeibullReliabilityModel.from_dict(
+                merge_shared_inputs(self.options["tech_config"]["model_inputs"], "reliability"),
+                additional_cls_name=self.__class__.__name__,
+            )
+        self.use_reliability = use_reliability
 
         # Add natural gas consumed output
         self.add_output(
@@ -153,6 +161,8 @@ def compute(self, inputs, outputs):
             inputs["electricity_command_value"],
         )
         natural_gas_demand = electricity_command_value * heat_rate_mmbtu_per_mwh
+        if self.use_reliability:
+            natural_gas_demand * self.reliability_model.availability
 
         # available feedstock, saturated at maximum system feedstock consumption
         natural_gas_available = np.where(
diff --git a/h2integrate/reliability/__init__.py b/h2integrate/reliability/__init__.py
new file mode 100644
index 000000000..e69de29bb
diff --git a/h2integrate/reliability/models.py b/h2integrate/reliability/models.py
new file mode 100644
index 000000000..6d7eb32a5
--- /dev/null
+++ b/h2integrate/reliability/models.py
@@ -0,0 +1,349 @@
+from abc import ABC, abstractmethod
+
+import numpy as np
+from attrs import field, define, validators
+from numpy.typing import ArrayLike
+
+from h2integrate.core.utilities import BaseConfig
+
+
+# generated from np.random.SeedSequence().entropy
+rng = np.random.default_rng(279299947538423226929715083173412195503)
+
+
+N_TIMESTEPS = 8760
+
+
+def create_failure_model(config: dict):
+    """Retrieves and initializes a matching reliability model."""
+    name = config["failure_model"]
+    fail_config = config["failure_parameters"]
+    match name:
+        case "WeibullReliability":
+            return WeibullReliability.from_dict(fail_config)
+        case "FixedIntervalReliability":
+            return FixedIntervalReliability.from_dict(fail_config)
+        case _:
+            raise NotImplementedError(f"{name} is not a valid model name")
+
+
+def create_maintenance_model(config: dict):
+    """Retrieves and initializes a matching reliability model."""
+    name = config["maintenance_model"]
+    maintenance_config = config["maintenance_parameters"]
+    match name:
+        case "WeibullReliability":
+            return WeibullReliability.from_dict(maintenance_config)
+        case "FixedIntervalReliability":
+            return FixedIntervalReliability.from_dict(maintenance_config)
+        case _:
+            raise NotImplementedError(f"{name} is not a valid model name")
+
+
+def generate_downtime_model(config: dict | int):
+    if not isinstance(config, dict):
+        return FixedDowntime(hours=config)
+    name = config["model"]
+    parameters = {k: v for k, v in config.items() if k != "model"}
+    match name:
+        case "LogNormalDowntime":
+            return LogNormalDowntime.from_dict(parameters)
+        case "FixedDowntime":
+            return FixedDowntime.from_dict(parameters)
+        case _:
+            raise NotImplementedError(f"{name} is not a valid model name")
+
+
+def update_dimensions(n_components: int, *args: np.ndarray):
+    """Update the dimensionality of a series of arrays and the value of :py:attr:`n_components`
+    to match. If the size of an array passed :py:attr:`args` is 1 and :py:attr:`n_components`
+    is greater than 1, all arrays passed to :py:attr:`args` will be broadcast to an array shaped
+    (:py:attr:`n_components`, 1). If the arrays are already larger than 1, then
+    :py:attr:`n_components will be updated to the size of the arrays.
+
+    Args:
+        n_components (int): Number of components in the model
+        args (np.ndarray): NumPy array of attribute values used to create a model.
+
+    Returns:
+        n_components: The value as passed or updated to match the size of arrays in :py:attr:`args`.
+        args: The arrays as passed or the arrays reshaped to shape (:py:attr:`n_components`, 1).
+    """
+    if args[0].size == 1 and n_components > 1:
+        shape = (n_components, 1)
+        for arg in args:
+            arg = np.broadcast_to(arg, shape)
+    n_components = args[0].size
+    return n_components, *args
+
+
+def float_array_converter(val: int | float | ArrayLike):
+    return np.array(val).astype(float).reshape(-1, 1)
+
+
+def int_array_converter(val: int | float | ArrayLike):
+    return np.array(val).astype(int).reshape(-1, 1)
+
+
+# NOTE: yes, I see the irony in creating a custom converter after removing all the attrs duplicates
+# NOTE: yes, I'll also be removing the above comment, but I'm sure someone will appreciate it
+def array_ge(val):
+    """Validates that all values of an array are greater than or equal to :py:attr:`val`."""
+
+    def validator(instance, attribute, value):
+        if any(value < val):
+            raise ValueError(f"'{attribute.name}' must have all values of at least {val}.")
+
+    return validator
+
+
+def array_gt(val):
+    """Validates that all values of an array are greater than or equal to :py:attr:`val`."""
+
+    def validator(instance, attribute, value):
+        if any(value <= val):
+            raise ValueError(f"'{attribute.name}' must have all values greater than {val}.")
+
+    return validator
+
+
+def match_shape(other):
+    """Validates the shape of an array matches the shape of :py:attr:`other`."""
+
+    def validator(instance, attribute, value):
+        other_arr = getattr(instance, other, None)
+        if other_arr is None:
+            raise ValueError(f"'{other}' does not exist.")
+        if other_arr.shape != value.shape:
+            msg = (
+                f"Shape of '{attribute.name}' {value.shape} does not match"
+                f" '{other}' {other_arr.shape}."
+            )
+            raise ValueError(msg)
+
+    return validator
+
+
+@define(kw_only=True)
+class BaseDowntime(ABC, BaseConfig):
+    @abstractmethod
+    def sample_downtime(self) -> np.ndarray: ...
+
+
+@define(kw_only=True)
+class BaseReliability(ABC, BaseConfig):
+    """Base reliability class responsible for common definitions and functionality.
+
+    Args:
+        burn_in (float): Number of years into the simulation to use as the starting point of the
+            the simulation's availability record. Defaults to 0.
+        n_components (int): Number of identical components to sample to avoid defining an array of
+            mean and sigma values when they are the same. After initialization this value changes
+            to align with the number of components being simulated, regardless of identicality.
+            Defaults to 1.
+        downtime (dict): Configuration for a downtime model.
+
+    Attributes:
+        time_to_failures (np.ndarray): Array of the hours to the next failure for each modeled
+            component. Generated by each subclass' :py:meth:`sample_events`.
+        downtime_per_event (np.ndarray): Number of hours of downtime corresponding to each downtime
+            event in :py:attr:`time_to_failures`.
+        availability (np.ndarray): Operational ratio of each modeled component for every time step
+            of the simulation.
+        system_availability (np.ndarray): Minimum operational ratio of all components for every time
+            time step of the simulation.
+    """
+
+    burn_in: float = field(default=0, converter=float, validator=validators.ge(0))
+    n_components: int = field(default=1, validator=(validators.instance_of(int), validators.ge(1)))
+    downtime: dict | BaseDowntime = field(converter=generate_downtime_model)
+
+    time_to_failures: np.ndarray = field(init=False, validator=validators.instance_of(np.ndarray))
+    downtime_per_event: np.ndarray = field(init=False, validator=validators.instance_of(np.ndarray))
+    availability: np.ndarray = field(init=False, validator=validators.instance_of(np.ndarray))
+    system_availability: np.ndarray = field(
+        init=False, validator=validators.instance_of(np.ndarray)
+    )
+
+    def sample_events(self) -> np.ndarray:
+        raise NotImplementedError("Failed to successfully subclass, please implement me.")
+
+    def create_downtime_events(self):
+        """Creates the ``time_to_failure`` and ``downtime_per_event`` arrays."""
+        self.time_to_failures = self.sample_events()
+        self.downtime_per_event = self.downtime.sample_downtime()
+
+    def calculate_availability(self):
+        """Determine the timing and duration of outages for a single year of simulation time."""
+        burn_in_time = np.ceil(self.burn_in * N_TIMESTEPS).astype(int)
+        simulation_end = burn_in_time + N_TIMESTEPS
+
+        accumulated = np.zeros((self.n_components, 1), dtype=int)
+        availability = np.ones((self.n_components, simulation_end), dtype=float)
+
+        while any(accumulated < simulation_end):
+            if not self.time_to_failures.size == 0:
+                self.create_downtime_events()
+            event = self.time_to_failures[:, 0].reshape(-1, 1)
+            self.time_to_failures = self.time_to_failures[:, 1:]
+            duration = self.downtime_per_event[:, 0].reshape(-1, 1)
+            self.downtime_per_event = self.downtime_per_event[:, 1:]
+
+            if all(event + accumulated) > simulation_end:
+                break
+
+            start = accumulated + event
+            end = start + duration
+            for i, (s, e) in enumerate(zip(start.flatten(), end.flatten())):
+                if s < simulation_end:
+                    e = min(simulation_end, e)
+                    availability[i, s:e] = 0
+            accumulated = end
+
+        self.availability = availability[:, burn_in_time:simulation_end]
+        self.system_availability = np.min(self.availability, axis=0)
+
+
+@define(kw_only=True)
+class FixedDowntime(BaseDowntime):
+    """Basic log-normal downtime model for generating the length of downtime for a given event.
+
+    Args:
+        hours (int | array-like): Length of downtime per event, in hours. Must be at least 1 hour.
+        n_components (int): Number of identical components to sample to avoid defining an array of
+            :py:attr:`hours` values when they are the same. Defaults to 1.
+    """
+
+    hours: int | ArrayLike = field(
+        converter=int_array_converter,
+        validator=(validators.instance_of(np.ndarray), array_ge(1)),
+    )
+    n_components: int = field(default=1, validator=(validators.instance_of(int), validators.ge(1)))
+
+    def __attrs_post_init__(self):
+        self.n_components, self.hours = update_dimensions(self.n_components, self.hours)
+
+    def sample_downtime(self):
+        """Return an array of 100 :py:attr:`hours`."""
+        return np.ones((1, 100), dtype=int) * self.hours
+
+
+@define(kw_only=True)
+class LogNormalDowntime(BaseDowntime):
+    """Basic log-normal downtime model for generating the length of downtime for a given event.
+
+    Args:
+        mean (float | array-like): Average length of downtime per event, in hours.
+        sigma (float | array-like): Standard deviation of the distribution(s), in hours.
+    """
+
+    mean: float = field(
+        converter=float_array_converter,
+        validator=(validators.instance_of(np.ndarray), array_ge(0)),
+    )
+    sigma: float = field(
+        converter=float_array_converter,
+        validator=(validators.instance_of(np.ndarray), match_shape("mean"), array_ge(0)),
+    )
+    n_components: int = field(default=1, validator=(validators.instance_of(int), validators.ge(1)))
+
+    def __attrs_post_init__(self):
+        self.n_components, self.mean, self.sigma = update_dimensions(
+            self.n_components, self.mean, self.sigma
+        )
+
+    def sample_downtime(self) -> np.ndarray:
+        """Return an array of 100 samples of each lognormal distribution."""
+        return np.ceil(rng.lognormal(self.mean, self.sigma, size=(self.mean.shape[0], 100))).astype(
+            int
+        )
+
+
+@define(kw_only=True)
+class WeibullReliability(BaseReliability):
+    r"""Basic reliability model for operating/not operating statuses.
+
+    Assumes a full operational shutdown with zero ramping of production for an hourly, 1 year
+    simulation.
+
+    Args:
+        scale (float): Also referred to as :math:`\lambda` or :math:`\alpha`. Determines
+            the scale of distribution, and is equivalent to the mean time
+            between failure in years (MTBF), or 1 / annual failure rate.
+        shape (float): Also referred to as ``k`` or :math:`\beta`. A value less than 1
+            corresponds to a decreasing hazard rate over time (break-in period failures);
+            a value greater than 1 corresponds to an increasing hazard rate over time (
+            aging/wear-out failures); and a value of 1 corresponds to a constant hazard
+            rate over time (exponential distribution).
+        downtime (float): Average amount of downtime per failure.
+
+    Attributes:
+        rng (np.random._generator.Generator): NumPy random generator object.
+
+    TODO:
+        - how to pass n_timesteps through from plant?
+        - burn-in
+    """
+
+    scale: float = field(
+        converter=float_array_converter,
+        validator=validators.instance_of(np.ndarray),
+    )
+    shape: float = field(
+        converter=float_array_converter,
+        validator=(validators.instance_of(np.ndarray), match_shape("scale")),
+    )
+
+    def __attrs_post_init__(self):
+        self.n_components, self.scale, self.shape = update_dimensions(
+            self.n_components, self.scale, self.shape
+        )
+
+        self.create_downtime_events()
+        self.calculate_availability()
+
+    def sample_events(self):
+        """Samples 100 events for each simulated component, rounding up to the nearest timestep."""
+        return np.ceil(
+            rng.weibull(self.shape, size=(self.shape.size, 100)) * self.scale * N_TIMESTEPS
+        ).astype(int)
+
+
+@define(kw_only=True)
+class FixedIntervalReliability(BaseReliability):
+    """Basic fixed interval downtime reliability model.
+
+    Args:
+        frequency (int | float | array-like): The annual frequency of events, e.g., 4 is equivalent
+            to a quarterly downtime event and 0.25 is equivalent to an every 4 years downtime event.
+            For all events the timing of the first event will be sampled within the first year or
+            interval period to offset events from being based on January 1st in an 8760.
+        downtime (int | float | dict): Either fixed length of each downtime, in hours, or a
+            configuration dictionary for a downtime length model.
+    """
+
+    frequency: float = field(
+        converter=float_array_converter,
+        validator=(validators.instance_of(np.ndarray), array_gt(0)),
+    )
+
+    def __attrs_post_init__(self):
+        self.n_components, self.frequency = update_dimensions(self.n_components, self.frequency)
+
+        self.create_downtime_events()
+        self.calculate_availability()
+
+    def sample_events(self):
+        """Creates the time to next failure array for each event's modality with the first event
+        occurring randomly either in the first year or first failure interval, whichever
+        comes first.
+
+        Returns:
+            time_to_failures (np.ndarray): An array of the next 100 events' time to next failure.
+        """
+        interval = np.ceil(8760 / (1 / self.frequency)).astype(int)
+        first_occurrence = rng.integers(0, np.where(interval > 8760, 8760, interval))
+        time_to_failures = np.hstack(
+            (first_occurrence, np.broadcast_to(interval, (interval.size, 99)))
+        )
+        return time_to_failures