diff --git a/NEWS.md b/NEWS.md
index 9dd53b3e..5099c893 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,5 +1,11 @@
 # Release notes
 
+## Unversioned
+
+* New functions (`variables_flow_resource()`, `constraints_resource()`,  `constraints_couple_resource()`) that dispatch on resource types, which allow for creation of new resource-specific variables and constraints in extension packages.
+* New function to indentify the unique resource types of a vector of resources
+* New function that segments the vector of resources into sub-vectors based on each resource type
+
 ## Version 0.9.4 (2025-11-26)
 
 ### Bugfixes
diff --git a/Project.toml b/Project.toml
index a6366609..1f2b997f 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "EnergyModelsBase"
 uuid = "5d7e687e-f956-46f3-9045-6f5a5fd49f50"
 authors = ["Lars Hellemo <Lars.Hellemo@sintef.no>, Julian Straus <Julian.Straus@sintef.no>"]
-version = "0.9.4"
+version = "0.9.5"
 
 [deps]
 JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
diff --git a/docs/make.jl b/docs/make.jl
index f09693c4..9beed0a6 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -58,6 +58,7 @@ makedocs(
         "How to" => Any[
             "Create a new element"=>"how-to/create_new_element.md",
             "Create a new node"=>"how-to/create-new-node.md",
+            "Extend resource functionality"=>"how-to/extend-resource-functionality.md",
             "Utilize TimeStruct"=>"how-to/utilize-timestruct.md",
             "Update models"=>"how-to/update-models.md",
             "Contribute to EnergyModelsBase"=>"how-to/contribute.md",
diff --git a/docs/src/how-to/extend-resource-functionality.md b/docs/src/how-to/extend-resource-functionality.md
new file mode 100644
index 00000000..9c341fdf
--- /dev/null
+++ b/docs/src/how-to/extend-resource-functionality.md
@@ -0,0 +1,207 @@
+# [Extend resource functionality](@id how_to-res_funct)
+
+```@meta
+CurrentModule = EMB
+```
+
+## [Concept](@id how_to-res_funct-concept)
+
+This guide shows how to extend resource functionality by adding a custom resource type and connecting it to custom variables and constraints through resource-dispatch functions.
+This is useful for modelling more complex resource behavior that cannot be captured by the default resource types where the standard behavior is built around energy or mass flow.
+
+The pattern follows the same structure as the resource dispatch test in `test/test_resource.jl`:
+
+1. Define a resource subtype with extra parameters.
+2. (Optionally) create a custom node subtype that uses the resource.
+3. Add resource-specific variables with [`variables_flow_resource`](@ref).
+4. Add resource-specific constraints with [`constraints_resource`](@ref).
+5. Couple node and link resource variables with [`constraints_couple_resource`](@ref).
+
+## [Example](@id how_to-res_funct-example)
+
+The following example illustrates the different steps that are required for creating a new resource with additional properties.
+It defines a `PotentialPower` resource which has as property a potential with upper and lower bounds in addition to its energy flow.
+The flow of this potential in and out of junctions follows equality constraints, as opposed to the energy and mass flow which follow sum constraints.
+
+The notation below follows the same conventions as the implementation and tests:
+
+- `𝒩` for nodes,
+- `ℒ` for links,
+- `𝒫` for resources,
+- `𝒯` for the time structure,
+- `ℒᶠʳᵒᵐ`, `ℒᵗᵒ` for outgoing and incoming links of a node, and
+- `𝒫ᵒᵘᵗ`, `𝒫ⁱⁿ`, `𝒫ˡⁱⁿᵏ` for resource subsets on outputs, inputs, and links.
+
+### 1. Define a special resource
+
+Create a subtype of [`Resource`](@ref) and keep `co2_int` as the second field for consistency with existing resource structures.
+Alternatively, you can create a new method for the internal function [`co2_int`](@ref).
+
+```julia
+struct PotentialPower <: Resource
+    id::String
+    co2_int::Float64
+    potential_lower::Float64
+    potential_upper::Float64
+end
+
+EMB.is_resource_emit(::PotentialPower) = false
+lower_limit(p::PotentialPower) = p.potential_lower
+upper_limit(p::PotentialPower) = p.potential_upper
+```
+
+### 2. Define a custom node (optional)
+
+If your resource needs dedicated node behavior, create a custom node subtype.
+If the node subtype is parametrized, it can handle different types of resources in different ways without defining multiple node types.
+In the dispatch test, the custom node is an intermediate `NetworkNode` with a potential loss, but without a loss in energy flow.
+
+```julia
+struct PotentialLossNode{T<:PotentialPower} <: NetworkNode
+    id::Any
+    cap::TimeProfile
+    opex_var::TimeProfile
+    opex_fixed::TimeProfile
+    resource::T
+    input::Dict{<:Resource,<:Real}
+    output::Dict{<:Resource,<:Real}
+    data::Vector{<:ExtensionData}
+    loss_factor::Float64
+end
+
+function PotentialLossNode(
+    id,
+    cap::TimeProfile,
+    opex_var::TimeProfile,
+    opex_fixed::TimeProfile,
+    resource::T,
+    loss_factor::Float64,
+) where {T<:PotentialPower}
+    return PotentialLossNode{T}(
+        id,
+        cap,
+        opex_var,
+        opex_fixed,
+        resource,
+        Dict(resource => 1.0),
+        Dict(resource => 1.0),
+        ExtensionData[],
+        loss_factor,
+    )
+end
+```
+
+### 3. Declare resource-specific variables
+
+Use [`variables_flow_resource`](@ref) to create resource variables.
+
+Important:
+
+- Declare each variable name once.
+- Filter `𝒩` and `ℒ` down to the subsets that actually use the special resource.
+- You can create resource dependent bounds as well.
+
+```julia
+function EMB.variables_flow_resource(
+    m,
+    𝒩::Vector{<:EMB.Node},
+    𝒫::Vector{<:PotentialPower},
+    𝒯,
+    modeltype::EnergyModel,
+)
+    𝒩ᵒᵘᵗ = filter(n -> any(p ∈ 𝒫 for p ∈ outputs(n)), 𝒩)
+    𝒩ⁱⁿ = filter(n -> any(p ∈ 𝒫 for p ∈ inputs(n)), 𝒩)
+
+    @variable(m,
+        lower_limit(p) ≤
+            energy_potential_node_out[n ∈ 𝒩ᵒᵘᵗ, 𝒯, p ∈ intersect(outputs(n), 𝒫)] ≤
+        upper_limit(p)
+    )
+    @variable(m,
+        lower_limit(p) ≤
+            energy_potential_node_in[n ∈ 𝒩ⁱⁿ, 𝒯, p ∈ intersect(inputs(n), 𝒫)] ≤
+        upper_limit(p)
+    )
+end
+
+function EMB.variables_flow_resource(
+    m,
+    ℒ::Vector{<:Link},
+    𝒫::Vector{<:PotentialPower},
+    𝒯,
+    modeltype::EnergyModel,
+)
+    ℒᵉᵖ = filter(l -> any(p ∈ 𝒫 for p ∈ EMB.link_res(l)), ℒ)
+    @variable(m, energy_potential_link_in[ℒᵉᵖ, 𝒯, 𝒫])
+    @variable(m, energy_potential_link_out[ℒᵉᵖ, 𝒯, 𝒫])
+end
+```
+
+### 4. Add resource-specific constraints
+
+Create a new method [`constraints_resource`](@ref) for custom node or link behavior.
+These methods can be either for the complete set of [`Node`](@ref EnergyModelsBase.Node) and [`Link`](@ref)s or alternatively for only a specified subset of nodes.
+If you only specify it for a subset of nodes, it is important that the new resource is only an `input` or `output` of this subset.
+
+```julia
+function EMB.constraints_resource(
+    m,
+    n::PotentialLossNode,
+    𝒯,
+    𝒫::Vector{<:PotentialPower},
+    modeltype::EnergyModel,
+)
+    𝒫ᵒᵘᵗ = filter(p -> p ∈ 𝒫, outputs(n))
+    𝒫ⁱⁿ = filter(p -> p ∈ 𝒫, inputs(n))
+
+    @constraint(m, [t ∈ 𝒯, p ∈ 𝒫ᵒᵘᵗ],
+        m[:energy_potential_node_out][n, t, p] ==
+            n.loss_factor * m[:energy_potential_node_in][n, t, p]
+    )
+end
+
+function EMB.constraints_resource(
+    m,
+    l::Link,
+    𝒯,
+    𝒫::Vector{<:PotentialPower},
+    modeltype::EnergyModel,
+)
+    𝒫ˡⁱⁿᵏ = filter(p -> p ∈ 𝒫, EMB.link_res(l))
+    @constraint(m, [t ∈ 𝒯, p ∈ 𝒫ˡⁱⁿᵏ],
+        m[:energy_potential_link_in][l, t, p] ==
+            m[:energy_potential_link_out][l, t, p]
+    )
+end
+```
+
+### 5. Couple node and link variables
+
+Use [`constraints_couple_resource`](@ref) to connect node and link resource variables.
+
+```julia
+function EMB.constraints_couple_resource(
+    m,
+    𝒩::Vector{<:EMB.Node},
+    ℒ::Vector{<:Link},
+    𝒫::Vector{<:PotentialPower},
+    𝒯,
+    modeltype::EnergyModel,
+)
+    for n ∈ 𝒩
+        ℒᶠʳᵒᵐ, ℒᵗᵒ = EMB.link_sub(ℒ, n)
+        𝒫ᵒᵘᵗ = filter(p -> p ∈ 𝒫, outputs(n))
+        𝒫ⁱⁿ = filter(p -> p ∈ 𝒫, inputs(n))
+
+        @constraint(m, [t ∈ 𝒯, p ∈ 𝒫ᵒᵘᵗ, l ∈ ℒᶠʳᵒᵐ],
+            m[:energy_potential_node_out][n, t, p] ==
+                m[:energy_potential_link_in][l, t, p]
+        )
+
+        @constraint(m, [t ∈ 𝒯, p ∈ 𝒫ⁱⁿ, l ∈ ℒᵗᵒ],
+            m[:energy_potential_link_out][l, t, p] ==
+                m[:energy_potential_node_in][n, t, p]
+        )
+    end
+end
+```
diff --git a/docs/src/how-to/utilize-timestruct.md b/docs/src/how-to/utilize-timestruct.md
index 85c5dc50..79c253da 100644
--- a/docs/src/how-to/utilize-timestruct.md
+++ b/docs/src/how-to/utilize-timestruct.md
@@ -39,7 +39,7 @@ op_number = length(op_duration)
 operational_periods = SimpleTimes(op_number, op_duration)
 
 # output
-SimpleTimes{Int64}(11, [4, 2, 1, 1, 2, 4, 2, 1, 1, 2, 4])
+SimpleTimes{Int64}(11, [4, 2, 1, 1, 2, 4, 2, 1, 1, 2, 4], 24)
 ```
 
 In this case, we model the day not with hourly resolution, but only have hourly resolution in the morning and afternoon.
@@ -60,7 +60,7 @@ Instead, one can also write
 operational_periods = SimpleTimes(op_duration)
 
 # output
-SimpleTimes{Int64}(11, [4, 2, 1, 1, 2, 4, 2, 1, 1, 2, 4])
+SimpleTimes{Int64}(11, [4, 2, 1, 1, 2, 4, 2, 1, 1, 2, 4], 24)
 ```
 
 and a constructor will automatically deduce that there have to be 11 operational periods.
diff --git a/docs/src/index.md b/docs/src/index.md
index b58b721f..6a964788 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -61,6 +61,7 @@ Depth = 1
 Pages = [
     "how-to/create_new_element.md",
     "how-to/create-new-node.md",
+    "how-to/extend-resource-functionality.md",
     "how-to/utilize-timestruct.md",
     "how-to/update-models.md",
     "how-to/contribute.md",
diff --git a/docs/src/library/internals/functions.md b/docs/src/library/internals/functions.md
index 09cf9c18..3c9d1883 100644
--- a/docs/src/library/internals/functions.md
+++ b/docs/src/library/internals/functions.md
@@ -28,6 +28,8 @@ emissions_operational
 constraints_emissions
 constraints_elements
 constraints_couple
+constraints_couple_resource
+constraints_resource
 constraints_level_iterate
 constraints_level_rp
 constraints_level_scp
@@ -39,6 +41,7 @@ constraints_level_bounds
 ```@docs
 variables_capacity
 variables_flow
+variables_flow_resource
 variables_opex
 variables_capex
 variables_emission
@@ -96,4 +99,6 @@ res_sub
 ```@docs
 collect_types
 sort_types
+res_types
+res_types_vec
 ```
diff --git a/src/model.jl b/src/model.jl
index 894f5243..c6e5431c 100644
--- a/src/model.jl
+++ b/src/model.jl
@@ -60,7 +60,7 @@ function create_model(
     # Declaration of element variables and constraints of the problem
     for 𝒳 ∈ 𝒳ᵛᵉᶜ
         variables_capacity(m, 𝒳, 𝒳ᵛᵉᶜ, 𝒯, modeltype)
-        variables_flow(m, 𝒳, 𝒳ᵛᵉᶜ, 𝒯, modeltype)
+        variables_flow(m, 𝒳, 𝒳ᵛᵉᶜ, 𝒫, 𝒯, modeltype)
         variables_opex(m, 𝒳, 𝒳ᵛᵉᶜ, 𝒯, modeltype)
         variables_capex(m, 𝒳, 𝒳ᵛᵉᶜ, 𝒯, modeltype)
         variables_emission(m, 𝒳, 𝒳ᵛᵉᶜ, 𝒫, 𝒯, modeltype)
@@ -203,12 +203,13 @@ function variables_capacity(m, ℒ::Vector{<:Link}, 𝒳ᵛᵉᶜ, 𝒯, modelty
 end
 
 """
+    variables_flow(m, 𝒳::Vector{<:AbstractElement}, 𝒳ᵛᵉᶜ, 𝒫, 𝒯, modeltype::EnergyModel)
     variables_flow(m, 𝒩::Vector{<:Node}, 𝒳ᵛᵉᶜ, 𝒯, modeltype::EnergyModel)
     variables_flow(m, ℒ::Vector{<:Link}, 𝒳ᵛᵉᶜ, 𝒯, modeltype::EnergyModel)
 
-Declaration of flow OPEX variables for the element types introduced in
-`EnergyModelsBase`. `EnergyModelsBase` introduces two elements for an energy system, and
-hence, provides the user with two individual methods:
+Declaration of flow variables for the element types introduced in `EnergyModelsBase`.
+`EnergyModelsBase` introduces two elements for an energy system, and hence, provides the
+user with two individual methods:
 
 !!! note "Node variables"
     - `flow_in[n, t, p]` is the flow _**into**_ node `n` in operational period `t` for
@@ -217,6 +218,8 @@ hence, provides the user with two individual methods:
     - `flow_out[n, t, p]` is the flow _**from**_ node `n` in operational period `t`
       for resource `p`. The outflow resources of node `n` are extracted using the
       function [`outputs`](@ref).
+    - call of the function [`variables_flow_resource`](@ref) for introducing resource
+      specific flow variables.
 
 !!! tip "Link variables"
     - `link_in[l, t, p]` is the flow _**into**_ link `l` in operational period `t` for
@@ -225,12 +228,17 @@ hence, provides the user with two individual methods:
     - `link_out[l, t, p]` is the flow _**from**_ link `l` in operational period `t`
       for resource `p`. The outflow resources of link `l` are extracted using the
       function [`outputs`](@ref).
+    - call of the function [`variables_flow_resource`](@ref) for introducing resource
+      specific flow variables.
 
 By default, all nodes `𝒩` and links `ℒ` only allow for unidirectional flow. You can specify
 bidirectional flow through providing a method to the function [`is_unidirectional`](@ref)
 for new link/node types.
+
+The fallback solution for `𝒳::Vector{<:AbstractElement}` is in the current stage included
+to maintain backwards compatibility for packages that introduce additional [`AbstractElement`](@ref)s.
 """
-function variables_flow(m, 𝒩::Vector{<:Node}, 𝒳ᵛᵉᶜ, 𝒯, modeltype::EnergyModel)
+function variables_flow(m, 𝒩::Vector{<:Node}, 𝒳ᵛᵉᶜ, 𝒫, 𝒯, modeltype::EnergyModel)
     # Extract the nodes with inputs and outputs
     𝒩ⁱⁿ = filter(has_input, 𝒩)
     𝒩ᵒᵘᵗ = filter(has_output, 𝒩)
@@ -249,8 +257,14 @@ function variables_flow(m, 𝒩::Vector{<:Node}, 𝒳ᵛᵉᶜ, 𝒯, modeltype:
     for n_out ∈ 𝒩ᵒᵘᵗ⁻ᵘⁿⁱ, t ∈ 𝒯, p ∈ outputs(n_out)
         set_lower_bound(m[:flow_out][n_out, t, p], 0)
     end
+
+    # Create new flow variables for specific resource types
+    for p_sub ∈ res_types_vec(𝒫)
+        variables_flow_resource(m, 𝒩, p_sub, 𝒯, modeltype)
+    end
+
 end
-function variables_flow(m, ℒ::Vector{<:Link}, 𝒳ᵛᵉᶜ, 𝒯, modeltype::EnergyModel)
+function variables_flow(m, ℒ::Vector{<:Link}, 𝒳ᵛᵉᶜ, 𝒫, 𝒯, modeltype::EnergyModel)
     # Create the link flow variables
     @variable(m, link_in[l ∈ ℒ, 𝒯, inputs(l)])
     @variable(m, link_out[l ∈ ℒ, 𝒯, outputs(l)])
@@ -266,7 +280,30 @@ function variables_flow(m, ℒ::Vector{<:Link}, 𝒳ᵛᵉᶜ, 𝒯, modeltype::
             set_lower_bound(m[:link_out][l, t, p], 0)
         end
     end
+
+    # Create new flow variables for specific resource types
+    for p_sub ∈ res_types_vec(𝒫)
+        variables_flow_resource(m, ℒ, p_sub, 𝒯, modeltype)
+    end
 end
+variables_flow(m, 𝒳::Vector{<:AbstractElement}, 𝒳ᵛᵉᶜ, 𝒫, 𝒯, modeltype::EnergyModel) =
+    variables_flow(m, 𝒳, 𝒳ᵛᵉᶜ, 𝒯, modeltype)
+
+"""
+    variables_flow_resource(m, ℒ::Vector{<:Link}, 𝒫::Vector{<:Resource}, 𝒯, modeltype::EnergyModel)
+    variables_flow_resource(m, 𝒩::Vector{<:Node}, 𝒫::Vector{<:Resource}, 𝒯, modeltype::EnergyModel)
+
+Create resource-specific flow variables for links or nodes.
+
+This function is called from [`variables_flow`](@ref) for each subset of resources
+sharing the same type. It can be used to add variables and bounds for specialized
+resource classes while keeping the default flow variables unchanged.
+
+The default methods are empty and intended to be implemented in extension packages.
+"""
+function variables_flow_resource(m, ℒ::Vector{<:Link}, 𝒫::Vector{<:Resource}, 𝒯, modeltype::EnergyModel) end
+function variables_flow_resource(m, 𝒩::Vector{<:Node}, 𝒫::Vector{<:Resource}, 𝒯, modeltype::EnergyModel) end
+
 
 """
     variables_opex(m, 𝒩::Vector{<:Node}, 𝒳ᵛᵉᶜ, 𝒯, modeltype::EnergyModel)
@@ -562,9 +599,9 @@ end
     create_element(m, n::Node, 𝒯, 𝒫, modeltype::EnergyModel)
     create_element(m, l::Link, 𝒯, 𝒫, modeltype::EnergyModel)
 
-Default fallback method for an element type if no other method is defined for a given type.
-This function calls subfunctions to maintain backwards compatibility and simplify the
-differentiation in extension packages.
+Calls the create functions for the specific elements to add element specific constraints (by
+calling individual subfunctions) and add resource specific constraints by calling
+[`constraints_resource`](@ref).
 
 `EnergyModelsBase` provides the user with two element types, [`Link`](@ref) and
 [`Node`](@ref EnergyModelsBase.Node):
@@ -572,11 +609,41 @@ differentiation in extension packages.
 - `Node` - the subfunction is [`create_node`](@ref).
 - `Link` - the subfunction is [`create_link`](@ref).
 """
-create_element(m, n::Node, 𝒯, 𝒫, modeltype::EnergyModel) =
+function create_element(m, n::Node, 𝒯, 𝒫, modeltype::EnergyModel)
+
     create_node(m, n, 𝒯, 𝒫, modeltype)
-create_element(m, l::Link, 𝒯, 𝒫, modeltype::EnergyModel) =
+
+    # Constraints based on the resource types
+    node_resources = Vector{Resource}(unique(vcat(inputs(n), outputs(n))))
+    for 𝒫ˢᵘᵇ ∈ res_types_vec(node_resources)
+        constraints_resource(m, n, 𝒯, 𝒫ˢᵘᵇ, modeltype)
+    end
+end
+
+function create_element(m, l::Link, 𝒯, 𝒫, modeltype::EnergyModel)
+
     create_link(m, l, 𝒯, 𝒫, modeltype)
 
+    # Constraints based on the resource types
+    for 𝒫ˢᵘᵇ ∈ res_types_vec(link_res(l))
+        constraints_resource(m, l, 𝒯, 𝒫ˢᵘᵇ, modeltype)
+    end
+end
+
+"""
+    constraints_resource(m, n::Node, 𝒯, 𝒫::Vector{<:Resource}, modeltype::EnergyModel)
+    constraints_resource(m, l::Link, 𝒯, 𝒫::Vector{<:Resource}, modeltype::EnergyModel)
+
+Create constraints for the flow of resources through an [`AbstractElement`](@ref) for
+specific resource types. In `EnergyModelsBase`, this method is provided for
+[`Node`](@ref EnergyModelsBase.Node) and [`Link`](@ref).
+
+The function is empty by default and can be implemented in extension packages.
+"""
+function constraints_resource(m, n::Node, 𝒯, 𝒫::Vector{<:Resource}, modeltype::EnergyModel) end
+
+function constraints_resource(m, l::Link, 𝒯, 𝒫::Vector{<:Resource}, modeltype::EnergyModel) end
+
 """
     constraints_couple(m, 𝒩::Vector{<:Node}, ℒ::Vector{<:Link}, 𝒫, 𝒯, modeltype::EnergyModel)
     constraints_couple(m, ℒ::Vector{<:Link}, 𝒩::Vector{<:Node}, 𝒫, 𝒯, modeltype::EnergyModel)
@@ -590,6 +657,7 @@ for the coupling between a [`Link`](@ref) and a [`Node`](@ref EnergyModelsBase.N
 function constraints_couple(m, 𝒩::Vector{<:Node}, ℒ::Vector{<:Link}, 𝒫, 𝒯, modeltype::EnergyModel)
     for n ∈ 𝒩
         ℒᶠʳᵒᵐ, ℒᵗᵒ = link_sub(ℒ, n)
+
         # Constraint for output flowrate and input links.
         if has_output(n)
             @constraint(m, [t ∈ 𝒯, p ∈ outputs(n)],
@@ -597,6 +665,7 @@ function constraints_couple(m, 𝒩::Vector{<:Node}, ℒ::Vector{<:Link}, 𝒫,
                 sum(m[:link_in][l, t, p] for l ∈ ℒᶠʳᵒᵐ if p ∈ inputs(l))
             )
         end
+
         # Constraint for input flowrate and output links.
         if has_input(n)
             @constraint(m, [t ∈ 𝒯, p ∈ inputs(n)],
@@ -605,11 +674,29 @@ function constraints_couple(m, 𝒩::Vector{<:Node}, ℒ::Vector{<:Link}, 𝒫,
             )
         end
     end
+
+    # Create new constraints for specific resource types
+    for p_sub ∈ res_types_vec(𝒫)
+        constraints_couple_resource(m, 𝒩, ℒ, p_sub, 𝒯, modeltype)
+    end
 end
 function constraints_couple(m, ℒ::Vector{<:Link}, 𝒩::Vector{<:Node}, 𝒫, 𝒯, modeltype::EnergyModel)
     return constraints_couple(m, 𝒩, ℒ, 𝒫, 𝒯, modeltype)
 end
 
+"""
+    constraints_couple_resource(m, 𝒩::Vector{<:Node}, ℒ::Vector{<:Link}, 𝒫::Vector{<:Resource}, 𝒯, modeltype::EnergyModel)
+
+Create resource-specific coupling constraints between nodes and links.
+
+This function is called from [`constraints_couple`](@ref) for each subset of resources
+sharing the same type. It can be used to add additional coupling constraints for
+specialized resource classes while keeping the default node-link flow balance unchanged.
+
+The default method is empty and intended to be implemented in extension packages.
+"""
+function constraints_couple_resource(m, 𝒩::Vector{<:Node}, ℒ::Vector{<:Link}, 𝒫::Vector{<:Resource}, 𝒯, modeltype::EnergyModel) end
+
 """
     constraints_emissions(m, 𝒳ᵛᵉᶜ, 𝒫, 𝒯, modeltype::EnergyModel)
 
@@ -956,7 +1043,6 @@ function create_link(m, l::Direct, 𝒯, 𝒫, modeltype::EnergyModel)
     )
 end
 function create_link(m, 𝒯, 𝒫, l::Link, modeltype::EnergyModel, formulation::Formulation)
-
     # Generic link in which each output corresponds to the input
     @constraint(m, [t ∈ 𝒯, p ∈ link_res(l)],
         m[:link_out][l, t, p] == m[:link_in][l, t, p]
diff --git a/src/structures/resource.jl b/src/structures/resource.jl
index 14fa0197..b60fa43a 100644
--- a/src/structures/resource.jl
+++ b/src/structures/resource.jl
@@ -86,3 +86,17 @@ Returns all emission resources for a
 """
 res_em(𝒫::Array{<:Resource}) = filter(is_resource_emit, 𝒫)
 res_em(𝒫::Dict) = filter(p -> is_resource_emit(first(p)), 𝒫)
+
+"""
+    res_types(𝒫::Vector{<:Resource})
+
+Return the unique resource types in an Vector of resources `𝒫`.
+"""
+res_types(𝒫::Vector{<:Resource}) = unique(map(x -> typeof(x), 𝒫))
+
+"""
+    res_types_vec(𝒫::Vector{<:Resource})
+
+Return a Vector-of-Vectors of resources by the concrete sub-types, if the input is empty it returns an empty Vector.
+"""
+res_types_vec(𝒫::Vector{<:Resource}) = [Vector{rt}(filter(x -> isa(x, rt), 𝒫)) for rt in res_types(𝒫)]
\ No newline at end of file
diff --git a/test/runtests.jl b/test/runtests.jl
index b12417d6..898f27de 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -21,6 +21,10 @@ ENV["EMB_TEST"] = true # Set flag for example scripts to check if they are run a
         include("test_data.jl")
     end
 
+    @testset "Base | Resource" begin
+        include("test_resource.jl")
+    end
+
     @testset "Base | Node" begin
         include("test_nodes.jl")
     end
diff --git a/test/test_resource.jl b/test/test_resource.jl
new file mode 100644
index 00000000..25a1bf1c
--- /dev/null
+++ b/test/test_resource.jl
@@ -0,0 +1,270 @@
+@testset "Resource - utilities" begin
+    # Declare the resources
+    Power = ResourceCarrier("Power", 0.0)
+    Heat = ResourceCarrier("Heat", 0.0)
+    CO2 = ResourceEmit("CO2", 1.0)
+
+    𝒫 = [Power, Heat, CO2]
+    @testset "General" begin
+        # returns a Vector of DataTypes
+        @test EMB.res_types(𝒫) isa Vector{DataType}
+
+        # returns the correct number of unique resource types
+        @test length(EMB.res_types(𝒫)) == 2
+
+        # returns a Vector
+        @test EMB.res_types_vec(𝒫) isa Vector{Vector}
+
+        # returns the correct number of segments
+        @test length(EMB.res_types_vec(𝒫)) == 2
+
+        # the length of the first segment should be 2 (2 ResourceCarriers)
+        @test length(EMB.res_types_vec(𝒫)[1]) == 2
+
+        # the length of the second segment should be 1 (1 ResourceEmit)
+        @test length(EMB.res_types_vec(𝒫)[2]) == 1
+
+        # returns an empty vector when given an empty resource vector
+        @test isempty(EMB.res_types_vec(Resource[]))
+    end
+    @testset "Resource with parameters" begin
+        struct TestResource <: Resource
+            id::String
+            a::Float64
+            b::Int64
+        end
+
+        # Add a new resource of type TestResource to the resource vector
+        push!(𝒫, TestResource("Test", 0.5, 1))
+
+        # returns a Vector of DataTypes (now including TestResource)
+        @test isa(EMB.res_types(𝒫), Vector{DataType})
+
+        # returns the correct number of unique resource types (now 3)
+        @test length(EMB.res_types(𝒫)) == 3
+
+        # returns the correct number of segments (now 3)
+        @test length(EMB.res_types_vec(𝒫)) == 3
+
+        # the length of the first segment should be 2 (2 ResourceCarriers)
+        @test length(EMB.res_types_vec(𝒫)[1]) == 2
+
+        # the length of the second segment should be 1 (1 ResourceEmit)
+        @test length(EMB.res_types_vec(𝒫)[2]) == 1
+
+        # the length of the third segment should be 1 (1 TestResource)
+        @test length(EMB.res_types_vec(𝒫)[3]) == 1
+    end
+end
+
+# Implement a custom resource type and check that it is correctly handled in the model via dispatch
+@testset "Resource - implementation" begin
+        struct PotentialPower <: Resource
+            id::String
+            co2_int::Float64
+            potential_lower::Float64
+            potential_upper::Float64
+        end
+        EMB.is_resource_emit(::PotentialPower) = false
+        lower_limit(p::PotentialPower) = p.potential_lower
+        upper_limit(p::PotentialPower) = p.potential_upper
+
+        # A costum node type that represents a potential loss node
+        # which has an input and output resource and a loss factor that determines how much
+        # of the input potential is lost in the node but there is no loss in energy
+        struct PotentialLossNode{T <: PotentialPower} <: NetworkNode
+            id::Any
+            cap::TimeProfile
+            opex_var::TimeProfile
+            opex_fixed::TimeProfile
+            resource::T
+            input::Dict{<:Resource,<:Real}
+            output::Dict{<:Resource,<:Real}
+            data::Vector{<:ExtensionData}
+            loss_factor::Float64
+        end
+        function PotentialLossNode(
+            id,
+            cap::TimeProfile,
+            opex_var::TimeProfile,
+            opex_fixed::TimeProfile,
+            resource::T,
+            loss_factor::Float64,
+        )  where {T <: PotentialPower}
+            return PotentialLossNode{T}(id, cap, opex_var, opex_fixed, resource, Dict(resource=>1.0), Dict(resource=>1.0), ExtensionData[], loss_factor)
+        end
+
+
+    function res_test_case(loss_factor::Float64)
+        pp = PotentialPower("PotentialPower", 0.0, 0.9, 1.1)
+        CO2 = ResourceEmit("CO2", 1.0)
+        source = RefSource(
+            "pp_source",
+            FixedProfile(4),
+            FixedProfile(10),
+            FixedProfile(0),
+            Dict(pp => 1),
+        )
+        loss_node = PotentialLossNode(
+            "pp_loss",
+            FixedProfile(4),
+            FixedProfile(0),
+            FixedProfile(0),
+            pp,
+            loss_factor,
+        )
+        sink = RefSink(
+            "pp_sink",
+            FixedProfile(3),
+            Dict(:surplus => FixedProfile(4), :deficit => FixedProfile(100)),
+            Dict(pp => 1),
+        )
+
+        𝒯 = TwoLevel(2, 2, SimpleTimes(5, 2); op_per_strat = 10)
+        𝒩 = [source, loss_node, sink]
+        ℒ = [
+            Direct("src-loss", source, loss_node, Linear())
+            Direct("loss-snk", loss_node, sink, Linear())
+        ]
+        modeltype = OperationalModel(
+            Dict(CO2 => FixedProfile(100)),
+            Dict(CO2 => FixedProfile(0)),
+            CO2,
+        )
+        case = Case(𝒯, [pp, CO2], [𝒩, ℒ])
+
+        return case, modeltype
+    end
+
+    # Delcare new variables for the potential power resource
+    function EMB.variables_flow_resource(
+        m, 𝒩::Vector{<:EMB.Node}, 𝒫::Vector{<:PotentialPower}, 𝒯, modeltype::EnergyModel
+    )
+        𝒩ᵒᵘᵗ = filter(n -> any(p ∈ 𝒫 for p ∈ outputs(n)), 𝒩)
+        𝒩ⁱⁿ = filter(n -> any(p ∈ 𝒫 for p ∈ inputs(n)), 𝒩)
+
+        @variable(m,
+            lower_limit(p) ≤
+                energy_potential_node_out[n ∈ 𝒩ᵒᵘᵗ, 𝒯, p ∈ intersect(outputs(n), 𝒫)] ≤
+            upper_limit(p)
+        )
+        @variable(m,
+            lower_limit(p) ≤
+                energy_potential_node_in[n ∈ 𝒩ⁱⁿ, 𝒯, p ∈ intersect(inputs(n), 𝒫)] ≤
+            upper_limit(p)
+        )
+    end
+
+    function EMB.variables_flow_resource(
+        m, ℒ::Vector{<:Link}, 𝒫::Vector{<:PotentialPower}, 𝒯, modeltype::EnergyModel
+    )
+        ℒᵉᵖ = filter(l -> any(p ∈ 𝒫 for p ∈ EMB.link_res(l)), ℒ)
+        @variable(m, energy_potential_link_in[ℒᵉᵖ, 𝒯, 𝒫])
+        @variable(m, energy_potential_link_out[ℒᵉᵖ, 𝒯, 𝒫])
+    end
+
+    # Declare new constraints for the potential power resource using the newly declared variables
+    function EMB.constraints_resource(
+        m, n::PotentialLossNode, 𝒯, 𝒫::Vector{<:PotentialPower}, modeltype::EnergyModel
+    )
+        𝒫ᵒᵘᵗ = filter(p -> p ∈ 𝒫, outputs(n))
+
+        @constraint(m, [t ∈ 𝒯, p ∈ 𝒫ᵒᵘᵗ],
+            m[:energy_potential_node_out][n, t, p] == n.loss_factor * m[:energy_potential_node_in][n, t, p]
+        )
+    end
+
+    function EMB.constraints_resource(
+        m, l::Link, 𝒯, 𝒫::Vector{<:PotentialPower}, modeltype::EnergyModel
+    )
+        𝒫ˡⁱⁿᵏ = filter(p -> p ∈ 𝒫, EMB.link_res(l))
+        @constraint(m, [t ∈ 𝒯, p ∈ 𝒫ˡⁱⁿᵏ],
+            m[:energy_potential_link_in][l, t, p] == m[:energy_potential_link_out][l, t, p]
+        )
+    end
+
+    function EMB.constraints_couple_resource(
+        m, 𝒩::Vector{<:EMB.Node}, ℒ::Vector{<:Link},
+        𝒫::Vector{<:PotentialPower}, 𝒯, modeltype::EnergyModel
+    )
+        for n ∈ 𝒩
+            ℒᶠʳᵒᵐ, ℒᵗᵒ = EMB.link_sub(ℒ, n)
+            𝒫ᵒᵘᵗ = filter(p -> p ∈ 𝒫, outputs(n))
+            𝒫ⁱⁿ = filter(p -> p ∈ 𝒫, inputs(n))
+
+            @constraint(m, [t ∈ 𝒯, p ∈ 𝒫ᵒᵘᵗ, l ∈ ℒᶠʳᵒᵐ],
+                m[:energy_potential_node_out][n, t, p] == m[:energy_potential_link_in][l, t, p]
+            )
+
+            @constraint(m, [t ∈ 𝒯, p ∈ 𝒫ⁱⁿ, l ∈ ℒᵗᵒ],
+                m[:energy_potential_link_out][l, t, p] == m[:energy_potential_node_in][n, t, p]
+            )
+        end
+    end
+
+
+    case, modeltype = res_test_case(0.9)
+    pp, co2 = get_products(case)
+    source, loss_node, sink = get_nodes(case)
+
+    m = run_model(case, modeltype, HiGHS.Optimizer)
+    𝒯 = get_time_struct(case)
+    ℒ = get_links(case)
+    n_t = length(𝒯)
+
+    # Variable testing (calling of the correct function)
+    # - variables_flow(m, 𝒩::Vector{<:Node}, 𝒳ᵛᵉᶜ, 𝒫, 𝒯, modeltype::EnergyModel)
+    # - variables_flow(m, ℒ::Vector{<:Link}, 𝒳ᵛᵉᶜ, 𝒫, 𝒯, modeltype::EnergyModel)
+    # Check that the variables are created
+    @test haskey(m, :energy_potential_node_in)
+    @test haskey(m, :energy_potential_node_out)
+    @test haskey(m, :energy_potential_link_in)
+    @test haskey(m, :energy_potential_link_out)
+
+    ## Check that the variables have the correct length
+    @test length(m[:energy_potential_node_in]) == 2 * n_t
+    @test length(m[:energy_potential_node_out]) == 2 * n_t
+    @test length(m[:energy_potential_link_in]) == length(ℒ) * n_t
+    @test length(m[:energy_potential_link_out]) == length(ℒ) * n_t
+
+    ## Check that the bounds of the variables are enforced
+    @test all(value(m[:energy_potential_node_out][source, t, pp]) ≥ lower_limit(pp) for t ∈ 𝒯)
+    @test all(value(m[:energy_potential_node_out][source, t, pp]) ≤ upper_limit(pp) for t ∈ 𝒯)
+    @test all(value(m[:energy_potential_node_in][sink, t, pp]) ≥ lower_limit(pp) for t ∈ 𝒯)
+    @test all(value(m[:energy_potential_node_in][sink, t, pp]) ≤ upper_limit(pp) for t ∈ 𝒯)
+
+    # Test that the coupling constraints are correctly enforced
+    # - EMB.constraints_couple_resource
+    @test all(
+        value(m[:energy_potential_node_out][source, t, pp]) ≈
+            value(m[:energy_potential_link_in][ℒ[1], t, pp])
+    for t ∈ 𝒯)
+    @test all(
+        value(m[:energy_potential_link_out][ℒ[1], t, pp]) ≈
+            value(m[:energy_potential_node_in][loss_node, t, pp])
+    for t ∈ 𝒯)
+    @test all(
+        value(m[:energy_potential_node_out][loss_node, t, pp]) ≈
+            loss_node.loss_factor * value(m[:energy_potential_node_in][loss_node, t, pp])
+    for t ∈ 𝒯)
+    @test all(
+        value(m[:energy_potential_node_out][loss_node, t, pp]) ≈
+            value(m[:energy_potential_link_in][ℒ[2], t, pp])
+    for t ∈ 𝒯)
+    @test all(
+        value(m[:energy_potential_link_out][ℒ[2], t, pp]) ≈
+            value(m[:energy_potential_node_in][sink, t, pp])
+    for t ∈ 𝒯)
+    @test all(
+        value(m[:energy_potential_node_out][loss_node, t, pp]) <
+            value(m[:energy_potential_node_in][loss_node, t, pp])
+    for t ∈ 𝒯)
+    @test all(
+        value(m[:energy_potential_node_out][source, t, pp]) <
+            value(m[:flow_out][source, t, pp])
+    for t ∈ 𝒯)
+    @test all(
+        value(m[:energy_potential_node_in][sink, t, pp]) <
+            value(m[:flow_in][sink, t, pp])
+    for t ∈ 𝒯)
+end