Benchmark study for the simulation of Underground Hydrogen Storage operations

While the share of renewable energy sources increased within the last years with an ongoing upward trend, the energy sector is facing the problem of storing large amounts of electrical energy properly. To compensate daily and seasonal fluctuations, a sufficient storage system has to be developed. The storage of hydrogen in the subsurface, referred to as Underground Hydrogen Storage (UHS), shows potential to be a solution for this problem. Hydrogen, produced from excess energy via electrolysis, is injected into a subsurface reservoir and withdrawn when required. As hydrogen possesses unique thermodynamic properties, many commonly used correlations can not be simply transferred to a system with a high hydrogen content. Mixing processes with the present fluids are essential to be understood to achieve high storage efficiencies. Additionally, in the past, microbial activity, e.g. by methanogenic archaea, was observed, leading to a changing fluid composition over time. To evaluate the capability of reservoir simulators to cover these processes, the present study establishes a benchmark scenario of an exemplary underground hydrogen storage scenario. The benchmark comprises of a generic sandstone gas reservoir and a typical gas storage schedule is defined. Based on this benchmark, the present study assesses the capabilities of the commercial simulator Schlumberger ECLIPSE and the open-source simulator DuMux to mimic UHS related processes such as hydrodynamics but also microbial activity. While ECLIPSE offers a reasonable mix of user-friendliness and computation time, DuMux allows for a better adjustment of correlations and the implementation of biochemical reactions. The corresponding input data (ECLIPSE format) and relevant results are provided in a repository to allow this simulation study’s reproduction and extension.

While natural gas storage with a history of more than 100 years can be stated as well examined, hydrogen storage is less assessed, and its unique properties require further investigations. As an example of hydrogen's unique properties, the combination of a lower energy density in terms of volumes and worse compressibility than natural gas leads to a reduced storage capacity with hydrogen. Furthermore, thermodynamic properties such as density and viscosity lead to differences in fluid dynamics and may strengthen phenomena such as viscous fingering and density segregation. However, also, hydrogen-induced chemical reactions were proven in the last years [1][2][3]. Biochemical reactions are governed by microorganisms consuming the stored hydrogen. Additionally, geochemical reactions related to the surrounding rocks may lead to a loss in hydrogen and even lead to impurities, e.g. harmful gases such as hydrogen sulfide. With a focus on UHS operation, it is relevant to predict and model these processes properly. The present paper focuses on this modeling part and proposes a benchmark scenario to allow for a comparison of reservoir simulators to mimic these unique properties.

Fluid dynamics in porous UHS
Hydrogen is in supercritical state in the range of typical subsurface storage conditions (pressure and temperature). However, its fluid dynamic behavior is expected to be gaslike. Considering the initially present brine phase, UHS becomes a complex two-phase system and capillary effects and the relative permeability behavior control the flow. Different zones in the storage can contain gases with different composition so that gas mixing by dispersion plays a significant role. Mass exchange between the phases occurs due to dissolution and vaporization.
Due to the small molecular mass and size of hydrogen, its properties deviate strongly from other gases. The compressibility of hydrogen is almost ideal so that its density is predicted almost accurately by using the ideal gas law. For mixtures between hydrogen and other gases typically a correction by using a cubic equation of state is required. Under the initial storage conditions in this study (T = 60 • C and P = 82 bar) the density of pure hydrogen is around 6.3 kg m 3 . The viscosity (9.4 μPa · s under the initial storage conditions) is also smaller than for natural gas. Both, the low density and the low viscosity, can lead to an unstable displacement of brine and consequently gravity overriding or viscous fingering might occur during the development of an UHS in an aquifer or at the gas-water contact in an UHS [4,5].
Mixing between hydrogen and other gases in the subsurface becomes relevant when a former natural gas storage or a depleted gas reservoir is transformed into an UHS or when an alternative cushion gas (e.g. N 2 or CO 2 [4,6,7]) is used. In these cases, it is obvious that impure hydrogen could be produced during withdrawal periods. The intensity of mixing is influenced by mobility ratios, density differences, molecular diffusion and mechanical dispersion [4]. Similar to the interactions between hydrogen and brine, the displacement of other gases by hydrogen might become unstable due to high density differences and unfavorable mobility ratios. Hence, gravity overriding might occur during the injection of hydrogen. However, the miscibility between different gases due to a high dispersion at the displacement front acts as a stabilizing force so that this effect is expected to be much less pronounced [8].
Molecular diffusion is generally a very slow process but the molecular diffusion coefficient of gaseous hydrogen is relatively high (in the order of 1 · 10 −6 m 2 s ). As molecular diffusion is independent of any flow, it results in a continuous mixing also during idle periods. In contrast, during injection and withdrawal periods a more intense mixing is expected due to mechanical dispersion.
Mechanical dispersion results from variations in the flow velocity, which can occur on different scales [4]. Assuming flow velocities of several meters per day, which are common in gas storages, the longitudinal mechanical dispersion coefficient is estimated to be around 5·10 −4 m 2 s and thus 2 orders of magnitude higher than the molecular diffusion [9].
Studies about the hydrogen solubility in water and brine by Chabab et al. [10] have shown that up to 80 moles of hydrogen dissolve in 1 m 3 of water (at 100 bar and 50 • C). The solubility increases with increasing pressure and decreases with increasing salt content. Consequently, once-only losses of hydrogen by dissolution into the connate water and diffusion into over-and underlaying water-saturated layers have to be considered when hydrogen is injected for the first time [11].

