# Gas production in the Barnett Shale obeys a simple scaling theory

1. aDepartment of Petroleum and Geosystems Engineering and
2. bCenter for Nonlinear Dynamics and Department of Physics, The University of Texas at Austin, Austin, TX 78712
1. Edited by Michael Celia, Princeton University, Princeton, NJ, and accepted by the Editorial Board October 2, 2013 (received for review July 17, 2013)

## Significance

Ten years ago, US natural gas cost 50% more than that from Russia. Now, it is threefold less. US gas prices plummeted because of the shale gas revolution. However, a key question remains: At what rate will the new hydrofractured horizontal wells in shales continue to produce gas? We analyze the simplest model of gas production consistent with basic physics of the extraction process. Its exact solution produces a nearly universal scaling law for gas wells in each shale play, where production first declines as 1 over the square root of time and then exponentially. The result is a surprisingly accurate description of gas extraction from thousands of wells in the United States’ oldest shale play, the Barnett Shale.

## Abstract

Natural gas from tight shale formations will provide the United States with a major source of energy over the next several decades. Estimates of gas production from these formations have mainly relied on formulas designed for wells with a different geometry. We consider the simplest model of gas production consistent with the basic physics and geometry of the extraction process. In principle, solutions of the model depend upon many parameters, but in practice and within a given gas field, all but two can be fixed at typical values, leading to a nonlinear diffusion problem we solve exactly with a scaling curve. The scaling curve production rate declines as 1 over the square root of time early on, and it later declines exponentially. This simple model provides a surprisingly accurate description of gas extraction from 8,294 wells in the United States’ oldest shale play, the Barnett Shale. There is good agreement with the scaling theory for 2,057 horizontal wells in which production started to decline exponentially in less than 10 y. The remaining 6,237 horizontal wells in our analysis are too young for us to predict when exponential decline will set in, but the model can nevertheless be used to establish lower and upper bounds on well lifetime. Finally, we obtain upper and lower bounds on the gas that will be produced by the wells in our sample, individually and in total. The estimated ultimate recovery from our sample of 8,294 wells is between 10 and 20 trillion standard cubic feet.

The fast progress of hydraulic fracturing technology (SI Text, Figs. S1 and S2) has led to the extraction of natural gas and oil from tens of thousands of wells drilled into mudrock (commonly called shale) formations. The wells are mainly in the United States, although there is significant potential on all continents (1). The “fracking” technology has generated considerable concern about environmental consequences (2, 3) and about whether hydrocarbon extraction from mudrocks will ultimately be profitable (4). The cumulative gas obtained from the hydrofractured horizontal wells and the profits to be made depend upon production rate. Because large-scale use of hydraulic fracturing in mudrocks is relatively new, data on the behavior of hydrofractured wells on the scale of 10 y or more are only now becoming available.

There is more than a century of experience describing how petroleum and gas production declines over time for vertical wells. The vocabulary used to discuss this problem comes from a seminal paper by Arps (5), who discussed exponential, hyperbolic, harmonic, and geometric declines. Initially, these types of decline emerged as simple functions providing good fits to empirical data. Thirty-six years later, Fetkovich (6) showed how they arise from physical reasoning when liquid or gas flows radially inward from a large region to a vertical perforated tubing, where it is collected. For specialists in this area, the simplicity and familiarity of hyperbolic decline make it easy to overlook that this functional form reasonably arises only when specific physical conditions are met. For example, all early decline curves were proposed for unfractured vertical wells or vertical wells with vertical, noninteracting hydrofractures. Physics-based descriptions of such wells are readily available in textbooks, such as those by Kelkar (7) and Dake (8). However, these books focus on radial or 1D flow from infinite or semiinfinite regions to pipes or planes. The resulting decline curves do not apply to the wells we describe here.

The geometry of horizontal wells in gas-rich mudrocks is quite different from the configuration that has guided intuition for the past century. The mudrock formations are thin layers, on the order of 30–90 m thick, lying at characteristic depths of 2 km or more and extending over areas of thousands of square kilometers. Wells that access these deposits drop vertically from the surface of the earth and then turn so as to extend horizontally within the mudrock for 1–8 km. The mudrock layers have such low natural permeability that they have trapped gas for millions of years, and this gas becomes accessible only after an elaborate process that involves drilling horizontal wells, fracturing the rock with pressurized water, and propping the fractures open with sand. Gas seeps from the region between each two consecutive fractures into the highly permeable fracture planes and into the wellbore, and it is rapidly produced from there.

