Skip to content

attenuation

Provides models that return quality factors and anelastic seismic velocities.

AttenuationModelGoes

Bases: object

This class implements a mantle seismic attenuation model [Goes et al., 2004, Maguire et al., 2016].

Optionally, different \(Q\) models can be used that correspond to different mantle materials. A mixing model function should be passed to the constructor that takes pressure and temperature as inputs and returns the fractions of the different materials.

This class has an anelastic_properties method to calculate anelastic \(V_P\) and \(V_S\), \(Q_S\), and \(Q_K\). The effective \(Q_S\), \(Q_K\) and \(\alpha\) (frequency dependence of the quality factor) are given by the linearly weighted sum of the \(Q_S\), \(Q_K\) and \(\alpha\) calculated for each material.

Goes et al., 2004

Goes S., Cammarano F. and Hansen U., 2004. Synthetic seismic signature of thermal mantle plumes. Earth and Planetary Science Letters. 218, pp.403-419. 10.1016/S0012-821X(03)00680-0
Maguire et al., 2016

Maguire R., Ritsema J., van Keken P.E., Fichtner A. and Goes S., 2016. P- and S-wave delays caused by thermal plumes. Geophysical Journal International. 206, pp.1169-1178. 10.1093/gji/ggw187
Source code in terratools/properties/attenuation.py 
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
class AttenuationModelGoes(object):
    """
    This class implements a mantle seismic attenuation model
    [@Goes2004, @Maguire2016].

    Optionally, different $Q$ models can be used that correspond to
    different mantle materials. A mixing model function should be passed
    to the constructor that takes pressure and temperature as inputs and
    returns the fractions of the different materials.

    This class has an anelastic_properties method to calculate
    anelastic $V_P$ and $V_S$, $Q_S$, and $Q_K$.
    The effective $Q_S$, $Q_K$ and $\\alpha$
    (frequency dependence of the quality factor)
    are given by the linearly weighted sum
    of the $Q_S$, $Q_K$ and $\\alpha$ calculated for each material.
    """

    def __init__(self, T_solidus_function, model_mixing_function, Q_models):
        """
        Constructor for the AttenuationModelGoes class.

        :param function T_solidus_function:
            A function returning the temperature of the solidus
            as a function of pressure.
        :param function model_mixing_function:
            A function returning the amounts of different materials
            as a function of pressure and temperature.
        :param list Q_models:
            Parameter dictionaries for the attenuation models - one for each material.
        """
        self.T_solidus_function = T_solidus_function
        self.model_mixing_function = model_mixing_function
        self.Q_models = Q_models

    def anelastic_properties(
        self,
        elastic_Vp,
        elastic_Vs,
        pressure,
        temperature,
        frequency,
        dT_Q_constant_above_solidus=0,
    ):
        """
        Calculates the anelastic $V_P$ and $V_S$, $Q_S$, and $Q_K$
        according to a published model [@Maguire2016].

        The effects of anelasticity on shear wave velocity are incorporated
        using a model for the S-wave quality factor $Q_S$ that varies with
        pressure $P$ and temperature $T$ as

        $Q_S(\\omega,z,T) = Q_0 \\omega \\alpha \\exp(\\alpha g T_m(z) / T)$

        where $\\omega$ is frequency,
        $\\alpha$ is exponential frequency dependence,
        $g$ is a scaling factor and
        $T_m$ is the dry solidus melting temperature.

        $Q_K$ is chosen to be temperature independent.

        The anelastic seismic velocities are calculated as follows:

        $\\lambda = 4/3 (V_{S,\\text{el}}/V_{P,\\text{el}})^2$

        $1/Q_P = (1 - \\lambda)/Q_K + \\lambda/Q_S

        If $1/Q_P$ is negative, it is set to 0.

        $V_{P,\\text{an}} = V_{P,\\text{el}} (1 - Q_P^{-1}/(2 \\tan ( \\pi \\alpha/2)))$

        $V_{S,\\text{an}} = V_{S,\\text{el}} (1 - Q_S^{-1}/(2 \\tan ( \\pi \\alpha/2)))$

        :param elastic_Vp: The elastic P-wave velocity
        :type elastic_Vp: float or numpy array

        :param elastic_Vs: The elastic S-wave velocity
        :type elastic_Vs: float or numpy array

        :param pressure: The pressure in Pa
        :type pressure: float or numpy array

        :param temperature: The temperature in K
        :type temperature: float or numpy array

        :param frequency: The frequency of the seismic waves in Hz
        :type frequency: float

        :param dT_Q_constant_above_solidus: if the temperature > (solidus temperature + dT),
            the value of QS, QK and a are frozen at the values
            corresponding to (solidus temperature + dT).
        :type dT_Q_constant_above_solidus: float

        :return: An instance of an AnelasticProperties named tuple.
        Has the following attributes:
        V_P, V_S, Q_S, Q_K, Q_P
        """

        fractions = self.model_mixing_function(pressure, temperature)

        try:
            pressure = float(pressure)
            Tm = self.T_solidus_function(pressure)
            # Freezes QS if above a certain temperature
            if dT_Q_constant_above_solidus < temperature - Tm:
                Q_temperature = Tm + dT_Q_constant_above_solidus
            else:
                Q_temperature = deepcopy(temperature)

            QS = 0.0
            QK = 0.0
            alpha = 0.0
            for i, f in enumerate(fractions):
                Q_mod = self.Q_models[i]
                QS += f * (
                    Q_mod["Q0"]
                    * np.power(frequency, Q_mod["a"])
                    * np.exp(Q_mod["a"] * Q_mod["g"] * Tm / Q_temperature)
                )
                QK += f * Q_mod["QK"]
                alpha += f * Q_mod["a"]

        except TypeError:
            Q_temperature = deepcopy(temperature)
            Tm = self.T_solidus_function(pressure)
            idx = np.argwhere(temperature > Tm + dT_Q_constant_above_solidus)
            Q_temperature[idx] = Tm[idx] + dT_Q_constant_above_solidus

            QS = np.zeros_like(temperature)
            QK = np.zeros_like(temperature)
            alpha = np.zeros_like(temperature)
            for i, f in enumerate(fractions.T):
                Q_mod = self.Q_models[i]
                QS += f * (
                    Q_mod["Q0"]
                    * np.power(frequency, Q_mod["a"])
                    * np.exp(Q_mod["a"] * Q_mod["g"] * Tm / Q_temperature)
                )
                QK += f * Q_mod["QK"]
                alpha += f * Q_mod["a"]

        invQS = 1.0 / QS
        invQK = 1.0 / QK

        lmda = 4.0 / 3.0 * np.power(elastic_Vs / elastic_Vp, 2.0)
        invQP = (1.0 - lmda) * invQK + lmda * invQS

        try:
            if invQP < 0.0:
                invQP = 0.0
                QP = np.inf
            else:
                QP = 1.0 / invQP
        except ValueError:
            QP = np.zeros_like(temperature)
            idx = np.argwhere(invQP <= 0.0)
            invQP[idx] = 0.0
            QP[idx] = np.inf
            idx = np.argwhere(invQP > 0.0)
            QP[idx] = 1.0 / invQP[idx]

        anelastic_Vp = elastic_Vp * (1.0 - invQP / (2.0 * np.tan(np.pi * alpha / 2.0)))
        anelastic_Vs = elastic_Vs * (1.0 - invQS / (2.0 * np.tan(np.pi * alpha / 2.0)))

        return AnelasticProperties(
            V_P=anelastic_Vp, V_S=anelastic_Vs, Q_S=QS, Q_K=QK, Q_P=QP, T_solidus=Tm
        )

