agd.Metrics.Seismic.thomsen_data

  1# Copyright 2020 Jean-Marie Mirebeau, University Paris-Sud, CNRS, University Paris-Saclay
  2# Distributed WITHOUT ANY WARRANTY. Licensed under the Apache License, Version 2.0, see http://www.apache.org/licenses/LICENSE-2.0
  3
  4"""
  5We reproduce in this file list of VTI examples taken from "Weak elastic anisotropy" (Thomsen, 1986), 
  6Units are : 
  7- Vp,Vs : m/s
  8- 'ε','η','δ','γ' : dimensionless. (Originally, η is listed as δ^star.)
  9- ρ : g/cm^3
 10
 11
 12We also provide conversion utilities
 13"""
 14
 15
 16import numpy as np
 17from collections import OrderedDict,namedtuple
 18
 19# ------ Thomsen Elastica materials -----
 20
 21# At the time of writing, namedtuple requires normalized unicode characters.
 22# This is why the 'ε' is used instead of 'ϵ'.
 23# In jupyter notebooks : ε = \varepsilon, and ϵ = \epsilon. On the Greek keyboard, e -> ε.
 24# import unicodedata; 'ε' == unicodedata.normalize('NFKC','ϵ')
 25# https://stackoverflow.com/a/30505623/12508258
 26
 27
 28HexagonalMaterial = namedtuple('HexagonalMaterial',['c11','c12','c13','c33','c44'])
 29
 30def get_δ(Vp,Vs,ϵ,η):
 31	"""
 32	Reconstructs Thomsen parameter δ, based on the other parameters, 
 33	and the formula in the Thomsen paper.
 34	(Note : some of the published ThomsenData appears to be inconsistent in this regard.)
 35	"""
 36	return 0.5*(ϵ + η/(1.-Vs**2/Vp**2))
 37
 38TEM_data = namedtuple('ThomsenElasticMaterial_data',['Vp','Vs','ε','η','δ','γ','ρ'])
 39class ThomsenElasticMaterial(TEM_data):
 40	@classmethod
 41	def units(cls): return cls('m/s','m/s',1,1,1,1,'g/cm^3')
 42
 43	@classmethod
 44	def from_hexagonal(cls,hex,ρ):
 45		"""
 46		Produces the Thomsen parameters, from the coefficients 'hex' of 
 47		a Hooke tensor with hexagonal symmetry, and the density 'ρ'.
 48		"""
 49		c11,c12,c13,c33,c44 = hex
 50		c66 = (c11-c12)/2
 51
 52		Vp = np.sqrt(c33)
 53		Vs = np.sqrt(c44)
 54
 55		ε = (c11-c33)/(2*c33)
 56		γ = (c66-c44)/(2*c44)
 57		η = (2*(c13+c44)**2 - (c33-c44)*(c11+c33-2*c44))/(2*c33**2)
 58		δ = get_δ(Vp,Vs,ϵ,η)
 59
 60		return cls(Vp,Vs,ε,η,δ,γ,ρ)
 61
 62	def to_hexagonal(self):
 63		"""
 64		Returns the coefficients of the reduced Hooke tensor, with hexagonal symmetry, and the density.
 65		Units. Reduced hooke tensor : (m/s)^2, density : g/cm^3
 66		"""
 67		Vp,Vs,ε,η,δ,γ,ρ = self
 68		c33 = Vp**2
 69		c44 = Vs**2
 70
 71		c11 = 2*c33*ε+c33
 72		c66 = 2*c44*γ+c44
 73		c12 = c11-2*c66
 74		x = η*c33**2 + 0.5*(c33-c44)*(c11+c33-c44)
 75		c13 = np.sqrt(x) - c44 # Assuming c13+c44 >= 0
 76		return HexagonalMaterial(c11,c12,c13,c33,c44),ρ
 77
 78# ----- Thomsen geometric materials -----
 79
 80# If one is only interested in the speed of the pressure wave in the medium, 
 81# which defines the TTI geometric metric, then some of the Thomsen coefficients can be discarded.
 82# We refer to the geometric data as a Thomsen geometric material.
 83
 84
 85TGM_data=namedtuple('ThomsenGeometricMaterial_data',['Vp','Vs','ε','δ'])
 86class ThomsenGeometricMaterial(TGM_data):
 87	@classmethod
 88	def units(cls): return cls('m/s','m/s',1,1)
 89	
 90	@classmethod
 91	def from_c(cls,c11,c13,c33,c44,ρ=1):
 92		"Thomsen coefficients transformation (c11,c13,c33,c44) -> (Vp,Vs,ε,δ)"
 93		Vp = np.sqrt(c33/ρ)
 94		Vs = np.sqrt(c44/ρ)
 95		ε = (c11-c33)/(2*c33)
 96		δ = ((c13+c44)**2-(c33-c44)**2)/(2*c33*(c33-c44))
 97		return cls(Vp,Vs,ε,δ)
 98	
 99	def to_c(self):
