Merge pull request #15 from anowacki/SC-REM

Add SC-REM models of Kemper et al. (2023)

NEWS.md

@@ -1,3 +1,9 @@
+# SeisModels v1.6.0 release notes
+
+## New features
+- Added the SC-REM models (Kemper et al., 2023), giving the confidence
+  bounds.
+
 # SeisModels v1.5.0 release notes
 
 ## New features

data/Earth/SC_REM/plot_SC-REM.jl

@@ -0,0 +1,42 @@
+# Test the implementation of SC_REM_* models by repdroducing Figure 4
+# of Kemper et al. (2023)
+
+using SeisModels
+import Plots
+
+function ribbon(xs, y1s, y2s)
+    [xs; reverse(xs); xs[begin]], [y1s; reverse(y2s); y1s[begin]]
+end
+
+function ribbon_plot!(p, xs, y1s, y2s; kwargs...)
+    x, y = ribbon(xs, y1s, y2s)
+    Plots.plot!(p, Plots.Shape(y, x); kwargs...)
+end
+
+function ribbon_plot!(p, f, m1::SeisModels.SeisModel, m2::SeisModels.SeisModel, rs; kwargs...)
+    vals1 = f.(m1, rs)
+    vals2 = f.(m2, rs)
+    ribbon_plot!(p, r, vals1, vals2; kwargs...)
+end
+
+ribbon_plot(args...; kwargs...) = ribbon_plot!(Plots.plot(), args...; kwargs...)
+
+rs = 0:6371
+prem_colour = :red
+
+ps = Plots.Plot[]
+for f in (density, vp, vs, Qμ)
+    p = f == density ? Plots.plot() :
+        f == Qμ ? Plots.plot(xaxis=:log10, xlim=(80, 1e5), legend=false) :
+        Plots.plot(yticks=nothing, legend=false)
+    push!(ps, p)
+    for (low, high, color) in ((SC_REM_75_LOW, SC_REM_75_HIGH, :lightblue),
+                               (SC_REM_50_LOW, SC_REM_50_HIGH, :blue),
+                               (SC_REM_25_LOW, SC_REM_25_HIGH, :darkblue))
+        ribbon_plot!(ps[end], f, low, high, rs;
+            label="", fillcolor=color, title=f)
+    end
+    # Prem
+    Plots.plot!(ps[end], f.(PREM, rs), rs, lc=prem_colour, ls=:dash, label="PREM")
+end
+Plots.plot(ps...; layout=(1, length(ps)), size=(800,400)) |> display

docs/src/inbuilt_models.md

@@ -11,6 +11,12 @@ EK137
 IASP91
 PREM
 PREM_NOOCEAN