__init__(T_solidus_function, model_mixing_function, Q_models)

Constructor for the AttenuationModelGoes class.

Parameters:

Name Type Description Default
T_solidus_function function

A function returning the temperature of the solidus as a function of pressure.

 required
model_mixing_function function

A function returning the amounts of different materials as a function of pressure and temperature.

 required
Q_models list

Parameter dictionaries for the attenuation models - one for each material.

 required
Source code in terratools/properties/attenuation.py 
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
def __init__(self, T_solidus_function, model_mixing_function, Q_models):
    """
    Constructor for the AttenuationModelGoes class.

    :param function T_solidus_function:
        A function returning the temperature of the solidus
        as a function of pressure.
    :param function model_mixing_function:
        A function returning the amounts of different materials
        as a function of pressure and temperature.
    :param list Q_models:
        Parameter dictionaries for the attenuation models - one for each material.
    """
    self.T_solidus_function = T_solidus_function
    self.model_mixing_function = model_mixing_function
    self.Q_models = Q_models

anelastic_properties(elastic_Vp, elastic_Vs, pressure, temperature, frequency, dT_Q_constant_above_solidus=0)

Calculates the anelastic \(V_P\) and \(V_S\), \(Q_S\), and \(Q_K\) according to a published model [Maguire et al., 2016].

