A rock physics model for the characterization of organic-rich shale from elastic properties

We construct a rock physics model for organic-rich shales to quantify the effect of kerogen content and kerogen porosity using the Kuster and Toksöz theory and the self-consistent approximation method. Rock physics modeling results show that with the increase of kerogen content and kerogen-related porosity, the velocity and density of shales decrease, and the effect of kerogen porosity becomes more obvious only for higher kerogen content.

Abstract

Kerogen content and kerogen porosity play a significant role in elastic properties of organic-rich shales. We construct a rock physics model for organic-rich shales to quantify the effect of kerogen content and kerogen porosity using the Kuster and Toksöz theory and the self-consistent approximation method. Rock physics modeling results show that with the increase of kerogen content and kerogen-related porosity, the velocity and density of shales decrease, and the effect of kerogen porosity becomes more obvious only for higher kerogen content. We also find that the Poisson’s ratio of the shale is not sensitive to kerogen porosity for the case of gas saturation.

Authors:

Ying Li¹ • Zhi-Qi Guo² • Cai Liu² • Xiang-Yang Li^3,4 • Gang Wang¹

¹College of Computer Science and Technology, Jilin University, Changchun 130012, Jilin, China. ²College of Geo-Exploration Science and Technology, Jilin University, Changchun 130026, Jilin, China. ³British Geological Survey, Edinburgh EH9 3LA, UK. ⁴CNPC Key Laboratory of Geophysical Prospecting, China. University of Petroleum, Beijing 102249, China.

Received: 24 July 2014 / Published online: 12 April 2015

Edited by Jie Hao

Finally, for the seismic reflection responses of an organic-rich shale layer, forward modeling results indicate the fifth type AVO responses which correspond to a negative intercept and a positive gradient. The absolute values of intercept and gradient increase with kerogen content and kerogen porosity, and present predictable variations associated with velocities and density.

1 Introduction

Recently, research on shale reservoirs has become more important in China (Liu et al. 2009; Liu et al. 2011; Zou et al. 2011; Cheng et al. 2013; Chen et al. 2013a; Chen et al. 2013b; Diao 2013; Wang et al. 2013). However, exploration and exploitation of shale reservoirs requires deeper understanding of rock physics relations between reservoir properties and elastic parameters in shales. Compared to rock physics models for conventional reservoirs (Bai et al. 2012; Zhang et al. 2013; Chen et al. 2014), rock physics methods for shales are primarily focusing on the modeling of shale anisotropy resulting from preferred orientations of clay minerals and pores in shale rock matrix (Vernik and Nur 1992; Hornby et al. 1994; Carcione 2000; Wenk et al. 2007; Bobko and Ulm 2008; Lonardelli et al. 2007; Ortega et al. 2009; Dewhurst et al. 2011; Slatt and Abousleiman 2011). Recently, the impact of kerogen on elastic properties and seismic responses of shales has received more interest. Carcione et al. (2011) and Vernik and Milovac (2011) researched and compared rock physics models for organic-rich shales. Zou et al. (2011) analyzed the capacity of reservoirs in terms of nano-porosity in kerogen in shales. Zhu et al. (2012) improved rock physics models for shale gas reservoirs. Sayers (2013b) calculated the effect of kerogen on anisotropy of shales. Zhu et al. (2011) and Sayers (2013a) studied seismic AVO responses based on the proposed rock physics models for shale. Meanwhile, Sun et al. (2013) modeled seismic responses of shale gas reservoirs based on anisotropic wave equations. However, the effect of kerogen content and kerogen porosity on elastic properties and seismic responses of shale has not been well investigated.

According to the observations from SEM images, kerogen plays a more significant role in elastic properties than ever expected due to its low density (Passey et al. 2010); moreover, hydrocarbon-filled pores in such organic matter generated during different maturity levels will enhance the effect of kerogen. Therefore, key parameters for the characterization of organic-rich shales should include both kerogen content and kerogen porosity. In this study, our objective is to build a rock physics workflow for the organic-rich shales to quantify the effects of such two factors on elastic properties, and therefore on the seismic responses.

2 Microstructural model of shale