Microbial activity
While the surface offers favorable conditions for living organisms, the subsurface can be considered as the absolute opposite. Factors such as high temperature, low to no oxygen concentration, pressure or unfavorable pH-values lead to a harsh domain for those organisms. However, some microorganisms have adapted themselves to these conditions. With a size of a few micrometers (<10 μm), archaea and bacteria consume resources present in the subsurface in order to maintain their metabolism. This metabolism leads to the process of cell division, where one cell divides into two independent ones [12]. Overall, this process leads to a dynamic size of populations where the number of bacteria or archaea is mainly driven by the availability of substrates, but also space and other environmental influences have an impact.
In UHS, various species of bacteria and archaea may be present in the subsurface. These types of organisms highly vary in their characteristics: while some are adapted to higher temperatures, others may be resistant to high pHvalues. The bacteria present in the subsurface have to be evaluated for each reservoir separately by performing laboratory experiments [13]. Nevertheless, the microorganisms can be classified by the source and products of the metabolism. During the storage of highly hydrogen containing town gas the conversion of hydrogen into methane was observed [1]. Later on it was accounted for methanogenic archaea, which consume carbon dioxide and hydrogen and produce methane and water [14]. The metabolism takes place according to the Sabatier-reaction [15]: These methanogenic archaea will lead to the beneficial side effect of converting carbon dioxide into methane and play a crucial role in the development of a so-called Underground Methanation Reactor. During first field tests within the project of Underground Sun Storage [2] by RAG Austria, the microbial growth of methnaogenic microorganisms was stated to be sufficient in the formation so that its selective utilization was investigated in the follow-up project Underground Sun Conversion [3]. Beside the already mentioned type of microorganisms, minor chemical reactions caused by sulfate-reducing bacteria, homoacetogenic archaea or iron(III)-reducing bacteria may be present. The microorganisms can be present simultaneously; hence, competitive behavior between the different species is expected.

Benchmark problem description
In the following, the general benchmark case is defined and afterwards based on this model Eclipse and DuMu x are compared as an initial contribution. The corresponding simulation files including the geological model, initialization, and schedule are provided as ECLIPSE formats in [16,17] to allow the comparison of other simulators.

Field characteristics/static model
A corner-point grid based on a semi-artificial geological structure is used for the spatial discretization of the simulation. Overall, it consists of 44652 (61 × 61 × 12) grid cells with a dimension of 50 m × 50 m (x-and y-direction) and a varying thickness (Fig. 1). The petrophysical properties (porosity and permeability) are distributed heterogeneously, and the permeability is additionally defined as anisotropic.
The average porosity is 15%, and the mean horizontal permeability is 143 mD (k v ∼ 3 mD), which can be observed in some sandstone formations in Northern Germany (permeability distribution based on a modified poro-perm-correlation).
The statistical distribution of the petrophysic properties can be found in Table 1.
The system is initialized with a pressure of P GWC = 81.6 bar at the gas-water-contact at a depth of 1210 m. Further, a transition zone is established by the capillary pressure (Brooks-Corey parameter λ = 2.0, P e = 0.1 bar [18]) separating the gas and water zone. The initial gas composition in the gas zone is defined as natural gas (see Table 2) [19].