The effects of anelasticity on shear wave velocity are incorporated using a model for the S-wave quality factor \(Q_S\) that varies with pressure \(P\) and temperature \(T\) as

\(Q_S(\omega,z,T) = Q_0 \omega \alpha \exp(\alpha g T_m(z) / T)\)

where \(\omega\) is frequency, \(\alpha\) is exponential frequency dependence, \(g\) is a scaling factor and \(T_m\) is the dry solidus melting temperature.

\(Q_K\) is chosen to be temperature independent.

The anelastic seismic velocities are calculated as follows:

\(\lambda = 4/3 (V_{S,\text{el}}/V_{P,\text{el}})^2\)

$1/Q_P = (1 - \lambda)/Q_K + \lambda/Q_S

If \(1/Q_P\) is negative, it is set to 0.

\(V_{P,\text{an}} = V_{P,\text{el}} (1 - Q_P^{-1}/(2 \tan ( \pi \alpha/2)))\)

\(V_{S,\text{an}} = V_{S,\text{el}} (1 - Q_S^{-1}/(2 \tan ( \pi \alpha/2)))\)

Maguire et al., 2016

Maguire R., Ritsema J., van Keken P.E., Fichtner A. and Goes S., 2016. P- and S-wave delays caused by thermal plumes. Geophysical Journal International. 206, pp.1169-1178. 10.1093/gji/ggw187

Parameters:

Name Type Description Default
elastic_Vp float | numpy array

The elastic P-wave velocity

 required
elastic_Vs float | numpy array

The elastic S-wave velocity

 required
pressure float | numpy array

The pressure in Pa

 required
temperature float | numpy array

The temperature in K

 required
frequency float

The frequency of the seismic waves in Hz

 required
dT_Q_constant_above_solidus float

if the temperature > (solidus temperature + dT), the value of QS, QK and a are frozen at the values corresponding to (solidus temperature + dT).

 0

Returns:

Type Description