SEM images reveal the presence of a nano-porosity system contained within the organic matter (Passey et al. 2010). It is estimated that these pores comprise the majority of the total porosity space and may be hydrocarbon saturated during most of the thermal maturation processes. The schematics in Fig. 1 illustrate how 5 wt% (weight percentage) solid kerogen in Fig. 1a may correspond to 10 vol% (volume percentage) solid kerogen in Fig. 1b because the density of kerogen is about half that of average rock minerals, and then may eventually take up 20 vol% in a shale rock in Fig. 1c if 50 % of the organic matter volume is occupied by hydrocarbon-saturated pores.

$Fig. 1 Schematics demonstrate 5 wt% solid kerogen in (a) corresponds to 10 vol% kerogen in (b) due to the lower density of kerogen compared to that of average minerals in the rock matrix, and the generation of hydrocarbon-filled pores during different maturity levels results in the expansion of the volumetric fraction of the mixture of kerogen and hydrocarbons and may take up 20 vol% if half of the spaces are occupied by pores in kerogen as shown in (c). (After Passey et al. 2010)$

Fig. 1 Schematics demonstrate 5 wt% solid kerogen in (a) corresponds to 10 vol% kerogen in (b) due to the lower density of kerogen compared to that of average minerals in the rock matrix, and the generation of hydrocarbon-filled pores during different maturity levels results in the expansion of the volumetric fraction of the mixture of kerogen and hydrocarbons and may take up 20 vol% if half of the spaces are occupied by pores in kerogen as shown in (c). (After Passey et al. 2010).

Therefore, a small amount of kerogen in weight percentage corresponds to a relevant higher volume percentage which has a significant impact on elastic properties and seismic responses of shales. The objective of this study is to build a rock physics model for organic-rich shales to quantify the effect of kerogen content and kerogen porosity.

Figure 2 illustrates the evolution of solid and fluid constituents in a source rock associated with various maturity levels. For an immature source rock, kerogen has no pores generated inside (Scenario A). As the source rock becomes mature, hydrocarbon-filled pores are generated in organic matter (Scenario B). In addition, the kerogen content in source rock may be different (Scenario C), and as the source rock becomes more mature, more hydrocarbon-filled pores will be generated in organic matter (Scenario D). According to the above analysis, kerogen content and porosity associated with hydrocarbon-filled pores in kerogen are two important factors for the characterization of a source rock, so we define the two factors as K and φ_k, respectively, in this work.

Fig. 2 Schematics illustrate the evolution of solid and fluid constituents in a source rock associated with different maturity levels. Scenario A illustrates an immature source rock with no pores generated in kerogen. Scenario B demonstrates a mature source rock with hydrocarbon-filled pores generated in kerogen. In addition, the kerogen content K in a source rock may be different as shown in Scenario C, and hydrocarbon-filled porosity φ_k may also vary as shown in Scenario D.

3 Rock physics workflow for organic-rich shales

Figure 3 illustrates a rock physics workflow we built to qualify the effect of K and φ_k on elastic properties of the organic-rich shale.

In the rock physics model, the mixture of kerogen and hydrocarbon-filled pores has non-zero shear modulus (Passey et al. 1990), and the elastic properties of such a mixture are calculated using the method proposed by Carcione (2000) based on the theory of Kuster and Toksöz (1974). After that, the elastic properties of the shale consisting of clay and non-clay minerals, the mixture of kerogen and hydrocarbons, and fluid-saturated pores in the shale rock matrix are calculated using the self-consistent approximation (SCA) method given by Berryman (1980).

The constructed rock physics model is based on isotropic assumption and focuses on quantifying the effect of kerogen on elastic properties of shales. The results model the elastic parameters of shales in the vertical direction (one-dimensional problem) and can be applied to analyze the relation between kerogen and elastic properties in a vertical well.

Fig. 3 A schematic illustrates the shale rock physics workflow to model the effect of kerogen content and kerogen porosity. Effective medium theories are employed in the workflow.

As indicated in Fig. 4, we consider that the kerogen/oil/gas/water mixture is composed of oil/gas/water bubbles embedded in a kerogen matrix. The saturations are S_k, S_o, S_g, and S_w for kerogen, oil, gas, and water, respectively, and we have S_k + S_o + S_g + S_w = 1. We define the kerogen-related porosity as φ_k = S_o + S_g + S_w, which is the proportion of the oil/gas/water bubbles in the mixture. The bulk modulus K_ogw of the oil/gas/water mixture is obtained by using Wood’s equation (Wood 1955):

