Merge pull request #89 from APiccolo89/ap-mckenzie-temp

Slab-internal temperature and linear thermal averaging

Author: boriskaus

docs/src/man/geodynamic_setups.md

@@ -20,4 +20,6 @@ GeophysicalModelGenerator.LithosphericTemp
 GeophysicalModelGenerator.ConstantPhase
 GeophysicalModelGenerator.Compute_Phase
 GeophysicalModelGenerator.LithosphericPhases
+GeophysicalModelGenerator.McKenzie_subducting_slab
+GeophysicalModelGenerator.LinearWeightedTemperature
 ```

src/Setup_geometry.jl

@@ -13,7 +13,8 @@ export  AddBox!, AddSphere!, AddEllipsoid!, AddCylinder!, AddLayer!,
         makeVolcTopo,
         ConstantTemp, LinearTemp, HalfspaceCoolingTemp, SpreadingRateTemp, LithosphericTemp,
         ConstantPhase, LithosphericPhases,
-        Compute_ThermalStructure, Compute_Phase
+        Compute_ThermalStructure, Compute_Phase,
+        McKenzie_subducting_slab, LinearWeightedTemperature
 
 
 """
@@ -101,10 +102,10 @@ function AddBox!(Phase, Temp, Grid::AbstractGeneralGrid;                 # requi
 
 
     # Set phase number & thermal structure in the full domain
-    ztop = zlim[2] - Origin[3]
-    zbot = zlim[1] - Origin[3]
-    ind = findall(  (Xrot .>= (xlim[1] - Origin[1])) .& (Xrot .<= (xlim[2] - Origin[1])) .&
-                    (Yrot .>= (ylim[1] - Origin[2])) .& (Yrot .<= (ylim[2] - Origin[2])) .&
+    ztop = maximum(zlim) - Origin[3]
+    zbot = minimum(zlim) - Origin[3]
+    ind = findall(  (Xrot .>= (minimum(xlim) - Origin[1])) .& (Xrot .<= (maximum(xlim) - Origin[1])) .&
+                    (Yrot .>= (minimum(ylim) - Origin[2])) .& (Yrot .<= (maximum(ylim) - Origin[2])) .&
                     (Zrot .>= zbot) .& (Zrot .<= ztop)  )
 
     # Compute thermal structure accordingly. See routines below for different options
@@ -1032,3 +1033,147 @@ end
 
 # allow AbstractGeneralGrid instead of Z and Ztop
 Compute_Phase(Phase, Temp, Grid::LaMEM_grid, s::LithosphericPhases) = Compute_Phase(Phase, Temp, Grid.X, Grid.Y, Grid.Z, s::LithosphericPhases, Ztop=maximum(Grid.coord_z))
+
+
+
+"""
+    McKenzie_subducting_slab
+
+Thermal structure by McKenzie for a subducted slab that is fully embedded in the mantle.
+
+Parameters
+===
+- Tsurface:     Top T [C]
+- Tmantle:      Bottom T [C]
+- Adiabat:      Adiabatic gradient in K/km
+- v_cm_yr:      Subduction velocity [cm/yr]
+- κ:            Thermal diffusivity [m2/s]
+- it:           Number iterations employed in the harmonic summation
+
+"""
+@with_kw_noshow mutable struct McKenzie_subducting_slab <: AbstractThermalStructure
+    Tsurface::Float64 = 20.0       # top T
+    Tmantle::Float64  = 1350.0     # bottom T
+    Adiabat::Float64  = 0.4        # Adiabatic gradient in K/km
+    v_cm_yr::Float64  = 2.0        # velocity of subduction [cm/yr]
+    κ::Float64        = 1e-6       # Thermal diffusivity [m2/s]
+    it::Int64         = 36         # number of harmonic summation (look Mckenzie formula)
+end
+
+""" 
+    Compute_ThermalStructure(Temp, X, Y, Z, s::McKenzie_subducting_slab)
+
+Compute the temperature field of a `McKenzie_subducting_slab`. Uses the analytical solution
+of McKenzie (1969) ["Speculations on the consequences and causes of plate motions"]. The functions assumes
+that the bottom of the slab is the coordinate Z=0. Internally the function shifts the coordinate. 
+Parameters
+=============================
+Temp Temperature array
+X    X Array 
+Y    Y Array 
+Z    Z Array 
+Phase Phase array 
+s    McKenzie_subducting_slab
+
+"""
+function Compute_ThermalStructure(Temp, X, Y, Z,Phase, s::McKenzie_subducting_slab)
+    @unpack Tsurface, Tmantle, Adiabat, v_cm_yr, κ, it = s
+
+    # Thickness of the layer: 
+    D0          =   (maximum(Z)-minimum(Z));
+    Zshift      =   Z .- Z[end]       # McKenzie model is defined with Z = 0 at the bottom of the slab
+
+    # Convert subduction velocity from cm/yr -> m/s; 
+    convert_velocity = 1/(100.0*365.25*60.0*60.0*24.0);
+    v_s = v_cm_yr*convert_velocity;
+    
+    # calculate the thermal Reynolds number
+    Re = (v_s*D0*1000)/2/κ;     # factor 1000 to transfer D0 from km to m
+    
+    # McKenzie model
+    sc = 1/D0
+    σ  = ones(size(Temp));
+    # Dividi et impera
+    for i=1:it
+        a   = (-1.0).^(i)./(i.*pi)
+        b   = (Re .- (Re.^2 .+ i^2.0 .* pi^2.0).^(0.5)) .*X .*sc;
+        c   = sin.(i .*pi .*(1 .- abs.(Zshift .*sc))) ;
+        e   = exp.(b);
+        σ .+= 2*a.*e.*c 
+    end
+
+    Temp           .= Tsurface .+ (Tmantle-Tsurface).*σ;
+    Temp           .= Temp + (Adiabat*abs.(Z))
+    
+    return Temp
+end
+
+"""
+    LinearWeightedTemperature
+
+Structure that defined a linear average temperature between two temperature fields as a function of distance
+
+Parameters
+===
+- w_min:        Minimum weight
+- w_max:        Maximum weight
+- crit_dist:    Critical distance
+- dir:          Direction of the averaging (`:X`, `:Y` or `:Z`)
+- F1:           First temperature field
+- F2:           Second temperature field
+
+"""
+@with_kw_noshow mutable struct LinearWeightedTemperature <: AbstractThermalStructure 
+    w_min::Float64 = 0.0; 
+    w_max::Float64 = 1.0; 
+    crit_dist::Float64 = 100.0;
+    dir::Symbol =:X; 
+    F1::AbstractThermalStructure = ConstantTemp();
+    F2::AbstractThermalStructure = ConstantTemp();
+end
+
+"""
+    Weight average along distance
+    Do a weight average between two field along a specified direction 
+    Given a distance {could be any array, from X,Y} -> it increase from the origin the weight of 
+    F1, while F2 decreases. 
+    This function has been conceived for averaging the solution of Mckenzie and half space cooling model, but in 
+    can be used to smooth the temperature field from continent ocean: 
+    -> Select the boundary to apply; 
+    -> transform the coordinate such that dist represent the perpendicular direction along which you want to apply
+    this smoothening and in a such way that 0.0 is the point in which the weight of F1 is equal to 0.0; 
+    -> Select the points that belongs to this area -> compute the thermal fields {F1} {F2} -> then modify F. 
+"""
+function Compute_ThermalStructure(Temp, X, Y, Z, Phase, s::LinearWeightedTemperature)
+    @unpack w_min, w_max, crit_dist,dir = s; 
+    @unpack F1, F2 = s; 
+    
+    if dir === :X
+        dist = X; 
+    elseif dir ===:Y 
+        dist = Y; 
+    else
+        dist = Z; 
+    end
+  
+    # compute the 1D thermal structures
+    Temp1 = zeros(size(Temp));
+    Temp2 = zeros(size(Temp));
+    Temp1 = Compute_ThermalStructure(Temp1, X, Y, Z, Phase, F1);
+    Temp2 = Compute_ThermalStructure(Temp2, X, Y, Z, Phase, F2);
+
+    # Compute the weights
+    weight = w_min .+(w_max-w_min) ./(crit_dist) .*(dist)
+
+    ind_1 = findall(weight .>w_max);
+    ind_2 = findall(weight .<w_min);
+
+    # Change the weight 
+    weight[ind_1] .= w_max; 
+    weight[ind_2] .= w_min;
+    
+    # Average temperature
+    Temp .= Temp1 .*(1.0 .- weight) + Temp2 .* weight; 
+
+    return Temp
+end

test/test_setup_geometry.jl

@@ -161,4 +161,53 @@ AddBox!(Phases,Temp,Grid, xlim=(-100,100), zlim=(-100,0),
     phase=LithosphericPhases(Layers=[20 15 65], Phases = [1 2 3], Tlab=nothing), 
     DipAngle=30.0, T=LithosphericTemp(lbound="flux",nz=201))
 
-@test sum(Temp[Int64(nel/2),Int64(nel/2),:]) ≈ 37359.648604722104
+@test sum(Temp[Int64(nel/2),Int64(nel/2),:]) ≈ 37359.648604722104
+
+
+# Test the McKenzie thermal structure
+
+# Create CartGrid struct
+x        = LinRange(0.0,1200.0,64);
+y        = LinRange(0.0,1200.0,64);
+z        = LinRange(-660,50,64);
+Cart     = CartData(XYZGrid(x, y, z));
+
+# initialize phase and temperature matrix
+Phase   = ones(Int32,size(Cart));
+Temp    = ones(Float64,size(Cart))*1350;
+
+# Create thermal structures
+TsHC = HalfspaceCoolingTemp(Tsurface=20.0, Tmantle=1350, Age=120, Adiabat=0.4)
+TsMK = McKenzie_subducting_slab(Tsurface = 20.0, Tmantle = 1350.0, v_cm_yr = 4.0, Adiabat = 0.0)
+
+@test TsMK.v_cm_yr == 4.0
+@test TsMK.it == 36
+
+# Add a box with a McKenzie thermal structure
+
+# horizontal 
+Temp    = ones(Float64,size(Cart))*1350;
+AddBox!(Phase, Temp, Cart; xlim=(0.0,600.0),ylim=(0.0,600.0), zlim=(-80.0, 0.0), phase = ConstantPhase(5), T=TsMK);
+@test sum(Temp)  ≈ 3.518172093383281e8
+
+# inclined slab
+Temp    = ones(Float64,size(Cart))*1350;
+AddBox!(Phase, Temp, Cart; xlim=(0.0,600.0),ylim=(0.0,600.0), zlim=(-80.0,0),StrikeAngle=0, DipAngle=45, phase = ConstantPhase(5), T=TsMK);
+@test sum(Temp)  ≈ 3.5125017626287365e8
+
+
+
+# horizontal slab, constant T
+T_slab  = LinearWeightedTemperature(0,1,600.0,:X,ConstantTemp(1000), ConstantTemp(2000));
+Temp    = ones(Float64,size(Cart))*1350;
+AddBox!(Phase, Temp, Cart; xlim=(0.0,600.0),ylim=(0.0,600.0), zlim=(-80.0, 0.0), phase = ConstantPhase(5), T=T_slab);
+@test sum(Temp)  ≈ 3.549127111111111e8
+
+# horizontal slab, halfspace and McKenzie
+T_slab = LinearWeightedTemperature(crit_dist=600, F1=TsHC, F2=TsMK);
+Temp    = ones(Float64,size(Cart))*1350;
+AddBox!(Phase, Temp, Cart; xlim=(0.0,600.0),ylim=(0.0,600.0), zlim=(-80.0, 0.0), phase = ConstantPhase(5), T=T_slab);
+@test sum(Temp)  ≈ 3.499457641038468e8
+
+
+Data_Final =   AddField(Cart,"Temp",Temp)