Schedule
In the present study, the schedule comprises of two sections: 1) Conversion from a natural gas storage into UHS and 2) Regular storage operation. In both sections, the injection/production occurs along a single well located in the center of the structure. Further, the injected gas compositions remain constant (see Table 2). With respect on the stepwise development of UHS, low hydrogen concentration of 10% may be interesting in first field projects (e.g. [2]), but also higher concentrations may become attractive later on. The first section, the conversion cycle, is characterized by a bottom-hole-pressure controlled injection which is incrementally increased from 90 bar to 102 bar (step size: 4 bar) to increase the reservoir pressure and raise the hydrogen content in the storage. In total, four cycles consisting of 60 days with one month idle time between each cycle are performed. After the conversion cycles, the regular storage cycles are conducted. The regular operation consists of alternating injection and production with a constant rate of q = 293.13 mol/s = 6 · 10 5 Sm 3 /day. The duration of the injection is identical to the withdrawal (90 days) to equalize the cumulative volumes. Similar to the conversion cycles, the regular storage cycles are separated by idle times.

Implementation
To mimic the flow in the subsurface on continuum scale, various simulators are available on the market. While commercial simulators are focused on user-friendly usability with reasonable outcomes, open-source simulators often allow better adjustments but as a consequence typically leading to a lack in simplicity. In the present study, the potential of two different simulators to mimic typical phenomena in UHS is evaluated: Firstly, the commonly used simulator ECLIPSE, developed by Schlumberger [20] with the focus on hydrocarbon fields, and secondly, DuMu x [21], which is developed by the University of Stuttgart on open-source base for scientific purposes. A brief overview of common and differing features is listed in Table 3.
Four different cases, which vary in the activity of microorganisms and the hydrogen content in the injection stream, are defined (see Fig. 2). While DuMu x can be used in all cases, ECLIPSE is not capable of simulating complex biochemical reactions and is therefore only assessed in the first two cases.

Schlumberger ECLIPSE
ECLIPSE [20] is a commercial reservoir simulator developed and distributed by Schlumberger and widespread in the oil and gas industry.

Fluid model
The fluid model in ECLIPSE (E300) consists of two phases (gas and water) with a composition of eight components: water (H 2 O), methane (CH 4 ), hydrogen (H 2 ), carbon  [25]. More precisely, the EoS is used for the determination of the fluid density and in a modified version (Søreide and Whitson [22]) for the evaluation of the phase equilibria. The dynamic viscosity of the gaseous phase is calculated by a combination of correlations of Stiel and Thodos [23] (low-pressure) and Lohrenz, Bray and Clark [24] (high-pressure correction).
Remarkable is the usage of a simplified version of Stiel and Thodos in which the additional term for the component hydrogen is not included. This exception was published within the initial paper of Stiel and Thodos [23] and is also valid for the component of helium. The simplified form may have an impact on the accuracy of the viscosity in cases with high hydrogen containing systems.
As the implementation of complex biochemical reactions is not possible in ECLIPSE, the simulation is limited on the prediction of the two-phase flow.

DuMu x
DuMu x is an open-source simulator which is in development by the University of Stuttgart (Institute of Modeling Hydraulic and Environmental Systems) since 2007 [21]. It is based on DUNE [26] and is provided as an additional  Comb. of original Stiel and Thodos [23] and LBC [24] Comb. of simplified Stiel Thodos [23] and LBC [24] module in order to mimic the fluid flow in porous media including chemical reactions. DuMu x is the abbreviation for 'DUNE for Multi-{Phase, Component, Scale, Physics, ...} flow and transport in porous media'. In recent publications, biochemical reactions by microbial activities were implemented by Hagemann [27] and underground hydrogen storage operations were simulated. For the present study, DuMu x 3.2 is used.

