The Onset of Convection in a Two-Layered Porous Medium with Anisotropic Permeability

We consider convection in a horizontal porous layer of uniform thickness which is heated from below and which is composed of two anisotropic sublayers with principal axes lying in the three coordinate directions. The aim is to determine criteria for the onset of convection by finding the critical Rayleigh number, wavenumber and roll orientation relative to the coordinate axes. The full set of nondimensional parameters has at least six members even when the sublayers are considered to be thermally isotropic, and therefore, we select some special cases in order to illuminate the type of qualitative behaviour which may be expected. One such case is where the anisotropic sublayers are identical except that one sublayer is rotated by an angle of 90∘\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$90^\circ $$\end{document} to the other. In this situation, the most unstable roll is found to lie at an angle of ±45∘\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\pm \,45^\circ $$\end{document} to the principal axes. It is also found that fluid particles exhibit a mean longitudinal flow as they circulate about the vortex axis. This drift along the vortex is balanced by an equal and opposite drift in the two neighbouring vortices. Convection in each sublayer is shown to be two dimensional even though the full flow field is three dimensional.


Introduction
The great majority of published theoretical and experimental investigations on convective instabilities in porous media has considered isotropic media. However, it is often the case that porous media, especially those which arise in geothermally active zones, are not only anisotropic, but layered. Many authors have considered the effect on the onset and subsequent development of these naturally occurring properties, but the combination of the two has not yet been considered. Our aim here is to provide some first steps in this direction, and especially to uncover new phenomena which arise only when both of these effects are present.
With regard to layering, there is now a modest presence in the literature of works which have considered the onset of convection when the layer is heated from below. A very early example is that of Donaldson (1962), who used a finite-difference method to compute the nonlinear flow and temperature fields in a two-layer system, the lower layer of which was impermeable (i.e. solid), but finitely conducting. Later, papers have studied the onset problem in more detail. For example, Mojtabi and Rees (2011) undertook a comprehensive study of the effect of conducting boundaries above and below the porous layer and their work showed that the solid boundaries may be used in practice to mimic a very wide range of boundary conditions for the porous layer itself when considered in isolation-this is a very useful practical result. Weakly nonlinear studies by Riahi (1983) and Rees and Mojtabi (2011) give detailed criteria for whether the convection pattern which is selected takes the form of a roll pattern or has a square planform. Masuoka et al. (1979) derived criteria for the onset of convection in a porous system with two sublayers, and they also computed two-dimensional nonlinear flows. Rana et al. (1979) also modelled the Pahoa reservoir using a three-layer model of porous strata. A comprehensive analysis of the onset of convection and the post-critical heat transfer were then presented in two papers by O'Sullivan 1980, 1981). This was extended to include the effects of thin, highly impermeable 'sheets' (McKibbin and Tyvand 1983), and thin highly permeable cracks (McKibbin and Tyvand 1983) within the layer. A threedimensional weakly nonlinear stability analysis of convection in a layered porous medium was performed by Rees and Riley (1990). They found that, in some circumstances, threedimensional convection with the form of cells of square planform is sometimes favoured over rolls. In some restricted ranges of parameter values, the neutral curve has a double minimum.
In many practical situations, the porous materials are anisotropic in their mechanical and thermal properties. Anisotropy is generally a consequence of preferential orientation or asymmetric geometry of grain or fibres. An important example is geological systems with anisotropic sediments and rocks. Another example of such a medium is loft insulation which usually has lower permeability across the insulating layer than it has in the perpendicular directions. Convection in anisotropic porous media has attracted the interest of many researchers over the last 40 years. In the pioneering works by Castinel and Combarnous (1974) and Epherre (1975), it was shown that anisotropy in the mechanical and thermal properties affects the marginal stability condition as well as the preferred width of the convection cells. On the other hand, Kvernvold and Tyvand (1979) showed that even three-dimensional anisotropy does not lead to any new mathematical difficulties or essential new flow patterns at convection onset compared with isotropy. This is true only as far as one of the principal axes of the anisotropic medium is vertical. This requirement was maintained in almost all former works in the field, see the review articles by McKibbin (1984) and Storesletten (1998), well as McKibbin (1986 and Nilsen and Storesletten (1990). The first papers, where none of the principal axes is vertical, seem to be Tyvand and Storesletten (1991) and Storesletten (1993). They considered a horizontal porous layer with anisotropy in the permeability or in the diffusivity. This was sufficient to achieve qualitative new flow patterns with a tilted plane of motion or tilted lateral cell walls. Other earlier papers considering onset of convection in anisotropic porous media are: Trew and McKibbin (1994), Rees and Storesletten (1995) and Rees et al. (2006), to mention a few, while recent research has gravitated towards determining how the combination of anisotropy with other effects modifies onset criteria; see Hill and Morad (2014) Gaikwad and Dhanraj (2016) and Raghunatha et al. (2018).
The configuration considered in the present paper combines anisotropy with layering, and this appears to be the first time where these have been considered. Such a combination also has competing mechanisms for deciding on the favoured orientation of the convection rolls in a fairly obvious way. When the layer consists of two sublayers with different anisotropies, then the individual sublayers may favour different roll orientations. However, such a system could, in principle, be governed by 17 nondimensional parameters: five permeability ratios, five diffusivity ratios, six rotation angles for the principle axes of the two permeability tensors (assuming that the principle axes of the diffusivity tensors coincide with the respective permeability tensor) and the height of sublayer 1 relative to that of the layer. This unwieldy generality shall be stripped down to having only two parameters in order to simplify the analysis and the presentation of results, but this will nevertheless reveal a novel aspect to the flow which arises from the simultaneous presence of anisotropy and layering. Thus, 123 the principle axes for the permeability will coincide with coordinate axes, the diffusivity is assumed to be isotropic and the two layers will have identical properties where the two principal permeabilities in the horizontal directions are different. In addition, the principle axes for sublayer 2 are obtained by rotating those for sublayer 1 by 90 • about the vertical axis. The two parameters which are retained are a permeability ratio, γ , and the relative depth of sublayer 1, δ.

Governing Equations
We consider a horizontal fluid-saturated porous layer of uniform vertical thickness which is heated from below and cooled from above. The layer consists of two anisotropic sublayers not necessarily of the same thickness and not necessarily with the same principal axes. The anisotropy is manifest in terms of the permeability of the sublayers but not in terms of the thermal diffusivity. We choose the x and y axes to be horizontal and the z-axis to be vertical. The lower sublayer is designated layer 1, and the upper sublayer is layer 2. The permeabilities in the three coordinate directions of layer i are given by K i,x , K i,y , K i,z where i = 1, 2. Subject to the Oberbeck-Boussinesq approximation and on assuming that Darcy's law is valid, the governing equations in layer i may be written in the form, In these equations, the coefficients take their usual meanings: p is the pressure, ρ the density, μ the dynamic viscosity, g gravity, α thermal diffusivity, T the temperature. The temperature of the lower and upper surfaces is T h and T c , respectively, where T h > T c . The velocities in the x, y and z directions are given by u, v and w, respectively. The values λ i are the heat capacity ratios for the fluid relative to that of each saturated sublayer.
As we are considering three-dimensional flows, we may eliminate all three velocities using Eqs. 2 to 4. Therefore, we obtain the pair of governing equations, The boundary conditions at the lower surface are that the temperature is T h and that the vertical velocity is zero; in terms of pressure and temperature, these become while the equivalent conditions on the upper surface are The specification of the problem is completed by the conditions at the interface which is located at z = δ d. Here, we require the continuity of temperature, heat flux, pressure and vertical velocity, and therefore, the following four quantities need to be continuous at the interface. We note that the continuity of the vertical velocity and the temperature causes the heat flux condition to be identical to that of a conduction condition.
Equations (6) and (7) and the boundary and interface conditions (8) to (10) may be nondimensionalised using the physical properties of layer 1 using the following substitutions, On dropping the overbars, we obtain where the permeability ratios are given by and where Therefore, we have the following fixed values, Equations 12 and 13 are to be solved subject to the following boundary and interface conditions, 123 The basic conducting state is given by where p 0 is the pressure at the upper boundary, and where the subscript, b, denotes the basic state. In Eqs. (12) and (13), the Rayleigh number is based upon the vertical permeability of layer 1, and is, therefore,

Linear Stability Analysis
The equations governing the evolution of small-amplitude disturbances to the basic conducting state are obtained by subtracting out the conducting state and by linearising the resulting equations. Specifically, we set where i = 1, 2, and neglect all nonlinear terms inp i andθ i . In this manner, we obtain the following disturbance equations The principle of the exchange of stabilities applies here, and therefore, the criterion for the onset of convection corresponds to a zero exponential growth rate which is equivalent to neglecting the time-derivative term in (24). We may now introduce roll solutions using the following substitution, into Eqs. (23) and (24) where φ is the angle between the roll axis and the y-axis in the horizontal plane. Thus, φ = 0 corresponds to a roll with its axis pointing in the y-direction. These equations then reduce to the ordinary differential eigenvalue form, where the boundary and interface conditions are identical to those given in (17) to (19), and where the primes denote ordinary derivatives with respect to z. We regard the Rayleigh number as the eigenvalue, and it depends on the values of the five independent ξ -values, the wavenumber k, the roll angle φ and the thickness of layer 1, δ. In general, the neutral curve, i.e. the variation of Ra with k, exhibits a single well-defined minimum value which is termed the critical Rayleigh number and is denoted by Ra c . The corresponding value of the wavenumber is denoted by k c . In some cases, such as those exemplified by the layered configurations of McKibbin and O'Sullivan (1980) and Rees and Riley (1990), it is possible to have a neutral curve which consists of two minima; given that their configurations consisted of isotropic sublayers, we would expect by continuity that double minima are generally possible for anisotropic sublayers. However, the simplification of the general problem which was mentioned in the Introduction serves to eliminate that possibility in the present work. Numerically, the most convenient manner in which critical Rayleigh number and wavenumber might be determined is to solve Eqs. (26) and (27) together with the following system, This latter pair of equations was obtained by differentiating (26) and 27 with respect to k, whereP andˆ are the respective k-derivatives of P and . We have set because it is equivalent to the condition that the value of k corresponds to a minimum (or maximum) in the neutral curve. Equations (26) to (29) form a double eigenvalue problem with Ra and k as the eigenvalues; as far as we are aware, this type of extended system was first used in Rees and Bassom (2000) in such contexts.
In the present configuration, the critical Rayleigh number also depends on the roll angle, φ, and therefore, it is essential to minimise Ra c with respect to φ. On definingP and¯ as the respective φ-derivatives of P and , the result of taking the φ-derivative of Eqs. (26) and (27) is to yield the system, where ∂Ra ∂φ = 0 (33) has been imposed. This triple eigenvalue problem has eigenvalues, Ra, k and φ.

Numerical Method
Equations (26) and (27) in order to obtain a fifth-order system of first-order ODEs. The reason is straightforward: Eqs. (26), 27 and the boundary conditions are homogeneous and the only manner in which an initial value problem solver such as the Runge-Kutta method might find nonzero solutions is to impose a normalisation condition such as This brings the total number of boundary conditions to five, and therefore, the full system may now be solved numerically. We note that any nonzero value of (0) may be used as the computed value of Ra is independent of this choice. For System 2, a similar idea was used and the equation supplemented the equations already quoted. The extra boundary condition needed was taken to beˆ although any value may be used. For System 3, we used as the extra equation and¯ as the extra boundary condition. For each case, the codes were organised so that any of the parameters might be changed slowly in order to determine how the critical Rayleigh number varies with the chosen parameter. We chose to use 100 intervals with 50 intervals in each of the ranges 0 ≤ z ≤ δ and δ ≤ z ≤ 1, i.e. in each sublayer. Given the order of accuracy of the method, the computed Rayleigh numbers were found to have more than five significant figures of accuracy when compared with analytical solutions or published data for layered systems. For example, Rees and Riley (1990), using an eigenfunction expansion technique, quotes the critical Rayleigh number Ra c = 6.917323 × 4π 2 for the isotropic layer with δ = 0.323889, k c = 1.540927π, ξ 1 = 1 and ξ 2 = 0.0657659 (a special case where the minimum of the neutral curve is a quartic minimum rather than the usual parabolic minimum), and this critical value was realised to the given number of figures.
We have restricted severely our choice of governing parameters in order not to overload the paper with merely numerical information. Here, we shall concentrate on cases where the two sublayers are identical in properties except that the principle axes of the permeability tensor for the upper layer are rotated by 90 • compared with those of the lower layer. Thus, we take the following set of permeability ratios: where γ will take various values.

Some Analytical Results and Axial Drift
We begin by noting that there is one instance when it is possible to write down relatively easily the solution to the perturbation equations given in Eqs. (26) and (27). When the roll angle is 45 • , the perturbation equations are identical in the two sublayers. In such an instance, both reduce to the form The solution may then be written as, where Ra = (π 2 + k 2 ) π 2 + k 2 2 (γ + 1) The Rayleigh number is minimised when and hence, The expressions given in Eqs. (45) and (46) reduce to π and 4π 2 in the isotropic case, γ = 1.
In the extreme case when γ = 0, all the fluid which occupies sublayer 1 is constrained to move within the y-z plane, given that the permeability in the x-direction is zero. Likewise, in sublayer 2, fluid may only move in the x-z plane. This observation provides a mechanism for the fluid to undergo an overall drift in the direction of the axis of a convective roll. This may be demonstrated as follows. Given that we have taken the roll direction to be φ = 45 • , we may use the solutions given in (43) to write down the three velocity components in the form, The fact that u and v are proportional to one another (which is true in general) means that the direction of the horizontal component of the velocity of a fluid particle is determined solely by the ratio, ξ i,x /ξ i,y , and the sign of cos π z. Once the fluid particle passes through the interface, then the direction of the horizontal component changes suddenly. This is illustrated in Fig. 1 where the depiction corresponds to γ = 3 and δ = 0.5, and it shows a particle which is moving along the top and bottom surfaces, as opposed to one within the interior of the roll. Blue arrows indicate the particle travelling along the cold upper surface. When it reaches the grey line, which is the roll boundary, it descends to the hot lower surface and then travels along a red arrow. From the above, the particle describes a simple zigzag pattern with the mean motion of the particle being in the direction of the roll axis. This means that there is a fluid drift along that axis. However, this is not a mechanism for causing an overall horizontal movement of fluid because the two neighbouring rolls are subject to an equal and opposite drift. The situation is likely to be very considerably different were the fluid to be contained within a cuboidal cavity, but with the same anisotropy and layering profiles.  Figure 2 is designed to show that the speed of drift depends on how close the fluid particle gets to the upper and lower surfaces. Particle paths have been determined numerically using a simple fourth-order Runge-Kutta scheme to solve where u, v and w are given in Eq. (47). The initial conditions correspond to x and y being located precisely equidistant between the cell boundaries, but with y = 0. The initial values of z are indicated in the figure. In this figure, the sublayers have equal width, so δ = 0.5. For each of the four values of z(0), two subfigures are displayed; the left-hand one shows how x, y and z vary in time, while the one on the right displays the plan view. When z(0) = 0.01, the particle spends most of its time close to either the upper or the lower surface with a brief transition when it is either ascending or descending very close to a vertical boundary of the roll; during this time, its horizontal velocity is small. The associated plan view looks very much like Fig. 1, and the particle resides close to the corners of the zigzag for a substantial period of time. As z(0) increases, the particle becomes more confined to a region surrounding the central core of the roll. Initially, this means that the residence time near the corners of the roll reduces, and this allows the particle to drift more rapidly along the roll axis, although the time taken to execute one circuit about the axis decreases. Eventually, as z(0) approaches the interface at δ = 0.5, the overall speed of movement decreases. In this region, u and v, which are proportional to cos π z, take small values, and this also has the  Fig. 2. The values of δ are as shown effect of decreasing the drift velocity. Therefore, the drift velocity of a particle increases at first as z(0) increases from zero, reaches a maximum, and then decreases towards zero when z(0) = 0.5. Figure 3 shows the analogous situation where z(0) = 0.19 has been chosen, but where δ takes four different values. When δ = 0.1, the interface is below where the chosen particle is closest to the lower boundary. Therefore, the particle will always reside in sublayer 2, and its path will always be in exactly the same horizontal direction. This particle will not drift, as may be seen in the right-hand subfigure. In the second case, δ = 0.2, the interface is just above where the particle is closest to the lower surface. Therefore, it spends little time in layer 1 where its direction of movement could be described as being ESE (East South East). Very quickly, the particle enters layer 2 and switches direction to SSE. After a while, it reverses direction, NNW, continuing to rise to its maximum height. Thereafter, it sinks, reverses direction to SSE, and then briefly enters layer 1 once more, travelling in the ESE direction. Thus, the particle executes most of a closed loop such as one would have in an isotropic layer. As δ increases towards 0.5, the distance over which the particle passes in opposite directions decreases until we obtain no reversal when δ = 0.

Critical Values
It is important to note that restricted set of parameters which we have chosen for this study are such that the neutral curves will always have a single minimum value for all choices of γ , δ and φ. Therefore, we shall present only those values which have been minimised over k. Figures 4 and 5 show how Ra c and k c vary with the roll angle for various choices of γ . In addition, the two sublayers have the same depth, and therefore, δ = 0.5. In each case, the vertical axes have been scaled with respect to the critical values for the isotropic single layer (i.e. for γ = 1), namely 4π 2 and π, respectively. When γ = 1, the sublayers are identical and isotropic, the stability problem reduces to that of the classical Darcy-Bénard problem, and therefore, we see that Ra c = 4π 2 and k c = π independently of the roll angle, φ.
For other values of γ , the sublayers are anisotropic, and therefore, Ra c varies with φ. We present the curves for the range 0 ≤ φ ≤ 45 • since the curves are even about both φ = 0 • and φ = 45 • . It is readily seen that φ = 45 • represents the smallest value of Ra c for all values of γ . In the lower sublayer, the x-z plane (i.e. φ = 0) is the preferred plane of motion given the permeability tensor, while the y-z plane (φ = 90 • ) is preferred in the upper layer. Our results suggest that the compromise direction, φ = 45 • , is preferred within the composite layer. Of course, φ = −45 • also minimises the Rayleigh number, by symmetry.
When γ > 1, we obtain the familiar qualitative result that the critical Rayleigh number is reduced compared with the isotropic case and that the critical wavenumber is also smaller. Of interest is the extreme case where γ = 0. This case represents porous sublayers where convection can only take place in the x-z plane in the lower layer and the y-z plane in the upper sublayer. Clearly, then, a fluid particle will trace out a path along one convection roll where its horizontal direction of movement alternates between the y-direction and the xdirection, and it will therefore undertake a spiral-like drift along the roll, unlike what happens in isotropic media where the particle paths form closed loops. It is to be noted that while drift takes place in one direction along any chosen roll, the drift will be in the opposite direction for the neighbouring rolls. Although we do not prove the point here, the drift of particles and their two-dimensional motion in each sublayer also apply in general since the properties of each sublayer are homogeneous. Figures 4 and 5 represent a symmetrical case where the sublayers have equal depth. We now relax that assumption by allowing δ to vary, and we illustrate the effect of that variation for the case γ = 2, which is representative. Figure 6 shows how Ra c varies with φ for a variety of values of δ. Now that the layer is no longer symmetric, it is essential to allow φ to vary between 0 • and 90 • . We present curves for which δ ≥ 0.5 and for which, therefore, the anisotropy of layer 1 is increasingly dominant as δ increases. Indeed, this increasing dominance is seen easily by virtue of the fact that the value of φ at which Ra c is minimised reduces towards φ = 0 • as δ increases. We find that the minimising value of φ reaches 0 • when δ 0.770308 (obtained using an extended version of System 2 where it is insisted that ∂ 2 Ra/∂φ 2 = 0). By symmetry, then, we can predict that the critical value of φ lies between 0 • and 90 • when 0.229692 < δ < 0.770308.
Concentrating now on the these values of φ which minimise the critical value of Ra c , Figure 7 shows how φ c varies with δ for a variety of values of γ . We see that the transition range for φ c becomes smaller as γ approaches unity, for when γ = 1 all roll directions are equally likely. We find numerically that when γ = 1.001, then φ lies between 0 • and 90 • for 0.460522 < δ < 0.539478. Figure 7 shows the situation for positive values of γ , but Figure 8 summarises the full behaviour of the stability characteristics by displaying the regions in γ -δ space where φ c takes the values 0 • or 90 • or else is transitional. It is interesting to note that γ = 1 marks the boundary where φ c changes suddenly from 0 • to 90 • or vice versa. The reason for such a change may be traced to the preferred roll orientation for uniform anisotropic layers. Thus, when γ > 1, and when δ = 1, sublayer 1 occupies the whole layer and φ c = 0 is the preferred roll orientation. At the opposite extreme where δ = 0, sublayer 2 occupies the whole layer, and then, φ c = 90 • . In both cases, there is a range of values of δ (i.e. greater than 0 and less than 1, respectively) within which the preferred roll direction does not change, but there is also an intermediate range of values of δ where the preferred roll direction varies smoothly from φ c = 90 • to φ c = 0 • as δ increases. Precisely, the opposite behaviour in terms of φ c is found when γ < 1. We also note that φ c = 45 • for all values of γ when δ = 0.5, i.e. the sublayers have equal heights.

Conclusions
We have studied the onset of convection in a porous layer consisting of two sublayers with anisotropic permeability. The principal axes for both sublayers lie in the three coordinate directions of which one is vertical. The permeability coefficients in the x, y, and z directions of the lower layer are (γ , 1, 1) and of the upper layer are (1, γ, 1). This means that the sublayers are identical in terms of their anisotropy, but the principal axes are rotated through 90 • relative to one other. The thicknesses are different in general, where d 1 and d 2 are the thicknesses of the lower and upper layers, respectively, and δ = d 1 /(d 1 + d 2 ) is the thickness ratio.
The following is a summary of our conclusions.
• When the thicknesses are identical (δ = 0.5), the preferred roll at convection onset occurs at an angle 45 • to the horizontal principal axes. This means that the critical Rayleigh number Ra c obtains its smallest value at this angle; the resulting values of Ra c are illustrated in Fig. 4, while the corresponding wavenumbers are shown in Fig. 5. • When the thicknesses of the sublayers are different (δ = 0.5 ), the variation in the value of φ at which Ra c is minimised is not straightforward. This is illustrated in Fig. 6 where the variation with φ of Ra c (minimised with respect to the wavenumber only) is shown for values of δ varying between 0.5 and 0.8 and for γ = 2. For this value of γ , the critical roll orientation is 90 • when 0 ≤ δ ≤ 0.229692, is 0 • when 0.770308 ≤ δ ≤ 1, and makes its transition between these two extremes in the intermediate range of values of δ.
The transitions for this value of γ and others are depicted in Fig. 7. • Figure 8 summarises the full behaviour of the stability characteristics by displaying the region in γ -δ space where φ c takes the values 0 • or 90 • , or else is transitional. It is interesting to note that γ = 1, i.e. the isotropic case, marks the boundary where φ c changes suddenly from 0 • to 90 • or vice versa. • It has also been found that the projection of the particle paths onto the x-y plane yields straight lines within each sublayer with a discontinuous change in direction at the interface. In addition, we have determined that individual particles may exhibit a mean drift along the roll axes, although there is a zero mean flow in general; this local drift, which is caused by the combination of anisotropy and layering, is a novel phenomenon.
Finally, it is important to consider what the consequences would be of lifting some of the restrictions on the layer which we have imposed. Specifically, we could retain the fact that the sublayers are identical in their properties, but they could each be allowed to have three different principal permeabilities, three different principal diffusivities, and that the rotation angle (about the vertical axis) between the sublayers may be different from 90 • . In such a case, we believe that the qualitative nature of the stability problem described here would be retained, given that the preferred roll direction for the porous medium which forms the upper sublayer will be different from that in the lower sublayer. Thus, the sole difference would be quantitative in nature given that this case will still display the drifting phenomenon. Further lifting of the restrictions would then involve the sublayers being completely different in their properties, and while exactly the same arguments still apply it then becomes possible to have neutral curves with multiple minima; see McKibbin and O'Sullivan (1980) and Rees and Riley (1990).
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.