diff --git a/src/jop_prop2DAcoIsoDenQ_DEO2_FDTD.jl b/src/jop_prop2DAcoIsoDenQ_DEO2_FDTD.jl
index 894a368..6baebae 100644
--- a/src/jop_prop2DAcoIsoDenQ_DEO2_FDTD.jl
+++ b/src/jop_prop2DAcoIsoDenQ_DEO2_FDTD.jl
@@ -33,6 +33,7 @@ function JetProp2DAcoIsoDenQ_DEO2_FDTD(;
         wavelet = WaveletCausalRicker(f=5.0),
         freesurface = false,
         imgcondition = "standard",
+        RTM_weight = 0.5,
         nthreads = Sys.CPU_THREADS,
         reportinterval = 500)
     # active and passive earth model properties.  The active set is in the model-space
@@ -126,10 +127,11 @@ function JetProp2DAcoIsoDenQ_DEO2_FDTD(;
     icdict = Dict(
         lowercase("standard") => WaveFD.ImagingConditionStandard(),
         lowercase("FWI") => WaveFD.ImagingConditionWaveFieldSeparationFWI(),
-        lowercase("RTM") => WaveFD.ImagingConditionWaveFieldSeparationRTM())
+        lowercase("RTM") => WaveFD.ImagingConditionWaveFieldSeparationRTM(),
+        lowercase("MIX") => WaveFD.ImagingConditionWaveFieldSeparationMIX())
 
     if lowercase(imgcondition) ∉ keys(icdict)