An instance of an AnelasticProperties named tuple. Has the following attributes: V_P, V_S, Q_S, Q_K, Q_P
Source code in terratools/properties/attenuation.py 
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
def anelastic_properties(
    self,
    elastic_Vp,
    elastic_Vs,
    pressure,
    temperature,
    frequency,
    dT_Q_constant_above_solidus=0,
):
    """
    Calculates the anelastic $V_P$ and $V_S$, $Q_S$, and $Q_K$
    according to a published model [@Maguire2016].

    The effects of anelasticity on shear wave velocity are incorporated
    using a model for the S-wave quality factor $Q_S$ that varies with
    pressure $P$ and temperature $T$ as

    $Q_S(\\omega,z,T) = Q_0 \\omega \\alpha \\exp(\\alpha g T_m(z) / T)$

    where $\\omega$ is frequency,
    $\\alpha$ is exponential frequency dependence,
    $g$ is a scaling factor and
    $T_m$ is the dry solidus melting temperature.

    $Q_K$ is chosen to be temperature independent.

    The anelastic seismic velocities are calculated as follows:

    $\\lambda = 4/3 (V_{S,\\text{el}}/V_{P,\\text{el}})^2$

    $1/Q_P = (1 - \\lambda)/Q_K + \\lambda/Q_S

    If $1/Q_P$ is negative, it is set to 0.

    $V_{P,\\text{an}} = V_{P,\\text{el}} (1 - Q_P^{-1}/(2 \\tan ( \\pi \\alpha/2)))$

    $V_{S,\\text{an}} = V_{S,\\text{el}} (1 - Q_S^{-1}/(2 \\tan ( \\pi \\alpha/2)))$

    :param elastic_Vp: The elastic P-wave velocity
    :type elastic_Vp: float or numpy array

    :param elastic_Vs: The elastic S-wave velocity
    :type elastic_Vs: float or numpy array

    :param pressure: The pressure in Pa
    :type pressure: float or numpy array

    :param temperature: The temperature in K
    :type temperature: float or numpy array

    :param frequency: The frequency of the seismic waves in Hz
    :type frequency: float

    :param dT_Q_constant_above_solidus: if the temperature > (solidus temperature + dT),
        the value of QS, QK and a are frozen at the values
        corresponding to (solidus temperature + dT).
    :type dT_Q_constant_above_solidus: float

    :return: An instance of an AnelasticProperties named tuple.
    Has the following attributes:
    V_P, V_S, Q_S, Q_K, Q_P
    """

    fractions = self.model_mixing_function(pressure, temperature)

    try:
        pressure = float(pressure)
        Tm = self.T_solidus_function(pressure)
        # Freezes QS if above a certain temperature
        if dT_Q_constant_above_solidus < temperature - Tm:
            Q_temperature = Tm + dT_Q_constant_above_solidus
        else:
            Q_temperature = deepcopy(temperature)

        QS = 0.0
        QK = 0.0
        alpha = 0.0
        for i, f in enumerate(fractions):
            Q_mod = self.Q_models[i]
            QS += f * (
                Q_mod["Q0"]
                * np.power(frequency, Q_mod["a"])
                * np.exp(Q_mod["a"] * Q_mod["g"] * Tm / Q_temperature)
            )
            QK += f * Q_mod["QK"]
            alpha += f * Q_mod["a"]

    except TypeError:
        Q_temperature = deepcopy(temperature)
        Tm = self.T_solidus_function(pressure)
        idx = np.argwhere(temperature > Tm + dT_Q_constant_above_solidus)
        Q_temperature[idx] = Tm[idx] + dT_Q_constant_above_solidus

        QS = np.zeros_like(temperature)
        QK = np.zeros_like(temperature)
        alpha = np.zeros_like(temperature)
        for i, f in enumerate(fractions.T):
            Q_mod = self.Q_models[i]
            QS += f * (
                Q_mod["Q0"]
                * np.power(frequency, Q_mod["a"])
                * np.exp(Q_mod["a"] * Q_mod["g"] * Tm / Q_temperature)
            )
            QK += f * Q_mod["QK"]
            alpha += f * Q_mod["a"]

    invQS = 1.0 / QS
    invQK = 1.0 / QK

    lmda = 4.0 / 3.0 * np.power(elastic_Vs / elastic_Vp, 2.0)
    invQP = (1.0 - lmda) * invQK + lmda * invQS

    try:
        if invQP < 0.0:
            invQP = 0.0
            QP = np.inf
        else:
            QP = 1.0 / invQP
    except ValueError:
        QP = np.zeros_like(temperature)
        idx = np.argwhere(invQP <= 0.0)
        invQP[idx] = 0.0
        QP[idx] = np.inf
        idx = np.argwhere(invQP > 0.0)
        QP[idx] = 1.0 / invQP[idx]

    anelastic_Vp = elastic_Vp * (1.0 - invQP / (2.0 * np.tan(np.pi * alpha / 2.0)))
    anelastic_Vs = elastic_Vs * (1.0 - invQS / (2.0 * np.tan(np.pi * alpha / 2.0)))

    return AnelasticProperties(
        V_P=anelastic_Vp, V_S=anelastic_Vs, Q_S=QS, Q_K=QK, Q_P=QP, T_solidus=Tm
    )

