Impact of pertinent parameters on foam behavior in the entrance region of porous media: mathematical modeling

Foam injection is a promising solution for control of mobility in oil and gas field exploration and development, including enhanced oil recovery, matrix-acidization treatments, contaminated-aquifer remediation and gas leakage prevention. This study presents a numerical investigation of foam behavior in a porous medium. Fractional flow method is applied to describe steady-state foam displacement in the entrance region. In this model, foam flow for the cases of excluding and including capillary pressure and for two types of gas, nitrogen (N2) and carbon dioxide (CO2) are investigated. Effects of pertinent parameters are also verified. Results indicate that the foam texture strongly governs foam flow in porous media. Required entrance region may be quite different for foam texture to accede local equilibrium, depending on the case and physical properties that are used. According to the fact that the aim of foaming of injected gas is to reduce gas mobility, results show that CO2 is a more proper injecting gas than N2. There are also some ideas presented here on improvement in foam displacement process. This study will provide an insight into future laboratory research and development of full-field foam flow in a porous medium.

Gas relative permeability in the absence of foam k rg Gas relative permeability in the presence of foam k rw Water relative permeability MRF Mobility reduction factor n D Dimensionless foam texture ( n D = n f ∕n max ) n LE D Dimensionless local-equilibrium foam texture n f Foam texture (number of lamellae per unit volume), m −3 n LE f Local-equilibrium foam texture, m −3 n max Maximum foam texture, m −3 P c Gas-water capillary pressure, Pa P *

