Assessment of Paleocene to lower Oligocene formations and basement to estimate the potential hydrocarbon reservoirs using seismic inversion: a case study in the Upper Assam Shelf, India

The upper Assam shelf is a self-slope basin in north-east India, filled with nearly 7 km of sedimentary rocks of tertiary period with the granite basement and various oil fields along the border of the Naga thrust. The major producing fields are structural and strati-structural. The study area is placed in between the Mikir hills and Naga thrust. The objective of the study is to identify potential hydrocarbon reservoir zones in the geologically complex south upper of the Assam shelf using estimates of acoustic impedance and porosity derived by 3D post-stack seismic inversion. Well data, such as sonic velocity and density logs, from two wells (namely, KA and TE) are used in the inversion and validation of results. Inversion results are used to build a geological model in the form of acoustic impedance from which we derive 3D porosity cube which are used for hydrocarbon potential in the Paleocene to lower Oligocene sands, and the Precambrian basement. Although the amplitude maps provide an indication of potential reservoirs, the extent of these zones are much better identified in the inverted impedance maps and the corresponding estimated high-porosity zones. The analysis predicted the potential reservoir rocks in the Sylhet, Kopili and Barail formations, in which the Sylhet and Kopili appear to have good potential zones. Near the vicinity of the Naga thrust belt, the proximity of potential reservoir is predicted in the Kopili, Sylhet formation and in the fractured basement, respectively.


Introduction
Derivation of elastic parameters such as acoustic and shear impedance and density is the primary task of seismic inversion (e.g. Sen 2006). It involves iterative updating of elastic properties while reducing the misfit between the observed and synthetic data assuming that the source wavelet is known (normally derived from the seismic trace) (Latimer et al. 2000). Synthetic angle stack traces are computed by convolving the angle-dependent reflection coefficients with the source wavelet. Seismic inversion is also known as iterative forward modelling because it relies on an initial model and perturbs the model parameters until it generates output with the minimum misfit between the original and synthetic one (Singha et al. 2019;Veeken and Silva 2004). However, different variants of seismic inversion exist. These include band-limited impedance inversion, maximum likelihood inversion, sparse spike inversion, Bayesian regularization, model based inversion, etc. (Berteussen and Ursin 1983;Hampson and Russell 1985;Russell 1988;Sacchi and Ulrych 1996;Ferguson and Margrave 1996). Compared to other methods, model-based seismic inversion is advantageous because of its higher correlation and the capability to generate the subsurface images beyond the seismic band limit and within the seismic band with the estimation of low-frequency model (Pendrel 2006;Russell and Hampson 1991;Lorenzen 2000;Latimer et al. 2000). In addition to deriving elastic parameter values, seismic inversion generally results in a higher-resolution image of the subsurface since the wavelet effect is deconvolved out of the seismic data (Das et al. 2017;Kumar et al. 2016). Acoustic impedances (it is the product of density and velocity -Russell 1988;Morozov and Ma 2009;Othman et al. 2015) are closely related to rock properties such as porosity; therefore, seismic inversion results can be used to obtain lateral information on petrophysical properties (Latimer et al. 2000;Dolberg et al. 2000;Helgesen et al. 2000;Marion and Jizba 1997). Prospect evaluation and determination of the quality of a hydrocarbon reservoir require estimation of petrophysical properties of a reservoir (Pandey et al. 2020;Chatterjee et al. 2016;Das and Chatterjee 2018). This is generally achieved by first establishing a relationship between the petrophysical properties and inverted seismic attributes at the well locations-this relationship is later used to estimate petrophysical properties in the entire volume (Rijks and Jauffred 1991;Lefeuvre et al. 1995;Russell 2004;Sancevero et al. 2005;Soubotcheva 2006). This commonly used method known as 'seismic reservoir characterization' or 'quantitative seismic interpretation' requires careful implementation of each step in areas of complex subsurface geology such as the upper Assam shelf, India-the area of our study.
The prime objective of this paper is to evaluate the hydrocarbon potential in the upper Assam basin by accurate estimation of acoustic impedance and porosity for Sylhet, Kopili and Barail formations and the granite basement in the study area (Zamanek et al. 1970;Hussain et al. 2017a,b). The shelf part of the upper Assam basin which is primarily in southeast (SE) dipping shelf (Gogoi and Chatterjee 2019;Kent and Dasgupta 2004) has spread over on the Brahmaputra and Dhansiri valley, and the slope part lies in Mikir hills and Naga thrust zone. The subsurface is structurally complex characterized by numerous faults. In this work, the 3D poststack seismic volume in the south upper Assam shelf has been interpreted for better assessments for the underlying low acoustic impedance sand in the formations of Paleocene to lower Oligocene age, and the Precambrian basement rock as well as well data. The well logs recorded parameters such as sonic velocity and density log, from each used wells (KA and TE) which are used in the inversion and to validate the results. Both the wells are drilled to the basement depth and placed within the seismic volume. The wells are addressed the depth of the formations and used to extract the wavelet from each of them. Both the wavelets are used in the seismic-tie correlation. To understand rock properties away from the well location, the lateral and vertical variation can be inferred via the quantitative interpretation of the seismic data. Here, we combine the seismic inversion analysis with well data to successfully resolve the problem of prospect evaluation. A workflow is established here to infer the potential zones of reservoir rock in the basement and the sedimentary formations of Paleocene to lower Oligocene ages. Within these formations, low-impedance and highporosity properties are used as diagnostic of the presence of hydrocarbon.