where K_o, K_g, and K_w are the moduli of oil, gas, and water, and are calculated by

Fig. 4 A schematic indicates the mixture of kerogen, oil, gas, and water, and the relationship of the associated saturations.

Then, the stiffnesses of the mixtures of kerogen and gas/oil/water can be calculated by using the method of Kuster and Toksöz (1974):

where the superscript m indicates the mixture, K_k and μ_k are the bulk and shear moduli of kerogen, and φ_k is the kerogen porosity. Equations (3) and (4) are used for the solutions of c₁₃ and c₅₅, which are two components in the stiffness matrix c_ij (i, j = 1, 2,…, 6), and c₁₃ = c₁₁−2c₅₅, c₁₁ for P-wave in the horizontal direction, and c₅₅ for S-wave in the vertical direction.

By solving Eqs. (3) and (4), bulk and shear moduli of the mixture of kerogen and gas/oil/water can be obtained by

The superscript m denotes the mixture of kerogen/oil/gas/water. Given the density of kerogen (ρ_k), oil (ρ_o), gas (ρ_g), and water (ρ_w), the density of the mixture is ρ^(m) = S_kρ_k + S_oρ_o + S_gρ_g + S_wρ_w.

After obtaining the elastic properties of the kerogen and gas/oil/water mixture, we employ the self-consistent approximation (SCA) method to calculate the bulk and shear moduli of the shale rock. The shale is composed of quartz, clay, carbonate, the mixture of kerogen and hydrocarbons, and fluid-saturated pores in the shale rock matrix. Elastic properties of the constituents and oil, gas, and water are given in Table 1 (Carcione et al. 2011; Mavko et al. 2009).

The shale is composed of quartz, clay, carbonate, the mixture of kerogen and hydrocarbons, and fluid-saturated pores in the shale rock matrix. Elastic properties of the constituents and oil, gas, and water are given in Table 1 (Carcione et al. 2011; Mavko et al. 2009).

Table 1 Material properties (Carcione 2000; Mavko et al. 2009).

The SCA method estimates the self-consistent elastic modulus K_SC and μ_SC of rocks for a given n phases of mineralogy and pore space:

$Each j indicates a phase of mineralogy or pore space with a corresponding volume fraction fj, and bulk (Kj) and shear (μj) modulus. The factors β*j and ζ*j describe the geometry of an inclusion made of phase j within a background medium (denoted with subscript b)$

Each j indicates a phase of mineralogy or pore space with a corresponding volume fraction f_j, and bulk (K_j) and shear (μ_j) modulus. The factors β^*j and ζ^*j describe the geometry of an inclusion made of phase j within a background medium (denoted with subscript b), and are given as follows:

Equations (7) and (8) represent the geometries of pores with aspect ratio smaller than 1 (Mavko et al. 2009).

The implementation of Eqs. (5)–(8) requires volumetric percentage, bulk and shear modulus, and aspect ratio for each inclusion, and the outputs are bulk and shear modulus of rocks. For each inclusion j, background medium b corresponds to the inclusion other than j. Equations (5) and (6) are coupled, and therefore need to be solved by a simultaneous iteration method. Equations (7) and (8) represent the geometries of pores with aspect ratio smaller than 1 (Mavko et al. 2009).

4 The effect of K and φ_k on elastic and geo-mechanical properties

According to the observations and analysis by Passey et al. 1990, constituents of the organic-rich shale consist of clay, quartz, dolomite, hydrocarbon-saturated kerogen, and fluid-saturated pores in shale matrix. In the rock physics modeling, we set the mineralogical volumetric contents of clay and dolomite to 0.3 and 0.15, respectively, and the porosity associated with shale matrix is 0.05. We assume the contents of kerogen and quartz have a negative correlation, that is, an increase in kerogen content K corresponds to a decrease in quartz content. So, if K increases from 0 to 0.3, quartz content decreases from 0.5 to 0.2 correspondingly. Associated elastic properties are given in Table 1. For the constituents in shales, in this study, n in Eqs. (5)–(8) corresponds to n = 5.

We assume the pores in kerogen are totally gas saturated to imitate high maturity level, and the pores in shale matrix are water saturated. Another key factor is the aspect ratio of the pores in the shale matrix, and we set it as 0.1 according to the study on the Barnett Shale (Guo et al. 2013a, b). According to the work on shale rock physics by Carcione (2000), we assume spherical geometry of pores in kerogen.