The simplest model of horizontal wells consistent with this setting is a cuboid region within which gas can diffuse to a set of parallel planar boundaries. Fig. 1 illustrates the well as 10–20 hydrofractures that are high and long, spaced at distances of around . The fact that this is the right starting point for these wells was recognized by Al-Ahmadi et al. (9), and the diffusion problem in this setting has been studied by both Silin and Kneafsy (10), and Nobakht et al. (11).

Examining Fig. 1 helps one to understand how gas production evolves. When a well is drilled and completed, the flow of gas is complicated and difficult to predict, particularly because the water used to create it is back-produced. In practice, the resulting initial transients last around 3 mo. After that time, gas should enter a phase where it flows into the fracture planes as if coming from a semiinfinite region.

Gas flows according to Darcy’s law through a system of microfractures, cracks, reopened natural fractures, faults, and failed rock. This multiscale and loosely connected flow system is created by the high-rate hydrofracturing of shale rock. It is fed by the rock matrix, where gas is stored (adsorbed) in very small pores.

It turns out that the gas effectively flows along paths that are straight lines (hence, the setting is sometimes called “linear flow”) perpendicular to the fracture planes. During the flow, the initially high gas pressure diffuses toward the hydrofractures, which are kept at a low pressure. This gas pressure diffusion creates a gas production rate proportional to the inverse of the square root of time on production.

At some point in time, gas flow causes the pressure along the midplane between the hydrofractures to drop below the original reservoir pressure, and gas production slows down relative to the square-root-of-time behavior. We call the time when this happens the interference time. Eventually, the gas is so depleted that the amount coming out per time is proportional to the amount of gas remaining. This is the classic condition for exponential decay. Thus, after a long enough time, the rate of production declines exponentially. The pressure-dependent coefficient describing the diffusion of gas pressure is called the “hydraulic diffusivity of gas.” Physically, it is unrelated to the molecular diffusion coefficients.

The more closely spaced the hydrofractures are, the higher will be the initial rate of gas production but the more quickly will the interference time be reached. These intuitive considerations are consistent with the mathematical results of Silin and Kneafsy (10) and Nobakht et al. (11), and with the analysis that follows.

## Results

### Model.

Hydraulic fracture in horizontal wells creates a network of cracks in rock that was previously impermeable, allowing gas to move. The true geometry is very complicated. Here, we explore the possibility that it suffices to treat the rock surrounding a well as a cuboid region in which the permeability is greatly enhanced over the surrounding area but is uniform. We focus first on the single region depicted in Fig. 1 (Lower). The solution of a problem with N hydrofractures is obtained trivially by multiplication once the problem for a single region bounded by a pair of hydrofractures has been solved.

In our initial treatments of the problem, we neglected variations in the hydraulic gas diffusivity as gas pressure changes. This made it possible to solve the problem with closed-form expressions, but the results did not match well with experimental data. The thermodynamic properties of natural gas must be treated properly, as recognized by, for example, Kelkar (7). Natural gas is not an ideal gas; its compressibility and viscosity depend upon its molar composition and vary strongly with temperature and pressure. Failing to take variable gas properties into account led to errors on the order of 50%.

We remark on some additional effects that we do not include. Injecting water into the gas-bearing rock leads to interactions between gas and water termed spontaneous imbibition. Although these effects are amenable to precise analysis (12) and have experimentally measurable consequences (13), we are able to neglect them because we discarded the first 3 mo of gas production, when these and other transient effects are most pronounced. In addition, as the pressure falls in a reservoir, gas adsorbed in the rock may escape, producing additional contributions to the gas flow. This process is described by the Langmuir desorption isotherm. In the particular field studied in this paper, we have carried out a detailed analysis and found this effect to be negligible, although it might not be so in other cases. Finally, desorption and flow of gas in unfractured shale have nonlinear properties at the microscopic level (14). We can neglect this phenomenon because gas transport is dominated by the effective properties of a fracture network, and the empirical evidence presented below shows that the net effect is pressure diffusion at an enhanced rate in a homogeneous medium.

Thus, we arrive at a specific nonlinear pressure diffusion problem to solve: It involves gas alone, permeability is uniform but enhanced in a cuboid volume, the experimental equation of state for gas is treated exactly, and spontaneous imbibition and desorption are neglected. The precise formulation and exact solution of the diffusion problem are contained in Methods, and we describe only the main results here.