Geological background of the study area
The Upper Assam shelf is an area of Upper Assam basin located in north-eastern India, and it is a self-slope basin shown in Fig. 1a, b (Bhandari et al. 1973;Kent and Dasgupta 2004). The basement rock is granite, and the sedimentation of ~7-km-thick sedimentary rocks has been reported in the shelf part of tertiary period (Kent and Dasgupta 2004). The basin evolved through the process of rifting and drifting of the Indian plate from super continent Gondwana. As mentioned in erstwhile studies, during the late Permian to early Mesozoic, a large number of structural features, such as Horst and Graben, were created on the granitic crust in the southern slope of the Dhansiri valley which is located in the south Assam Shelf Evans 1964). South Assam Shelf is separated from the North Assam shelf by the east-west trending Jorhat fault. It has been also observed that most of the oil field in the Assam Shelf lies along the border of the Naga thrust, and major producing fields are basically structural and strati-structural (Bastia et al. 1993;Murty 1984;Ranga 1984;Singh et al. 2010). The block of the study area is placed nearby area of the Naga thrust fault. In the erstwhile research, the discovery of oil production has been reported and tested within the fractures of the basement, Basal sandstone and sand within Kopili and Sylhet formation; the steepening slope of the basement is identified towards SE direction (Saha et al. 2006;Pahari et al. 2008). The deposition of the sediments overlying transgressive Tura sandstone, Sylhet formation, and Kopili formation is tilted over on the fractured basement and created a favourable structural entrapment condition for hydrocarbons in the above-mentioned formations as shown in Fig. 2 (Singh et al. 2010;Ranga 1984). The block of the seismic survey area is located within Latitude 26° 14 ′ N to 26° 22 ′ N and Longitude 93° 56 ′ E to 94° 2 ′ E, and bounded by Mikir hills in the west and north-east (NE)-south-west (SW) trending Naga Schuppen belt which makes its eastern boundary (Fig. 1b).
The study deciphers that the structural evolution in the upper Assam basin exerted a major influence on hydrocarbon occurrence. Most likely, the source characteristic and trapping condition are identified in the Tura sandstone of lower Paleocene age, Sylhet formation of late Paleocene to Eocene age and Kopili formations of late Eocene age, whereas Barail formations of Oligocene age are found to be highly organic-rich among all other sequences (Fig. 2). In the Dhansiri valley, the effective source rock characteristics are identified in the shallow marine thick carbonaceous shale of Kopili formation of late Eocene age .  (Wandrey 2004). The study area includes the presence of gas and oil production within the basement, Sylhet, Kopili and Barail formations, respectively, which are considered for our study These studies demonstrated the presence of reservoir rocks not only in the fractured basement, but also in the overlying transgressive Sylhet limestone-Kopili interbedded sandstone of Eocene age and Tura sandstones of Paleocene age in the upper Assam shelf (Mandal and Dasgupta 2013; Gogoi and Chatterjee 2019).
The processed crosslines 205 and 214 seismic section are shown Fig. 3a, b as examples of subsurface structural complexity. These demonstrate the basement steepening slope direction to be NE-SW; the overlying deposition of Tura sandstone, Sylhet limestone, and basal part of the Kopili formation are indicated as anticlinal features as shown in Fig. 3a, b.