Then based on the rock physics model in Fig. 3 and properties in Table 1, we calculate elastic and geo-mechanical properties of the organic-rich shale along with the variations in kerogen content K and kerogen-related porosity φ_k. Figure 5a, b, and c illustrate the obtained P- and S-wave velocity V_p, and V_s, and density ρ of the shale for varied kerogen content K and kerogen porosity φ_k. We can see that the values of V_p, V_s, and ρ decrease significantly as K increases from 0 to 0.3 and φ_k from 0 to 0.5. For example, V_p has a value about 3370 km/s for K = 0.1 and φ_k = 0.1, and this value decreases to 2783 km/s for K = 0.2 and φ_k = 0.2.

Fig. 5 Calculated templates showing a P-wave velocity V_p, b S-wave velocity V_s, and c Density ρ for the variations in kerogen content K and kerogen-related porosity φ_k in the shale.

We also find the correlation between the effect of K and that of φ_k on elastic properties, which impact elastic properties of shale simultaneously. Due to low velocities and density, the increase in K decreases velocities and density, and the presence of φ_k enhances such impacts (as shown in Fig. 1). For K < 0.05, the effect of φ_k on the elastic properties is not obvious; for K > 0.05, the increase in φ_k significantly decreases velocities and density of shales.

Figure 6a and b illustrate the calculated Young’s modulus E and Poisson’s ratio v for the same range of K and φ_k as those in Fig. 5 (Zong et al. 2012). The two geo-mechanical properties are critical for the evaluation of the brittleness index of shales, as higher E and lower v indicate higher potential for shales to be fractured during hydraulic stimulation (Rickman et al. 2008). In Fig. 6b, we find that in the range of 0 < K < 0.2, for each φ_k, the Poisson’s ratio v increases with K, while for each K, v nearly keeps constant for varied φ_k. This may imply that v is not sensitive to the gas-saturated porosity in kerogen.

Fig. 6 Calculated templates showing geo-mechanical properties a Young’s modulus E and b Poisson’s ratio v for the variations in K and φ_k in the shale.

5 The effect of K and φ_k on seismic AVO responses

In Fig. 7, we demonstrate the seismic AVO responses resulting from the variations in K and φ_k. For the geological model shown in Fig. 7a, b illustrates how the variations in K and φ_k affect the acoustic impedance I_p of the shale based on the rock physics model. We consider two cases as indicated by blue and red solid circles in Fig. 7b. The blue solid circle represents values of K and φ_k of 0.05 and 0.1, respectively, while the red solid circle corresponds to the values of 0.15 and 0.3.

The properties of the overlaid chalk are V_p = 4100 m/s, V_s = 2500 m/s, and ρ = 2800 kg/m³. We calculate AVO responses using numerical solutions of the Zoeppritz equations and display the results in Fig. 7c by the same colors for the cases in Fig. 7b.

Fig. 7 a The model in which a chalk layer overlays on a shale layer. b Calculated acoustic impedance I_p of the shale varying with K and φ_k. Two cases with different properties are highlighted. c Seismic AVO responses corresponding to the two cases in (b), with the same color scheme.

We observe significant differences in AVO responses resulting from the variations in K and φ_k of the shale. Reflection amplitudes are negative for normal incidence, and their absolute values decrease with increasing incidence angles, which corresponds to the fifth type AVO responses (Roger and Robert 2010). Meanwhile, normal incidence amplitudes and AVO gradients have relatively higher values for higher kerogen content and kerogen porosity for low wave impedance shale (in red).

Figure 8 illustrates the calculated AVO intercept P and gradient G corresponding to the model in Fig. 7a for the variations in K and φ_k. In Fig. 8a, we can see that the values of intercept P are negative and its absolute values increase with increasing K and φ_k, which has similar trends as the velocities and density as shown in Fig. 5. In Fig. 8b, the AVO gradient G represents positive and increases with increasing K and φ_k. These predictable trends of AVO responses can be used for the characterization of the variations in K and φ_k.

Fig. 8 Templates illustrating a AVO intercept P and b gradient G for the variations in K and φ_k, corresponding to the geological model and acoustic impedance variations in Fig. 7a and b.