Two planar hydrofractures in a well separated by distance 2d interfere with one another after a characteristic interference time (15), which we define aswhere is called the hydraulic diffusivity. It is related to the permeability of the rock k bywhere ? is porosity, Sg is the fraction of pore space occupied by gas, is the gas viscosity, and cg is the gas compressibility. In Eq. 1, τ is a constant defined at the initial state of the reservoir. It does not depend on the instantaneous gas pressure that varies in space and time as the reservoir is depleted. This does not mean that our solution relies on any approximation where quantities are fixed at reservoir values. It simply means that we adopt a time unit that is defined in terms of the initial reservoir properties; the final results take a particularly simple form when we do so.

We measure time in units of τ, defining a dimensionless time asNext, let m be the cumulative production of gas mass from a horizontal well with N hydrofractures (Fig. 1), and let ? be the original mass of gas contained in the reservoir volume drained by this well. The exact solution of the model for cumulative gas production is given by a dimensionless recovery factor (RF):We compute the RF by solving a particular boundary value problem (Methods) and plot it in Figs. 2 and 3. Although this RF is obtained from a numerical solution, it is a scaling function that, for practical purposes, provides the benefits of insight and convenience commonly associated with closed-form analytical solutions. To describe essentially all wells in the Barnett Shale, one has only to rescale this function in the time and gas production coordinates. In SI Text, we provide a spreadsheet (Dataset S1) in which the function is tabulated for convenience.

View larger version:
Fig. 2.

Cumulative production and production rate from scaling theory. (A) Dimensionless RF vs. dimensionless time computed from the scaling solution (black) compared with five typical wells (burnt orange). The fracture pressure pf is 500 psi, and the initial reservoir pressure pi is 3,500 psi. (B) Dimensionless well production rate vs. dimensionless time (black) under the same conditions compared with the same five typical wells (burnt orange). Production rates of individual wells are noisy, although cumulative production matches the scaling function well. Because the production rate becomes linear on a semilog plot, production decline is exponential for .

View larger version:
Fig. 3.

Comparison of 8,294 wells with scaling function. (A) Time history of 2,057 wells in the Barnett Shale, scaled so as to fit our scaling function (initial reservoir pressure of 3,500 psi and well flowing pressure of 500 psi), for which the dimensionless time starts below 0.25 and reaches 0.64 or more. The burnt orange curves give the scaled production of each well, and the black curve is the scaling function. Overall agreement is satisfactory. (B) Time history of 6,237 wells in the Barnett Shale for which the scaled maximum time comes out as (burnt orange). These wells are too young to trust our estimate of the interference time τ; therefore, we simply compare them with a square root function (black line). Time is scaled by the maximum time reached for each well, and production m is scaled by .

Our nearly universal solution to the boundary value problem depends, in principle, upon the initial state of the reservoir, , the well flowing pressure, pf, and gas composition, y, although it is independent of the details of the well geometry and the hydraulic diffusivity. In practice, pi, pf, T, and y can be set to typical values within a given shale gas play. For dimensionless times much less than 1, we show in SI Text that the solution takes a particularly simple intermediate asymptotic form:where is the dimensionless time necessary to extinguish the initial transients in gas flow.

The constant κ depends on the gas composition and temperature, as well as on the limiting pressures, pi and pf. For the wells we present here from the Barnett Shale, we set it equal to a typical value of 0.645. Table S1 shows that it varies rather little as the limiting pressures range over realistic values. Once the scaled time reaches 1, the growth in gas recovery slows, and it eventually reaches a plateau, which describes the maximum recovery possible for the given problem parameters. The way this slowing down occurs depends in detail upon the thermodynamics of gas expansion, the reservoir permeability, and the initial and final pressures in the reservoir. Eventually, as also shown in SI Text, production declines exponentially.

As a first illustration of Eqs. 4 and 5, suppose one knows the original gas in place, ?. After transients of the first few months of production have subsided, cumulative production takes the form . The constant is obtained by fitting a curve of this form to the measured cumulative production. Then,Therefore, to estimate the time τ after which well production declines exponentially, measure from the first year of production, estimate ? from the well geometry, and insert , and τ follows from Eq. 6.

However, the practical difficulty we face with gas production from the hydraulically fractured horizontal wells is greater than this example indicates. Neither the total mass of gas in place nor the time scale for interference to begin is known with any precision. The original mass of gas in place is uncertain, mainly because the effective hydrofracture length, 2L, and the number of active hydrofractures are uncertain. The time to interference is uncertain because the hydrofracturing process greatly increases the effective permeability k of the rock in the vicinity of the well; laboratory values of k obtained from core samples are on the order of nanodarcies (16), whereas accounting for observed well production requires effective values of k on the order of 100-fold greater.