Data availability
In Fig. 1b, an 81 sq. km block is shown in the south upper Assam shelf in which the 3D seismic survey covered 34 sq. km which is used in this work (Fig. 1c). The study area is located in the east of Mikir hills and the west of Naga Schuppen belt. The general strike direction of the study area is in NE to SW direction. In the paper, processed 3D poststack seismic data are used for the geological interpretation, which include 650-inlines and 621-crosslines in NE-SW and NW-SE direction, respectively, as shown in Fig. 1d. The orientation of the inline and crossline is 29.986 and 119.988 degree from north, respectively (Fig. 1d). The 3D seismic survey of sampling rate 2 ms are recorded for the time length of 5 s. The frequency content of the seismic data ranges from ~8 to 70 Hz. The study included the information of transit time from sonic log and density from density log from two wells KA and TE from within the 3D seismic survey area. The location of well KA is within the seismic inline-285 and crossline-358, and drilled to the depth of 2338.3 m at the basement. The basement rock of the well KA is granitic rock overlying the deposition of sediments of Sylhet formation, Kopili formation, Barail formation, Bokabil formation, Tipam formation, respectively, and the well TE is located within the seismic inline-353 and crossline-267, and drilled to the depth of 2690.4 m, and contained similar stratigraphic succession with varying thickness as discussed in well KA.

Methodology
The workflow for inversion and rock property estimation is given below: • Generation of 3D acoustic impedance (AI) model using model-based seismic inversion. • Estimation of porosity volume for the formations using the extracted acoustic impedance model. • Mapping of the top of Barail, Kopili, Sylhet and Basement to study the steepening slope of the formations using the inverted volume. • Mapping of unconformable Isochronopach trend to study the deposition of sediments. • Mapping of low acoustic impedance model, porosity model for the assessment of potential reservoir zone. Thus, the basic workflow comprises two distinct steps: acoustic impedance inversion and porosity mapping. We describe below each of these steps, in brief, for completeness.

Seismic inversion
The post-stack seismic data are modelled by convolution of a wavelet with a series of reflection events: where S(t) is the seismic trace amplitude as a function of two way reflection time t, W(t) is the wavelet, r(t) is the reflectivity (normal incidence reflection coefficient), and n(t) is the noise.
Acoustic impedance is defined as the product of density and velocity given by equation (Lindseth 1979): where Z i is the acoustic impedance of ith layer, i is the density of ith layer, and v i is the velocity of ith layer.
The reflectivity at the boundary between two homogeneous layers is given by (Lindseth 1979): where r i is the reflectivity of ith layer, Z(i + 1) is the acoustic impedance of (i + 1)th layer (Eq. 2), and Z(i) is the acoustic impedance of ith layer (Eq. 2).
Given an acoustic impedance model, the reflection coefficient for each layer boundary is computed from Eq. 3. This results in a reflection coefficient time series, which is convolved with the given wavelet to generate a synthetic seismic trace. The model-based seismic inversion iteratively updates the impedance model while improving the fitness between the observed and the synthetic seismograms. At the convergence, the final impedance model is considered to be the subsurface impedance model at the trace location. The inversion is capable of generating a model of beyond and within the limit of seismic data by the addition of well log data (Russell and Hampson 1991;Austin et al. 2018). The suitable parameters are estimated in order to obtain minimum residual between synthetic inverted data and input seismic data by least square methods at well site, and those parameters are applied on the whole initial acoustic impedance model with the iteration based on generalized linear inversion (Cooke and Schneider 1983). (1)

Mapping of density-porosity model in 3D seismic data
The main steps used to estimate the porosity are as follows: 1. In Eq. 2, it is required to obtain an accurate velocity model to get a density cube from the impedance. The general relationship is used to estimate the velocity model from the acoustic impedance model using Lindseth's approach is shown (Atat et al. 2020): where a = 0.272926, and b = 0.821244. 2. Using the velocity model, the density model estimated from the acoustic impedance is calculated from Eq. 2. 3. The estimated density model is then used to estimate the porosity model. The general equation used to estimate the porosity from the density model is represented (Serra 1984): where is the porosity, d m is the density of matrix, and taken the value ~2.71 g/cc, d is the density of the density model, d f is the density of fluid, and taken value ~1.01 g/cc.