Q4Goes

Bases: AttenuationModelGoes

Implements the weak T dependence attenuation model after [Goes et al., 2004].

The model uses the peridotite_solidus and mantle_domain_fractions functions to determine the attenuation and proportions of upper mantle, transition zone and lower mantle materials as a function of pressure and temperature.

The parameter values for the Q4 attenuation models are given in the following table:

Layer \(Q_0\) \(g\) \(\alpha\) \(Q_K\)
upper mantle 0.1 38.0 0.15 1000.0
transition zone 3.5 20.0 0.15 1000.0
lower mantle 35.0 10.0 0.15 1000.0
Goes et al., 2004

Goes S., Cammarano F. and Hansen U., 2004. Synthetic seismic signature of thermal mantle plumes. Earth and Planetary Science Letters. 218, pp.403-419. 10.1016/S0012-821X(03)00680-0
Source code in terratools/properties/attenuation.py 
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
class Q4Goes(AttenuationModelGoes):
    """
    Implements the weak T dependence attenuation model
    after [@Goes2004].

    The model uses the
    [peridotite_solidus][terratools.properties.profiles.peridotite_solidus] and
    [mantle_domain_fractions][terratools.properties.attenuation.mantle_domain_fractions]
    functions to determine the attenuation
    and proportions of upper mantle, transition zone and lower mantle
    materials as a function of pressure and temperature.

    The parameter values for the Q4 attenuation models are given in the
    following table:

    | Layer           | $Q_0$ | $g$  | $\\alpha$ | $Q_K$  |
    | --------------- | ----- | ---- | --------- | ------ |
    | upper mantle    | 0.1   | 38.0 | 0.15      | 1000.0 |
    | transition zone | 3.5   | 20.0 | 0.15      | 1000.0 |
    | lower mantle    | 35.0  | 10.0 | 0.15      | 1000.0 |

    """

    def __init__(self):
        super().__init__(
            peridotite_solidus,
            mantle_domain_fractions,
            Q_models=[
                {"Q0": 0.1, "g": 38.0, "a": 0.15, "QK": 1000.0},
                {"Q0": 3.5, "g": 20.0, "a": 0.15, "QK": 1000.0},
                {"Q0": 35.0, "g": 10.0, "a": 0.15, "QK": 1000.0},
            ],
        )

Q6Goes

Bases: AttenuationModelGoes

Implements the strong T dependence attenuation model [Goes et al., 2004].

The model uses the peridotite_solidus and mantle_domain_fractions functions to determine the attenuation and proportions of upper mantle, transition zone and lower mantle materials as a function of pressure and temperature.

The parameter values for the Q6 attenuation models are given in the following table:

Layer \(Q_0\) \(g\) \(\alpha\) \(Q_K\)
upper mantle 0.1 38.0 0.15 1000.0
transition zone 0.5 30.0 0.15 1000.0
lower mantle 3.5 20.0 0.15 1000.0
Goes et al., 2004