Thus, we arrive at the following question: Can we extract enough information from existing field production data to estimate both the interference time τ and the original gas in place ? at the same time? In the early stages of gas production, when , the production rate declines purely as and τ and ? are impossible to determine separately. Wells delivering a small ultimate amount of gas at a relatively high rate cannot be distinguished from those where lower permeability rock or a small number of hydrofractures deliver ultimately larger quantities of gas at a relatively lower rate. Only the onset of interference between adjacent hydrofractures makes it possible to disentangle the two scenarios.

### Comparison with Field Data.

We display the dimensionless RF in Fig. 2. To illustrate its correspondence to data, we begin with a sample of 66 wells hand-selected by an experienced reservoir engineer as examples of good wells. In 5 of them, we find evidence of interference, meaning that cumulative production is not acceptably fit simply by . They do, however, fit the full scaling curve well, as we show with a graph of the cumulative production and production rate of these five wells in Fig. 2.

We then proceed to a more comprehensive study. We obtained data for 16,533 wells in the Barnett Shale, and from them, we selected the 8,807 horizontal wells that had operated continuously for 18 mo or more and had not been recompleted (i.e., the hydrofracturing process was not repeated to increase production). We allow ourselves only two fitting parameters on a per-well basis, horizontal and vertical scale factors, which correspond physically to the interference time, τ, and the original mass of gas in place, ?. Details of the fitting process are contained in SI Text. We find 2,057 horizontal wells for which the dimensionless time starts with a value less than 0.25 and reaches a value greater than 0.64. These are the wells for which interference is sufficiently advanced that it can be detected with an average uncertainty in parameters of less than 20%. We plot the RF of all these wells vs. scaled time and compare the results with the predicted scaling function in Fig. 5 and Fig. S3. Most of the wells show interference because the interference time τ is around 5 y, but a few of them have interference times of 10 y or more (Fig. 4). The fact that production from these more than 2,000 wells falls so well on the predicted curve provides evidence that the simple model we adopted is sufficiently realistic to estimate gas production in the future. We note that upper bounds on the total mass of gas in place are available from measurements of well geometry. As a check on our results, we show in Fig. S4 that the volumes of gas ? we calculate with our theory from production data are indeed less than these upper bounds.

View larger version:
Fig. 4.

Values of interference time τ and gas in place ? for the 2,057 wells in Fig. 3A. Error bars indicate two standard uncertainties. Maximum interference times here are around 10 y due to the fact that wells more than 10 y old are still rare; interference times of, say, 30 y will only be reliably detected when wells are 19 y old or more.

View larger version:
Fig. 5.

Upper and lower bounds on cumulative production from 8,294 wells in our sample. Vertical wells are excluded from the analysis, whereas twofold more wells will ultimately be drilled; thus, the upper bound is not an upper bound on the whole field.

We acknowledge that for any given well at particular points in time, production is noisy for many reasons (Fig. 2B). However, the cumulative production of individual wells falls remarkably well on our scaling curve (Fig. 2A), as does the expected production of thousands of wells (Fig. 3A).

There are an additional 6,237 wells for which interference is not yet visible, and which we say are in the square root decline phase. We cannot calculate τ and ? for these wells, but we can make use of our theory to put upper and lower bounds on them. We present these bounds in Fig. S5.

Summing up the production of the 8,294 wells in our sample, we obtain the lower and upper bounds on cumulative production over time, as shown in Fig. 5.

## Discussion

• i) We have found the minimal ingredients that suffice to model thousands of wells in the Barnett Shale with acceptable accuracy. The geometry of each well is a cuboid volume with a uniform array of absorbing boundaries. Between those boundaries, rock permeability is enhanced above laboratory values but is constant. Spontaneous imbibition can be neglected, but the gas equation of state must be treated realistically. Gas desorption is also negligible in the Barnett Shale but not elsewhere (e.g., in the Fayetteville shale). The scaling curve we find as a result provides surprisingly good agreement with all wells that can reasonably be analyzed in the Barnett Shale.

• ii) Inserting characteristic values into Eqs. 1 and 2, one deduces rock permeability k of 50 nanodarcies for τ of 50 y and 500 nanodarcies for τ of 5 y. These values of permeability are 20- to 200-fold larger than the values of a few nanodarcies found for shale core samples in laboratory experiments (16). This enhanced permeability must result from the hydrofracturing process. Many processes could be involved, including the reopening of preexisting fracture networks.