Fluid model
The fluid model in DuMu x comprises of two phases (gas and water) which are composed of seven components: water (H 2 O), methane (CH 4 ), hydrogen (H 2 ), carbon dioxide (CO 2 ), ethane (C 2 H 6 ), one pseudo component (C 3+ ), and nitrogen (N 2 ). Apart from the chemical components, the methanogenic archaea are implemented as an additional pseudo component, which does not affect the fluids' viscosity and density. Analogous to the implementation in ECLIPSE, the Peng-Robinson EoS is used to describe thermodynamic relations. However, in DuMu x , the EoS is merely used to determine the density of the gaseous phase. Further, the phase equilibria are defined by a combination of Raoult's law and Henry's law (vaporization and dissolution in the phases), leading to a different implementation in DuMu x and ECLIPSE. The impact of temperature and pressure on the fluid mixture's viscosity is modeled using two correlations. First, the correlation of Stiel and Thodos in its full extended form which is not covered by ECLIPSE, is used to determine the low-pressure viscosity. Further, the Lohrenz, Bray and Clark correlation is applied to determine the corrected high pressure viscosity.

Microbial activity
In contrast to the simulation in ECLIPSE, the fluid system is extended by an additional component representing the methanogenic organisms. The dynamics of this pseudo component are mainly governed by the bacterial growth (double Monod model) and the decay/maintenance as expressed in Eq. 2.
The mathematical model of the reactive transport in porous media in relation to UHS was implemented first by Hagemann [28]. The implementation in DuMu x is based on this model with minor adaptions: φ ∂ ρ g c k g S g + ρ w c k w S w ∂t In accordance with Strobel et al. [29], the characteristic growth parameters are based on recent literature values and are defined as listed in Table 4. However, these parameters may strongly vary due to the presence of different strains of microorganisms coming with diverse favorable conditions (pH, temperature and salinity). Therefore, an independent assessment has to be conducted for different formations.

Fluid model
While most of the fluid properties are identical in both simulators, merely the dynamic viscosity is divergent due to varying correlations. In DuMu x , the extended form of Stiel and Thodos [23] is implemented for this study, while ECLIPSE uses the simplified version, which may be stated as sufficient for systems containing less than 10% of hydrogen. Figure 3 depicts the behavior of the dynamic viscosity in dependency of pressure for study-relevant fluid compositions (see Table 2). For cases with a minor content of hydrogen a satisfying match is observable. However, during UHS higher percentages of hydrogen may be interesting and thus, the difference between both correlations increases. Here, the extended correlation is stated as more accurate [23] and leads to an approximate 6% higher viscosity than

Simulation
Overall, six simulation runs (2x ECLIPSE, 4x DuMu x ) varying in microbial activity and hydrogen concentration were performed. In the first step, the fluid dynamic processes (convection and diffusion) are compared to classify the similarity of both simulators. Afterwards, the impact of metabolism maintained by methanogenic microorganisms on the reproduction of hydrogen is investigated.

Fluid dynamics
To evaluate the fluid dynamic behavior of both simulators, the field rates, the average field pressure, and the spatial spreading of the component hydrogen are evaluated. Here, the focus is placed on Case I and II as these instances neglect the effects of bio-reactions maintained by microorganisms. Figures 4, 5 and 6 depict the general field parameters consisting of average reservoir pressure and injection/production rates of the different cases. During the conversion period, the injection rates are pressure controlled with an incremental increase of the bottomhole-flowing pressure. Almost no differences are remarkable between both cases (high and low injection concentrations) according to Figs. 5 and 6).
Altogether, it can be said: the higher the hydrogen content, the lower the density, yielding a higher average reservoir pressure. Nevertheless, the reservoir pressures behave With a focus on the consistency of both simulators, an almost ideal match in terms of pressure and field rate is achieved for Case I. However, differences in the results of both simulators are apparent in the second case. They are mainly influenced by the varying viscosity models of both simulators, yielding a deviation of 0.4 bar (≈0.4%) after four storage cycles.
In the next step, the spatial distribution of the component hydrogen in the gaseous phase is evaluated (see Figs. 7 and 8). As expected, the highest hydrogen content remains close to the wellbore and decreases with growing distance. However, further dissimilarities can also be seen here: DuMu x leads to a smaller cushion containing higher concentrations of hydrogen, which can be explained by the mixing behavior. ECLIPSE shows a higher mixing potential with the initial gas, leading to more extensive spreading of the component hydrogen with a lower concentration around the wellbore. This trend can be observed for all cases.
In general, the efficiency of an UHS can be defined by the amount of reproduced hydrogen previously injected into the storage formation. As the production rates are equal for all simulations, the H 2 -fraction in the production stream (see Figs. 9 and 10) can be used as a sensitive parameter.
Generally, it is reasonable that the share of hydrogen in the production stream is close to the injection concentration at the beginning of each cycle. This fraction decreases with time as the gas is recovered from more distant regions of the formation, which is more likely to be mixed with the initial gas. Overall, this mixing can be counted as a loss of hydrogen and is therefore mainly responsible for the efficiency of the energy storage. The observation of a higher mixing potential of hydrogen with the cushion gas in ECLIPSE is also indicated in the evaluation of the H 2fraction in the withdrawal stream. The mixing behavior during idle times is mainly governed by molecular diffusion. Therefore, stronger mixing can be explained by stronger (molecular and/or numerical) diffusion in ECLIPSE. Overall, the results obtained with DuMu x in the first two cases show always a higher concentration of hydrogen than ECLIPSE (after 4th cycle ≈1.5%).
Concerning the modeling of fluid dynamics, these results can be stated as a desirable match. Both simulators show only minor differences in the range of less than five percent. Even the selection of varying viscosity models seems to have a lower impact than expected.