Model-based seismic inversion of 3D seismic volume
Derivation of acoustic impedance from seismic data by seismic inversion requires a wavelet and a low-frequency starting model. The steps for achieving these are described below: Step 1 Seismic inversion as used here and described above requires that the wavelet be known (Cui and Margrave 2014;Edgar and Baan 2011;Yi et al. 2013). The wavelet is extracted using both the wells KA and TE keeping the phase constant and using the window length of 150 ms, and taper length of 20 ms. Given these parameters, an optimal wavelet is chosen for which the correlation between synthetic and real seismic trace is the highest. The first step, however, is to convert the well logs to time domain. This is achieved by using the following equation (Karim et al. 2016): where t(z) is the two-way travel time at depth z , t(0) is the vertical two-way travel time at z(0) , and V(z) is the two-way p-wave velocity, taken from log data.
Step 2 The wavelet extracted from the well KA that gives a high correlation value between the synthetic trace and the (4) log(velocity) = a + b * log (impedance) original trace is shown in Fig. 4. The variation of amplitude of the wavelet with time and frequency is shown in Fig. 4a and b. In Fig. 4b, the red line is showing the average value of the minimum phase-89°-for each of the frequency component.
Step 3 Figure 5a, b displays the cross-correlation between the synthetic seismogram and the original seismic traces. The traces of two different colours of blue and red are shown on the right panel of the image; the blue traces represent the synthetic seismogram, and the red traces represent the original seismogram; the correlation value obtained between these two traces at each well KA and TE locations are ~0.732 and ~0.754, respectively (Fig. 5a, b).
Step 4 The depositional environments separated by imaginary lines are marked from the information given in well log data (Faraklioti and Petrou 2004). The depth of the formations Barail, Kopili, Sylhet and Basement is extracted from both the wells KA and TE, and interpolated between the wells within the 3D seismic volume. The illustrative 2D time slice of the top of each horizons Barail, Kopili, Sylhet and Basement is shown in Fig. 6a, b, c, d.
Step 5 The current step is used to compensate for the loss of low-frequency information from the seismic data which gives more realistic geological model (Oldenburg et al. 1983). The high cut filter of range ~10/15 Hz is used to build initial acoustic impedance along all the four horizons Barail, Kopili, Sylhet, and Basement with the p-wave velocity extracted from both the wells (Fig. 5). The acoustic impedance value extracted for each of the formations is shown in Fig. 7. Step 6 The step is used for the analysis of misfit between synthetic acoustic impedance and the original acoustic impedance computed from the well data. All three formations and the basement are considered in the analysis window for each of the wells KA and TE as shown in Fig. 8. The analysis of minimum misfit is verified for the well KA as shown in Fig. 8a, and the analysis of minimum misfit is verified for the well TE as shown in Fig. 8b. Figure 8a, b displays the analysis of model-based seismic inversion at each well. There are two panels: the left panel of the display shows three impedance curves, namely the blue curve which shows original impedance, the black curve which shows initial guess curve, and the red curve which shows the final inversion result, and the error between the real log impedances and inverted impedances is shown above the panel. At the location of well KA, the estimated error is ~1287.22 (m/s)*(g/cc), and at the location of well TE, the estimated error is ~1209.97 (m/s)*(g/cc). The second panel of the display is shown with two colours-red and black traces-the red traces are the synthetic traces calculated from the inversion results, and the black traces are the seismic traces. The correlation value between the synthetic trace and the seismic b basement, extracted along the horizons from an initial acoustic impedance model. The locations of the wells are marked with asterisks traces is mentioned above the panel which is estimated for the well KA and TE to be ~0.961 and ~0.971, respectively. The third panel of the display shows the error between the second panel traces, and the errors estimated for well KA and TE are ~0.275 and 0.238, respectively, which are closer to zero, indicating the consistency of acoustic impedance trace with the wavelet and the input seismic trace (Fig. 8). In this study, we have used two wells simultaneously instead of one well for inversion, and the results are remarkably improved in the distributions of impedance.