• iii) Cumulative gas production follows a nearly universal function scaled by two parameters, interference time τ and mass of gas in place ?.

• iv) For 2,057 of the horizontal wells in the Barnett Shale, interference is far enough advanced for us to verify that wells behave as predicted by the scaling form. The typical interference time in these wells is around 5 y.

• v) For 6,237 additional horizontal wells, no significant deviation from cumulative production growing as the square root of time is observed; these wells are too young to show evidence of interference. We provide upper and lower bounds on time to interference and original gas in place for each of these wells. The median lower bound on time to interference is 5 y, and the median upper bound is 100 y. The bounds on gas in place are somewhat tighter; the mean of the lower bounds is 1 billion standard cubic feet (Bscf), and the mean of the upper bounds is 7 Bscf. The lower bound on cumulative production from the wells we analyzed is 10 trillion standard cubic feet (Tscf) extracted over the next 10 y, whereas the upper bound is more than 20 Tscf that will continue to be recovered, at declining rates, over the next 50 y. By way of comparison, a recent estimate of the total gas production from all wells to be drilled in the Barnett Shale by 2050 is 40 Tscf (17, 18).

• vi) The contributions of shale gas to the US economy are so enormous (SI Text) that even small corrections to production estimates are of great practical significance.

Gas released by hydraulic fracturing can only be extracted from the finite volume where permeability is enhanced. Exponential decline of production once the interference time has been reached is inevitable, and extrapolations based upon the power law that prevails earlier are inaccurate. The majority of wells are too young to be displaying interference yet. The precise amount of gas they produce, and therefore their ultimate profitability, will depend upon when interference sets in.

For the moment, it is necessary to live with some uncertainty. Upper and lower bounds on gas in place are still far apart, even in the Barnett Shale with the longest history of production. Pessimists (4) see only the lower bounds, whereas optimists (19) look beyond the upper bounds. A detailed economic analysis based on the model presented here is possible, however, and is being published elsewhere (17, 18, 20, 21). The theoretical tools we are providing should make it possible to detect the onset of interference at the earliest possible date, provide increasingly accurate production forecasts as data become available, and assist with rational decisions about how hydraulic fracturing should proceed in light of its impact on the US environment and economy.

## Methods

We begin with an expression for mass balance of gas flowing in a porous rock:

where ug is the Darcy (superficial) velocity of gas, is gas saturation (with Swc being the connate water saturation), is the free gas density, is the adsorbed gas density (kilograms of gas per cubic meter of solid), and ? is the rock porosity.

By applying Darcy’s law to the linear, horizontal flow of gas, we can substitute

and obtain the following nonlinear partial differential equation:

The gas density is related to its pressure and temperature through an equation of state for real gases:

where is the compressibility factor of gas, Mg is the pseudomolecular mass of gas, J/kmol-K is the universal gas constant, and T is a constant temperature of the reservoir.

The isothermal compressibility of gas is defined as

We define as the differential equilibrium partitioning coefficient of gas at a constant temperature (e.g., ref. 22):

By inserting Eqs. 11 and 12 into Eq. 9, the general nonlinear equation of transient, linear, and horizontal flow of gas is obtained:

This nonlinear differential Eq. 13 can be simplified by introducing the Kirchhoff integral transform of gas pressure after Al-Hussainy et al. (23), which, in the present context, is also called “the real gas pseudopressure”:

Here, is a reference pressure that will be set to pf. After differentiation of Eq. 14 and cancelation of terms, one obtains the following nonlinear diffusion equation for gas pseudopressure:

with

The initial condition for Eq. 15 is

Note that mi is a constant only in a virgin reservoir. During refracturing, it will vary with the distance to the old hydrofractures.

We apply this equation to a finite region between two fractures, as shown in Fig. 1 (Lower):

where the hydrofracture pseudopressure, mf, might be a constant or a slow function of time. At the midpoint between two fractures, one has by symmetry

Eq. 15 is most useful in a scaled form. We define dimensionless time, distance, and pseudopressure by

Here, the subscript i refers to the quantities at the initial reservoir pressure pi and temperature T.

Consider the linear flow of gas into a transverse planar hydrofracture of height H and length 2L, and separated by distance 2d from the next hydrofracture planes, as depicted in Fig. 1. The scaled transport equation isand