Introduction
The worldwide economy, technology advancement, and population growth have led to a significant increase in energy demand. Hydrocarbon reservoirs are predicted to stay as the main energy resources for the upcoming decades. Oil reservoir initial pressure is adequate to induce just a small volume of oil flow to production wells at the end of primary recovery (Al Adasani and Bai 2011;Radler and Bell 2009). Therefore, large and considerable amount of hydrocarbon remains in reservoirs at the end of this stage. Enhanced oil recovery (EOR) methods are used as a usual practice in the oil industry to extract the residual hydrocarbon in porous media (Bedrikovetsky 2013;Farajzadeh et al. 2012;Lake 1989). In EOR, external forces and materials are used to handle physico-chemical interactions in the reservoirs for promoting favorable recovery conditions (Al Adasani and Bai 2011). Gas flooding has been applied widely as one of the most productive EOR methods in hydrocarbon reservoirs through various injection modes like water-alternating-gas (WAG) process (Sheng 2015). Currently, nitrogen (N 2 ), carbon dioxide (CO 2 ), air and hydrocarbon gasses are injected into the reservoir to displace remaining hydrocarbon in place Boersma and Hagoort 1994;Bruining and Marchesin 2007;Green and Willhite 1998;Orr 2007;van Batenburg et al. 2010). Gas-based EOR methods are classified into immiscible and miscible flooding (Al Adasani and Bai 2011). In the former, injected gas dissolves partly in oil while in the latter it mixes completely with oil, causing oil to swell, and a considerable reduction in the interfacial tension between oil and displacing phases. Thus, the saturation of residual oil is reduced and oil is mobilized ). However, low volumetric sweep efficiency of the injected gas is major challenge associated with the most of gas-based EOR techniques, because the gas cannot touch a huge quantity of oil and, consequently, the ultimate recovery remains poor (Rossen and van Duijn 2004;Wellington and Vinegar 1988). This occurs as the results of channeling, viscous fingering and gravity override (Chang et al. 1994;Joekar-Niasar and Hassanizadeh 2011;Koval 1963;Waggoner et al. 1992). Foam, as a mobility-control agent, may be injected to reservoirs to improve the volumetric sweep efficiency (Du et al. 2008;Farajzadeh et al. 2009;Friedmann et al. 1991;Gauglitz et al. 2002;Schramm 1994). Foam is described as a dispersion of gas throughout a continuous liquid phase (containing surfactant) where gas bubbles are segregated by thin liquid films named lamellae Farajzadeh et al. 2009;Hiraski 1989). The ability of foam to improve oil recovery is due to: (1) it can control viscous fingering and gravity override by increasing viscosity of gas phase , (2) block the high-permeable swept zones and divert the displacing fluid into lower permeability zones (Yan et al. 2006) and (3) decrease the capillary forces by reducing interfacial tension in the attendance of a surfactant (Derikvand and Riazi 2016;Panahpoori et al. 2018).
Bond and Halbrook (Boud and Holbrook 1958) were the first ones who stated that in situ foam generation in the hydrocarbon reservoirs would increase the volumetric sweep efficiency. Foam as a mixture of gas, water and a surface active material consists of lots of gas/liquid interfaces (lamellae) that divide gas bubbles (Kam and Rossen 2003;Li et al. 2014;Schramm and Wassmuth 1994;Vikingstad et al. 2006). Plateau border is where three lamellae are jointed, at an angle of 120° (Al-Attar 2011; Al-Mossawy et al. 2011). Film thickness could be in a range of nanometers while the area is usually up to a few square meters. These films become stable in the attendance of molecules of surfactant at the interface of a gas/surfactant solution (Bergeron et al. 1993;Farajzadeh et al. 2008). Based on foam quality (the ratio of gas flow rate to the summation of gas and liquid flow rates), they can be classified into "dry foams" and "wet foams" which are referred to high quality (strong foams) and low quality (weak foams), respectively (Al-Attar 2011). Foams can also be categorized in terms of their bubble size as "coarse" and "fine" for large and small bubble sizes, sequentially (Gauglitz et al. 2002). Foam as a metastable system coarsens inherently and finally crumbles completely owing to capillary suction and liquid drainage from the lamella to the plateau border. This results in the fracture of foam films and foam destruction . The strength of foam relies on the magnitude of limiting capillary pressure ( P * c ). In a favorable condition, foam would have a finer bubble size and gas would possess a lower mobility up to the rate of foam destruction which is equal to the rate of foam formation. Water comes to a lower saturation, as the gas approaches to a lower mobility. Consequently, a higher capillary pressure of gas/water will be created. As the capillary pressure reaches the limiting value, the deformation rate will be raised. The limiting capillary pressure depends on the formulation and concentration of foam, gas velocity, oil existence and permeability of the porous medium. Foam would be unstable in water saturation range of lower than the corresponding saturation at the limiting capillary pressure ( S * c ) (Khatib et al. 1988). Population-balance and local-equilibrium models are two main techniques for modeling foam displacement in porous media. Population-balance approach is utilized to simulate the steady state and transient flow of foam in porous media (Kovscek et al. 1995). Power of population-balance model to predict foam behavior and successively lowering of the gas mobility relies upon quantifying the foam bubble size variation (Kovscek et al. 1995;Kam and Rossen 2003). The population-balance provides effective gas viscosity and relative permeability in terms of flow rates, bubble size and differential equations for foam texture (n f ) estimation (Kam and Rossen 2003). In the population-balance model, the foam texture is developed explicitly, writing balance equation on lamella like to the mass balances for water or surfactant. Considered mechanisms for lamella formation and destruction are described in the balance on n f , explicitly. Also, the gas mobility is presented as a function of some factors including foam texture n f and water saturation S w . Above all, the population-balance model supplies a general framework to express all the pertinent physics of foam generation and transport (Kovscek et al. 1995).
A substitute for population-balance modeling with assumption of local equilibrium at all formation locations is a local-equilibrium approximation (Kovscek et al. 2010). Under more-restricted conditions (some additional assumptions in one-dimensional displacement of local-equilibrium model), the method of characteristics, or fractional flow theory, can be applied. Fractional flow methods are applied to investigating foam and other improved oil recovery processes (Shan and Rossen 2002). The governing hyperbolic partial differential equations derived by fractional flow theory are solved by a mathematical method known as the method of characteristics analytically. The assumptions consist (Dholkawala et al. 2007): • Incompressible phases • Newtonian mobilities • One-dimensional displacement • No dispersion, viscous fingering and capillary pressure gradients • Quick achievement of local steady state Despite the fractional flow method requires simplifying the above-mentioned hypothesizes, even when the hypothesizes are not satisfied completely, this method is still capable of describing complicated mechanisms involved in foam flow in porous media (Dholkawala et al. 2007).
There are several experimental (Kumar and Mandal 2017;Lv et al. 2017;Panahpoori et al. 2018) and theoretical studies (Du et al. 2017;Gassara et al. 2017;Shojaei et al. 2018) to investigate the foam flow behavior in porous media. As the foam texture determines the strength and mobility of foam, theoretical models have been developed to examine and modify gas mobility in the presence of foam. The prediction of foam performance was investigated using population-balance models (Du et al. 2017;Falls et al. 1988;Kovscek et al. 1995;Schramm 1994;Kovscek et al. 1997;Patzek 1986) and fractional flow models (Gassara et al. 2017;Rossen and Zhou 1995;Zhou and Rossen 1994) as the principal foam-modeling methods. Nevertheless, other numerical works have been accomplished such as using digital rock physics (DRP) technology to examine foam flow behavior in a complex pore geometry (Yang et al. 2016;Shojaei et al. 2018;Du et al. 2014). Dholkawala et al. (2007) suggested that increasing the total injection velocity would change the form of fractional flow curves, which successively would affect volumetric sweep efficiency and pressure distribution. Ashoori et al. (2010) investigated water saturation and foam texture in the entrance region (the area near the well, where gas, water and surfactant make foam) of porous media. Minssieux (1974) reported a considerable water saturation entrance region about 10 cm. Also, Myers and Radke (2000) observed an entrance region as long as 15 cm in their examinations.
Although various studies have been accomplished on the topic of foam displacement in porous media, less attention has been paid to investigation of the impact of pertinent parameters such as water fractional flow, residual gas saturation, foam constant and gas viscosity on foam behavior. In the present study, first-order kinetic (FOK) foam model (Ashoori et al. 2010(Ashoori et al. , 2011 is validated and implemented to describe steady-state behavior of foam flow in the entrance region. Then, foam flow for two different cases (i.e., excluding and including capillary pressure) and for two types of gases (i.e., N 2 , CO 2 ) are investigated and compared for the entrance region. All pertinent parameters that affect foam texture to reach local equilibrium have been also examined in this study.

Mathematical modeling
The following assumptions were considered in this work: absence of mobile oil in the entrance region, uniform surfactant concentration, steady-state incompressible displacement of an immiscible gas-water rectilinear flow through a porous medium. Water mass balance equation, considering the above assumptions, is governed by Rapoport-Leas equation (Barenblatt et al. 1991) as follows: Foam texture is expressed by the population-balance equation as bellow (Barenblatt et al. 1991): n D is the dimensionless foam texture and can be expressed as n D = n f ∕n max .
It should be stated that in foam processes, flow properties depend on water saturation and bubble size or foam texture ( n f ) which is introduced as the number of lamellae per unit volume of gas phase. Hence, mathematical modeling of foam processes needs, in addition to the water mass balance, a numerical representation for foam texture, named population-balance equation (Eq. 2). As previously mentioned, in this model, foam texture is developed explicitly through writing balance equation on lamella. Important mechanisms for lamella formation and destruction are included in the balance on n f . In present work, for solving Eqs. (1) and (2), first-order kinetic (FOK) foam model (Ashoori et al. 2010) was applied. This model describes foam flow behavior in porous media at steady-state condition.

Foam in the entrance region
The fractional flow method is applied as follows: the entrance region at location x = 0 is indicated by E. There is no bubble texture at this point. As foam is generated, x approaches toward the end point of the entrance region (J) in which the bubble texture reaches to local equilibrium (LE) with S w = S J w and r g − r c = 0 . Following boundary conditions at locations of E and J are defined: In the present study, the effect of various parameters on length of the entrance region (from point E to state J) will be determined and discussed in detail. By simplifying Eqs. (1) and (2) along with some required correlations listed in Table 1, following results are achieved for foam displacement in the entrance region at the presence and the absence of capillary pressure.

Case 1: absence of capillary pressure
Water mass balance excluding capillary pressure is obtained as below (Barenblatt et al. 1991): The foam texture is obtained based on water fractional flow: The boundary condition is:

Case 2: presence of capillary pressure
Water saturation including capillary pressure and foam texture is obtained using Eq. (8)-(10) (Barenblatt et al. 1991): The boundary conditions are: Water fractional flow Gas mobility

Auxiliary correlations
Auxiliary correlations should be added to these equations to complement the modeling. In this study, the correlations of local-equilibrium texture, foam model, foam fractional flow, capillary pressure function, gas mobility reduction as well as foam-free relative permeabilities should be considered. The equations are listed in Table 2.

Numerical method
Governing equations consist of a set of ordinary differential equations along with auxiliary correlations. Algebraic equations in the model incorporate aforementioned correlations listed in Table 2. These equations create a set of nonlinear algebraic equations. The equations are solved in MATLAB programming environment using Gauss-Newton approach.

Results and discussion
Present study verifies foam injection into a core (linear displacement) at the presence of gas and surfactant and at constant temperature. The FOK model was applied to describe steadystate behavior of foam flow in porous media. In this model, foam displacement in two different cases (i.e., excluding and including capillary pressure) and for two types of gas (i.e., N 2 , CO 2 ) were investigated and compared for the entrance region.

Model validation
To examine the capability and accuracy of the model used in this study, the reported experimental data by Persoff et al. in 1991(Persoff et al. 1989) and Ma et al. in 2013(Ma et al. 2013 were modeled. In the study done by Persoff et al. in 1991, the steady-state water saturation was measured over a range of injected water fractional flow from 0.008 to 0.20 in a sandstone sample with porosity and permeability of 0.25 and 1 Darcy, respectively (Persoff et al. 1989). Ma et al. in 2013 measured foam apparent viscosity, as a representation of mobility, at a fixed superficial velocity for injected water Local-equilibrium texture  Zhou and Rossen (1994) Gas relative permeability in the absence of foam Zhou and Rossen (1994) Water fractional flow function Ashoori et al. (2010) Gas-water capillary pressure P c = 15000 × It should be noted that in the modeling of the foam experimental data, generally some of the model parameters should be adjusted (Rossen and Boeije 2015). Hence, S * w , A and the coefficient of n D in the mobility reduction factor correlation (Table 2) were tuned to 0.362, 70 and 54,330 for the data of Persoff et al. in 1991(Persoff et al. 1989) and 0.098, 207 and 60,000 for the data of Ma et al. (2013), respectively, using pattern search algorithm (Abramson 2002;Vasant 2013).

Case 1: absence of capillary pressure
The values of the required parameters for the FOK model which were considered for this study are presented in Table 3. In the absence of capillary pressure, Eqs. (4)-(6) are applied to predict foam behavior in the entrance region. Before discussing the results, it seems to be necessary to explain foam generation mechanisms.
Foam generation is responsible for the formation of liquid films (lamellas), that segregates gas bubbles (Risal et al. 2019). In situ lamella generation is as the result of different mechanisms, such as snap-off, leave-behind and lamella division. Snap-off takes place when there is a capillary pressure difference between the pore bodies and pore throats. Leave-behind mechanism happens during injection of a gas phase to a continuous liquid phase in a porous medium (Kovscek and Radke 1993;Ransohoff and Radke 1986). Lamella division as the primary mechanism of foam generation in steady gas-liquid flow (Gauglitz et al. 2002;Qichong Li 2006) increases the number of lamella or bubbles by splitting a moving lamella into two, one in each   Fig. 3 Water saturation and dimensionless foam texture in the entrance region for K c = 1 s −1 branch of the paths when confronts a branch in the flow path (Dicksen et al. 2002). Therefore, in the absence of capillary pressure, foam generation is dominantly caused by lamella division mechanisms. Water saturation and dimensionless foam texture are depicted in Fig. 3 when the foam constant is equal to one (i.e. K c = 1 s −1 ). The entrance region begins with no bubble texture and then with foam generation, the dimensionless foam texture increases to approach the local equilibrium (LE), state J. Since at the core inlet, no foam is presented, the nitrogen gas mobility remains unchanged, therefore water saturation has a higher value than that in the rest of the core. If the bubble sizes are equal, with increasing the number of foam bubbles, the gas saturation increases and consequently, the water saturation reduces. As can be seen, the foam texture is a strong function of water saturation because a small change in S w leads to a significant augmentation of n D (in the range of 20-40 cm length). The entrance region in this case is of order of cm (43 cm).

The effect of water fractional flow
In the absence of capillary pressure, using Eqs. (1) and (2) along with the required correlations presented in Table 1, the following relationship is achieved: The effects of injected water fractional flow ( f J w ) on foam texture, water saturation and gas relative permeability are shown in Fig. 4. On an adequately large scale, where the entrance region could be neglected, the fractional flow at the end point J of the entrance region indicates the so-called injection condition for the displacement. As can been seen from Fig. 4a, an increment of water fractional flow increases the foam texture. According to the definition of water fractional flow (Eq. 12), the gas mobility reduces when the fractional flow increases. This means that the number of foam bubbles increases. On the other hand, at a constant gas viscosity (despite the fact that gas relative permeability and gas viscosity are both responsible for dynamic behavior of foam displacement, in the present model, reduction in gas relative permeability is assumed as the only affecting mechanism for simplicity, and the impact of viscosity is not taken into account) the gas relative permeability reduces as the gas mobility decreases (see Eq. 13), consequently the gas saturation decays and as a result the water saturation increases. This behavior can be seen in Fig. 4b and c.  Figure 5 demonstrates the foam texture and water saturation as a function of residual gas saturation. The residual gas saturation was considered as a representation of the injected (13) g = k rg k 0 g gas-water ratio to investigate the effect of this parameter on foam behavior in the entrance region. As clearly can be seen, the length of the entrance region increases as the residual gas saturation reduces. It means that gas and water have more time to direct contact, hence more foam texture with smaller size is generated.

The effect of foam constant in the FOK model
The influence of foam constant (which appears in Table 2 in the foam model) on foam texture is investigated as shown in Fig. 6 and Table 4. The length of the entrance region reduces as the foam constant increases. It means that at higher values of foam constant, the foam is formed faster. Foam constant does not affect the final foam texture value and only influences the length of the entrance region. For faster kinetic rates (e.g., K c = 1 s −1 ) the entrance region has an extremely short length. According to Table 4, by decreasing K c one order of magnitude, the length of the entrance region increases about ten times. Therefore, the length of the entrance region is inversely proportional to K c .

The effect of gas viscosity (gas type)
Figure 7 compares the foam texture and mobility of nitrogen and carbon dioxide gases. Nitrogen has a lower viscosity than carbon dioxide (its viscosity is almost onethird of that of carbon dioxide), therefore more distance will travel in order to reach equilibrium. According to Eq. (6), the foam texture is an inverse function of gas viscosity. Since the carbon dioxide viscosity is almost three times greater than the nitrogen viscosity, the final foam texture of carbon dioxide is almost one-third of nitrogen (see Fig. 7a). As demonstrated in Fig. 7b),the carbon dioxide mobility is lower than that of nitrogen. It shows that the effect of viscosity reduction for nitrogen is more than the influence of decreasing gas relative permeability caused by foam generation (see Eq. 13). Since the aim of foaming of the injected gas is reduction in gas mobility, from these figures, it can be understood that carbon dioxide is a more proper displacing gas than nitrogen. It should also be noted that although the nitrogen relative permeability is lower than that of carbon dioxide, more foam texture in the presence of nitrogen covers the main disadvantages of nitrogen (low viscosity is considered as the main disadvantage in foaming processes). Consequently, the nitrogen mobility approaches that of carbon dioxide.  Fig. 8 The effect of foam constant on a relative gas permeability and b foam texture in the entrance region 4.2.5 The effect of foam constant on gas relative permeability and foam texture Figure 8a shows gas relative permeability as a function of foam constant and the length of the entrance region. Gas relative permeability reduces as the foam is generated in the entrance region and once it reached an equilibrium state, it would follow a constant value. Foam constant alters from 0.01 to 100. With a reduction in foam constant, the length of the entrance region increases, and more time is required for foam to be generated. Figure 8b also confirms that with a augmentation of foam constant, the equilibrium is achieved faster.

Case 2: Presence of capillary pressure
To predict foam flow behavior in the entrance area including capillary pressure, Eqs. (8) and (9) are used. In this case, in addition to lamella division, leave-behind and snap-off mechanisms as by cause of the difference in the capillary pressure, also play a role in foam generation. In Fig. 9a, the capillary pressure as a function of water saturation is illustrated. Foam stability in porous media strongly relies on capillary pressure. In the favorable condition of foam generation, the foam would have less mobility with finer texture until the foam destruction rate becomes equal to the foam generation rate. Consequently, lower water saturation will be attained. As the water saturation is reduced, the capillary pressure and consequently the gas fractional flow (f g ) increase. As the capillary pressure reaches around the limiting capillary pressure ( P * c ), the coalescence rate will be raised and in the case of exceeding the limiting value, unstable foam appears. Below the corresponding liquid saturation ( S * w ) corresponds to P * c , foam is unstable. At the conditions of high stability of lamella during a foam displacement, the magnitude of capillary pressure cannot exceed a limiting value, P * c , then the water saturation must be greater than or equal to the corresponding liquid saturation, S * w owing to the capillary pressure is a monotonic function of saturation. Figure 9b shows the gas fractional flow as a function of water saturation and foam texture. By increasing the gas fractional flow, the water saturation declines until to approach the limiting value S * w (at P * c ). In the limiting capillary pressure, foam bubble sizes are changed as needed to prevent the capillary pressure rising above the limiting value. In this region, the water saturation is constant and equal to S * w . The capillary pressure stays at its limiting value while the foam texture becomes coarser, and a strong foam converts slightly to a weak foam, with more increment in gas fractional flow. Therefore, in porous media, the coalescence of all lamellas do not take place at once. Figure 10 shows the water saturation and the dimensionless foam texture in the case of considering capillary pressure. As illustrated, the length of the entrance region is of the order of mm (about 20 mm). Figure 11 compares the foam textures for two different cases including and excluding capillary pressure for K c = 100 s −1 . As discussed already, at the presence of P c , snap-off and leave-behind mechanisms are also responsible in the foam generation, therefore foam is formed faster in a smaller length of the entrance region compared to the case of absence of capillary pressure in which only the lamella division mechanism causes foam generation. The length of the entrance region for the case of including capillary pressure is almost one order of magnitude less than the case of excluding capillary pressure. It should be stated that when leave-behind, snap-off, and lamellae division mechanisms are responsible for foam generation simultaneously (at the presence of P c ), they could affect each other according to the system conditions and specifications such as pore throat geometry, aspect ratio distribution, pore-network topology, flow rate, wettability of rock surface, etc. (Shojaei et al. 2018). Generally, one of these mechanisms controls the foam formation rate. Nevertheless, a pore level investigation is required to understand which one is dominant. Model presented in this study is not able to determine the dominant mechanism due to its macro-scale feature. However, the existence of three mechanisms and their effects on foam formation rate can be observed from the outcomes presented in this section. Figure 12a shows that although in the entrance region, water fractional flow for carbon dioxide is higher than that for nitrogen, it finally reaches a lower value. Water fractional  Figure 12b also compares water fractional flow in terms of water saturation for nitrogen and carbon dioxide gases.

Conclusions
This study presents the influence of pertinent parameters on foam behavior in the entrance region of porous media. Fractional flow method was implemented for analyzing foam behavior in the entrance area. In this model, the ability of nitrogen and carbon dioxide gases in reduction in gas mobility were also investigated and compared. The results show that the required length of the entrance region for foam texture to approach local equilibrium is a strong function of water fractional flow, residual gas saturation, foam constant and gas viscosity as the pertinent properties in foaming the injected gases in coreflooding experiments and also scale-up process for field applications. The key findings of this study are summarized as below: • The length of the entrance region increases as the residual gas saturation reduces and more foam texture is generated. • The length of the entrance region reduces as the foam constant increases; consequently, the foam is formed faster. • Foam constant does not affect the final foam texture value. However, it influences the length of the entrance region. • Carbon dioxide is a more proper gas than nitrogen due to more reduction in its gas mobility. However, more foam texture in the presence of nitrogen nearly could cover the main disadvantages of nitrogen (low viscosity). • In the presence of capillary pressure, foam is formed faster in a smaller length of the entrance region, because in addition to lamella division, leave-behind and snap-off mechanisms are also responsible for foam generation. • Capillary pressure has a significant role in the stability of foam in porous medium. Hence, neglecting this important factor, as in many simulation and modeling studies, could eventuate the inaccurate foam behavior in the entrance area. Nevertheless, several variables such as surfactant type and concentration, gas velocity and permeability of formation impress on limiting capillary pressure ( P * c ), which should be screened and examined in different cases in future works.
The present study proves that foam can be a successful technology for improving oil recovery, when there is an exhaustively understanding of all the possible factors that could influence on foam behavior in porous media.
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/.