diff --git a/docs/explanations/components/property_package/general/eos/cubic.rst b/docs/explanations/components/property_package/general/eos/cubic.rst
index 1893289ff0..115896bb98 100644
--- a/docs/explanations/components/property_package/general/eos/cubic.rst
+++ b/docs/explanations/components/property_package/general/eos/cubic.rst
@@ -67,6 +67,7 @@ Cubic equations of state require the following parameters to be defined:
 
 1. `omega` (Pitzer acentricity factor) needs to be defined for each component (in the `parameter_data` for each component).
 2. `kappa` (binary interaction parameters) needs to be defined for each component pair in the system. This parameter needs to be defined in the general `parameter_data` argument for the overall property package (as it can be used in multiple phases).
+3. Component critical properties (`compress_fact_crit`, `dens_mol_crit`, `pressure_crit`, and `temperature_crit`) are often required as well.
 
 Calculation of Properties
 -------------------------
@@ -77,6 +78,26 @@ Many thermophysical properties are calculated using an ideal and residual term,
 
 The residual term is derived from the partial derivatives of the cubic equation of state, whilst the ideal term is determined using pure component properties for the ideal gas phase defined for each component.
 
+Critical Properties of Mixtures
+-------------------------------
+
+The critical properties of the mixture are calculated by solving the cubic equation of state at the critical point using the following constraints.
+
+.. math:: \Omega_A R^2 T_c^2 = a_{c,m}P_c
+.. math:: \Omega_BRT_c = b_{c,m}P_c
+.. math:: 0 = Z_c^3 - (1+B_c-uB_c)Z_c^2 + (A_c-uB_c-(u-w)B_c^2)Z_c - A_cB_c-wB_c^2-wB_c^3
+.. math:: P_c = Z_cRT_c\rho_{c}
+
+
+with the following expressions:
+
+.. math:: a_{c,j} = \frac{\Omega_AR^2T_{c,j}^2}{P_{c, j}}\alpha_j(T_c)
+.. math:: a_{c,m} = \sum_i{\sum_j{y_iy_j(a_{c,i}a_{c,j})^{1/2}(1-\kappa_{ij})}}
+.. math:: b_{c,m} = \sum_i{y_ib_i}
+.. math:: A_c = \frac{a_{c,m}P_c}{R^2T_c^2}
+.. math:: B_c = \frac{b_{c,m}P_c}{RT_c}
+
+
 Mass Density by Phase
 ---------------------
 
@@ -109,14 +130,14 @@ Component Molar Enthalpy by Phase
 
 Component molar enthalpies by phase are calculated using the pure component method provided by the users in the property package configuration arguments.
 
-Molar Isobaric Heat Capcity (:math:`C_p`)
------------------------------------------
+Molar Isobaric Heat Capacity (:math:`C_p`)
+------------------------------------------
 
-The ideal molar isobaric heat capcity term is calculated from the weighted sum of the (ideal) component molar isobaric heat capacity:
+The ideal molar isobaric heat capacity term is calculated from the weighted sum of the (ideal) component molar isobaric heat capacity:
 
 .. math:: C_{p, ig}^0 = \sum_j y_j C_{p, ig, j}
 
-The residual molar isobaric heat capcity term is given by:
+The residual molar isobaric heat capacity term is given by:
 
 .. math:: C_p^r = R \left[ T \left(\frac{\partial Z}{\partial T}\right)_P + Z - 1 \right] +  \frac{ T \frac{d^2a_m}{dT^2}}{\sqrt{u^2 - 4w} \cdot b_m} \ln \left[ \frac{2Z + uB + \sqrt{u^2 - 4w} B}{2Z + uB - \sqrt{u^2 - 4w} B} \right]
 .. math:: + \left(a_m - T \frac{da_m}{dT}\right) \cdot \frac{B}{b_m} \cdot \frac{\left(\frac{\partial Z}{\partial T}\right)_P + \frac{Z}{T}}{Z^2 + Z uB + wB^2}
diff --git a/docs/explanations/components/property_package/general/eos/ideal.rst b/docs/explanations/components/property_package/general/eos/ideal.rst
index f4b1a37d36..38fc9cd32b 100644
--- a/docs/explanations/components/property_package/general/eos/ideal.rst
+++ b/docs/explanations/components/property_package/general/eos/ideal.rst
@@ -9,6 +9,11 @@ Introduction
 
 Ideal behavior represents the simplest possible equation of state that ensures thermodynamic consistency between different properties.
 
+Critical Properties of Mixtures
+-------------------------------
+
+The Ideal Equation of State module does not support the calculation of mixture critical properties as the ideal equation of state cannot represent the critical point.
+
 Mass Density by Phase
 ---------------------
 
diff --git a/docs/explanations/components/property_package/general/phase_def.rst b/docs/explanations/components/property_package/general/phase_def.rst
index 22d64e9f17..d7df89afb3 100644
--- a/docs/explanations/components/property_package/general/phase_def.rst
+++ b/docs/explanations/components/property_package/general/phase_def.rst
@@ -34,7 +34,7 @@ Equations of state (or equivalent methods) describe the relationship between dif
 
 A wide range of equations of states are available in literature for different applications and levels of rigor, and the IDAES Generic Property Package Framework provides a number of prebuilt modules for users, which are listed below.
 