100		"Thomsen coefficients transformation (Vp,Vs,ε,δ) -> (c11,c13,c33,c44)"
101		Vp,Vs,ε,δ = self
102		c33 = Vp**2
103		c44 = Vs**2
104		c11 = c33*(2*ϵ+1)
105		c13 = np.sqrt(2*c33*(c33-c44)*δ+(c33-c44)**2)-c44
106		return c11,c13,c33,c44
107
108	@classmethod
109	def from_Elastic(cls,thomsenElasticMaterial):
110		"Extracts the geometric coefficients (relevant for the pressure wave velocity)"
111		Vp,Vs,ε,_,δ,_,_ = thomsenElasticMaterial
112		return cls(Vp,Vs,ε,δ)
113
114# ----- Some routines based on the minimal description c11,c13,c33,c44 -----
115
116
117def is_definite(c11,c13,c33,c44): 
118	"""Wether the coefficients define a positive definite hooke tensor, hence a convex inner sheet."""
119	return (c11>0) & (c33>0) & (c11*c33>c13**2) & (c44>0) 
120
121def is_separable(c11,c13,c33,c44): 
122	"""Wether the inner and outer sheets are well separated"""
123	return (c11>c44) & (c33>c44) & (c11*c33>c13**2) & (c44>0) & (c13+c44>0)
124
125def is_second_sheet_convex(c11,c13,c33,c44):
126	"""
127	Wether the second sheet defined by the coefficients is convex.
128	Formula obtained using formal computing software. 
129	No proof, but checked on Thomsen's data and some examples by hand.
130
131	"""
132	det = c13**2 - c11*c33 + 4*c13*c44 + 4*c44**2
133	slope_z = -c13**2 + c11*c33 - 2*c13*c44 - c44*(c33 + c44)
134	slope_x = -c13**2 + c11*(c33 - c44) - 2*c13*c44 - c44**2
135	inflexion = -c13**6 + 4*c11**3*c33**2*(c33 - c44) - 6*c13**5*c44 + 2*c11*c33*(7*c33 - 9*c44)*c44**3 - c33**2*c44**4 + c13**4*(6*c11*c33 - 2*c11*c44 - 2*c33*c44 - 9*c44**2) + 4*c13**3*c44*(6*c11*c33 - 2*c11*c44 - 2*c33*c44 + c44**2) - c11**2*c44*(4*c33**3 + 13*c33**2*c44 - 14*c33*c44**2 + c44**3) - 2*c13*c44*(c33*(c33 - 2*c44)*c44**2 + c11**2*(-3*c33 + c44)**2 - 2*c11*c44*(3*c33**2 - 4*c33*c44 + c44**2)) -  c13**2*(c11**2*(-3*c33 + c44)**2 + c44**2*(c33**2 + 6*c33*c44 - 12*c44**2) + 2*c11*c44*(-3*c33**2 - 8*c33*c44 + 3*c44**2))
136#	print(f"{det=}, {slope_z=}, {slope_x=}, {inflexion=}")
137	return np.where(det>0, (slope_x>0) & (slope_z>0), inflexion<0)
138
139
140# ----- Thomsen tabulated data -----
141
142TEM = ThomsenElasticMaterial
143
144ThomsenData = OrderedDict([
145	("Taylor sandstone",                     TEM(3368,1829,0.110,-0.127,-0.035,0.255,2.500)),
146	("Mesaverde (4903) mudshale",            TEM(4529,2703,0.034,0.250,0.211,0.046,2.520)),
147	("Mesaverde (4912) immature sandstone",  TEM(4476,2814,0.097,0.051,0.091,0.051,2.500)),
148	("Mesaverde (4946) immature sandstone",  TEM(4099,2346,0.077,-0.039,0.010,0.066,2.450)),
149	("Mesaverde (5469.5) silty sandstone",   TEM(4972,2899,0.056,-0.041,-0.003,0.067,2.630)),
150	("Mesaverde (5481.3) immature sandstone",TEM(4349,2571,0.091,0.134,0.148,0.105,2.460)),
151	("Mesaverde (5501) clayshale",           TEM(3928,2055,0.334,0.818,0.730,0.575,2.590)),
152	("Mesaverde (5555.5) immature sandstone",TEM(4539,2706,0.060,0.147,0.143,0.045,2.480)),
153	("Mesaverde (5566.3) laminated siltstone",TEM(4449,2585,0.091,0.688,0.565,0.046,2.570)),
154	("Mesaverde (5837.5) immature sandstone",TEM(4672,2833,0.023,-0.013,0.002,0.013,2.470)),
155	("Mesaverde (5858.6) clayshale",         TEM(3794,2074,0.189,0.154,0.204,0.175,2.560)),
156	("Mesaverde (6423.6) calcareous sandstone",TEM(5460,3219,0.000,-0.345,-0.264,-0.007,2.690)),
157	("Mesaverde (6455.1) immature sandstone",TEM(4418,2587,0.053,0.173,0.158,0.133,2.450)),
158	("Mesaverde (6542.6) immature sandstone",TEM(4405,2542,0.080,-0.057,-0.003,0.093,2.510)),
159	("Mesaverde (6563.7) mudshale",          TEM(5073,2998,0.010,0.009,0.012,-0.005,2.680)),
160	("Mesaverde (7888.4) sandstone",         TEM(4869,2911,0.033,0.030,0.040,-0.019,2500)),
161	("Mesaverde (7939.5) mudshale",          TEM(4296,2471,0.081,0.118,0.129,0.048,2.660)),
162
163	("Mesaverde shale (350)",                TEM(3383,2438,0.065,-0.003,0.059,0.071,2.35)),
164	("Mesaverde sandstone (1582)",           TEM(3688,2774,0.081,0.010,0.057,0.000,2.73)),
165	("Mesaverde shale (1599)",               TEM(3901,2682,0.137,-0.078,-0.012,0.026,2.64)),
166	("Mesaverde sandstone (1958)",           TEM(4237,3018,0.036,-0.037,-0.039,0.030,2.69)),
167	("Mesaverde shale (1968)",               TEM(4846,3170,0.063,-0.031,0.008,0.028,2.69)),
168	("Mesaverde sandstone (3512)",           TEM(4633,3231,-0.026,-0.004,-0.033,0.035,2.71)),
169	("Mesaverde shale (3511)",               TEM(4359,3048,0.172,-0.088,0.000,0.157,2.81)),
170	("Mesaverde sandstone (3805)",           TEM(3962,2926,0.055,-0.066,-0.089,0.041,2.87)),
171	("Mesaverde shale (3883)",               TEM(3749,2621,0.128,-0.025,0.078,0.100,2.92)),
172	("Dog Creek shale",                      TEM(1875,826,0.225,-0.020,0.100,0.345,2.000)),
173	("Wills Point shale - 1",                TEM(1058,387,0.215,0.359,0.315,0.280,1.800)),
174	("Wills Point shale - 2",                TEM(4130,2380,0.085,0.104,0.120,0.185,2640)),
175	("Cotton Valley shale",                  TEM(4721,2890,0.135,0.172,0.205,0.180,2.640)),
176	("Pierre shale - 1",                     TEM(2074,869,0.110,0.058,0.090,0.165,2.25)), # rho ?
177	("Pierre shale - 2",                     TEM(2106,887,0.195,0.128,0.175,0.300,2.25)), # rho ?
178	("Pierre shale - 3",                     TEM(2202,969,0.015,0.085,0.060,0.030,2.25)), # rho ?
179	("shale (5000) - 1",                     TEM(3048,1490,0.255,-0.270,-0.050,0.480,2.420)),
180
181	("shale (5000) - 2",                     TEM(3377,1490,0.200,-0.282,-0.075,0.510,2.420)),
182	("Oil Shale",                            TEM(4231,2539,0.200,0.000,0.100,0.145,2.370)),
183	("Green River shale - 1",                TEM(4167,2432,0.040,-0.013,0.010,0.145,2.370)),
184	("Green River shale - 2",                TEM(4404,2582,0.025,0.056,0.055,0.020,2.310)),
185	("Berea sandstone - 1",                  TEM(4206,2664,0.002,0.023,0.020,0.005,2.140)),
186	("Berea sandstone - 2",                  TEM(3810,2368,0.030,0.037,0.045,0.030,2.160)),
187	("Green River shale - 3",                TEM(3292,1768,0.195,-0.45,-0.220,0.180,2.075)),
188	("Lance sandstone",                      TEM(5029,2987,-0.005,-0.032,-0.015,0.005,2.430)),
189	("Ft. Union siltstone",                  TEM(4877,2941,0.045,-0.071,-0.045,0.040,2.600)),
190	("Timber Mtn tuff",                      TEM(4846,1856,0.020,-0.003,-0.030,0.105,2.330)),
191	("Muscovite crystal",                    TEM(4420,2091,1.12,-1.23,-0.235,2.28,2.79)),
192	("Quartz crystal (hexag. approx.)",      TEM(6096,4481,-0.096,0.169,0.273,-0.159,2.65)),
193	("Calcite crystal (hexag. approx.)",     TEM(5334,3353,0.369,0.127,0.579,0.169,2.71)),
194	("Biotite crystal",                      TEM(4054,1341,1.222,-1.437,-0.388,6.12,3.05)),
195	("Apatite crystal",                      TEM(6340,4389,0.097,0.257,0.586,0.079,3.218)),
196	("Ice I crystal",                        TEM(3627,1676,-0.038,-0.10,-0.164,0.031,1.064)),
197	("Aluminium-lucite composite",           TEM(2868,1350,0.97,-0.89,-0.09,1.30,1.86)),
198
199	("Sandstone-shale",                      TEM(3009,1654,0.013,-0.010,-0.001,0.035,2.34)),
200	("SS-anisotropic shale",                 TEM(3009,1654,0.059,-0.042,-0.001,0.163,2.34)),
201	("Limestone-shale",                      TEM(3306,1819,0.134,-0.094,0.000,0.156,2.44)),
202	("LS-anisotropic shale",                 TEM(3306,1819,0.169,-0.123,0.000,0.271,2.44)),
203	("Anisotropic shale",                    TEM(2745,1508,0.103,-0.073,-0.001,0.345,2.34)),
204	("Gas sand-water sand",                  TEM(1409,780,0.022,-0.002,0.018,0.004,2.03)),
205	("Gypsum-weathered material",            TEM(1911,795,1.161,-1.075,-0.140,2.781,2.35))
206])