Goes S., Cammarano F. and Hansen U., 2004. Synthetic seismic signature of thermal mantle plumes. Earth and Planetary Science Letters. 218, pp.403-419. 10.1016/S0012-821X(03)00680-0
Source code in terratools/properties/attenuation.py 
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
class Q6Goes(AttenuationModelGoes):
    """
    Implements the strong T dependence attenuation model
    [@Goes2004].

    The model uses the
    [peridotite_solidus][terratools.properties.profiles.peridotite_solidus] and
    [mantle_domain_fractions][terratools.properties.attenuation.mantle_domain_fractions]
    functions to determine the attenuation
    and proportions of upper mantle, transition zone and lower mantle
    materials as a function of pressure and temperature.

    The parameter values for the Q6 attenuation models are given in the
    following table:

    | Layer           | $Q_0$ | $g$  | $\\alpha$ | $Q_K$  |
    | --------------- | ----- | ---- | --------- | ------ |
    | upper mantle    | 0.1   | 38.0 | 0.15      | 1000.0 |
    | transition zone | 0.5   | 30.0 | 0.15      | 1000.0 |
    | lower mantle    | 3.5   | 20.0 | 0.15      | 1000.0 |

    """

    def __init__(self):
        super().__init__(
            peridotite_solidus,
            mantle_domain_fractions,
            Q_models=[
                {"Q0": 0.1, "g": 38.0, "a": 0.15, "QK": 1000.0},
                {"Q0": 0.5, "g": 30.0, "a": 0.15, "QK": 1000.0},
                {"Q0": 3.5, "g": 20.0, "a": 0.15, "QK": 1000.0},
            ],
        )

Q7Goes

Bases: AttenuationModelGoes

Implements the intermediate strength T dependence attenuation model [Goes et al., 2004]. This model is most consistent with observational data [Matas and Bukowinski, 2007].

The model uses the peridotite_solidus and mantle_domain_fractions functions to determine the attenuation and proportions of upper mantle, transition zone and lower mantle materials as a function of pressure and temperature.

The parameter values for the Q7 attenuation models are given in the following table:

Layer \(Q_0\) \(g\) \(\alpha\) \(Q_K\)
upper mantle 0.1 38.0 0.15 1000.0
transition zone 0.5 30.0 0.15 1000.0
lower mantle 1.5 26.0 0.15 1000.0
Goes et al., 2004

Goes S., Cammarano F. and Hansen U., 2004. Synthetic seismic signature of thermal mantle plumes. Earth and Planetary Science Letters. 218, pp.403-419. 10.1016/S0012-821X(03)00680-0
Matas and Bukowinski, 2007

Matas J. and Bukowinski M.S., 2007. On the anelastic contribution to the temperature dependence of lower mantle seismic velocities. Earth and Planetary Science Letters. 259, pp.51-65. 10.1016/j.epsl.2007.04.028
Source code in terratools/properties/attenuation.py 
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
class Q7Goes(AttenuationModelGoes):
    """
    Implements the intermediate strength T dependence attenuation model
    [@Goes2004]. This model is most consistent with observational data
    [@Matas2007].

    The model uses the
    [peridotite_solidus][terratools.properties.profiles.peridotite_solidus] and
    [mantle_domain_fractions][terratools.properties.attenuation.mantle_domain_fractions]
    functions to determine the attenuation
    and proportions of upper mantle, transition zone and lower mantle
    materials as a function of pressure and temperature.

    The parameter values for the Q7 attenuation models are given in the
    following table:

    | Layer           | $Q_0$ | $g$  | $\\alpha$ | $Q_K$  |
    | --------------- | ----- | ---- | --------- | ------ |
    | upper mantle    | 0.1   | 38.0 | 0.15      | 1000.0 |
    | transition zone | 0.5   | 30.0 | 0.15      | 1000.0 |
    | lower mantle    | 1.5   | 26.0 | 0.15      | 1000.0 |


    """

    def __init__(self):
        super().__init__(
            peridotite_solidus,
            mantle_domain_fractions,
            Q_models=[
                {"Q0": 0.1, "g": 38.0, "a": 0.15, "QK": 1000.0},
                {"Q0": 0.5, "g": 30.0, "a": 0.15, "QK": 1000.0},
                {"Q0": 1.5, "g": 26.0, "a": 0.15, "QK": 1000.0},
            ],
        )

mantle_domain_fractions(pressure, temperature)

This function defines the proportions of upper mantle, transition zone, and lower mantle domains as a function of pressure and temperature.

To avoid step-changes in fractions at the top and base of the mantle transition zone, transition regions 2.2 GPa wide are implemented. At a reference temperature of 750K, the center of the ol-wd transition is at 11.1 GPa. At the same reference temperature, the center of the postspinel transition is at 26.1 GPa. Clapeyron slopes of 2.4e6 Pa/K and -2.2e6 Pa/K are applied.

