diff --git a/docs/src/assets/tnrkit.bib b/docs/src/assets/tnrkit.bib
index 98e6bb59..59abde24 100644
--- a/docs/src/assets/tnrkit.bib
+++ b/docs/src/assets/tnrkit.bib
@@ -101,3 +101,30 @@ @article{yangLoopOptimizationTensor2017
   doi       = {10.1103/PhysRevLett.118.110504},
   keywords  = {TNR Algorithm}
 }
+
+@article{kadohTensorNetworkAnalysis2019,
+  title = {Tensor network analysis of critical coupling in two dimensional $\varphi^{4}$ theory},
+  author = {Kadoh, Daisuke and Kuramashi, Yoshinobu and Nakamura, Yoshifumi and Sakai, Ryo and Takeda, Shinji and Yoshimura, Yusuke},
+  year = 2019,
+  month = may,
+  journal = {Journal of High Energy Physics},
+  volume = {2019},
+  number = {5},
+  eprint = {1811.12376},
+  primaryclass = {hep-lat},
+  pages = {184},
+  issn = {1029-8479},
+  doi = {10.1007/JHEP05(2019)184}
+}
+
+@article{moritaMultiimpurityMethodBondweighted2024,
+  title = {Multi-Impurity Method for the Bond-Weighted Tensor Renormalization Group},
+  author = {Morita, Satoshi and Kawashima, Naoki},
+  year = 2024,
+  month = nov,
+  number = {arXiv:2411.13998},
+  eprint = {2411.13998},
+  primaryclass = {cond-mat},
+  publisher = {arXiv},
+  doi = {10.48550/arXiv.2411.13998}
+}
diff --git a/src/TNRKit.jl b/src/TNRKit.jl
index a6c4b36f..180fe740 100644
--- a/src/TNRKit.jl
+++ b/src/TNRKit.jl
@@ -20,7 +20,6 @@ include("schemes/hotrg.jl")
 include("schemes/hotrg3d.jl")
 include("schemes/atrg.jl")
 include("schemes/atrg3d.jl")
-include("schemes/impurityhotrg.jl")
 # CTM methods
 include("schemes/ctm/utility.jl")
 include("schemes/ctm/c4ctm.jl")
@@ -30,6 +29,10 @@ include("schemes/ctm/ctm_hotrg.jl")
 include("schemes/ctm/onesite_ctm.jl")
 include("schemes/ctm/sublattice_ctm.jl")
 
+# Impurity methods
+include("schemes/impuritytrg.jl")
+include("schemes/impurityhotrg.jl")
+
 # Loop Methods
 include("schemes/looptnr.jl")
 include("schemes/symmetric_looptnr.jl")
@@ -43,7 +46,6 @@ export HOTRG
 export HOTRG_3D
 export ATRG
 export ATRG_3D
-export ImpurityHOTRG
 
 export CTM
 export Sublattice_CTM
@@ -53,6 +55,9 @@ export ctm_TRG
 export ctm_HOTRG
 export lnz
 
+export ImpurityTRG
+export ImpurityHOTRG
+
 export LoopTNR
 export SLoopTNR
 