Results and discussion
An acoustic impedance volume has been generated by seismic inversion using the 3D post-stack seismic data, the wavelet and the low-frequency acoustic impedance. In this section, we provide an interpretation of the inverted volume in terms of subsurface geology.
In our study area, the strike direction is marked in NE-SW direction based on geological information. Our goal is to identify potential reservoir zones within formations. The extracted horizons for Barail, Kopili, Sylhet and Basement shown in Fig. 6a, b, c, d, respectively, show the steepening slope directions of the formations to be NW-SE, towards the Naga Schuppen belt. The Isochronopach between the top of  Fig. 9a. Maximum Isochronopach in range ~462 ms to ~546 ms demonstrates the presence of a plunging syncline structure around the NE and as well as, Isochronopach of range ~41 ms to ~104 ms demarcates the presence of a plunging anticline structure around the area containing the TE well in the study area (Fig. 9a). Transgressive Sylhet deposition is demarcated with the presence of plunging anticline and syncline structures in Paleocene and lower Eocene ages, to which it prompted the pressure around the study area due to local and regional tectonic activity in the Paleocene and lower Eocene age. The Kopili formation unconformities of upper Eocene age overlays the Sylhet formation. The subsidence and sedimentation is more or less continuous towards SE direction, and the maximum areal coverage of the accommodation of SE sediment is thickening towards the north direction of the study area, as shown in Fig. 9b. The maximum and minimum values of Isochronopach in the proximal region draw the presence of accommodations of sedimentation space in plunging syncline and anticline structure, to which it prompted the occurrence of tectonic events, regional and local sediment transport, climatic controls and the effect of eustasy in the Upper Eocene era (Fig. 9b). The Isochronopach of Barail formation demarcates lower Oligocene age. The Isochronopach map indicates that the variation in accommodation throughout the volume of study is in the lower Oligocene era as shown in Fig. 9c.
Next, we examine the high reflectivity zone of the formations and the areal extent of low-impedance sand and the availability of maximum porosity in each formation for the assessment of reservoir in the formations. Fig. 10 Preparation of porosity cube from density section to velocity cube which is further derived from impedance cube. IVT: impedance-tovelocity transformation, VDT: velocity-to-density transformation, DPT: density-to-porosity transformation The porosity cube has been transformed from density cube which is extracted from velocity cube, as shown in Fig. 10. The velocity has been extracted from inverted impedance data using Eq. 4 (Fig. 10).
The high reflectivity zone encircled in seismic amplitude slice is shown in Fig. 10a which indicated the presence of low-impedance sand of ranges ~2000 to ~5000 Fig. 11 Illustrative diagram is the assessment of the presence of highreflectivity seismic amplitude in the basement of the study area: a the encircled zone indicates the high reflectivity zone in the acquired seismic data; b model-based seismic inversion is used to extract the acoustic impedance model, and low impedances ~3000 to ~4667 (m/s)*(g/cc) correspond to high reflectivity zones; c the porosity model: encircled porosity ranges from ~12 to ~20% at high reflectivity zone. The study is the assessment of the presence of low-impedance sand reservoir inside the fractured basement (m/s)*(g/cc) encircled in Fig. 11b to show the presence of reservoir rock in the Sylhet formation; the porosity extracted for the similar low-impedance zone ranges from ~28 to 38% which is encircled in Fig. 11c. In Fig. 11 b, c, we identify proximal variation of reservoir rock more or less in the NE and SE directions around the location of both the wells KA and TE. Larger areal extent of the reservoir rock is found in the SE part of the study area in the Sylhet formation. In Fig. 12b, the low-impedance sand of ranges ~2000 to 4000 (m/s)*(g/cc) is encircled to show the presence of reservoir rock in the high reflective zone (Fig. 12a) within Kopili formation; the porosity extracted for the similar low-impedance zone is ~31% to ~39% which is encircled in Fig. 12c. In Fig. 13b, the zone of the reservoir low-impedance sand ranges from ~2000 to 4000 (m/s)*(g/cc) which is encircled with the high reflectivity zone (Fig. 13a), and the porosity of ranges ~37% to 41% estimated in the Barail formation is shown in Fig. 13c. Figure 14b displays the top of the basement with 10 ms window, and the encircled low impedance of range ~3000 to ~4667 (m/s)*(g/cc) within the high reflectivity zone (Fig. 14a) demarcated the presence reservoir with the porosity of ranges ~12% to 21% which is encircled in Fig. 14b. In Fig. 12 Illustrative diagram is the assessment of the presence of highreflectivity seismic amplitude in Sylhet formation: a the encircled zone indicates the high reflectivity zone in the seismic data; b modelbased seismic inversion is used to extract the acoustic impedance model, and encircled the low impedances ~2000-~5000 (m/s)*(g/cc) at high reflectivity zone; c the porosity model encircled the porosity ~28-~38% at high reflectivity zone. The study is the assessment of the presence of low-impedance sand reservoir in the Sylhet formation Fig. 14b, the nonzero value of porosity indicates the presence of fractures in the basement.
Oligocene Barail shales, which act as a source rock in the Oligocene to Miocene of the petroleum system, have higher porosity as observed in time section of the porosity of Fig. 13c. This higher porosity in the Barail formation could be associated with the higher rate of sedimentation in the Middle Miocene period in fluvio-deltaic conditions. The Eocene Kopili formation has fine-grained sand during a marine transgression and contains huge amounts of lowerimpedance sand reservoir with higher porosity but lesser than Barail formation. The lower Eocene Sylhet which is made of limestone with interbedded fine shale and sand layers has lower porosity than the Kopili formation.

