Uw2 - (phase) functions

  • Explain how I handle physical parameters in my models 
  • Describe the parametrisation of thermodynamic phase transitions (Christenson and Yuen )*
  • Overview of the implementation of these in Underworld2

Christensen, U. R., Yuen, D. A., 1984. The Interaction of a Subducting Litho-
spheric Slab with a Chemical or Phase Boundary. Journal of Geophysical
Research

Uw2 - (phase) functions

  • how I handle physical parameters in my models
dp = edict({#Main physical paramters
           'depth':2*660e3, #Depth
           'rho':3300.,  #reference density
           'g':9.8, #surface gravity
           'eta0':1e20, #Dislocation creep at 250 km, 1573 K, 1e-15 s-1 
           'k':1e-6, #thermal diffusivity
           'a':3e-5, #surface thermal expansivity

           ...

           'low_mantle_visc_fac':1.
             })

Benefits of storing your parameters in dictionaries.

  • Easy to save entire parameter suite.
  • harder to overwrite variables (although at some cost to brevity)

         rho vs dp.rho

Uw2 - (phase) functions

  • how I handle physical parameters in my models
#Scaling factors, used to map the dimensional parameters to dimensionless
    
sf = edict({'stress':dp.LS**2/(dp.k*dp.eta0),
            'lith_grad':dp.rho*dp.g*(dp.LS)**3/(dp.eta0*dp.k) ,
            'vel':dp.LS/dp.k,
            'SR':dp.LS**2/dp.k,
            'W':(dp.rho*dp.g*dp.LS)/(dp.R*dp.deltaTa), #Including adiabatic compression, and deltaTa
            'E': 1./(dp.R*dp.deltaTa), #using deltaTa, the guesstimated adiabatic temp difference
            'Ads':1./((dp.eta0**(-1.*dp.n))*(dp.k**(1. - dp.n))*(dp.LS**(-2.+ (2.*dp.n)))),
            'Adf':dp.eta0,
...


#dimensionless parameters
ndp = edict({'RA':(dp.g*dp.rho*dp.a*dp.deltaTa*(dp.LS)**3)/(dp.k*dp.eta0),
             'Di': dp.a*dp.g*dp.LS/dp.Cp, #Dissipation number
             'H':0.,
             #Temperatures and reference depth

...

Uw2 - (phase) functions

  • how I handle physical parameters in my models

Benefits of storing your parameters in dictionaries.

  • scaling is transparent
  • easy to incorporate command line args
docker run -v $PWD:/workspace  -i -t --rm dansand/underworld2-dev \

mpirun -np 8 python kaplan_dev.py md.RES=48 md.ppc=50 dp.Edf*=0.7

Uw2 - (phase) functions

  • Describe the parametrisation of thermodynamic phase transitions (Yuen, etc. )
\nabla \sigma + \left( Ra T - \sum_k Rb_k \Gamma_k - \sum_n Rc_n C_n \right) \delta_{ir}= 0
∇σ+(RaT−∑kRbkΓk−∑nRcnCn)δir=0\nabla \sigma + \left( Ra T - \sum_k Rb_k \Gamma_k - \sum_n Rc_n C_n \right) \delta_{ir}= 0
Ra = \frac{\alpha \rho_0 g \Delta T L^3}{\eta_0 \kappa_0}
Ra=αρ0gΔTL3η0κ0 Ra = \frac{\alpha \rho_0 g \Delta T L^3}{\eta_0 \kappa_0}
Rb_k = \frac{X_k \rho_k g \ L^3}{\eta_0 \kappa_0} = Ra \frac{X_k \rho_k }{\rho_0 \alpha_0 \Delta T}
Rbk=Xkρkg L3η0κ0=RaXkρkρ0α0ΔTRb_k = \frac{X_k \rho_k g \ L^3}{\eta_0 \kappa_0} = Ra \frac{X_k \rho_k }{\rho_0 \alpha_0 \Delta T}
Rc_n = \frac{\rho_n g \ L^3}{\eta_0 \kappa_0} = Ra \frac{ \rho_n }{\rho_0 \alpha_0 \Delta T}
Rcn=ρng L3η0κ0=Raρnρ0α0ΔTRc_n = \frac{\rho_n g \ L^3}{\eta_0 \kappa_0} = Ra \frac{ \rho_n }{\rho_0 \alpha_0 \Delta T}

Uw2 - (phase) functions

  • Describe the parametrisation of thermodynamic phase transitions (Yuen and Christenson )
\nabla \sigma + \left( Ra T - \sum_k Rb_k \Gamma_k - \sum_n Rc_n C_n \right) \delta_{ir}= 0
∇σ+(RaT−∑kRbkΓk−∑nRcnCn)δir=0\nabla \sigma + \left( Ra T - \sum_k Rb_k \Gamma_k - \sum_n Rc_n C_n \right) \delta_{ir}= 0
\Gamma_k = \frac{1}{2} \left( 1 + \tanh(\frac{\delta y}{ w}) \right)
Γk=12(1+tanh(δyw))\Gamma_k = \frac{1}{2} \left( 1 + \tanh(\frac{\delta y}{ w}) \right)
\delta y = (1 - y -d_{ph}) - \gamma_{ph}(T - T_{ph})
δy=(1−y−dph)−γph(T−Tph)\delta y = (1 - y -d_{ph}) - \gamma_{ph}(T - T_{ph})

Uw2 - (phase) functions

\nabla \sigma + \left( Ra T - \sum_k Rb_k \Gamma_k - \sum_n Rc_n C_n \right) \delta_{ir}= 0
∇σ+(RaT−∑kRbkΓk−∑nRcnCn)δir=0\nabla \sigma + \left( Ra T - \sum_k Rb_k \Gamma_k - \sum_n Rc_n C_n \right) \delta_{ir}= 0
\Gamma_k = \frac{1}{2} \left(1 + \tanh(\frac{\delta y}{ w}) \right)
Γk=12(1+tanh(δyw))\Gamma_k = \frac{1}{2} \left(1 + \tanh(\frac{\delta y}{ w}) \right)
\delta y = (1 - y -d_{ph}) - \gamma_{ph}(T - T_{ph})
δy=(1−y−dph)−γph(T−Tph)\delta y = (1 - y -d_{ph}) - \gamma_{ph}(T - T_{ph})
Rb_k = \frac{X_k \rho_k g \ L^3}{\eta_0 \kappa_0} = Ra \frac{X_k \rho_k }{\rho_0 \alpha_0 \Delta T}
Rbk=Xkρkg L3η0κ0=RaXkρkρ0α0ΔTRb_k = \frac{X_k \rho_k g \ L^3}{\eta_0 \kappa_0} = Ra \frac{X_k \rho_k }{\rho_0 \alpha_0 \Delta T}

Uw2 - (phase) functions

class component_phases():
    """
    Class that allows you to create 'phase functions' for a mineral component
    """
    def __init__(self,name, depths,temps, widths,claps,densities):
        """
        Class initialiser.
        Parameter
        ---------
        name : str
            'Component', e.g olivine, pyroxene-garnet
        depths: list
            list of transition depths in kilometers
        widths: list
            list of transition widths in kilometers
        claps: list
            list of Clapeyron slopes in Pa/K
        densities: list
            list of density changes in kg/m3
        Returns
        -------
        dp : Dictionary storing the phase-transition vales

Uw2 - (phase) functions

def nd_phase(self, reduced_p, widthPh):
    """
    Creates an Underworld function, representing the phase function in the domain
    """
    return 0.5*(1. + fn.math.tanh(reduced_p/(widthPh)))

Uw2 - (phase) functions

 def buoyancy_sum(self, temperatureField, depthFn, gravityscale, lengthscale, diffusivityscale, viscosityscale):
        """
        Creates an Underworld function, representing the Sum of the individual phase functions...
        and the associated density changes:
        pf_sum = Sum_k{ (Ra*delRho_k*pf_k/rho_0*eta_0*delta_t)}
        -----------
        temperatureField : underworld.mesh._meshvariable.MeshVariable
        ...need to put warning in about running build_nd_dict first
        """
        bouyancy_factor = (gravityscale*lengthscale**3)/(viscosityscale*diffusivityscale)

        pf_sum = uw.function.misc.constant(0.)

        for phaseId in range(len(self.dp['depths'])):
            #build reduced pressure
            rp = self.nd_reduced_pressure(depthFn,
                                   temperatureField,
                                   self.ndp['depths'][phaseId ],
                                   self.ndp['claps'][phaseId ],
                                   self.ndp['temps'][phaseId ])
            #build phase function
            pf = self.nd_phase(rp, self.ndp['widths'][phaseId ])
            pf_sum += bouyancy_factor*pf*self.dp['densities'][phaseId ] #we want the dimensional densities here

Uw2 - (phase) functions

buoyancy_factor = (dp.g*dp.LS**3)/(dp.eta0*dp.k)

basalt_comp_buoyancy  = (dp.rho - 2940.)*buoyancy_factor
harz_comp_buoyancy = (dp.rho - 3235.)*buoyancy_factor
pyrolite_comp_buoyancy = (dp.rho - 3300.)*buoyancy_factor




if not md.compBuoyancy:
    pyrolitebuoyancyFn =  (ndp.RA*temperatureField*taFn)
    harzbuoyancyFn =      (ndp.RA*temperatureField*taFn) 
    basaltbuoyancyFn =    (ndp.RA*temperatureField*taFn)

else : 
    pyrolitebuoyancyFn =  (ndp.RA*temperatureField*taFn) -\
                          (0.6*olivine_phase_buoyancy + 0.4*garnet_phase_buoyancy) +\
                           pyrolite_comp_buoyancy
    harzbuoyancyFn =      (ndp.RA*temperatureField*taFn) -\
                          (0.8*olivine_phase_buoyancy + 0.2*garnet_phase_buoyancy) +\
                           harz_comp_buoyancy
    basaltbuoyancyFn =    (ndp.RA*temperatureField*taFn) -\
                          (1.*garnet_phase_buoyancy) +\
                           basalt_comp_buoyancy

Uw2 - (phase) functions

Uw2 - (phase) functions

Uw2 - (phase) functions

Underworld meetup

My thesis research involves three main themes:

  • faults in continuum models, particularly the subduction fault,
  • controls on the subduction dip, asymmetry,  and upscaling of subduction models to conve ction models
  • evolution of 3d collisional and congested subduction scenarios

Underworld meetup

This talk:

Subduction velocity patterns in young subduction zones

Background:

  • 50 % of known subduction zones initiated in the Cenozoic.
  • Much work on velocities for current SZs, what about changes over time?
  • Often initiation is followed by intense back arc-spreading.

 

First: some gratuitous uw2 

Even with comparable approaches (rheology), different researchers seem to  generate quite different behaviour

By the way, what do we like / dislike about these kind of models...

  • viscous rheology is consistent with experiments
  • models can be tuned to near equilibrium 
  • decoupling is still fudged 
  • Not clear whether these can be upscaled
  • Fairly expensive to run (3d?)

By the way, what do we like / dislike about these kind of models...

  • viscous rheology is consistent with experiments
  • models can be tuned to near equilibrium 
  • decoupling is still fudged 
  • Not clear whether these can be upscaled
  • Fairly expensive to run (3d?)

One thing these types of models do agree on is the a distinctive velocity pattern in the free sinking and transition-zone interaction phase

This is a pretty big signal - hard to disguise, if this is the way SZs really evolve

So what might be wrong with this model?

So what might be wrong with this models?

Discounting all effects due to 3d, and non- plate-scale mantle flow...

perhaps we can isolate three factors to investigate...

 

Discounting all effects due to 3d, and non- plate-scale mantle flow...

perhaps we can isolate three factors to investigate...

 

  • upper mantle rheology, particularly non-linearity
  • Subducting plate age-bouyancy structure
  • evolution of subduction interface

 

Questions

  • Do we see this is in nature?
    • Not sure yet, but I here's what I plan to do to find out
  • And if not why not?
    • How much does slab pull contribute to the overall force balance?
    • How to the forces evolve?

 

Questions

  • Do we see this is in nature?
    • Not sure yet, but I here's what I plan to do to find out
  • And if not why not?
    • How much does slab pull contribute to the overall force balance?
    • How to the forces evolve?

Underworld meetup

By Dan Sandiford

Made with Slides.com

Underworld meetup

  • 562

More from Dan Sandiford