Biochemical effects
In another step, the specific effects of hydrogenotrophic microorganisms (here: methanogenic archaea) on the storage operation are investigated. As mentioned previously, the microbial conversion process is only covered by DuMu x (Case III and Case IV). Figure 11 depicts the spatial distribution of the microbial population where the parameter n MG represents the ratio of number of microorganisms per m 3 to the initial number (so-called: microbial density).
As expected, the bacteria population follows the spatial distribution of hydrogen and achieves the highest numbers around the wellbore. Due to the fact that hydrogen invades further regions in the second case, the spatial spreading of the microorganisms is more extended. However, the maximum density of the microbes is 57.3 · 10 8 m −3 respectively, which strongly depends on the characteristic growth parameters and indicates that the growth and decay of the population are balanced. In the present study, this seems to be already achieved at low percentages of hydrogen (no substrate limitation) yielding a large cushion with the maximum microbial density.
Further field parameters also indicate the higher microbial activity in highly hydrogen-containing systems (Figs. 4-6). Due to the microbial conversion process, the average field pressure is lower than in the previous cases, and the injection rate increases during the initialization cycles (pressurecontrolled) with an amplifying tendency of higher H 2concentrations. Based on Figs. 9 and 10, the effect of microorganisms on the recovered hydrogen is visible. Here, two contrary effects can be identified: while Case III yields a 0.7% lower H 2 -fraction (7% of the injection concentration respectively) in the production stream, Case IV leads to a slightly higher concentration of H 2 during the withdrawal of gas. The first phenomenon is logical with the previously stated conclusions of the consumption of hydrogen by microorganisms. However, the second observation may be more influenced by the initialization/conversion cycles, where lower reservoir pressures increase the injection rate and therefore result in a minimal higher H 2 -content in the storage formation compared with Case II. Further, this yields a higher concentration in the production stream while this phenomenon vanishes with a higher number of storage cycles and may turn over to a reduced fraction of hydrogen.
Here, a remarkable drop in the hydrogen concentration is expected in subsequent cycles.
In the next step, it may be relevant to evaluate the absolute and relative quantities of hydrogen converted by the microorganisms (Fig. 12). It can be seen that the absolute amount of substance of hydrogen converted is approximately 50% higher in Case IV (n H 2 conv = 7.68 · 10 8 mol) than in Case III (n H 2 conv = 4.93·10 8 mol), indicating a nonlinear relationship of the injected hydrogen concentration to the conversion rate. Nevertheless, the relative share of the converted hydrogen of the cumulative injected hydrogen shows a contrary behavior. As observed earlier, the impact of the microbial activity is more sensitive to storage scenarios with hydrogen as a minor component. Additionally, the injection operation of the storage is visible with small drops indicating a higher injection than conversion rate.
As previously shown, the hydrogen concentration in the production stream is very sensitive and the microbial effects may be superimposed by other phenomena. Therefore, the CO 2 concentration is additionally analyzed (Figs. 13 and 14) as CO 2 acts as a carbon and energy source for the methanation.
Here, the impact of the microorganisms becomes more visible. In cases without microorganisms, the CO 2 -fraction  increases with time and the number of cycles, while the microorganisms lead to a reduction of CO 2 in the production stream. Overall, the microbial activity influences hydrogen recovery, and the efficiency in scenarios with hydrogen as a minor component is more substantial than in cases with higher concentrations within the first storage cycles. Further, the microbial activity does not behave linearly in dependency of the hydrogen injection concentration.
The percentages of hydrogen recovered after four regular storage cycles are presented in Table 5. Both simulators show congruent results with approximately 39% of total injected volume and almost 99% concerning the hydrogen injected during the regular storage operation. However, DuMu x achieves slightly higher percentages than ECLIPSE, which coincides with previous observations. Surprisingly, the consumption of 21% of hydrogen by microorganisms (Case III) leads only to a minor reduction of approximately 1.3% (incl. cushion gas) of the recovery. However, the efficiency of the actual storage (excl. cushion gas) seems to be superimposed by the hydrogen injection during the conversion from natural gas storage to a UHS. Therefore, it is expected that these values decrease within the following storage cycles.