Our approach is somewhat more general than that of Silin and Kneafsey (10) because we do not require any particular equation of state for natural gas and do not use the more limited p2 formulation (8). The and p2 solutions are equivalent only if is a linear function of pressure; however, generally, it is not (ref. 8, pp. 254–255). The price we pay is that our model must be solved numerically, but the cost is just a couple of seconds of delay before the full solution is computed on an average laptop. For the Barnett Shale, we use the values of well flowing pressure pf = 500 psi and initial reservoir pressure pi = 3,500 psi.

The superficial velocity of gas flowing into the right face of the hydrofracture at the origin is

The mass flow rate into this fracture is

Using Eq. 22,

Next, we replace the pressure with the real gas pseudopressure in Eq. 14:

The partial derivative can now, in turn, be rewritten with use of the scaled pseudopressure from Eq. 20, and the permeability, k, can be eliminated in favor of the gas diffusivity, , and the characteristic interference time, τ.

Let ? be the total mass of gas contained originally in the reservoir within the volume 4LHd between two consecutive hydrofractures, ; then the gas flow rate into the fracture plane at the origin takes the final form:

where the function depends only on gas composition, the initial and fracture pressures, and reservoir temperature. The scaling of has been devised so that this relation is exact.

The total flow into each pair of hydrofractures is twice that in Eq. 26. More generally, when there are N fracture stages, as depicted in Fig. 1, and we include the contribution to mass flow from the exterior faces of the left- and right-most hydrofractures, the original mass in place is

and the total mass transport out of the reservoir is given by

Here, we are treating the left- and right-most exterior hydrofracture faces approximately, as extensions of the wellbore length by d at each end. The reason it is not appropriate to treat the two ends as semiinfinite is that without the great enhancement of permeability brought about by the hydrofracturing process, gas transport is negligible. Our assumption is that volumetric rock damage extends beyond the ends of the two last fractures for characteristic distance d.

Integrating Eq. 28 with respect to the dimensionless time gives the final result

The initial boundary value problem (Eq. 21) is solved numerically with an efficient fully implicit solver and a sequential implicit solver. The first solver has been implemented in Python, and the second has been implemented in MATLAB (MathWorks). Accurate numerical solutions can be obtained in both cases within a few seconds on an average laptop. Essential properties of the result are revealed by exact solution of simplified equations, depicted in Fig. S6.

## Acknowledgments

We thank John Browning for help in estimating reserves, for his deep insights into well performance in the Barnett Shale, and for help in selection of well-behaved groups of wells. We thank D. Silin, S. Bhattacharaya, and R. Dombrowski for detailed comments on the manuscript. The gas production data were extracted from the IHS Cambridge Energy Research Associates database, licensed to the Bureau of Economic Geology. This paper was supported by the Shell Oil Company/University of Texas at Austin project “Physics of Hydrocarbon Recovery,” with T.W.P. and M.M. as coprincipal investigators, and the Bureau of Economic Geology’s Sloan Foundation-funded project “The Role of Shale Gas in the U.S. Energy Transition: Recoverable Resources, Production Rates, and Implications.” M.M. acknowledges partial support from the National Science Foundation Condensed Matter and Materials Theory program.

## Footnotes

• ?1To whom correspondence should be addressed. E-mail: patzek{at}mail.utexas.edu.
• Author contributions: T.W.P., F.M., and M.M. designed research, performed research, analyzed data, and wrote the paper.

• The authors declare no conflict of interest.

• See Commentary on page 19660.

Freely available online through the PNAS open access option.

## References

1. ?
2. ?
3. ?
4. ?
5. ?
6. ?
7. ?
8. ?
9. ?
10. ?
11. ?
12. ?
13. ?
14. ?
15. ?
16. ?
17. ?
18. ?
19. ?
20. ?
21. ?
22. ?
23. ?

#### Online Impact

• 956115858 2018-01-22
• 730379857 2018-01-22
• 346624856 2018-01-22
• 201609855 2018-01-22
• 72549854 2018-01-21
• 795928853 2018-01-21
• 752345852 2018-01-21
• 566508851 2018-01-21
• 615722850 2018-01-21
• 689612849 2018-01-21
• 846903848 2018-01-21
• 674896847 2018-01-21
• 11197846 2018-01-21
• 986896845 2018-01-21
• 667601844 2018-01-21
• 385442843 2018-01-21
• 496686842 2018-01-21
• 915288841 2018-01-21
• 885256840 2018-01-21
• 726268839 2018-01-21