Conclusions
Inversion of a 3D seismic volume in the upper Assam shelf was carried out, and the resulting impedance and porosity volumes were examined to identify potential reservoir zones. The data and the corresponding inverted volume clearly highlight the basement, and the formations of Paleocene to lower Oligocene age of the upper Assam shelf and the steepening slope of the basement. It also displays that each formation dips towards the SE direction, approaching the Naga thrust fault. The maximum and minimum values of Isochronopach in the proximal region demarcate the presence of accommodations of sedimentation in the plunging syncline and anticline structures, to which it prompted the Fig. 13 Illustrative diagram is the assessment of the presence of high-reflectivity seismic amplitude in Kopili formation: a the encircled zone is indicated the high reflectivity zone in the seismic data; b model-based seismic inversion is used to extract the acoustic impedance model, and encircled the low impedances ~2000-~4000 (m/s)*(g/cc) at high reflectivity zone; c the porosity model encircled the porosity ~31-~39% at high reflectivity zone. The study is the assessment of the presence of low-impedance sand reservoir in the Kopili formation occurrence of tectonic events, regional and local sediment transport, climatic controls and the effect of eustasy in Paleocene to lower Oligocene. Considering the thickening sediments displaying low impedance and high porosity, we identify the most prospective zones for reservoir rocks in the formations. The study identifies the potential reservoir rock in the Sylhet, Kopili and Barail formations, in which the Sylhet and Kopili show good potential zones of reservoir rock which are found to be consistent with the well log data. As the area approaches the vicinity of Naga thrust belt, the proximity of the potential reservoir is predicted in the Kopili formation, Sylhet formation and in the fractured basement. The analysis of porosity modelling reveals that the higher porosity values reach 41% in the Barail formation, and the distributions of porosity in the granite basement clearly indicate potential sand reservoirs. The fund was used for the data collection, research analysis and manpower for the preparation of the manuscript.

Declarations
Conflict of interest It is hereby informed that we do not have any conflict among authors and regarding data for possible publication the manuscript in the Journal of Petroleum Exploration and Production Technology. I would like to inform you that we have used one 3D seismic data and three well log data in the Assam upper shelf. The data Fig. 14 Illustrative diagram is the assessment of the presence of highreflectivity seismic amplitude in Barail formation: a the encircled zone indicates the high reflectivity zone in the seismic data; b modelbased seismic inversion is used to extract the acoustic impedance model, and encircled the low impedances ~2000-~4000 (m/s)*(g/cc) at high reflectivity zone; c the porosity model encircled the porosity ~37-~41% at high reflectivity zone. The study is the assessment of the presence of low-impedance sand reservoir in the Barail formation were collected from the Oil India with the help of Prof Rima Chatterjee who is one of the co-author of this manuscript. The results of the manuscript were not published to any other journal previously.

Ethical approval
We, three authors (Neha Rai, Dip Kumar Singha and Rima Chatterjee), have maintained all the ethical transparency and professional conducts during research analysis, data sharing and preparing the manuscript. We have mentioned the data source Oil India in the acknowledgement section of the manuscript.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.