6 Conclusions

We have developed a rock physics model based on the effective medium theories for the characterization of the effects of kerogen content K and kerogen porosity φ_k on elastic and geo-mechanical properties in the organic-rich shale. The parameter φ_k indicates the proportion of hydrocarbon-saturated pores generated within kerogen. Rock physics modeling results show that the increases in K and φ_k significantly reduce the values of V_p, V_s, and ρ and the effects of K and φ_k are correlated with each other. Due to low velocities and density, the increase in K decreases velocities and density and the presence of φ_k enhances such impacts. For the case of K < 0.05, in this study, the effect of φ_k on the elastic properties is not obvious, while for K > 0.05, the increase in φ_k significantly decreases velocities and density of shales. In the range of 0 < K < 0.2, Poisson’s ratio v is not sensitive to the gas-saturated porosity in kerogen. AVO modeling reveals that for the model in this work, reflection amplitudes for normal incidence are negative and the corresponding absolute values decrease with increasing incidence angles, which represents the fifth type AVO responses. For increasing K and φ_k, the absolute values of intercept P and gradient G increase. The predictable trends in AVO responses can be used for the characterization of the variations in K and φ_k.

Acknowledgments

This work is supported by the National Natural Science Foundation of China under Grants U1262208 and the National Natural Science Foundation of China under Grants 41404090. We thank two anonymous reviewers for their constructive suggestions to improve the paper.

References