Parameters:

Name Type Description Default
pressure float | numpy array

Pressure (Pa)

 required
temperature float | numpy array

Temperature (K)

 required

Returns:

Type Description

A 1D or 2D numpy array containing the effective fractions of upper mantle, transition zone and lower mantle material. If 2D, the fractions[i,j] corresponds to the ith P-T point and jth material.
Source code in terratools/properties/attenuation.py 
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
def mantle_domain_fractions(pressure, temperature):
    """
    This function defines the proportions of
    upper mantle, transition zone, and lower mantle
    domains as a function of pressure and temperature.

    To avoid step-changes in fractions at the top and base of
    the mantle transition zone, transition regions 2.2 GPa wide
    are implemented. At a reference temperature of 750K,
    the center of the ol-wd transition
    is at 11.1 GPa. At the same reference temperature, the center
    of the postspinel transition is at 26.1 GPa. Clapeyron slopes of
    2.4e6 Pa/K and -2.2e6 Pa/K are applied.

    :param pressure: Pressure (Pa)
    :type pressure: float or numpy array

    :param temperature: Temperature (K)
    :type temperature: float or numpy array


    :return:
        A 1D or 2D numpy array containing the effective fractions of
        upper mantle, transition zone and lower mantle material.
        If 2D, the fractions[i,j] corresponds to the ith
        P-T point and jth material.
    """

    P_smooth_halfwidth = 1.1e9
    T_ref = 750.0  # K
    pressure_tztop = 11.1e9 + 2.4e6 * (temperature - T_ref)
    pressure_tzbase = 26.1e9 - 2.2e6 * (temperature - T_ref)

    try:
        fractions = np.zeros(3)
        if pressure < pressure_tztop - P_smooth_halfwidth:
            fractions[0] = 1.0

        elif pressure < pressure_tztop + P_smooth_halfwidth:
            f = (pressure - (pressure_tztop - P_smooth_halfwidth)) / (
                2.0 * P_smooth_halfwidth
            )
            fractions[:2] = [1.0 - f, f]

        elif pressure < pressure_tzbase - P_smooth_halfwidth:
            fractions[1] = 1.0

        elif pressure < pressure_tzbase + P_smooth_halfwidth:
            f = (pressure - (pressure_tzbase - P_smooth_halfwidth)) / (
                2.0 * P_smooth_halfwidth
            )
            fractions[1:] = [1.0 - f, f]

        else:
            fractions[2] = 1.0
    except ValueError:
        fractions = np.zeros((len(pressure), 3))

        f_umtz = (pressure - (pressure_tztop - P_smooth_halfwidth)) / (
            2.0 * P_smooth_halfwidth
        )
        f_tzlm = (pressure - (pressure_tzbase - P_smooth_halfwidth)) / (
            2.0 * P_smooth_halfwidth
        )

        um_idx = np.argwhere(f_umtz <= 0.0)
        umtz_idx = np.argwhere(np.all([f_umtz >= 0.0, f_umtz <= 1.0], axis=0)).T[0]
        tz_idx = np.argwhere(np.all([f_umtz >= 1.0, f_tzlm <= 0.0], axis=0)).T[0]
        tzlm_idx = np.argwhere(np.all([f_tzlm >= 0.0, f_tzlm <= 1.0], axis=0)).T[0]
        lm_idx = np.argwhere(f_tzlm >= 1.0)

        fractions[um_idx, 0] = 1.0
        fractions[umtz_idx, :2] = np.array([1.0 - f_umtz[umtz_idx], f_umtz[umtz_idx]]).T
        fractions[tz_idx, 1] = 1.0
        fractions[tzlm_idx, 1:] = np.array([1.0 - f_tzlm[tzlm_idx], f_tzlm[tzlm_idx]]).T
        fractions[lm_idx, 2] = 1.0

    return fractions