diff --git a/src/schemes/impuritytrg.jl b/src/schemes/impuritytrg.jl
new file mode 100644
index 00000000..b9df3c9f
--- /dev/null
+++ b/src/schemes/impuritytrg.jl
@@ -0,0 +1,88 @@
+"""
+$(TYPEDEF)
+
+Impurity method for Tensor Renormalization Group
+
+### Constructors
+    $(FUNCTIONNAME)(T, T_imp1, T_imp2, T_imp3, T_imp4 [, finalize=finalize!])
+
+### Running the algorithm
+    run!(::ImpurityTRG, trunc::TensorKit.TruncationSheme, stop::Stopcrit[, finalize_beginning=true, verbosity=1])
+
+Each step rescales the lattice by a (linear) factor of 2
+
+!!! info "verbosity levels"
+    - 0: No output
+    - 1: Print information at start and end of the algorithm
+    - 2: Print information at each step
+
+### Fields
+
+$(TYPEDFIELDS)
+
+### References
+* [Morita et. al. 10.48550/arXiv.2411.13998 (2024)](@cite moritaMultiimpurityMethodBondweighted2024)
+* [Nakamoto et. al. 10.1007/JHEP05(2019)184 (2019)](@cite kadohTensorNetworkAnalysis2019)
+"""
+mutable struct ImpurityTRG <: TNRScheme
+    T::TensorMap
+    T_imp1::TensorMap
+    T_imp2::TensorMap
+    T_imp3::TensorMap
+    T_imp4::TensorMap
+
+    finalize!::Function
+    function ImpurityTRG(
+            T::TensorMap{E, S, 2, 2}, T_imp1::TensorMap{E, S, 2, 2}, T_imp2::TensorMap{E, S, 2, 2},
+            T_imp3::TensorMap{E, S, 2, 2}, T_imp4::TensorMap{E, S, 2, 2}; finalize = (finalize!)
+        ) where {E, S}
+
+
+        @assert space(T, 1) == space(T_imp1, 1) == space(T_imp2, 1) == space(T_imp3, 1) == space(T_imp4, 1) "First space of T, T_imp1, T_imp2, T_imp3 and T_imp4 must be the same"
+        @assert space(T, 2) == space(T_imp1, 2) == space(T_imp2, 2) == space(T_imp3, 2) == space(T_imp4, 2) "Second space of T, T_imp1, T_imp2, T_imp3 and T_imp4 must be the same"
+        @assert space(T, 3) == space(T_imp1, 3) == space(T_imp2, 3) == space(T_imp3, 3) == space(T_imp4, 3) "Third space of T, T_imp1, T_imp2, T_imp3 and T_imp4 must be the same"
+        @assert space(T, 4) == space(T_imp1, 4) == space(T_imp2, 4) == space(T_imp3, 4) == space(T_imp4, 4) "Fourth space of T, T_imp1, T_imp2, T_imp3 and T_imp4 must be the same"
+
+        return new(T, T_imp1, T_imp2, T_imp3, T_imp4, finalize)
+    end
+end
+
+function step!(scheme::ImpurityTRG, trunc::TensorKit.TruncationScheme)
+    # Tensor1
+    A1, B1 = SVD12(scheme.T_imp1, trunc)
+
+    # Tensor2
+    tensor2p = transpose(scheme.T_imp2, ((2, 4), (1, 3)))
+    C2, D2 = SVD12(tensor2p, trunc)
+
+    # Tensor3
+    A3, B3 = SVD12(scheme.T_imp3, trunc)
+
+    # Tensor4
+    tensor4p = transpose(scheme.T_imp4, ((2, 4), (1, 3)))
+    C4, D4 = SVD12(tensor4p, trunc)
+
+    # Tensor pure
+    Ap, Bp = SVD12(scheme.T, trunc)
+    tensorpurep = transpose(scheme.T, ((2, 4), (1, 3)))
+    Cp, Dp = SVD12(tensorpurep, trunc)
+
+    # Contract
+    @planar scheme.T[-1, -2; -3, -4] := Dp[-2; 1 2] * Bp[-1; 4 1] * Cp[4 3; -3] * Ap[3 2; -4]
+    @planar scheme.T_imp1[-1, -2; -3, -4] := Dp[-2; 1 2] * B1[-1; 4 1] * C2[4 3; -3] * Ap[3 2; -4]
+    @planar scheme.T_imp2[-1, -2; -3, -4] := D2[-2; 1 2] * B3[-1; 4 1] * Cp[4 3; -3] * Ap[3 2; -4]
+    @planar scheme.T_imp3[-1, -2; -3, -4] := D4[-2; 1 2] * Bp[-1; 4 1] * Cp[4 3; -3] * A3[3 2; -4]
+    @planar scheme.T_imp4[-1, -2; -3, -4] := Dp[-2; 1 2] * Bp[-1; 4 1] * C4[4 3; -3] * A1[3 2; -4]
+
+    return scheme
+end
+
+function Base.show(io::IO, scheme::ImpurityTRG)
+    println(io, "ImpurityTRG - Impurity TRG")
+    println(io, "  * T: $(summary(scheme.T))")
+    println(io, "  * T_imp1: $(summary(scheme.T_imp1))")
+    println(io, "  * T_imp2: $(summary(scheme.T_imp2))")
+    println(io, "  * T_imp3: $(summary(scheme.T_imp3))")
+    println(io, "  * T_imp4: $(summary(scheme.T_imp4))")
+    return nothing
+end
diff --git a/src/utility/finalize.jl b/src/utility/finalize.jl
index c9435758..78fb2f06 100644
--- a/src/utility/finalize.jl
+++ b/src/utility/finalize.jl
@@ -71,6 +71,22 @@ function finalize!(scheme::SLoopTNR)
     return tr_norm^0.25
 end
 
+# finalize! for ImpurityTRG
+function finalize!(scheme::ImpurityTRG)
+    # First normalize everything by the pure tensor
+    npure = norm(@tensor scheme.T[1 2; 2 1])
+    scheme.T_imp1 /= npure
+    scheme.T_imp2 /= npure
+    scheme.T_imp3 /= npure
+    scheme.T_imp4 /= npure
+    scheme.T /= npure
+
+    # Then calculate the contracted/traced 4 impurity tensors
+    nimp = norm(@tensoropt scheme.T_imp1[5 4;6 1] * scheme.T_imp2[1 2;7 5] * scheme.T_imp3[3 7;2 8] * scheme.T_imp4[8 6;4 3])
+
+    return npure, nimp
+end
+
 # finalize! for ImpurityHOTRG
 function finalize!(scheme::ImpurityHOTRG)
     n = norm(@tensor scheme.T[1 2; 2 1])
diff --git a/test/schemes.jl b/test/schemes.jl
index 61d9164e..b6644af4 100644
--- a/test/schemes.jl
+++ b/test/schemes.jl
@@ -250,3 +250,43 @@ end
     m2_lowT = data[end][4] / data[end][1]
     @test m2_lowT ≈ 1 rtol = 1.0e-2
 end
+
+# ImpurityTRG
+@testset "ImpurityTRG - Ising Model" begin
+    T = classical_ising()
+    T_imp = classical_ising_impurity()
+
+    scheme = ImpurityTRG(T, T_imp, T, T, T)
+
+    data = run!(scheme, truncdim(24), maxiter(25))
+
+    @test free_energy(getindex.(data, 1), ising_βc) ≈ f_onsager rtol = 2.0e-6
+end
+
+@testset "ImpurityTRG - Magnetisation" begin
+    # High T
+    β = 0.1
+
+    T = classical_ising(β)
+    T_imp = classical_ising_impurity(β)
+
+    scheme = ImpurityTRG(T, T_imp, T, T, T)
+
+    data = run!(scheme, truncdim(16), maxiter(25))
+
+    m_expection = data[end][2] / data[end][1]
+    @test m_expection ≈ 0.0 atol = 1.0e-4
+
+    # Low T
+    β = 2
+
+    T = classical_ising(β; h = 1.0e-6)
+    T_imp = classical_ising_impurity(β; h = 1.0e-6)
+
+    scheme = ImpurityTRG(T, T_imp, T, T, T)
+
+    data = run!(scheme, truncdim(16), maxiter(25))
+
+    m_expection = data[end][2] / data[end][1]
+    @test m_expection ≈ 1.0 rtol = 1.0e-4
+end