-        error("Supplied imaging condition 'imgcondition' is not in [standard, FWI, RTM]")
+        error("Supplied imaging condition 'imgcondition' is not in [standard, FWI, RTM, MIX]")
     end
         
     Jet(
@@ -170,6 +172,7 @@ function JetProp2DAcoIsoDenQ_DEO2_FDTD(;
             nbx_cache = nbx_cache,
             wavelet = wavelet,
             freesurface = freesurface,
+            RTM_weight,
             imgcondition = get(icdict, lowercase(imgcondition), WaveFD.ImagingConditionStandard()),
             nthreads = nthreads,
             reportinterval = reportinterval,
@@ -356,7 +359,11 @@ Defaults for arguments are shown inside square brackets.
 * `imgcondition` ["standard"] Selects the type of imaging condition used. Choose from "standard", "FWI", 
     and "RTM". "FWI" and "RTM" will perform Kz wavenumber filtering prior to the imaging condition
     in order to promote long wavelengths (for FWI), or remove long wavelength backscattered energy (for 
-    RTM). Note the true adjoint only exists for "standard" imaging condition currently.
+    RTM).  "MIX" mixes the FWI imaging condition using the parameter RTM_weight. Note the true adjoint 
+    only exists for "standard" imaging condition currently.
+* `RTM_weight` determines the balance of short wavelengths and long wavelengths in the imaging condition. 
+    A value of 0.0 is equivalent to the FWI imaging condition, a value of 0.5 is equivalent to the standard
+    imaging condtion, and a value of 1.0 is equivalent to the RTM imaging condition. 
 * `nthreads [Sys.CPU_THREADS]` The number of threads to use for OpenMP parallelization of the modeling.
 * `reportinterval [500]` The interval at which information about the propagtion is logged.
 
@@ -722,10 +729,21 @@ function JopProp2DAcoIsoDenQ_DEO2_FDTD_df′!(δm::AbstractArray, δd::AbstractA
     end
 
     δm_ginsu = Dict{String,Array{Float32,2}}()
+
+        # Temporary model updates if mixed imaging condition used. 
+        if kwargs[:imgcondition] ===  WaveFD.ImagingConditionWaveFieldSeparationMIX()
+        δm_RTM = Dict{String,Array{Float32,2}}()
+        δm_All = Dict{String,Array{Float32,2}}()
+    end
+
     for prop in keys(kwargs[:active_modelset])
         δm_ginsu[prop] = zeros(Float32, nz_ginsu, nx_ginsu)
+        if kwargs[:imgcondition] ===  WaveFD.ImagingConditionWaveFieldSeparationMIX()
+            δm_RTM[prop] = zeros(Float32, nz_ginsu, nx_ginsu)
+            δm_All[prop] = zeros(Float32, nz_ginsu, nx_ginsu)
+        end
     end
-
+    
     ginsu_interior_range = interior(kwargs[:ginsu])
 
     # Get receiver interpolation coefficients
@@ -797,7 +815,24 @@ function JopProp2DAcoIsoDenQ_DEO2_FDTD_df′!(δm::AbstractArray, δd::AbstractA
             end
 
             # born accumulation
-            cumtime_im += @elapsed WaveFD.adjointBornAccumulation!(p, kwargs[:modeltype], kwargs[:imgcondition], δm_ginsu, wavefields)
+            if kwargs[:imgcondition] !== WaveFD.ImagingConditionWaveFieldSeparationMIX()
+                # Standard born accumulation for FWI, RTM, or standard imaging condition 
+                cumtime_im += @elapsed WaveFD.adjointBornAccumulation!(p, kwargs[:modeltype], kwargs[:imgcondition], δm_ginsu, wavefields)
+            else
+                # Born accumulation for mix imaging condition 
+                imcon =  WaveFD.ImagingConditionStandard()
+                cumtime_im += @elapsed WaveFD.adjointBornAccumulation!(p, kwargs[:modeltype], imcon, δm_All, wavefields)
+
+                imcon = WaveFD.ImagingConditionWaveFieldSeparationRTM()
+                cumtime_im += @elapsed WaveFD.adjointBornAccumulation!(p, kwargs[:modeltype], imcon, δm_RTM, wavefields)
+
+                weightAll = 1.0f0 - kwargs[:RTM_weight]
+                weightShort = 2.0f0*kwargs[:RTM_weight] - 1.0f0
+
+                for prop in keys(kwargs[:active_modelset])
+                    δm_ginsu[prop] = (δm_All[prop]  .* weightAll) .+ (δm_RTM[prop]  .* weightShort)
+                end
+            end
         end
     end
     for prop in keys(kwargs[:active_modelset])
diff --git a/src/jop_prop2DAcoTTIDenQ_DEO2_FDTD.jl b/src/jop_prop2DAcoTTIDenQ_DEO2_FDTD.jl
index d14d95f..afd0162 100644
--- a/src/jop_prop2DAcoTTIDenQ_DEO2_FDTD.jl
+++ b/src/jop_prop2DAcoTTIDenQ_DEO2_FDTD.jl
@@ -35,6 +35,7 @@ function JetProp2DAcoTTIDenQ_DEO2_FDTD(;
         wavelet = WaveletCausalRicker(f=5.0),
         freesurface = false,
         imgcondition = "standard",
+        RTM_weight = 0.5, 
         nthreads = Sys.CPU_THREADS,
         reportinterval = 500)
     # active and passive earth model properties.  The active set is in the model-space
@@ -116,10 +117,11 @@ function JetProp2DAcoTTIDenQ_DEO2_FDTD(;
     icdict = Dict(
         lowercase("standard") => WaveFD.ImagingConditionStandard(),
         lowercase("FWI") => WaveFD.ImagingConditionWaveFieldSeparationFWI(),
-        lowercase("RTM") => WaveFD.ImagingConditionWaveFieldSeparationRTM())
+        lowercase("RTM") => WaveFD.ImagingConditionWaveFieldSeparationRTM(),
+        lowercase("MIX") => WaveFD.ImagingConditionWaveFieldSeparationMIX())
 
     if lowercase(imgcondition) ∉ keys(icdict)
-        error("Supplied imaging condition 'imgcondition' is not in [standard, FWI, RTM]")
+        error("Supplied imaging condition 'imgcondition' is not in [standard, FWI, RTM, MIX]")
     end
         
     # construct:
@@ -163,6 +165,7 @@ function JetProp2DAcoTTIDenQ_DEO2_FDTD(;
             wavelet = wavelet,
             freesurface = freesurface,
             imgcondition = get(icdict, lowercase(imgcondition), WaveFD.ImagingConditionStandard()),
+            RTM_weight = RTM_weight,
             nthreads = nthreads,
             reportinterval = reportinterval,
             stats = Dict{String,Float64}("MCells/s"=>0.0, "%io"=>0.0, "%inject/extract"=>0.0, "%imaging"=>0.0)))
@@ -401,7 +404,11 @@ Defaults for arguments are shown inside square brackets.
 * `imgcondition` ["standard"] Selects the type of imaging condition used. Choose from "standard", "FWI", 
     and "RTM". "FWI" and "RTM" will perform Kz wavenumber filtering prior to the imaging condition
     in order to promote long wavelengths (for FWI), or remove long wavelength backscattered energy (for 
-    RTM). Note the true adjoint only exists for "standard" imaging condition currently.
+    RTM). "MIX" mixes the FWI imaging condition using the parameter RTM_weight. Note the true adjoint 
+    only exists for "standard" imaging condition currently.
+* `RTM_weight` determines the balance of short wavelengths and long wavelengths in the imaging condition. 
+    A value of 0.0 is equivalent to the FWI imaging condition, a value of 0.5 is equivalent to the standard
+    imaging condtion, and a value of 1.0 is equivalent to the RTM imaging condition. 
 * `nthreads [Sys.CPU_THREADS]` The number of threads to use for OpenMP parallelization of the modeling.
 * `reportinterval [500]` The interval at which information about the propagtion is logged.
 
@@ -763,8 +770,18 @@ function JopProp2DAcoTTIDenQ_DEO2_FDTD_df′!(δm::AbstractArray, δd::AbstractA
     earth_ginsu["f"] .= kwargs[:f]
 
     δm_ginsu = Dict{String,Array{Float32,2}}()
+     # Temporary model updates if mixed imaging condition used. 
+    if kwargs[:imgcondition] ===  WaveFD.ImagingConditionWaveFieldSeparationMIX()
+        δm_RTM = Dict{String,Array{Float32,2}}()
+        δm_All = Dict{String,Array{Float32,2}}()
+    end
+
     for prop in keys(kwargs[:active_modelset])
         δm_ginsu[prop] = zeros(Float32, nz_ginsu, nx_ginsu)
+        if kwargs[:imgcondition] ===  WaveFD.ImagingConditionWaveFieldSeparationMIX()
+            δm_RTM[prop] = zeros(Float32, nz_ginsu, nx_ginsu)
+            δm_All[prop] = zeros(Float32, nz_ginsu, nx_ginsu)
+        end
     end
 
     ginsu_interior_range = interior(kwargs[:ginsu])
@@ -835,7 +852,24 @@ function JopProp2DAcoTTIDenQ_DEO2_FDTD_df′!(δm::AbstractArray, δd::AbstractA
             end
 
             # born accumulation
-            cumtime_im += @elapsed WaveFD.adjointBornAccumulation!(p, kwargs[:modeltype], kwargs[:imgcondition], δm_ginsu, wavefields)
+            if kwargs[:imgcondition] !== WaveFD.ImagingConditionWaveFieldSeparationMIX()
+                # Standard born accumulation for FWI, RTM, or standard imaging condition 
+                cumtime_im += @elapsed WaveFD.adjointBornAccumulation!(p, kwargs[:modeltype], kwargs[:imgcondition], δm_ginsu, wavefields)
+            else
+                # Born accumulation for mix imaging condition 
+                imcon =  WaveFD.ImagingConditionStandard()
+                cumtime_im += @elapsed WaveFD.adjointBornAccumulation!(p, kwargs[:modeltype], imcon, δm_All, wavefields)
+
+                imcon = WaveFD.ImagingConditionWaveFieldSeparationRTM()
+                cumtime_im += @elapsed WaveFD.adjointBornAccumulation!(p, kwargs[:modeltype], imcon, δm_RTM, wavefields)
+
+                weightAll = 1.0f0 - kwargs[:RTM_weight]
+                weightShort = 2.0f0*kwargs[:RTM_weight] - 1.0f0
+
+                for prop in keys(kwargs[:active_modelset])
+                    δm_ginsu[prop] = (δm_All[prop]  .* weightAll) .+ (δm_RTM[prop]  .* weightShort)
+                end
+            end
         end
     end
     for prop in keys(kwargs[:active_modelset])
diff --git a/src/jop_prop2DAcoVTIDenQ_DEO2_FDTD.jl b/src/jop_prop2DAcoVTIDenQ_DEO2_FDTD.jl
index 84e9cd3..28246f7 100644
--- a/src/jop_prop2DAcoVTIDenQ_DEO2_FDTD.jl
+++ b/src/jop_prop2DAcoVTIDenQ_DEO2_FDTD.jl
@@ -34,6 +34,7 @@ function JetProp2DAcoVTIDenQ_DEO2_FDTD(;
         wavelet = WaveletCausalRicker(f=5.0),
         freesurface = false,
         imgcondition = "standard",
+        RTM_weight = 0.5,
         nthreads = Sys.CPU_THREADS,
         reportinterval = 500)
     # active and passive earth model properties.  The active set is in the model-space
@@ -113,10 +114,11 @@ function JetProp2DAcoVTIDenQ_DEO2_FDTD(;
     icdict = Dict(
         lowercase("standard") => WaveFD.ImagingConditionStandard(),
         lowercase("FWI") => WaveFD.ImagingConditionWaveFieldSeparationFWI(),
-        lowercase("RTM") => WaveFD.ImagingConditionWaveFieldSeparationRTM())
+        lowercase("RTM") => WaveFD.ImagingConditionWaveFieldSeparationRTM(),
+        lowercase("MIX") => WaveFD.ImagingConditionWaveFieldSeparationMIX())
 
     if lowercase(imgcondition) ∉ keys(icdict)
-        error("Supplied imaging condition 'imgcondition' is not in [standard, FWI, RTM]")
+        error("Supplied imaging condition 'imgcondition' is not in [standard, FWI, RTM, MIX]")
     end
         
     Jet(
@@ -159,6 +161,7 @@ function JetProp2DAcoVTIDenQ_DEO2_FDTD(;
             wavelet = wavelet,
             freesurface = freesurface,
             imgcondition = get(icdict, lowercase(imgcondition), WaveFD.ImagingConditionStandard()),
+            RTM_weight = RTM_weight, 
             nthreads = nthreads,
             reportinterval = reportinterval,
             stats = Dict{String,Float64}("MCells/s"=>0.0, "%io"=>0.0, "%inject/extract"=>0.0, "%imaging"=>0.0)))
@@ -393,7 +396,11 @@ Defaults for arguments are shown inside square brackets.
 * `imgcondition` ["standard"] Selects the type of imaging condition used. Choose from "standard", "FWI", 
     and "RTM". "FWI" and "RTM" will perform Kz wavenumber filtering prior to the imaging condition
     in order to promote long wavelengths (for FWI), or remove long wavelength backscattered energy (for 
-    RTM). Note the true adjoint only exists for "standard" imaging condition currently.
+    RTM). "MIX" mixes the FWI imaging condition using the parameter RTM_weight. Note the true adjoint 
+    only exists for "standard" imaging condition currently.
+* `RTM_weight` determines the balance of short wavelengths and long wavelengths in the imaging condition. 
+    A value of 0.0 is equivalent to the FWI imaging condition, a value of 0.5 is equivalent to the standard
+    imaging condtion, and a value of 1.0 is equivalent to the RTM imaging condition. 
 * `nthreads [Sys.CPU_THREADS]` The number of threads to use for OpenMP parallelization of the modeling.
 * `reportinterval [500]` The interval at which information about the propagtion is logged.
 
@@ -762,8 +769,18 @@ function JopProp2DAcoVTIDenQ_DEO2_FDTD_df′!(δm::AbstractArray, δd::AbstractA
     earth_ginsu["f"] .= kwargs[:f]
 
     δm_ginsu = Dict{String,Array{Float32,2}}()
+     # Temporary model updates if mixed imaging condition used. 
+    if kwargs[:imgcondition] ===  WaveFD.ImagingConditionWaveFieldSeparationMIX()
+        δm_RTM = Dict{String,Array{Float32,2}}()
+        δm_All = Dict{String,Array{Float32,2}}()
+    end
+
     for prop in keys(kwargs[:active_modelset])
         δm_ginsu[prop] = zeros(Float32, nz_ginsu, nx_ginsu)
+        if kwargs[:imgcondition] ===  WaveFD.ImagingConditionWaveFieldSeparationMIX()
+            δm_RTM[prop] = zeros(Float32, nz_ginsu, nx_ginsu)
+            δm_All[prop] = zeros(Float32, nz_ginsu, nx_ginsu)
+        end
     end
 
     ginsu_interior_range = interior(kwargs[:ginsu])
@@ -836,7 +853,24 @@ function JopProp2DAcoVTIDenQ_DEO2_FDTD_df′!(δm::AbstractArray, δd::AbstractA
             end
 
             # born accumulation
-            cumtime_im += @elapsed WaveFD.adjointBornAccumulation!(p, kwargs[:modeltype], kwargs[:imgcondition], δm_ginsu, wavefields)
+            if kwargs[:imgcondition] !== WaveFD.ImagingConditionWaveFieldSeparationMIX()
+                # Standard born accumulation for FWI, RTM, or standard imaging condition 
+                cumtime_im += @elapsed WaveFD.adjointBornAccumulation!(p, kwargs[:modeltype], kwargs[:imgcondition], δm_ginsu, wavefields)
+            else
+                # Born accumulation for mix imaging condition 
+                imcon =  WaveFD.ImagingConditionStandard()
+                cumtime_im += @elapsed WaveFD.adjointBornAccumulation!(p, kwargs[:modeltype], imcon, δm_All, wavefields)
+
+                imcon = WaveFD.ImagingConditionWaveFieldSeparationRTM()
+                cumtime_im += @elapsed WaveFD.adjointBornAccumulation!(p, kwargs[:modeltype], imcon, δm_RTM, wavefields)
+
+                weightAll = 1.0f0 - kwargs[:RTM_weight]
+                weightShort = 2.0f0*kwargs[:RTM_weight] - 1.0f0
+
+                for prop in keys(kwargs[:active_modelset])
+                    δm_ginsu[prop] = (δm_All[prop]  .* weightAll) .+ (δm_RTM[prop]  .* weightShort)
+                end
+            end
         end
     end
     for prop in keys(kwargs[:active_modelset])
diff --git a/src/jop_prop3DAcoIsoDenQ_DEO2_FDTD.jl b/src/jop_prop3DAcoIsoDenQ_DEO2_FDTD.jl
index 929f6bd..309619a 100644
--- a/src/jop_prop3DAcoIsoDenQ_DEO2_FDTD.jl
+++ b/src/jop_prop3DAcoIsoDenQ_DEO2_FDTD.jl
@@ -41,6 +41,7 @@ function JetProp3DAcoIsoDenQ_DEO2_FDTD(;
         wavelet = WaveletCausalRicker(f=5.0),
         freesurface = false,
         imgcondition = "standard",
+        RTM_weight = 0.5, 
         nthreads = Sys.CPU_THREADS,
         reportinterval = 500)
 
@@ -137,10 +138,11 @@ function JetProp3DAcoIsoDenQ_DEO2_FDTD(;
     icdict = Dict(
         lowercase("standard") => WaveFD.ImagingConditionStandard(),
         lowercase("FWI") => WaveFD.ImagingConditionWaveFieldSeparationFWI(),
-        lowercase("RTM") => WaveFD.ImagingConditionWaveFieldSeparationRTM())
+        lowercase("RTM") => WaveFD.ImagingConditionWaveFieldSeparationRTM(),
+        lowercase("MIX") => WaveFD.ImagingConditionWaveFieldSeparationMIX())
 
     if lowercase(imgcondition) ∉ keys(icdict)
-        error("Supplied imaging condition 'imgcondition' is not in [standard, FWI, RTM]")
+        error("Supplied imaging condition 'imgcondition' is not in [standard, FWI, RTM, MIX]")
     end
         
     Jet(
@@ -187,6 +189,7 @@ function JetProp3DAcoIsoDenQ_DEO2_FDTD(;
             wavelet = wavelet,
             freesurface = freesurface,
             imgcondition = get(icdict, lowercase(imgcondition), WaveFD.ImagingConditionStandard()),
+            RTM_weight = RTM_weight,
             nthreads = nthreads,
             reportinterval = reportinterval,
             stats = Dict{String,Float64}("MCells/s"=>0.0, "%io"=>0.0, "%inject/extract"=>0.0, "%imaging"=>0.0)))
@@ -378,7 +381,11 @@ Defaults for arguments are shown inside square brackets.
 * `imgcondition` ["standard"] Selects the type of imaging condition used. Choose from "standard", "FWI", 
     and "RTM". "FWI" and "RTM" will perform Kz wavenumber filtering prior to the imaging condition
     in order to promote long wavelengths (for FWI), or remove long wavelength backscattered energy (for 
-    RTM). Note the true adjoint only exists for "standard" imaging condition currently.
+    RTM). "MIX" mixes the FWI imaging condition using the parameter RTM_weight. Note the true adjoint 
+    only exists for "standard" imaging condition currently.
+* `RTM_weight` determines the balance of short wavelengths and long wavelengths in the imaging condition. 
+    A value of 0.0 is equivalent to the FWI imaging condition, a value of 0.5 is equivalent to the standard
+    imaging condtion, and a value of 1.0 is equivalent to the RTM imaging condition.
 * `nthreads [Sys.CPU_THREADS]` The number of threads to use for OpenMP parallelization of the modeling.
 * `reportinterval [500]` The interval at which information about the propagtion is logged.
 
@@ -752,8 +759,18 @@ function JopProp3DAcoIsoDenQ_DEO2_FDTD_df′!(δm::AbstractArray, δd::AbstractA
     end
 
     δm_ginsu = Dict{String,Array{Float32,3}}()
+    # Temporary model updates if mixed imaging condition used. 
+    if kwargs[:imgcondition] ===  WaveFD.ImagingConditionWaveFieldSeparationMIX()
+        δm_RTM = Dict{String,Array{Float32,3}}()
+        δm_All = Dict{String,Array{Float32,3}}()
+    end
+
     for prop in keys(kwargs[:active_modelset])
         δm_ginsu[prop] = zeros(Float32, nz_ginsu, ny_ginsu, nx_ginsu)
+        if kwargs[:imgcondition] ===  WaveFD.ImagingConditionWaveFieldSeparationMIX()
+            δm_RTM[prop] = zeros(Float32, nz_ginsu, ny_ginsu, nx_ginsu)
+            δm_All[prop] = zeros(Float32, nz_ginsu, ny_ginsu, nx_ginsu)
+        end
     end
 
     ginsu_interior_range = interior(kwargs[:ginsu])
@@ -829,7 +846,24 @@ function JopProp3DAcoIsoDenQ_DEO2_FDTD_df′!(δm::AbstractArray, δd::AbstractA
             end
 
             # born accumulation
-            cumtime_im += @elapsed WaveFD.adjointBornAccumulation!(p, kwargs[:modeltype], kwargs[:imgcondition], δm_ginsu, wavefields)
+            if kwargs[:imgcondition] !== WaveFD.ImagingConditionWaveFieldSeparationMIX()
+                # Standard born accumulation for FWI, RTM, or standard imaging condition 
+                cumtime_im += @elapsed WaveFD.adjointBornAccumulation!(p, kwargs[:modeltype], kwargs[:imgcondition], δm_ginsu, wavefields)
+            else
+                # Born accumulation for mix imaging condition 
+                imcon =  WaveFD.ImagingConditionStandard()
+                cumtime_im += @elapsed WaveFD.adjointBornAccumulation!(p, kwargs[:modeltype], imcon, δm_All, wavefields)
+
+                imcon = WaveFD.ImagingConditionWaveFieldSeparationRTM()
+                cumtime_im += @elapsed WaveFD.adjointBornAccumulation!(p, kwargs[:modeltype], imcon, δm_RTM, wavefields)
+
+                weightAll = 1.0f0 - kwargs[:RTM_weight]
+                weightShort = 2.0f0*kwargs[:RTM_weight] - 1.0f0
+
+                for prop in keys(kwargs[:active_modelset])
+                    δm_ginsu[prop] = (δm_All[prop]  .* weightAll) .+ (δm_RTM[prop]  .* weightShort)
+                end
+            end
         end
     end
     for prop in keys(kwargs[:active_modelset])
diff --git a/src/jop_prop3DAcoTTIDenQ_DEO2_FDTD.jl b/src/jop_prop3DAcoTTIDenQ_DEO2_FDTD.jl
index 2a5900e..0d009fa 100644
--- a/src/jop_prop3DAcoTTIDenQ_DEO2_FDTD.jl
+++ b/src/jop_prop3DAcoTTIDenQ_DEO2_FDTD.jl
@@ -43,6 +43,7 @@ function JetProp3DAcoTTIDenQ_DEO2_FDTD(;
         wavelet = WaveletCausalRicker(f=5.0),
         freesurface = false,
         imgcondition = "standard",
+        RTM_weight = 0.5, 
         nthreads = Sys.CPU_THREADS,
         reportinterval = 500)
 
@@ -128,10 +129,11 @@ function JetProp3DAcoTTIDenQ_DEO2_FDTD(;
     icdict = Dict(
         lowercase("standard") => WaveFD.ImagingConditionStandard(),
         lowercase("FWI") => WaveFD.ImagingConditionWaveFieldSeparationFWI(),
-        lowercase("RTM") => WaveFD.ImagingConditionWaveFieldSeparationRTM())
+        lowercase("RTM") => WaveFD.ImagingConditionWaveFieldSeparationRTM(),
+        lowercase("MIX") => WaveFD.ImagingConditionWaveFieldSeparationMIX())
 
     if lowercase(imgcondition) ∉ keys(icdict)
-        error("Supplied imaging condition 'imgcondition' is not in [standard, FWI, RTM]")
+        error("Supplied imaging condition 'imgcondition' is not in [standard, FWI, RTM, MIX]")
     end
         
     # construct:
@@ -180,6 +182,7 @@ function JetProp3DAcoTTIDenQ_DEO2_FDTD(;
             wavelet = wavelet,
             freesurface = freesurface,
             imgcondition = get(icdict, lowercase(imgcondition), WaveFD.ImagingConditionStandard()),
+            RTM_weight = RTM_weight, 
             nthreads = nthreads,
             reportinterval = reportinterval,
             stats = Dict{String,Float64}("MCells/s"=>0.0, "%io"=>0.0, "%inject/extract"=>0.0, "%imaging"=>0.0)))
@@ -432,7 +435,11 @@ Defaults for arguments are shown inside square brackets.
 * `imgcondition` ["standard"] Selects the type of imaging condition used. Choose from "standard", "FWI", 
     and "RTM". "FWI" and "RTM" will perform Kz wavenumber filtering prior to the imaging condition
     in order to promote long wavelengths (for FWI), or remove long wavelength backscattered energy (for 
-    RTM). Note the true adjoint only exists for "standard" imaging condition currently.
+    RTM). "MIX" mixes the FWI imaging condition using the parameter RTM_weight. Note the true adjoint 
+    only exists for "standard" imaging condition currently.
+* `RTM_weight` determines the balance of short wavelengths and long wavelengths in the imaging condition. 
+    A value of 0.0 is equivalent to the FWI imaging condition, a value of 0.5 is equivalent to the standard
+    imaging condtion, and a value of 1.0 is equivalent to the RTM imaging condition. 
 * `nthreads [Sys.CPU_THREADS]` The number of threads to use for OpenMP parallelization of the modeling.
 * `reportinterval [500]` The interval at which information about the propagtion is logged.
 
@@ -808,8 +815,18 @@ function JopProp3DAcoTTIDenQ_DEO2_FDTD_df′!(δm::AbstractArray, δd::AbstractA
     earth_ginsu["f"] .= kwargs[:f]
 
     δm_ginsu = Dict{String,Array{Float32,3}}()
+    # Temporary model updates if mixed imaging condition used. 
+    if kwargs[:imgcondition] ===  WaveFD.ImagingConditionWaveFieldSeparationMIX()
+        δm_RTM = Dict{String,Array{Float32,3}}()
+        δm_All = Dict{String,Array{Float32,3}}()
+    end
+
     for prop in keys(kwargs[:active_modelset])
         δm_ginsu[prop] = zeros(Float32, nz_ginsu, ny_ginsu, nx_ginsu)
+        if kwargs[:imgcondition] ===  WaveFD.ImagingConditionWaveFieldSeparationMIX()
+            δm_RTM[prop] = zeros(Float32, nz_ginsu, ny_ginsu, nx_ginsu)
+            δm_All[prop] = zeros(Float32, nz_ginsu, ny_ginsu, nx_ginsu)
+        end
     end
 
     ginsu_interior_range = interior(kwargs[:ginsu])
@@ -884,7 +901,24 @@ function JopProp3DAcoTTIDenQ_DEO2_FDTD_df′!(δm::AbstractArray, δd::AbstractA
             end
 
             # born accumulation
-            cumtime_im += @elapsed WaveFD.adjointBornAccumulation!(p, kwargs[:modeltype], kwargs[:imgcondition], δm_ginsu, wavefields)
+            if kwargs[:imgcondition] !== WaveFD.ImagingConditionWaveFieldSeparationMIX()
+                # Standard born accumulation for FWI, RTM, or standard imaging condition 
+                cumtime_im += @elapsed WaveFD.adjointBornAccumulation!(p, kwargs[:modeltype], kwargs[:imgcondition], δm_ginsu, wavefields)
+            else
+                # Born accumulation for mix imaging condition 
+                imcon =  WaveFD.ImagingConditionStandard()
+                cumtime_im += @elapsed WaveFD.adjointBornAccumulation!(p, kwargs[:modeltype], imcon, δm_All, wavefields)
+
+                imcon = WaveFD.ImagingConditionWaveFieldSeparationRTM()
+                cumtime_im += @elapsed WaveFD.adjointBornAccumulation!(p, kwargs[:modeltype], imcon, δm_RTM, wavefields)
+
+                weightAll = 1.0f0 - kwargs[:RTM_weight]
+                weightShort = 2.0f0*kwargs[:RTM_weight] - 1.0f0
+
+                for prop in keys(kwargs[:active_modelset])
+                    δm_ginsu[prop] = (δm_All[prop]  .* weightAll) .+ (δm_RTM[prop]  .* weightShort)
+                end
+            end
         end
     end
     for prop in keys(kwargs[:active_modelset])
diff --git a/src/jop_prop3DAcoVTIDenQ_DEO2_FDTD.jl b/src/jop_prop3DAcoVTIDenQ_DEO2_FDTD.jl
index 223bb5d..76e7a38 100644
--- a/src/jop_prop3DAcoVTIDenQ_DEO2_FDTD.jl
+++ b/src/jop_prop3DAcoVTIDenQ_DEO2_FDTD.jl
@@ -41,6 +41,7 @@ function JetProp3DAcoVTIDenQ_DEO2_FDTD(;
         wavelet = WaveletCausalRicker(f=5.0),
         freesurface = false,
         imgcondition = "standard",
+        RTM_weight = 0.5,
         nthreads = Sys.CPU_THREADS,
         reportinterval = 500)
 
@@ -122,10 +123,11 @@ function JetProp3DAcoVTIDenQ_DEO2_FDTD(;
     icdict = Dict(
         lowercase("standard") => WaveFD.ImagingConditionStandard(),
         lowercase("FWI") => WaveFD.ImagingConditionWaveFieldSeparationFWI(),
-        lowercase("RTM") => WaveFD.ImagingConditionWaveFieldSeparationRTM())
+        lowercase("RTM") => WaveFD.ImagingConditionWaveFieldSeparationRTM(),
+        lowercase("MIX") => WaveFD.ImagingConditionWaveFieldSeparationMIX())
 
     if lowercase(imgcondition) ∉ keys(icdict)
-        error("Supplied imaging condition 'imgcondition' is not in [standard, FWI, RTM]")
+        error("Supplied imaging condition 'imgcondition' is not in [standard, FWI, RTM, MIX]")
     end
 
     # construct:
@@ -174,6 +176,7 @@ function JetProp3DAcoVTIDenQ_DEO2_FDTD(;
             wavelet = wavelet,
             freesurface = freesurface,
             imgcondition = get(icdict, lowercase(imgcondition), WaveFD.ImagingConditionStandard()),
+            RTM_weight = RTM_weight,
             nthreads = nthreads,
             reportinterval = reportinterval,
             stats = Dict{String,Float64}("MCells/s"=>0.0, "%io"=>0.0, "%inject/extract"=>0.0, "%imaging"=>0.0)))
@@ -418,7 +421,11 @@ Defaults for arguments are shown inside square brackets.
 * `imgcondition` ["standard"] Selects the type of imaging condition used. Choose from "standard", "FWI", 
     and "RTM". "FWI" and "RTM" will perform Kz wavenumber filtering prior to the imaging condition
     in order to promote long wavelengths (for FWI), or remove long wavelength backscattered energy (for 
-    RTM). Note the true adjoint only exists for "standard" imaging condition currently.
+    RTM). "MIX" mixes the FWI imaging condition using the parameter RTM_weight. Note the true adjoint 
+    only exists for "standard" imaging condition currently.
+* `RTM_weight` determines the balance of short wavelengths and long wavelengths in the imaging condition. 
+    A value of 0.0 is equivalent to the FWI imaging condition, a value of 0.5 is equivalent to the standard
+    imaging condtion, and a value of 1.0 is equivalent to the RTM imaging condition. 
 * `nthreads [Sys.CPU_THREADS]` The number of threads to use for OpenMP parallelization of the modeling.
 * `reportinterval [500]` The interval at which information about the propagtion is logged.
 
@@ -791,8 +798,18 @@ function JopProp3DAcoVTIDenQ_DEO2_FDTD_df′!(δm::AbstractArray, δd::AbstractA
     earth_ginsu["f"] .= kwargs[:f]
 
     δm_ginsu = Dict{String,Array{Float32,3}}()
+    # Temporary model updates if mixed imaging condition used. 
+    if kwargs[:imgcondition] ===  WaveFD.ImagingConditionWaveFieldSeparationMIX()
+        δm_RTM = Dict{String,Array{Float32,3}}()
+        δm_All = Dict{String,Array{Float32,3}}()
+    end
+
     for prop in keys(kwargs[:active_modelset])
         δm_ginsu[prop] = zeros(Float32, nz_ginsu, ny_ginsu, nx_ginsu)
+        if kwargs[:imgcondition] ===  WaveFD.ImagingConditionWaveFieldSeparationMIX()
+            δm_RTM[prop] = zeros(Float32, nz_ginsu, ny_ginsu, nx_ginsu)
+            δm_All[prop] = zeros(Float32, nz_ginsu, ny_ginsu, nx_ginsu)
+        end
     end
 
     ginsu_interior_range = interior(kwargs[:ginsu])
@@ -867,7 +884,25 @@ function JopProp3DAcoVTIDenQ_DEO2_FDTD_df′!(δm::AbstractArray, δd::AbstractA
             end
 
             # born accumulation
-            cumtime_im += @elapsed WaveFD.adjointBornAccumulation!(p, kwargs[:modeltype], kwargs[:imgcondition], δm_ginsu, wavefields)
+            if kwargs[:imgcondition] !== WaveFD.ImagingConditionWaveFieldSeparationMIX()
+                # Standard born accumulation for FWI, RTM, or standard imaging condition 
+                cumtime_im += @elapsed WaveFD.adjointBornAccumulation!(p, kwargs[:modeltype], kwargs[:imgcondition], δm_ginsu, wavefields)
+            else
+                # Born accumulation for mix imaging condition 
+                imcon =  WaveFD.ImagingConditionStandard()
+                cumtime_im += @elapsed WaveFD.adjointBornAccumulation!(p, kwargs[:modeltype], imcon, δm_All, wavefields)
+
+                imcon = WaveFD.ImagingConditionWaveFieldSeparationRTM()
+                cumtime_im += @elapsed WaveFD.adjointBornAccumulation!(p, kwargs[:modeltype], imcon, δm_RTM, wavefields)
+
+                weightAll = 1.0f0 - kwargs[:RTM_weight]
+                weightShort = 2.0f0*kwargs[:RTM_weight] - 1.0f0
+
+                for prop in keys(kwargs[:active_modelset])
+                    δm_ginsu[prop] = (δm_All[prop]  .* weightAll) .+ (δm_RTM[prop]  .* weightShort)
+                end
+            end
+
         end
     end
     for prop in keys(kwargs[:active_modelset])