+SC_REM_25_LOW
+SC_REM_25_HIGH
+SC_REM_50_LOW
+SC_REM_50_HIGH
+SC_REM_75_LOW
+SC_REM_75_HIGH
 STW105
 ```
 

docs/src/references.md

@@ -14,12 +14,16 @@
 - Kennett, B.L.N., 2020. Radial earth models revisited.
   Geophysical Journal International 222, 2189–2204.
   doi:[10.1093/gji/ggaa298](https://doi.org/10.1093/gji/ggaa298)
+- Kemper, J., Khan, A., Helffrich, G., van Driel, M., Giardini, D., 2023.
+  Self-consistent models of Earth’s mantle and core from long-period seismic and tidal constraints.
+  _Geophys J Int_ 235, 690–717.
+  doi:[10.1093/gji/ggad254](https://doi.org/10.1093/gji/ggad254)
 - Kustowski, B., Ekström, G., Dziewoński, A., 2008.
   Anisotropic shear-wave velocity structure of the Earth’s mantle: A global model.
   _J Geophys Res-Sol Ea_ 113, B06306.
   doi:[10.1029/2007JB005169](https://doi.org/10.1029/2007JB005169)
 - Nowacki, A., 2020. SeisModels.jl: A Julia package for models of the
-  Earth’s interior. Journal of Open Source Software 5, 2043. doi:[10.21105/joss.02043](https://doi.org/10.21105/joss.02043)
+  Earth’s interior. _Journal of Open Source Software_ 5, 2043. doi:[10.21105/joss.02043](https://doi.org/10.21105/joss.02043)
 - Weber, R.C.,, Lin, P.-Y.,, Garnero, E.J., Williams, Q., Lognonné, P., 2011.
   Seismic Detection of the Lunar Core. _Science_, 331, 309–312.doi:[10.1126/science.1199375](https://doi.org/10.1126/science.1199375)
 

src/Earth/Earth.jl

@@ -6,6 +6,7 @@ This module defines seismic models for the Earth.
 module Earth
 
 using ..SeisModels
+using DelimitedFiles: readdlm
 
 """
 	ALL_MODELS
@@ -14,12 +15,26 @@ List of Earth model instances as `Symbol`s.
 
 If adding new inbuilt Earth models, append the name of the new model to this list
 """
-const ALL_MODELS = (:AK135, :EK137, :IASP91, :PREM, :PREM_NOOCEAN, :STW105)
+const ALL_MODELS = (
+	:AK135,
+	:EK137,
+	:IASP91,
+	:PREM,
+	:PREM_NOOCEAN,
+	:SC_REM_25_LOW,
+	:SC_REM_25_HIGH,
+	:SC_REM_50_LOW,
+	:SC_REM_50_HIGH,
+	:SC_REM_75_LOW,
+	:SC_REM_75_HIGH,
+	:STW105,
+)
 
 include("ak135.jl")
 include("ek137.jl")
 include("iasp91.jl")
 include("prem.jl")
+include("sc_rem.jl")
 include("stw105.jl")
 
 end # module

src/Earth/sc_rem.jl

@@ -0,0 +1,56 @@
+# Generate the SC-REM models
+
+let data_path = joinpath(@__DIR__, "..", "..", "data", "Earth", "SC_REM")
+    for percent in (25, 50, 75), high_low in ("high", "low")
+        model_file = joinpath(data_path, "screm-$(percent)p-$(high_low).dat")
+        model_data = readdlm(model_file, ',', Float64; skipstart=1)
+
+        r = model_data[:,1]
+        ρ = model_data[:,3]
+        vp = model_data[:,4]
+        vs = model_data[:,5]
+        Qκ = model_data[:,6]
+        Qμ = model_data[:,7]
+
+        lower_upper = high_low == "high" ? "Upper" : "Lower"
+
+        model = Symbol(:SC_REM_, percent, :_, uppercase(high_low))
+        model_name = String(model)
+
+        @eval begin
+            """
+            # `$($model_name)`
+
+            $($lower_upper) bound of the $($percent)% credible interval of the SC-REM model of
+            Kemper et al. (2023).
+
+            See also:
+            [`SC_REM_25_LOW`](@ref),
+            [`SC_REM_25_HIGH`](@ref),
+            [`SC_REM_50_LOW`](@ref),
+            [`SC_REM_50_HIGH`](@ref),
+            [`SC_REM_75_LOW`](@ref),
+            [`SC_REM_75_HIGH`](@ref).
+
+            ## References
+
+            - Kemper, J., Khan, A., Helffrich, G., van Driel, M., Giardini, D., 2023.
+              Self-consistent models of Earth’s mantle and core from long-period seismic and tidal constraints.
+              Geophys J Int 235, 690–717.
+              https://doi.org/10.1093/gji/ggad254
+            - Kemper, J., Khan, A., Helffrich, G., Giardini, D., 2023.
+              Self-consistent models of Earth's mantle and core from long-period seismic and tidal constraints (0.2) [Data set].
+              Zenodo.
+              https://doi.org/10.5281/zenodo.10037465
+            """
+            const $model = LinearLayeredModel(
+                r = $r,
+                vp = $vp,
+                vs = $vs,
+                density = $ρ,
+                Qκ = $Qκ,
+                Qμ = $Qμ
+            )
+        end
+    end
+end

test/inbuilt_models.jl

@@ -62,6 +62,22 @@ using SeisModels
             @test vsv(PREM_NOOCEAN, 0, depth=true) == 3.2
         end
 
+        @testset "SC_REM" begin
+            sc_rem_models = (
+                SC_REM_25_LOW, SC_REM_25_HIGH, SC_REM_50_LOW, SC_REM_50_HIGH,
+                SC_REM_75_LOW, SC_REM_75_HIGH
+            )
+            for m in sc_rem_models
+                @test surface_radius(m) == 6371.0
+                @test !isanisotropic(m)
+                @test hasdensity(m)
+                @test hasattenuation(m)
+                I_MR2 = moment_of_inertia(m)/(mass(m)*surface_radius(m)^2*1e6)
+                # Eyeballed from Figure 4 of Kemper et al. (2023)
+                @test 0.3297 <= I_MR2 <= 0.331
+            end
+        end
+
         @testset "STW105" begin
             @test surface_radius(STW105) == 6371.0
             @test isanisotropic(STW105)