General simulation parameters
Besides factors such as a proper realization of physical behavior, also simulation parameters such as computation time and simplicity of implementation should be considered for the simulator selection. As previously mentioned, ECLIPSE offers a good keyword-based way to define simulation cases, especially as a plug-in in Schlumberger Petrel. However, the keyword library is limited to predefined functions and may not implement more complex processes. Here, DuMu x provides a better adjustment for specific needs but along with the downside of a more challenging implementation because of its open-source structure in C++.
Apart from user-friendliness, computation efficiency can play a significant role during the selection of a simulator. The simulation runtime of DuMu x is approximately 100 times the one of ECLIPSE, which can be explained by a less sophisticated numerical implementation in the case of the open-source simulator. Additionally, a simulator's accuracy and numerical correctness are relevant and have to be considered. While DuMu x uses the Finite-Volume-Method (FVM) to solve the numerical problems, ECLIPSE applies the Finite-Differences-Method (FDM). This different  implementation has an impact on the material balance error: In the case of FVM, the material balance is solved more conservatively leading to smaller errors in this relation.
Concerning the material balance error of the component hydrogen, the simulation in DuMu x achieved a maximum error of 9.8 · 10 −2 mol. Contrarily, ECLIPSE occasionally shows errors in the magnitude of 10 8 mol, leading to less accurate results. However, this error, in the range of a few percentages of the stored mass, may be influenced by a lower output accuracy, and these errors seem to equalize due to differing signs.

Conclusion
With increasing interest in the storage of energy in the subsurface, the prediction of the processes occurring in UHS is becoming more relevant. In this field the selection of a suitable reservoir simulator plays a crucial role. In the present study, both, the commercial software ECLIPSE and the open-source simulator DuMu x showed a good potential to cover these processes. While ECLIPSE offers a good mix of user-friendliness and efficiency, DuMu x allows an excellent potential to adapt it to specific needs. With respect to the simulation of fluid dynamics, both simulators lead to congruent results within a margin of less than 5%. However, the mixing behavior of the injected fluids with the initial cushion gas seemed to be more intensive in ECLIPSE. Concerning specific processes in UHS, it was possible to implement biochemical reactions maintained by microorganisms (here: methanogenic) into the simulation model in DuMu x . Although the metabolism of the methanogenic microorganisms became stronger with increasing hydrogen concentrations, the impact on the reproduction of hydrogen is more remarkable in cases with low hydrogen percentages. With a focus on the detection and monitoring of activity of methanogenic archaea, it was observed that a reduction in field pressure and CO 2fraction in the production stream are the most sensitive indicators. Surprisingly, the H 2 -fraction became a less reliable indicator as it may be superimposed by other effects. The simulation results indicated a good efficiency in terms of hydrogen reproduction within the first four storage cycles, which will most likely decrease with advancing time.
As the present study is based on a semi-artificial grid and no observed data are available, no statement about the physical correctness of both simulators can be made.
Overall, it can be said that both simulators are good options to simulate the fluid dynamics in UHS. ECLIPSE is a good choice if time is an important parameter leading to reasonable results. DuMu x on the other hand, may be interesting for cases, where unique effects as microbial activity play a crucial role. Here, the implementation of other bacteria (e.g. sulfate-reducing) or potential geochemical effects due to the presence of hydrogen may be good extensions for the simulation model in DuMu x .

Supplementary information Not applicable
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://creativecommons. org/licenses/by/4.0/.