Bai JY, Song ZX, Su L, et al. Error analysis of shear-velocity prediction by the Xu-White model. Chin J Geophys. 2012;55(2):589–95 (in Chinese).
Berryman JG. Long-wavelength propagation in composite elastic media. J Acoust Soc Am. 1980;68(6):1809–31.
Bobko C, Ulm FJ. The nano-mechanical morphology of shale. Mech Mater. 2008;40(4–5):318–37.
Carcione JM. A model for seismic velocity and attenuation in petroleum source rocks. Geophysics. 2000;65(4):1080–92.
Carcione JM, Helle HB, Avseth P. Source-rock seismic-velocity models: Gassmann versus Backus. Geophysics. 2011;76(5): N37–45.
Chen HZ, Yin XY, Gao CG, et al. AVAZ inversion for fluid factor based on fracture anisotropic rock physics theory. Chin J Geophys. 2014;57(3):968–78 (in Chinese).
Chen P, Zhang MQ, Xu YZ, et al. The shale reservoir characteristic of Dalong Formation Upper Permian in Chaohu-Jingxian, Lower Yangtze area. Acta Petrologica Sinica. 2013a;29(8):2925–35 (in Chinese).
Chen WL, Zhou W, Luo P, et al. Analysis of the shale gas reservoir in the Lower Silurian Longmaxi Formation, Changxin 1 well, Southeast Sichuan Basin, China. Acta Petrologica Sinica. 2013b;29(3):1073–86 (in Chinese).
Cheng LX, Wang YJ, Chen HD, et al. Sedimentary and burial environment of black shales of Sinian to early Palaeozoic in upper Yangtze region. Acta Petrologica Sinica. 2013;29(8): 2906–12 (in Chinese).
Dewhurst DN, Siggins AF, Sarout J, et al. Geomechanical and ultrasonic characterization of a Norwegian Sea shale. Geo- physics. 2011;76(3):WA101–11.
Diao HY. Rock mechanical properties and brittleness evaluation of shale reservoir. Acta Petrologica Sinica. 2013;29(9):3300–6 (in Chinese).
Guo ZQ, Li XY, Chapman M, et al. Microstructure characterization and S-wave velocity prediction in the Barnett Shale formation. In: 75th EAGE conference & exhibition, extended abstracts. 2013a. P06.
Guo ZQ, Li XY, Liu C, et al. A shale rock physics model for analysis of brittleness index, mineralogy and porosity in the Barnett Shale. J Geophys Eng. 2013b;10(2):025006. doi:10.1088/1742-2132/10/2/025006.
Hornby BE, Schwartz LM, Hudson JA. Anisotropic effective medium modeling of the elastic properties of shales. Geophysics. 1994;59(10):1570–83.
Kuster GT, Tokso¨ z MN. Velocity and attenuation of seismic waves in two-phase media, Part I: theoretical formulations. Geophysics. 1974;39(5):587–606.
Liu SG, Ma WX, Luba J, et al. Characteristics of the shale gas reservoir rocks in the Lower Silurian Longmaxi Formation, East Sichuan Basin, China. Acta Petrologica Sinica. 2011;27(8):2239–52 (in Chinese).
Liu ZJ, Meng QT, Liu R, et al. Geochemical characteristics of oil shale of Eocene Jijuntun Formation and its geological sig- nificance, Fushun Basin. Acta Petrologica Sinica. 2009;25(10):2340–50 (in Chinese).
Lonardelli I, Wenk HR, Ren Y. Preferred orientation and elastic anisotropy in shales. Geophysics. 2007;72(2):D33–40.
Mavko G, Mukerji T, Dovrkin J. The rock physics handbook: tools for analysis of porous media. Cambridge: Cambridge University Press; 2009. p. 185–90.
Ortega JA, Ulm FJ, Abousleiman Y. The nanogranular acoustic signature of shale. Geophysics. 2009;74(3):D65–84.
Passey QR, Bohacs KM, Esch WL, et al. From oil-prone source rock to gas-producing shale reservoir—geologic and petrophysical characterization of unconventional shale-gas reservoirs. SPE Conference Paper. 2010. p. 131350.
Passey QR, Creaney S, Kulla JB, et al. A practical model for organic richness from porosity and resistivity logs. AAPG Bull. 1990;74(12):1777–94.
Rickman R, Mullen M, Petre E, et al. A practical use of shale petrophysics for simulation design optimization: All shale plays are not clones of the Barnett Shale. SEP Conference Paper. 2008. p. 115258.
Roger AY, Robert DL. A comprehensive AVO classification. Lead Edge. 2010;22(10):1030–7.
Sayers CM. The effect of kerogen on the AVO response of organic- rich shales. Lead Edge. 2013a;32(12):1514–9.
Sayers CM. The effect of kerogen on the elastic anisotropy of organic-rich shales. Geophysics. 2013b;78(2):D65–74.
Slatt RM, Abousleiman Y. Merging sequence stratigraphy and geomechanics for unconventional gas shales. Lead Edge. 2011;30(3):274–82.
Sun WJ, Fu LT, Guan XZ, et al. A study on anisotropy of shale using seismic forward modeling in shale gas exploration. Chin J Geophys. 2013;56(3):961–70 (in Chinese).
Vernik L, Milovac J. Rock physics of organic shales. Lead Edge. 2011;30(3):318–23.
Vernik L, Nur A. Ultrasonic velocity and anisotropy of hydrocarbon source rocks. Geophysics. 1992;57(5):727–35.
Wang LB, Jiu K, Zeng WT, et al. Characteristics of Lower Cambrian marine black shales and evaluation of shale gas prospective area in Qianbei area, Upper Yangtze region. Acta Petrologica Sinica. 2013;29(9):3263–78 (in Chinese).
Wenk HR, Lonardelli I, Franz H, et al. Preferred orientation and elastic anisotropy of illite-rich shale. Geophysics. 2007;72(2): E69–75.
Wood AB. A textbook of sound. New York: Macmillan Co.; 1955. Zhang GZ, Chen HZ, Wang Q, et al. Estimation of S-wave velocity and anisotropic parameters using fractured carbonate rock physics model. Chin J Geophys. 2013;56(5):1707–15 (in Chinese).
Zhu Y, Liu E, Martinez A, et al. Understanding geophysical responses of shale-gas plays. Lead Edge. 2011;30(3):332–8.
Zhu Y, Xu SY, Payne M, et al. Improved rock-physics model for shale gas reservoirs. In: 82nd annual international meeting, SEG, expanded abstracts. 2012. pp. 1–5.
Zong ZY, Yin XY, Zhang F, et al. Reflection coefficient equation and pre-stack seismic inversion with Young’s modulus and Poisson ratio. Chin J Geophys. 2012;55(11):3786–94 (in Chinese).
Zou CN, Zhu RK, Bai B, et al. First discovery of nano-pore throat in oil and gas reservoir in China and its scientific value. Acta Petrologica Sinica. 2011;27(6):1857–64 (in Chinese).

Author:

Zhi-Qi Guo zhiqiguo@aliyun.com

Open AccessThis article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.