-Equation of state packages may allow for user options (e.g. choosing a specific type of cubic equation of state). The options are set using the `equation_of_state_options` argument, and the options available are described in the documentation of each equation of state module.
+Equation of state packages may allow for user options (e.g., choosing a specific type of cubic equation of state). The options are set using the `equation_of_state_options` argument, and the options available are described in the documentation of each equation of state module.
 
 Equation of State Libraries
 """""""""""""""""""""""""""
@@ -44,6 +44,20 @@ Equation of State Libraries
 
     eos/ideal
     eos/cubic
+    
+Critical Properties of Mixtures
+"""""""""""""""""""""""""""""""
+
+Calculation of the critical properties of mixtures depends on the equation of state being used. As the Modular Property Package framework allows users to specify different equations of state for each phase, the following logic is used to determine which equation of state to use for calculating critical properties.
+
+1. If a vapor-liquid equilibrium pair is defined in the `"phases_in_equilibrium"` configuration argument, then the liquid phase from this pair is used (IDAES generally assumes supercritical fluids are liquid-like).
+2. If no vapor-liquid equilibrium pair is defined, then the first liquid phase defined is used.
+3. If no liquid phases are defined then the vapor phase is used (it is assumed there will only be one vapor phase).
+4. If no vapor phase is defined, then a `PropertyPackageError` is returned as there is no suitable phase for calculating critical properties.
+
+Note that not all Equations of State are suitable for calculating critical properties as many cannot represent the critical conditions. The Ideal equation of State is one common example of an equation of state that DOES NOT support critical properties.
+
+During initialization, the critical properties of the mixture are approximated using the mole fraction weighted sum of the component critical properties (note that users must provide values for all component critical properties (i.e., `compress_fact_crit`, `dens_mol_crit`, `pressure_crit` and `temperature_crit`) for initialization if critical properties are to be calculated. Note that it is up to the user to ensure these values are consistent (if necessary).
 
 Phase-Specific Parameter
 ^^^^^^^^^^^^^^^^^^^^^^^^
@@ -53,6 +67,6 @@ In some cases, a property package may include parameters which are specific to a
 Phases with Partial Component Lists
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 
-In many applications a mixture will contain species that only appear in a single phase (either by nature or assumption). Common examples include crystalline solids and non-condensable gases. The IDAES Generic Property Package Framework provides support for these behaviors and allows users to specify phase-specific component lists (i.e. a list of components which appear in a given phase).
+In many applications a mixture will contain species that only appear in a single phase (either by nature or assumption). Common examples include crystalline solids and non-condensable gases. The IDAES Generic Property Package Framework provides support for these behaviors and allows users to specify phase-specific component lists (i.e., a list of components which appear in a given phase).
 
 This is done by providing a phase with a `component_list` argument, which provides a `list` of component names which appear in the phase. The framework automatically validates the `component_list` argument to ensure that it is a sub-set of the master component list for the property package, and will inform the user if an unrecognized component is included. If a phase is not provided with a `component_list` argument it is assumed that all components defined in the master component list may be present in the phase.
diff --git a/idaes/core/base/components.py b/idaes/core/base/components.py
index f0e8baa507..5f9bcbe091 100644
--- a/idaes/core/base/components.py
+++ b/idaes/core/base/components.py
@@ -246,6 +246,7 @@ def build(self):
 
         # Create Vars for common parameters
         var_dict = {
+            "compress_fact_crit": pyunits.dimensionless,
             "dens_mol_crit": base_units.DENSITY_MOLE,
             "omega": pyunits.dimensionless,
             "pressure_crit": base_units.PRESSURE,
diff --git a/idaes/core/base/property_meta.py b/idaes/core/base/property_meta.py
index b12aee115c..502e11b15d 100644
--- a/idaes/core/base/property_meta.py
+++ b/idaes/core/base/property_meta.py
@@ -257,9 +257,20 @@ def VOLUME(self):
         return self._length**3
 
     @property
-    def MOLAR_VOLUME(self):
+    def VOLUME_MASS(self):
+        return self._length**3 * self._mass**-1
+
+    @property
+    def VOLUME_MOLE(self):
         return self._length**3 * self._amount**-1
 
+    # Backward compatibility name
+    @property
+    def MOLAR_VOLUME(self):
+        msg = "The unit name MOLAR_VOLUME is being deprecated in favor of VOLUME_MOL."
+        deprecation_warning(msg=msg, logger=_log, version="2.3.0", remove_in="3.0.0")
+        return self.VOLUME_MOLE
+
     # Flows
     @property
     def FLOW_MASS(self):
diff --git a/idaes/core/base/property_set.py b/idaes/core/base/property_set.py
index 70af3f584f..e9b8a41629 100644
--- a/idaes/core/base/property_set.py
+++ b/idaes/core/base/property_set.py
@@ -627,6 +627,11 @@ class StandardPropertySet(PropertySetBase):
         doc="Compressiblity Factor",
         units=pyunits.dimensionless,
     )
+    compress_fact_crit = PropertyMetadata(
+        name="compress_fact_crit",
+        doc="Compressiblity Factor at Critical Point",
+        units=pyunits.dimensionless,
+    )
     conc_mass = PropertyMetadata(
         name="conc_mass",
         doc="Concentration on a Mass Basis",
@@ -890,12 +895,12 @@ class StandardPropertySet(PropertySetBase):
     vol_mol = PropertyMetadata(
         name="vol_mol",
         doc="Molar Volume",
-        units="MOLAR_VOLUME",
+        units="VOLUME_MOLE",
     )
     vol_mol_crit = PropertyMetadata(
         name="vol_mol",
         doc="Molar Volume at Critical Point",
-        units="MOLAR_VOLUME",
+        units="VOLUME_MOLE",
     )
     # Log terms
     log_act = PropertyMetadata(
diff --git a/idaes/core/base/tests/test_components.py b/idaes/core/base/tests/test_components.py
index a0fda6ce6a..31d378d24f 100644
--- a/idaes/core/base/tests/test_components.py
+++ b/idaes/core/base/tests/test_components.py
@@ -185,6 +185,8 @@ def get_metadata(self):
         m.comp = Component(
             parameter_data={
                 "mw": 10,
+                "compress_fact_crit": 1,
+                "dens_mol_crit": 55,
                 "pressure_crit": 1e5,
                 "temperature_crit": 500,
             }
@@ -193,6 +195,12 @@ def get_metadata(self):
         assert isinstance(m.comp.mw, Param)
         assert m.comp.mw.value == 10
 
+        assert isinstance(m.comp.compress_fact_crit, Var)
+        assert m.comp.compress_fact_crit.value == 1
+
+        assert isinstance(m.comp.dens_mol_crit, Var)
+        assert m.comp.dens_mol_crit.value == 55
+
         assert isinstance(m.comp.pressure_crit, Var)
         assert m.comp.pressure_crit.value == 1e5
 
diff --git a/idaes/core/base/tests/test_property_meta.py b/idaes/core/base/tests/test_property_meta.py
index 4e4e30b6d2..3ed90a2086 100644
--- a/idaes/core/base/tests/test_property_meta.py
+++ b/idaes/core/base/tests/test_property_meta.py
@@ -175,6 +175,8 @@ def test_luminous_intensity(unit_set):
 derived_quantities = {
     "area": units.m**2,
     "volume": units.m**3,
+    "volume_mass": units.m**3 / units.kg,
+    "volume_mole": units.m**3 / units.mol,
     "flow_mass": units.kg * units.s**-1,
     "flow_mole": units.mol * units.s**-1,
     "flow_vol": units.m**3 * units.s**-1,
diff --git a/idaes/models/properties/modular_properties/base/generic_property.py b/idaes/models/properties/modular_properties/base/generic_property.py
index f90609b0fb..1a21f239c2 100644
--- a/idaes/models/properties/modular_properties/base/generic_property.py
+++ b/idaes/models/properties/modular_properties/base/generic_property.py
@@ -1073,6 +1073,7 @@ def define_metadata(cls, obj):
                     "method": "_act_coeff_phase_comp_apparent"
                 },
                 "compress_fact_phase": {"method": "_compress_fact_phase"},
+                "compress_fact_crit": {"method": "_critical_props"},
                 "conc_mol_comp": {"method": "_conc_mol_comp"},
                 "conc_mol_phase_comp": {"method": "_conc_mol_phase_comp"},
                 "conc_mol_phase_comp_apparent": {
@@ -1094,6 +1095,7 @@ def define_metadata(cls, obj):
                 "dens_mass": {"method": "_dens_mass"},
                 "dens_mass_phase": {"method": "_dens_mass_phase"},
                 "dens_mol": {"method": "_dens_mol"},
+                "dens_mol_crit": {"method": "_critical_props"},
                 "dens_mol_phase": {"method": "_dens_mol_phase"},
                 "energy_internal_mol": {"method": "_energy_internal_mol"},
                 "energy_internal_mol_phase": {"method": "_energy_internal_mol_phase"},
@@ -1132,6 +1134,7 @@ def define_metadata(cls, obj):
                 "mw_comp": {"method": "_mw_comp"},
                 "mw_phase": {"method": "_mw_phase"},
                 "prandtl_number_phase": {"method": "_prandtl_number_phase"},
+                "pressure_crit": {"method": "_critical_props"},
                 "pressure_phase_comp": {"method": "_pressure_phase_comp"},
                 "pressure_phase_comp_true": {"method": "_pressure_phase_comp_true"},
                 "pressure_phase_comp_apparent": {
@@ -1142,6 +1145,7 @@ def define_metadata(cls, obj):
                 "pressure_osm_phase": {"method": "_pressure_osm_phase"},
                 "pressure_sat_comp": {"method": "_pressure_sat_comp"},
                 "surf_tens_phase": {"method": "_surf_tens_phase"},
+                "temperature_crit": {"method": "_critical_props"},
                 "temperature_bubble": {"method": "_temperature_bubble"},
                 "temperature_dew": {"method": "_temperature_dew"},
                 "therm_cond_phase": {"method": "_therm_cond_phase"},
@@ -1292,9 +1296,41 @@ def initialization_routine(
                     k.inherent_equilibrium_constraint.deactivate()
 
         # ---------------------------------------------------------------------
-        # If present, initialize bubble and dew point calculations
+        # If present, initialize bubble, dew, and critical point calculations
         for k in model.values():
             T_units = k.params.get_metadata().default_units.TEMPERATURE
+
+            # List of bubble and dew point constraints
+            cons_list = [
+                "eq_pressure_dew",
+                "eq_pressure_bubble",
+                "eq_temperature_dew",
+                "eq_temperature_bubble",
+                "eq_mole_frac_tbub",
+                "eq_mole_frac_tdew",
+                "eq_mole_frac_pbub",
+                "eq_mole_frac_pdew",
+                "log_mole_frac_tbub_eqn",
+                "log_mole_frac_tdew_eqn",
+                "log_mole_frac_pbub_eqn",
+                "log_mole_frac_pdew_eqn",
+                "mole_frac_comp_eq",
+                "log_mole_frac_comp_eqn",
+            ]
+
+            # Critical point
+            with k.lock_attribute_creation_context():
+                # Only need to look for one, as it is all-or-nothing
+                if hasattr(k, "pressure_crit"):
+                    # Initialize critical point properties
+                    _initialize_critical_props(k)
+                    # Add critical point constraints to cons_list
+                    ref_phase = k._get_critical_ref_phase()
+                    p_config = k.params.get_phase(ref_phase).config
+                    cons_list += (
+                        p_config.equation_of_state.list_critical_property_constraint_names()
+                    )
+
             # Bubble temperature initialization
             if hasattr(k, "_mole_frac_tbub"):
                 model._init_Tbub(k, T_units)
@@ -1311,36 +1347,22 @@ def initialization_routine(
             if hasattr(k, "_mole_frac_pdew"):
                 model._init_Pdew(k, T_units)
 
-            # Solve bubble and dew point constraints
+            # Solve bubble, dew, and critical point constraints
             for c in k.component_objects(Constraint):
                 # Deactivate all constraints not associated with bubble and dew
-                # points
-                if c.local_name not in (
-                    "eq_pressure_dew",
-                    "eq_pressure_bubble",
-                    "eq_temperature_dew",
-                    "eq_temperature_bubble",
-                    "eq_mole_frac_tbub",
-                    "eq_mole_frac_tdew",
-                    "eq_mole_frac_pbub",
-                    "eq_mole_frac_pdew",
-                    "log_mole_frac_tbub_eqn",
-                    "log_mole_frac_tdew_eqn",
-                    "log_mole_frac_pbub_eqn",
-                    "log_mole_frac_pdew_eqn",
-                    "mole_frac_comp_eq",
-                    "log_mole_frac_comp_eqn",
-                ):
+                # points or critical points
+                if c.local_name not in cons_list:
                     c.deactivate()
 
-        # If StateBlock has active constraints (i.e. has bubble and/or dew
+        # If StateBlock has active constraints (i.e. has bubble, dew, or critical
         # point calculations), solve the block to converge these
         for b in model.values():
             if number_activated_constraints(b) > 0:
                 if not degrees_of_freedom(b) == 0:
                     raise InitializationError(
                         f"{b.name} Unexpected degrees of freedom during "
-                        f"initialization at bubble and dew point step: {degrees_of_freedom(b)}."
+                        f"initialization at bubble, dew, and critical point step: "
+                        f"{degrees_of_freedom(b)}."
                     )
                 with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc:
                     solve_strongly_connected_components(
@@ -1349,7 +1371,7 @@ def initialization_routine(
                         solve_kwds={"tee": slc.tee},
                         calc_var_kwds=self.config.calculate_variable_options,
                     )
-        init_log.info("Dew and bubble point initialization completed.")
+        init_log.info("Bubble, dew, and critical point initialization completed.")
 
         # ---------------------------------------------------------------------
         # Calculate _teq if required
@@ -1729,9 +1751,37 @@ def initialize(
         opt = get_solver(solver, optarg)
 
         # ---------------------------------------------------------------------
-        # If present, initialize bubble and dew point calculations
+        # If present, initialize bubble, dew , and critical point calculations
         for k in blk.values():
             T_units = k.params.get_metadata().default_units.TEMPERATURE
+
+            # List of bubble and dew point constraints
+            cons_list = [
+                "eq_pressure_dew",
+                "eq_pressure_bubble",
+                "eq_temperature_dew",
+                "eq_temperature_bubble",
+                "eq_mole_frac_tbub",
+                "eq_mole_frac_tdew",
+                "eq_mole_frac_pbub",
+                "eq_mole_frac_pdew",
+                "log_mole_frac_tbub_eqn",
+                "log_mole_frac_tdew_eqn",
+                "log_mole_frac_pbub_eqn",
+                "log_mole_frac_pdew_eqn",
+                "mole_frac_comp_eq",
+                "log_mole_frac_comp_eqn",
+            ]
+
+            # Critical point
+            with k.lock_attribute_creation_context():
+                # Only need to look for one, as it is all-or-nothing
+                if hasattr(k, "pressure_crit"):
+                    # Initialize critical point properties
+                    _initialize_critical_props(k)
+                    # Add critical point constraints to cons_list
+                    cons_list += k.list_critical_property_constraint_names()
+
             # Bubble temperature initialization
             if hasattr(k, "_mole_frac_tbub"):
                 blk._init_Tbub(k, T_units)
@@ -1748,29 +1798,14 @@ def initialize(
             if hasattr(k, "_mole_frac_pdew"):
                 blk._init_Pdew(k, T_units)
 
-            # Solve bubble and dew point constraints
+            # Solve bubble, dew, and critical point constraints
             for c in k.component_objects(Constraint):
-                # Deactivate all constraints not associated with bubble and dew
-                # points
-                if c.local_name not in (
-                    "eq_pressure_dew",
-                    "eq_pressure_bubble",
-                    "eq_temperature_dew",
-                    "eq_temperature_bubble",
-                    "eq_mole_frac_tbub",
-                    "eq_mole_frac_tdew",
-                    "eq_mole_frac_pbub",
-                    "eq_mole_frac_pdew",
-                    "log_mole_frac_tbub_eqn",
-                    "log_mole_frac_tdew_eqn",
-                    "log_mole_frac_pbub_eqn",
-                    "log_mole_frac_pdew_eqn",
-                    "mole_frac_comp_eq",
-                    "log_mole_frac_comp_eqn",
-                ):
+                # Deactivate all constraints not associated with bubble, dew,
+                # or critical points
+                if c.local_name not in cons_list:
                     c.deactivate()
 
-        # If StateBlock has active constraints (i.e. has bubble and/or dew
+        # If StateBlock has active constraints (i.e. has bubble, dew, or critical
         # point calculations), solve the block to converge these
         n_cons = 0
         dof = 0
@@ -1781,12 +1816,12 @@ def initialize(
             if dof > 0:
                 raise InitializationError(
                     f"{blk.name} Unexpected degrees of freedom during "
-                    f"initialization at bubble and dew point step: {dof}."
+                    f"initialization at bubble, dew, and critical point step: {dof}."
                 )
             with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc:
                 res = solve_indexed_blocks(opt, [blk], tee=slc.tee)
             init_log.info(
-                "Dew and bubble point initialization: {}.".format(
+                "Bubble, dew, and critical point initialization: {}.".format(
                     idaeslog.condition(res)
                 )
             )
@@ -2920,6 +2955,78 @@ def get_mole_frac(self, phase=None):
                     return self.mole_frac_phase_comp_apparent
             return self.mole_frac_phase_comp_true
 
+    # -------------------------------------------------------------------------
+    # Mixture critical point
+    # Critical properties will be based off liquid phase if present, as we assume
+    # supercritical fluid is liquid-like.
+    def _get_critical_ref_phase(self):
+        ref_phase = None
+        vap_phase = None
+        # First, look for VLE pair and use liquid phase if present
+        if self.params.config.phases_in_equilibrium is not None:
+            for pp in self.params.config.phases_in_equilibrium:
+                p1 = self.params.get_phase(pp[0])
+                p2 = self.params.get_phase(pp[1])
+
+                if p1.is_liquid_phase() and p2.is_vapor_phase():
+                    return pp[0]
+                elif p2.is_liquid_phase() and p1.is_vapor_phase():
+                    return pp[1]
+
+        # Next, iterate all phases and return either the first liquid phase
+        # Also collect vapor phase for final fall back
+        for p in self.params.phase_list:
+            if self.params.get_phase(p).is_liquid_phase():
+                return p
+            elif self.params.get_phase(p).is_vapor_phase():
+                vap_phase = p
+
+        if ref_phase is None and vap_phase is not None:
+            # Use vapor phase for reference phase
+            return vap_phase
+        if ref_phase is None:
+            # If still no reference phase, raise an exception
+            raise PropertyPackageError(
+                "No liquid or vapor phase found to use as reference phase "
+                "for calculating critical properties."
+            )
+
+    def _critical_props(self):
+        ref_phase = self._get_critical_ref_phase()
+
+        try:
+            base_units = self.params.get_metadata().default_units
+
+            self.compress_fact_crit = Var(
+                doc="Critical compressibility factor of mixture",
+                units=pyunits.dimensionless,
+            )
+
+            self.dens_mol_crit = Var(
+                doc="Critical molar density of mixture",
+                units=base_units.DENSITY_MOLE,
+            )
+
+            self.pressure_crit = Var(
+                doc="Critical pressure of mixture",
+                units=base_units.PRESSURE,
+            )
+
+            self.temperature_crit = Var(
+                doc="Critical temperature of mixture",
+                units=base_units.TEMPERATURE,
+            )
+
+            p_config = self.params.get_phase(ref_phase).config
+            p_config.equation_of_state.build_critical_properties(self, ref_phase)
+
+        except AttributeError:
+            self.del_component(self.compress_fact_crit)
+            self.del_component(self.dens_mol_crit)
+            self.del_component(self.pressure_crit)
+            self.del_component(self.temperature_crit)
+            raise
+
     # -------------------------------------------------------------------------
     # Bubble and Dew Points
 
@@ -5045,3 +5152,30 @@ def rule_log_mole_frac(b, p1, p2, j):
             eqn = getattr(b, "log_mole_frac_" + abbrv + "_eqn")
             b.del_component(eqn)
         raise
+
+
+def _initialize_critical_props(state_data):
+    params = state_data.params
+    # Use mole weighted sum of component critical properties
+    crit_props = [
+        "compress_fact_crit",
+        "dens_mol_crit",
+        "pressure_crit",
+        "temperature_crit",
+    ]
+
+    for prop in crit_props:
+        try:
+            getattr(state_data, prop).set_value(
+                sum(
+                    state_data.mole_frac_comp[j]
+                    * getattr(params.get_component(j), prop)
+                    for j in state_data.component_list
+                )
+            )
+        except AttributeError:
+            raise AttributeError(
+                f"Missing attribute found when initializing {prop}. "
+                f"Make sure you have provided values for {prop} in all Component "
+                "declarations."
+            )
diff --git a/idaes/models/properties/modular_properties/base/tests/test_generic_property.py b/idaes/models/properties/modular_properties/base/tests/test_generic_property.py
index 85b4950afb..8f00f5c26b 100644
--- a/idaes/models/properties/modular_properties/base/tests/test_generic_property.py
+++ b/idaes/models/properties/modular_properties/base/tests/test_generic_property.py
@@ -25,6 +25,7 @@
 from idaes.models.properties.modular_properties.base.generic_property import (
     GenericParameterData,
     GenericStateBlock,
+    _initialize_critical_props,
 )
 from idaes.models.properties.modular_properties.base.tests.dummy_eos import DummyEoS
 
@@ -1125,6 +1126,7 @@ def test_default_scaling_convert(self):
 # Dummy methods for testing build calls to sub-modules
 def define_state(b):
     b.state_defined = True
+    b.temperature = Var(initialize=100)
 
 
 def get_material_flow_basis(self, *args, **kwargs):
@@ -1276,6 +1278,13 @@ def test_build_all(self, frame):
                         frame.params.get_metadata().properties[p.name].method,
                     )
                 continue
+            elif p.name.endswith(("_crit")):
+                # Critical properties will be tested elsewhere
+                assert hasattr(
+                    frame.props[1],
+                    frame.params.get_metadata().properties[p.name].method,
+                )
+                continue
             elif p.name.endswith(("bubble", "bub", "dew")):
                 # Bubble and dew properties require phase equilibria, which are
                 # not tested here
@@ -1461,3 +1470,478 @@ def test_visc_d_phase(self, frame):
 
         for p in frame.props[1].phase_list:
             assert value(frame.props[1].visc_d_phase[p]) == 4
+
+
+class TestCriticalProps:
+    @pytest.mark.unit
+    def test_get_critical_ref_phase_exception(self):
+        # No vapor or liquid phases
+        m = ConcreteModel()
+        m.params = DummyParameterBlock(
+            components={
+                "a": {},
+            },
+            phases={
+                "p1": {
+                    "equation_of_state": DummyEoS,
+                },
+                "p2": {
+                    "equation_of_state": DummyEoS,
+                },
+            },
+            state_definition=modules[__name__],
+            pressure_ref=100000.0,
+            temperature_ref=300,
+            base_units=base_units,
+        )
+
+        m.props = m.params.build_state_block([1], defined_state=False)
+
+        with pytest.raises(
+            AttributeError,
+            match="No liquid or vapor phase found to use as reference phase "
+            "for calculating critical properties.",
+        ):
+            m.props[1]._get_critical_ref_phase()
+
+    @pytest.mark.unit
+    def test_get_critical_ref_phase_exception(self):
+        # No vapor or liquid phases
+        m = ConcreteModel()
+        m.params = DummyParameterBlock(
+            components={
+                "a": {},
+            },
+            phases={
+                "p1": {
+                    "equation_of_state": DummyEoS,
+                },
+                "p2": {
+                    "equation_of_state": DummyEoS,
+                },
+            },
+            state_definition=modules[__name__],
+            pressure_ref=100000.0,
+            temperature_ref=300,
+            base_units=base_units,
+        )
+
+        m.props = m.params.build_state_block([1], defined_state=False)
+
+        with pytest.raises(
+            PropertyPackageError,
+            match="No liquid or vapor phase found to use as reference phase "
+            "for calculating critical properties.",
+        ):
+            m.props[1]._get_critical_ref_phase()
+
+    @pytest.mark.unit
+    def test_get_critical_ref_phase_VL(self):
+        # No vapor or liquid phases
+        m = ConcreteModel()
+        m.params = DummyParameterBlock(
+            components={
+                "a": {},
+            },
+            phases={
+                "p1": {
+                    "type": LiquidPhase,
+                    "equation_of_state": DummyEoS,
+                },
+                "p2": {
+                    "type": VaporPhase,
+                    "equation_of_state": DummyEoS,
+                },
+            },
+            state_definition=modules[__name__],
+            pressure_ref=100000.0,
+            temperature_ref=300,
+            base_units=base_units,
+        )
+
+        m.props = m.params.build_state_block([1], defined_state=False)
+
+        assert m.props[1]._get_critical_ref_phase() == "p1"
+
+    @pytest.mark.unit
+    def test_get_critical_ref_phase_LL(self):
+        # No vapor or liquid phases
+        m = ConcreteModel()
+        m.params = DummyParameterBlock(
+            components={
+                "a": {},
+            },
+            phases={
+                "p1": {
+                    "type": LiquidPhase,
+                    "equation_of_state": DummyEoS,
+                },
+                "p2": {
+                    "type": LiquidPhase,
+                    "equation_of_state": DummyEoS,
+                },
+            },
+            state_definition=modules[__name__],
+            pressure_ref=100000.0,
+            temperature_ref=300,
+            base_units=base_units,
+        )
+
+        m.props = m.params.build_state_block([1], defined_state=False)
+
+        assert m.props[1]._get_critical_ref_phase() == "p1"
+
+    @pytest.mark.unit
+    def test_get_critical_ref_phase_VLL(self):
+        class DummyPE:
+            def phase_equil(self, *args, **kwargs):
+                pass
+
+            @staticmethod
+            def return_expression(b, *args, **kwargs):
+                return b.temperature == 42
+
+        # No vapor or liquid phases
+        m = ConcreteModel()
+        m.params = DummyParameterBlock(
+            components={
+                "a": {
+                    "phase_equilibrium_form": {("p2", "p3"): DummyPE},
+                },
+            },
+            phases={
+                "p1": {
+                    "type": LiquidPhase,
+                    "equation_of_state": DummyEoS,
+                },
+                "p2": {
+                    "type": VaporPhase,
+                    "equation_of_state": DummyEoS,
+                },
+                "p3": {
+                    "type": LiquidPhase,
+                    "equation_of_state": DummyEoS,
+                },
+            },
+            state_definition=modules[__name__],
+            pressure_ref=100000.0,
+            temperature_ref=300,
+            base_units=base_units,
+            phases_in_equilibrium=[("p2", "p3")],
+            phase_equilibrium_state={("p2", "p3"): DummyPE},
+        )
+
+        m.props = m.params.build_state_block([1], defined_state=False)
+
+        assert m.props[1]._get_critical_ref_phase() == "p3"
+
+    @pytest.mark.unit
+    def test_get_critical_ref_phase_VX(self):
+        # No vapor or liquid phases
+        m = ConcreteModel()
+        m.params = DummyParameterBlock(
+            components={
+                "a": {},
+            },
+            phases={
+                "p1": {
+                    "type": VaporPhase,
+                    "equation_of_state": DummyEoS,
+                },
+                "p2": {
+                    "equation_of_state": DummyEoS,
+                },
+            },
+            state_definition=modules[__name__],
+            pressure_ref=100000.0,
+            temperature_ref=300,
+            base_units=base_units,
+        )
+
+        m.props = m.params.build_state_block([1], defined_state=False)
+
+        assert m.props[1]._get_critical_ref_phase() == "p1"
+
+    @pytest.mark.unit
+    def test_critical_props_not_implemented(self):
+        m = ConcreteModel()
+        m.params = DummyParameterBlock(
+            components={
+                "a": {},
+            },
+            phases={
+                "p1": {
+                    "type": LiquidPhase,
+                    "equation_of_state": DummyEoS,
+                },
+            },
+            state_definition=modules[__name__],
+            pressure_ref=100000.0,
+            temperature_ref=300,
+            base_units=base_units,
+        )
+
+        m.props = m.params.build_state_block([1], defined_state=False)
+
+        with pytest.raises(
+            NotImplementedError,
+            match="props\[1\] Equation of State module has not implemented a method for "
+            "build_critical_properties. Please contact the EoS developer or use a "
+            "different module.",
+        ):
+            m.props[1]._critical_props()
+
+    @pytest.mark.unit
+    def test_critical_props_dummy_method(self):
+        class DummyEoS2(DummyEoS):
+            @staticmethod
+            def build_critical_properties(b, *args, **kwargs):
+                b._dummy_crit_executed = True
+
+        m = ConcreteModel()
+        m.params = DummyParameterBlock(
+            components={
+                "a": {},
+            },
+            phases={
+                "p1": {
+                    "type": LiquidPhase,
+                    "equation_of_state": DummyEoS2,
+                },
+            },
+            state_definition=modules[__name__],
+            pressure_ref=100000.0,
+            temperature_ref=300,
+            base_units=base_units,
+        )
+
+        m.props = m.params.build_state_block([1], defined_state=False)
+        m.props[1]._critical_props()
+
+        assert m.props[1]._dummy_crit_executed
+
+        assert isinstance(m.props[1].compress_fact_crit, Var)
+        assert isinstance(m.props[1].dens_mol_crit, Var)
+        assert isinstance(m.props[1].pressure_crit, Var)
+        assert isinstance(m.props[1].temperature_crit, Var)
+
+    @pytest.mark.unit
+    def test_critical_props_attribute_error(self):
+        class DummyEoS2(DummyEoS):
+            @staticmethod
+            def build_critical_properties(b, *args, **kwargs):
+                raise AttributeError()
+
+        m = ConcreteModel()
+        m.params = DummyParameterBlock(
+            components={
+                "a": {},
+            },
+            phases={
+                "p1": {
+                    "type": LiquidPhase,
+                    "equation_of_state": DummyEoS2,
+                },
+            },
+            state_definition=modules[__name__],
+            pressure_ref=100000.0,
+            temperature_ref=300,
+            base_units=base_units,
+        )
+
+        m.props = m.params.build_state_block([1], defined_state=False)
+        try:
+            m.props[1]._critical_props()
+        except AttributeError:
+            assert not hasattr(m.props[1], "compress_fact_crit")
+            assert not hasattr(m.props[1], "dens_mol_crit")
+            assert not hasattr(m.props[1], "pressure_crit")
+            assert not hasattr(m.props[1], "temperature_crit")
+
+    @pytest.mark.unit
+    def test_initialize_critical_props(self):
+        m = ConcreteModel()
+
+        class DummyEoS2(DummyEoS):
+            @staticmethod
+            def build_critical_properties(b, *args, **kwargs):
+                b._dummy_crit_executed = True
+
+        # Dummy params block
+        m.params = DummyParameterBlock(
+            components={
+                "a": {
+                    "parameter_data": {
+                        "compress_fact_crit": 0.1,
+                        "dens_mol_crit": 1,
+                        "pressure_crit": 1e5,
+                        "temperature_crit": 100,
+                    }
+                },
+                "b": {
+                    "parameter_data": {
+                        "compress_fact_crit": 0.2,
+                        "dens_mol_crit": 2,
+                        "pressure_crit": 2e5,
+                        "temperature_crit": 200,
+                    }
+                },
+                "c": {
+                    "parameter_data": {
+                        "compress_fact_crit": 0.3,
+                        "dens_mol_crit": 3,
+                        "pressure_crit": 3e5,
+                        "temperature_crit": 300,
+                    }
+                },
+            },
+            phases={
+                "Vap": {
+                    "type": LiquidPhase,
+                    "equation_of_state": DummyEoS2,
+                },
+                "Liq": {
+                    "type": LiquidPhase,
+                    "equation_of_state": DummyEoS2,
+                },
+            },
+            base_units={
+                "time": pyunits.s,
+                "length": pyunits.m,
+                "mass": pyunits.kg,
+                "amount": pyunits.mol,
+                "temperature": pyunits.K,
+            },
+            state_definition=modules[__name__],
+            pressure_ref=100000.0,
+            temperature_ref=300,
+            parameter_data={
+                "PR_kappa": {
+                    ("a", "a"): 0.0,
+                    ("a", "b"): 0.0,
+                    ("a", "c"): 0.0,
+                    ("b", "a"): 0.0,
+                    ("b", "b"): 0.0,
+                    ("b", "c"): 0.0,
+                    ("c", "a"): 0.0,
+                    ("c", "b"): 0.0,
+                    ("c", "c"): 0.0,
+                },
+            },
+        )
+
+        m.props = m.params.build_state_block([1], defined_state=False)
+
+        # Add common variables
+        m.props[1].mole_frac_comp = Var(
+            m.params.component_list, initialize=0.5, bounds=(1e-12, 1)
+        )
+        m.props[1].mole_frac_comp["a"].fix(0.4)
+        m.props[1].mole_frac_comp["b"].fix(0.6)
+        m.props[1].mole_frac_comp["c"].fix(1e-8)
+
+        # Build critical props
+        m.props[1]._critical_props()
+
+        # Initialize critical props
+        _initialize_critical_props(m.props[1])
+
+        assert (
+            value(m.props[1].compress_fact_crit) == 0.4 * 0.1 + 0.6 * 0.2 + 1e-8 * 0.3
+        )
+        assert value(m.props[1].dens_mol_crit) == 0.4 * 1 + 0.6 * 2 + 1e-8 * 3
+        assert value(m.props[1].pressure_crit) == 0.4 * 1e5 + 0.6 * 2e5 + 1e-8 * 3e5
+        assert value(m.props[1].temperature_crit) == 0.4 * 1e2 + 0.6 * 2e2 + 1e-8 * 3e2
+
+    @pytest.mark.unit
+    def test_initialize_critical_props_missing_value(self):
+        m = ConcreteModel()
+
+        class DummyEoS2(DummyEoS):
+            @staticmethod
+            def build_critical_properties(b, *args, **kwargs):
+                b._dummy_crit_executed = True
+
+        # Dummy params block
+        m.params = DummyParameterBlock(
+            components={
+                "a": {
+                    "parameter_data": {
+                        "dens_mol_crit": 1,  # Missing Z_crit
+                        "pressure_crit": 1e5,
+                        "temperature_crit": 100,
+                    }
+                },
+                "b": {
+                    "parameter_data": {
+                        "compress_fact_crit": 0.2,
+                        "dens_mol_crit": 2,
+                        "pressure_crit": 2e5,
+                        "temperature_crit": 200,
+                    }
+                },
+                "c": {
+                    "parameter_data": {
+                        "compress_fact_crit": 0.3,
+                        "dens_mol_crit": 3,
+                        "pressure_crit": 3e5,
+                        "temperature_crit": 300,
+                    }
+                },
+            },
+            phases={
+                "Vap": {
+                    "type": LiquidPhase,
+                    "equation_of_state": DummyEoS2,
+                },
+                "Liq": {
+                    "type": LiquidPhase,
+                    "equation_of_state": DummyEoS2,
+                },
+            },
+            base_units={
+                "time": pyunits.s,
+                "length": pyunits.m,
+                "mass": pyunits.kg,
+                "amount": pyunits.mol,
+                "temperature": pyunits.K,
+            },
+            state_definition=modules[__name__],
+            pressure_ref=100000.0,
+            temperature_ref=300,
+            parameter_data={
+                "PR_kappa": {
+                    ("a", "a"): 0.0,
+                    ("a", "b"): 0.0,
+                    ("a", "c"): 0.0,
+                    ("b", "a"): 0.0,
+                    ("b", "b"): 0.0,
+                    ("b", "c"): 0.0,
+                    ("c", "a"): 0.0,
+                    ("c", "b"): 0.0,
+                    ("c", "c"): 0.0,
+                },
+            },
+        )
+
+        m.props = m.params.build_state_block([1], defined_state=False)
+
+        # Add common variables
+        m.props[1].mole_frac_comp = Var(
+            m.params.component_list, initialize=0.5, bounds=(1e-12, 1)
+        )
+        m.props[1].mole_frac_comp["a"].fix(0.4)
+        m.props[1].mole_frac_comp["b"].fix(0.6)
+        m.props[1].mole_frac_comp["c"].fix(1e-8)
+
+        # Build critical props
+        m.props[1]._critical_props()
+
+        # Initialize critical props
+        with pytest.raises(
+            AttributeError,
+            match="Missing attribute found when initializing compress_fact_crit. "
+            "Make sure you have provided values for compress_fact_crit in all "
+            "Component declarations.",
+        ):
+            _initialize_critical_props(m.props[1])
diff --git a/idaes/models/properties/modular_properties/eos/ceos.py b/idaes/models/properties/modular_properties/eos/ceos.py
index 2157d78db5..d419ef149f 100644
--- a/idaes/models/properties/modular_properties/eos/ceos.py
+++ b/idaes/models/properties/modular_properties/eos/ceos.py
@@ -25,6 +25,7 @@
 from copy import deepcopy
 
 from pyomo.environ import (
+    Constraint,
     exp,
     Expression,
     log,
@@ -350,7 +351,7 @@ def rule_bm(m, p):
         else:
             raise ConfigurationError(
                 "{} Unrecognized option for Equation of State "
-                "mixing_rule_a: {}. Must be an instance of MixingRuleB "
+                "mixing_rule_b: {}. Must be an instance of MixingRuleB "
                 "Enum.".format(b.name, mixing_rule_b)
             )
 
@@ -901,6 +902,125 @@ def vol_mol_phase_comp(b, p, j):
             * (b.compress_fact_phase[p] + _N_dZ_dNj(b, p, j))
         )
 
+    @staticmethod
+    def build_critical_properties(m, ref_phase):
+        pobj = m.params.get_phase(ref_phase)
+        cname = pobj.config.equation_of_state_options["type"].name
+        ctype = pobj._cubic_type
+
+        # Add components for calculation of mixture critical properties
+        def func_a_crit(m, j):
+            cobj = m.params.get_component(j)
+            fw = getattr(m, cname + "_fw")
+            return (
+                EoS_param[ctype]["omegaA"]
+                * (
+                    (Cubic.gas_constant(m) * cobj.temperature_crit) ** 2
+                    / cobj.pressure_crit
+                )
+                * (
+                    (1 + fw[j] * (1 - sqrt(m.temperature_crit / cobj.temperature_crit)))
+                    ** 2
+                )
+            )
+
+        m.add_component(
+            "a_crit",
+            Expression(
+                m.component_list, rule=func_a_crit, doc="Component a coefficient"
+            ),
+        )
+
+        def rule_am_crit(m):
+            try:
+                rule = m.params.get_phase(ref_phase).config.equation_of_state_options[
+                    "mixing_rule_a"
+                ]
+            except (KeyError, TypeError):
+                rule = MixingRuleA.default
+
+            a_crit = getattr(m, "a_crit")
+            if rule == MixingRuleA.default:
+                return rule_am_crit_default(m, cname, a_crit)
+            else:
+                raise ConfigurationError(
+                    "{} Unrecognized option for Equation of State "
+                    "mixing_rule_a: {}. Must be an instance of MixingRuleA "
+                    "Enum.".format(m.name, rule)
+                )
+
+        m.add_component("am_crit", Expression(rule=rule_am_crit))
+
+        def rule_bm_crit(m):
+            try:
+                rule = m.params.get_phase(ref_phase).config.equation_of_state_options[
+                    "mixing_rule_b"
+                ]
+            except KeyError:
+                rule = MixingRuleB.default
+
+            b = getattr(m, cname + "_b")
+            if rule == MixingRuleB.default:
+                return rule_bm_crit_default(m, b)
+            else:
+                raise ConfigurationError(
+                    "{} Unrecognized option for Equation of State "
+                    "mixing_rule_b: {}. Must be an instance of MixingRuleB "
+                    "Enum.".format(m.name, rule)
+                )
+
+        m.add_component("bm_crit", Expression(rule=rule_bm_crit))
+
+        def rule_A_crit(m):
+            am_crit = getattr(m, "am_crit")
+            return EoS_param[ctype]["omegaA"] * (
+                Cubic.gas_constant(m) * m.temperature_crit
+            ) ** 2 == (am_crit * m.pressure_crit)
+
+        m.add_component("A_crit", Constraint(rule=rule_A_crit))
+
+        def rule_B_crit(m):
+            bm_crit = getattr(m, "bm_crit")
+            return EoS_param[ctype]["coeff_b"] * Cubic.gas_constant(
+                m
+            ) * m.temperature_crit == (bm_crit * m.pressure_crit)
+
+        m.add_component("B_crit", Constraint(rule=rule_B_crit))
+
+        Acrit = (
+            m.am_crit
+            * m.pressure_crit
+            / (Cubic.gas_constant(m) * m.temperature_crit) ** 2
+        )
+        Bcrit = (
+            m.bm_crit * m.pressure_crit / (Cubic.gas_constant(m) * m.temperature_crit)
+        )
+
+        expr_write = CubicThermoExpressions(m)
+        if pobj.is_vapor_phase():
+            Z_crit = expr_write.z_vap(eos=pobj._cubic_type, A=Acrit, B=Bcrit)
+        elif pobj.is_liquid_phase():
+            Z_crit = expr_write.z_liq(eos=pobj._cubic_type, A=Acrit, B=Bcrit)
+
+        m.compress_fact_crit_eq = Constraint(rule=m.compress_fact_crit == Z_crit)
+
+        m.dens_mol_crit_eq = Constraint(
+            rule=m.pressure_crit
+            == m.compress_fact_crit
+            * Cubic.gas_constant(m)
+            * m.temperature_crit
+            * m.dens_mol_crit
+        )
+
+    @staticmethod
+    def list_critical_property_constraint_names():
+        return [
+            "dens_mol_crit_eq",
+            "compress_fact_crit_eq",
+            "A_crit",
+            "B_crit",
+        ]
+
 
 def _dZ_dT(blk, p):
     pobj = blk.params.get_phase(p)
@@ -1277,5 +1397,23 @@ def rule_am_default(m, cname, a, p, pp=()):
     )
 
 
+def rule_am_crit_default(m, cname, a_crit):
+    k = getattr(m.params, cname + "_kappa")
+    return sum(
+        sum(
+            m.mole_frac_comp[i]
+            * m.mole_frac_comp[j]
+            * sqrt(a_crit[i] * a_crit[j])
+            * (1 - k[i, j])
+            for j in m.component_list
+        )
+        for i in m.component_list
+    )
+
+
 def rule_bm_default(m, b, p):
     return sum(m.mole_frac_phase_comp[p, i] * b[i] for i in m.components_in_phase(p))
+
+
+def rule_bm_crit_default(m, b):
+    return sum(m.mole_frac_comp[i] * b[i] for i in m.component_list)
diff --git a/idaes/models/properties/modular_properties/eos/eos_base.py b/idaes/models/properties/modular_properties/eos/eos_base.py
index 7cf4cd8b04..37db8a1f40 100644
--- a/idaes/models/properties/modular_properties/eos/eos_base.py
+++ b/idaes/models/properties/modular_properties/eos/eos_base.py
@@ -50,6 +50,30 @@ def calculate_scaling_factors(b, pobj):
     def build_parameters(b):
         raise NotImplementedError(_msg(b, "build_parameters"))
 
+    @staticmethod
+    def build_critical_properties(b, ref_phase):
+        """
+        This method is used to define the constraints used to calculate mixture
+        critical properties.
+        """
+        raise NotImplementedError(_msg(b, "build_critical_properties"))
+
+    @staticmethod
+    def list_critical_property_constraint_names():
+        """
+        If critical properties are supported, this method is required during
+        initialization to identify the constraints that are associated
+        with calculating the critical properties.
+
+        Returns:
+            list of constraint names
+        """
+        raise NotImplementedError(
+            "Equation of State module has not implemented a method for "
+            "list_critical_property_constraint_names. Please contact the "
+            "EoS developer or use a different module."
+        )
+
     @staticmethod
     def get_vol_mol_pure(b, phase, comp, temperature):
         try:
diff --git a/idaes/models/properties/modular_properties/eos/tests/test_ceos_PR.py b/idaes/models/properties/modular_properties/eos/tests/test_ceos_PR.py
index 7fbaffa440..2240c0160e 100644
--- a/idaes/models/properties/modular_properties/eos/tests/test_ceos_PR.py
+++ b/idaes/models/properties/modular_properties/eos/tests/test_ceos_PR.py
@@ -19,7 +19,9 @@
 from sys import modules
 
 from pyomo.environ import (
+    assert_optimal_termination,
     ConcreteModel,
+    Constraint,
     Expression,
     log,
     sqrt,
@@ -28,15 +30,22 @@
     units as pyunits,
 )
 from pyomo.core.expr.numeric_expr import ExternalFunctionExpression
+from pyomo.util.check_units import assert_units_consistent
 
 from idaes.core import declare_process_block_class, LiquidPhase, VaporPhase, SolidPhase
 from idaes.models.properties.modular_properties.eos.ceos import Cubic, CubicType
 from idaes.models.properties.modular_properties.base.generic_property import (
+    GenericParameterBlock,
     GenericParameterData,
+    ModularPropertiesInitializer,
 )
 from idaes.core.util.exceptions import PropertyNotSupportedError, ConfigurationError
 from idaes.core.util.constants import Constants as const
 from idaes.models.properties.modular_properties.eos.ceos import cubic_roots_available
+from idaes.core.solvers import get_solver
+from idaes.core.initialization.initializer_base import InitializationStatus
+from idaes.models.properties.modular_properties.state_definitions import FTPx
+
 
 # Dummy method for property method calls
 def dummy_call(b, j, T):
@@ -1065,3 +1074,310 @@ def test_vol_mol_phase_comp(m):
             for j in m.params.component_list
         )
     )
+
+
+class TestCEOSCriticalProps:
+    @pytest.fixture(scope="class")
+    def model(self):
+        m = ConcreteModel()
+
+        # Dummy params block
+        m.params = GenericParameterBlock(
+            components={
+                "a": {
+                    "parameter_data": {
+                        "omega": 0.1,
+                        "compress_fact_crit": 1,
+                        "dens_mol_crit": 1,
+                        "pressure_crit": 1e5,
+                        "temperature_crit": 100,
+                    }
+                },
+                "b": {
+                    "parameter_data": {
+                        "omega": 0.2,
+                        "compress_fact_crit": 1,
+                        "dens_mol_crit": 1,
+                        "pressure_crit": 2e5,
+                        "temperature_crit": 200,
+                    }
+                },
+                "c": {
+                    "parameter_data": {
+                        "omega": 0.3,
+                        "compress_fact_crit": 1,
+                        "dens_mol_crit": 1,
+                        "pressure_crit": 3e5,
+                        "temperature_crit": 300,
+                    }
+                },
+            },
+            phases={
+                "Liq": {
+                    "type": LiquidPhase,
+                    "equation_of_state": Cubic,
+                    "equation_of_state_options": {"type": CubicType.PR},
+                },
+            },
+            base_units={
+                "time": pyunits.s,
+                "length": pyunits.m,
+                "mass": pyunits.kg,
+                "amount": pyunits.mol,
+                "temperature": pyunits.K,
+            },
+            state_definition=FTPx,
+            pressure_ref=100000.0,
+            temperature_ref=300,
+            parameter_data={
+                "PR_kappa": {
+                    ("a", "a"): 0.0,
+                    ("a", "b"): 0.0,
+                    ("a", "c"): 0.0,
+                    ("b", "a"): 0.0,
+                    ("b", "b"): 0.0,
+                    ("b", "c"): 0.0,
+                    ("c", "a"): 0.0,
+                    ("c", "b"): 0.0,
+                    ("c", "c"): 0.0,
+                },
+            },
+        )
+
+        m.props = m.params.state_block_class(
+            [1], defined_state=True, parameters=m.params
+        )
+
+        # Fix state variables
+        m.props[1].flow_mol.fix(1)
+        m.props[1].pressure.fix(1e5)
+        m.props[1].temperature.fix(300)
+
+        m.props[1].mole_frac_comp["a"].fix(0.4)
+        m.props[1].mole_frac_comp["b"].fix(0.6)
+        m.props[1].mole_frac_comp["c"].fix(1e-8)
+
+        # Build critical properties
+        m.props[1]._critical_props()
+
+        return m
+
+    @pytest.mark.unit
+    def test_build_critical_props(self, model):
+        assert isinstance(model.props[1].a_crit, Expression)
+        assert isinstance(model.props[1].am_crit, Expression)
+        assert isinstance(model.props[1].bm_crit, Expression)
+
+        base_units = model.params.get_metadata().default_units
+
+        assert str(model.props[1].a_crit["a"].expr) == str(
+            0.45724
+            * (
+                (
+                    pyunits.convert(
+                        const.gas_constant, to_units=base_units.GAS_CONSTANT
+                    )
+                    * model.params.a.temperature_crit
+                )
+                ** 2
+                / model.params.a.pressure_crit
+            )
+            * (
+                1
+                + model.props[1].PR_fw["a"]
+                * (
+                    1
+                    - sqrt(
+                        model.props[1].temperature_crit
+                        / model.params.a.temperature_crit
+                    )
+                )
+            )
+            ** 2
+        )
+        assert str(model.props[1].a_crit["b"].expr) == str(
+            0.45724
+            * (
+                (
+                    pyunits.convert(
+                        const.gas_constant, to_units=base_units.GAS_CONSTANT
+                    )
+                    * model.params.b.temperature_crit
+                )
+                ** 2
+                / model.params.b.pressure_crit
+            )
+            * (
+                1
+                + model.props[1].PR_fw["b"]
+                * (
+                    1
+                    - sqrt(
+                        model.props[1].temperature_crit
+                        / model.params.b.temperature_crit
+                    )
+                )
+            )
+            ** 2
+        )
+        assert str(model.props[1].a_crit["c"].expr) == str(
+            0.45724
+            * (
+                (
+                    pyunits.convert(
+                        const.gas_constant, to_units=base_units.GAS_CONSTANT
+                    )
+                    * model.params.c.temperature_crit
+                )
+                ** 2
+                / model.params.c.pressure_crit
+            )
+            * (
+                1
+                + model.props[1].PR_fw["c"]
+                * (
+                    1
+                    - sqrt(
+                        model.props[1].temperature_crit
+                        / model.params.c.temperature_crit
+                    )
+                )
+            )
+            ** 2
+        )
+
+        assert str(model.props[1].am_crit.expr) == str(
+            (
+                model.props[1].mole_frac_comp["a"]
+                * model.props[1].mole_frac_comp["a"]
+                * sqrt(model.props[1].a_crit["a"] * model.props[1].a_crit["a"])
+                * (1 - model.params.PR_kappa["a", "a"])
+            )
+            + (
+                model.props[1].mole_frac_comp["a"]
+                * model.props[1].mole_frac_comp["b"]
+                * sqrt(model.props[1].a_crit["a"] * model.props[1].a_crit["b"])
+                * (1 - model.params.PR_kappa["a", "b"])
+            )
+            + (
+                model.props[1].mole_frac_comp["a"]
+                * model.props[1].mole_frac_comp["c"]
+                * sqrt(model.props[1].a_crit["a"] * model.props[1].a_crit["c"])
+                * (1 - model.params.PR_kappa["a", "c"])
+            )
+            + (
+                model.props[1].mole_frac_comp["b"]
+                * model.props[1].mole_frac_comp["a"]
+                * sqrt(model.props[1].a_crit["b"] * model.props[1].a_crit["a"])
+                * (1 - model.params.PR_kappa["b", "a"])
+            )
+            + (
+                model.props[1].mole_frac_comp["b"]
+                * model.props[1].mole_frac_comp["b"]
+                * sqrt(model.props[1].a_crit["b"] * model.props[1].a_crit["b"])
+                * (1 - model.params.PR_kappa["b", "b"])
+            )
+            + (
+                model.props[1].mole_frac_comp["b"]
+                * model.props[1].mole_frac_comp["c"]
+                * sqrt(model.props[1].a_crit["b"] * model.props[1].a_crit["c"])
+                * (1 - model.params.PR_kappa["b", "c"])
+            )
+            + (
+                model.props[1].mole_frac_comp["c"]
+                * model.props[1].mole_frac_comp["a"]
+                * sqrt(model.props[1].a_crit["c"] * model.props[1].a_crit["a"])
+                * (1 - model.params.PR_kappa["c", "a"])
+            )
+            + (
+                model.props[1].mole_frac_comp["c"]
+                * model.props[1].mole_frac_comp["b"]
+                * sqrt(model.props[1].a_crit["c"] * model.props[1].a_crit["b"])
+                * (1 - model.params.PR_kappa["c", "b"])
+            )
+            + (
+                model.props[1].mole_frac_comp["c"]
+                * model.props[1].mole_frac_comp["c"]
+                * sqrt(model.props[1].a_crit["c"] * model.props[1].a_crit["c"])
+                * (1 - model.params.PR_kappa["c", "c"])
+            )
+        )
+
+        assert str(model.props[1].bm_crit.expr) == str(
+            model.props[1].mole_frac_comp["a"] * model.props[1].PR_b["a"]
+            + model.props[1].mole_frac_comp["b"] * model.props[1].PR_b["b"]
+            + model.props[1].mole_frac_comp["c"] * model.props[1].PR_b["c"]
+        )
+
+        assert isinstance(model.props[1].A_crit, Constraint)
+        assert isinstance(model.props[1].B_crit, Constraint)
+
+        assert str(model.props[1].A_crit.expr) == str(
+            0.45724
+            * (
+                pyunits.convert(const.gas_constant, to_units=base_units.GAS_CONSTANT)
+                * model.props[1].temperature_crit
+            )
+            ** 2
+            == model.props[1].am_crit * model.props[1].pressure_crit
+        )
+
+        assert str(model.props[1].B_crit.expr) == str(
+            0.07780
+            * pyunits.convert(const.gas_constant, to_units=base_units.GAS_CONSTANT)
+            * model.props[1].temperature_crit
+            == model.props[1].bm_crit * model.props[1].pressure_crit
+        )
+
+        assert_units_consistent(model)
+
+        # Check list_critical_property_constraint_names
+        con_list = Cubic.list_critical_property_constraint_names()
+
+        assert con_list == [
+            "dens_mol_crit_eq",
+            "compress_fact_crit_eq",
+            "A_crit",
+            "B_crit",
+        ]
+
+        # Check that names match constraints that were built
+        for c in con_list:
+            assert isinstance(getattr(model.props[1], c), Constraint)
+
+    @pytest.mark.component
+    @pytest.mark.solver
+    def test_initialize_with_critical_props(self, model):
+        initializer = ModularPropertiesInitializer()
+
+        status = initializer.initialize(model.props)
+
+        assert status == InitializationStatus.Ok
+
+        # Confirm results
+        assert value(model.props[1].compress_fact_crit) == pytest.approx(
+            0.321379, rel=1e-5
+        )
+        assert value(model.props[1].pressure_crit) == pytest.approx(1.58138e5, rel=1e-5)
+        assert value(model.props[1].temperature_crit) == pytest.approx(
+            158.138, rel=1e-5
+        )
+        assert value(model.props[1].dens_mol_crit) == pytest.approx(374.238, rel=1e-5)
+
+    @pytest.mark.component
+    @pytest.mark.solver
+    def test_solve_critical_props(self, model):
+        # Solve model
+        solver = get_solver()
+        res = solver.solve(model, tee=True)
+        assert_optimal_termination(res)
+
+        # Confirm results
+        assert value(model.props[1].compress_fact_crit) == pytest.approx(
+            0.321379, rel=1e-5
+        )
+        assert value(model.props[1].pressure_crit) == pytest.approx(1.58138e5, rel=1e-5)
+        assert value(model.props[1].temperature_crit) == pytest.approx(
+            158.138, rel=1e-5
+        )
+        assert value(model.props[1].dens_mol_crit) == pytest.approx(374.238, rel=1e-5)
diff --git a/idaes/models/properties/modular_properties/eos/tests/test_ideal.py b/idaes/models/properties/modular_properties/eos/tests/test_ideal.py
index 61dd2082fa..087f8a2b3b 100644
--- a/idaes/models/properties/modular_properties/eos/tests/test_ideal.py
+++ b/idaes/models/properties/modular_properties/eos/tests/test_ideal.py
@@ -217,7 +217,6 @@ def test_cv_mol_phase_comp(m):
 
 @pytest.mark.unit
 def test_heat_capacity_ratio_phase(m):
-
     m.props[1].cp_mol_phase = Var(m.params.phase_list)
     m.props[1].cv_mol_phase = Var(m.params.phase_list)