Abstract:
Apparatus and methods for delivering and placing proppant to a subterranean formation fracture including identifying control variables and uncertain parameters of the proppant delivery and placement, optimizing a performance metric of the proppant delivery and placement under uncertainty, calculating sensitivity indices and ranking parameters according to a relative contribution in total variance for an optimized control variable, and updating a probability distribution for parameters, repeating optimizing comprising the probability distribution, and evaluating a risk profile of the optimized performance metric using a processor. Some embodiments may deliver proppant to the fracture using updated optimized values of control variables.

Description:
FIELD 
       [0001]    Embodiments herein relate to hydraulic fracturing including proppant placement. 
       BACKGROUND 
       [0002]    A standard approach to optimization under uncertainty is based on original Markovitz portfolio theory and more recently was tailored to oilfield applications with modified definition of efficient frontier (U.S. Pat. No. 6,775,578 B. Couet, R. Burridge, D. Wilkinson, Optimization of Oil Well Production with Deference to Reservoir and Financial Uncertainty, 2004) and Value of Information (Raghuraman, B., Couët, B., Savundararaj, P., Bailey, W. J. and Wilkinson, D.: “Valuation of Technology and Information for Reservoir Risk Management,” paper SPE 86568, SPE Reservoir Engineering, 6, No. 5, October 2003, pp. 307-316). However, these methods employ mean-variance approach and do not provide a much needed insight into the inherent uncertainty of the optimized model and, more importantly, any quantitative guidance on reducing this uncertainty, which is very desirable from the operational point of view. 
         [0003]    Application of Global Sensitivity Analysis to address various problems arising in oilfield industry has been described for reservoir performance evaluation, for measurement screening under uncertainty, for pressure transient test design and interpretation, for design and analysis of miscible fluid sampling clean-up, and for targeted survey design. However, these disclosures were focusing only on quantifying uncertainty in specific physical quantities and using that analysis to gain a new insight about the measurement program design and interpretation. The references did not look at optimization of the underlying physical processes. 
     
    
     
       BRIEF DESCRIPTION OF THE FIGURES 
         [0004]      FIG. 1  is a workflow summarizing adaptive GSA-optimization approach. 
           [0005]      FIG. 2  is a workflow summarizing the inputs and outputs for the example proppant placement and fracture conductivity calculation. 
           [0006]      FIG. 3  is a schematic diagram providing a definition of cycle phase shift and perforation spacing for two injectors from a vertical well into a vertical fracture. 
           [0007]      FIG. 4  is a schematic diagram illustrating the length of the cycle and length of the proppant-laden portion. 
           [0008]      FIG. 5  is a schematic diagram for one example considered. The final placed distribution of proppant is also influenced by mixing between the proppant-laden and clean fracturing fluid portions. The mixing process is characterized by a single mixing length. 
           [0009]      FIG. 6  is a workflow illustrating the inputs, outputs, and workflow for the example proppant placement and fracture conductivity calculation. 
           [0010]      FIG. 7  is a chart plotting three points of the efficient frontier from the optimization using initial ranges for uncertain variables. Lower values of the objective function (μ−λσ) for increasing values of λ illustrate the inherent penalty for risk. 
           [0011]      FIG. 8  is a chart plotting three points of the efficient frontier from the optimization using GSA-updated ranges for uncertain variables. Initial efficient frontier points are also included for comparison. Lower values of the objective function (μ−λσ) for increasing values of λ illustrate the inherent penalty for risk. 
       
    
    
     SUMMARY 
       [0012]    Embodiments herein relate to apparatus and methods for delivering and placing proppant to a subterranean formation fracture including identifying control variables and uncertain parameters of the proppant delivery and placement, optimizing a performance metric of the proppant delivery and placement under uncertainty, calculating sensitivity indices and ranking parameters according to a relative contribution in total variance for an optimized control variable, and updating a probability distribution for parameters, repeating optimizing comprising the updated probability distribution, and evaluating a risk profile of the optimized performance metric using a processor. Some embodiments may deliver proppant to the fracture using updated optimized values of control variables. 
       DETAILED DESCRIPTION 
       [0013]    This disclosed approach combines Global Sensitivity Analysis (Saltelli et al., 2004) with optimization under uncertainty in an adaptive workflow that results in guided uncertainty reduction of the optimized model predictions. Embodiments herein relate to a general area of optimization under uncertainty. The application of the disclosed method relates to well stimulation and hydraulic fracturing in particular. Heterogenous Proppant Placement (HPP) strategies seek to increase propped fracture conductivity by selectively placing the proppant such that the fracture is held open at discrete locations and the reservoir fluids can be transported through open channels between the proppant. Schlumberger Technology Corporation provides well services that include introducing proppant into the fractures in discrete slugs (Gillard, M. et al., 2010; Medvedev, A. et al., 2013). For the purposes of technology development and optimal implementation, tools must be developed for predicting the conductivity of the heterogeneously propped fractures during the increase in closure stress resulting from flow-back and subsequent production. In the presence of uncertainty in formation properties, optimal HPP strategies will result in inherently uncertain predictions of fracture conductivity. Herein, we describe a method to reduce uncertainty in predicted fracture conductivity and identify an optimal HPP operational strategy for an acceptable level of risk. 
         [0014]    Embodiments herein show how a predictive physics-based HPP model is used to estimate fracture conductivity under a given closure stress. The input parameters of the model are divided into control variables (operational controls may include dirty pulse fraction, injector spacing, proppant Young&#39;s modulus etc.) and uncertain variables (uncertain formation properties may include Poison ratio, Young&#39;s modulus, proppant diffusion rates etc.). The model is first optimized to obtain values of control variables maximizing mean fracture conductivity (for a given closure stress) under initial uncertainty of formation properties. An efficient frontier may be obtained at this step to characterize dependence between the optimized mean value of fracture conductivity and its uncertainty expressed by the standard deviation. Global sensitivity analysis (GSA) is then applied to quantify and rank contributions from uncertain input parameters to the standard deviation of the optimized values of fracture conductivity. Uncertain parameters are ranked according to their calculated sensitivity indices and additional measurements can be performed to reduce uncertainty in the high-ranking parameters. Constrained optimization of the model with reduced ranges of uncertain parameters is performed and a new efficient frontier is obtained. In most cases, the points of the updated efficient frontier will shift to the left indicating a reduction in the risk associated with achieving the desired fracture conductivity. The disclosed method provides an adaptive GSA-optimization approach that results in uncertainty reduction for optimized HPP performance. 
         [0015]    The workflow is applied for HPP optimization, which requires a capability for the prediction of the placement of proppant and the resultant conductivity within a potentially rough fracture under any prescribed closure stress. This capability receives inputs relating to the pumping schedule, proppant properties and formation properties and provides a prediction of the achieved fracture conductivity. For example, in our demonstration, we utilize the methods in U.S. Provisional Patent Application Ser. No. 61/870,901, filed Aug. 28, 2013 which is incorporated by reference herein in its entirety where the combination of fracture and proppant is represented by a collection of asperities arranged upon a regular grid attached to two deformable half-spaces. The deformation characteristics of the deformable half-spaces are pre-calculated, allowing for very efficient prediction of the deformation of the formation on either side of the fracture. The method automatically detects additional contact as the fracture closes during increasing closure stress (such as during flow-back and production). In addition, the asperity mechanical response is modified to account for the combined mechanical response of the rough fracture surface and any proppant that may be present in the fracture at that location. In this way, the deformation of any combination of fracture roughness and heterogeneous arrangement of proppant in the fracture can be evaluated. The deformed state is then converted into a pore network model which calculates the conductivity of the fracture during flow-back and subsequent production. Embodiments herein allow one to progressively reduce uncertainty in the performance of an optimized HPP operational strategy by iterative reduction of uncertainty in identified properties of the reservoir. 
       Optimization Under Uncertainty and Global Sensitivity Analysis 
       [0016]    Let us consider a general case when the underlying physical process is modeled by a function y=f(α, β), where α={α 1  . . . α N } and β={β 1  . . . β M } are two sets of parameters. Here, α represents the set of control parameters (to be used in optimization), and β denotes the set of uncertain parameters. Mathematically, β&#39;s are considered to be random variables represented by a joint probability density function (pdf). Therefore, for each vector of control variables α, the output of the model is itself a random variable with its own pdf (due to uncertainty in β). A mean-variance approach is commonly used for optimization, i.e. a function of the form 
         [0000]        F =μ(α,β)−λσ(α,β)
 
         [0000]    where μ, and σ are the mean and standard deviation of the output y of the numerical simulation, and λ is a non-negative parameter defining a tolerance to risk (uncertainty). The optimization problem can then be formulated as 
         [0000]    
       
         
           
             
               max 
               α 
             
              
             
                 
             
              
             
               F 
                
               
                 ( 
                 
                   α 
                   , 
                   β 
                 
                 ) 
               
             
           
         
       
     
         [0017]    For each optimization iteration, a sample of the random vector β is chosen, and the values of y(α, β) are first computed using this sample for a given a and then averaged over β. 
         [0018]    Various optimization algorithms can then be used to find the optimal value of α. The process of optimizing under uncertainty will lead to a set of parameters α opt  that provide the optimum of the objective function F. Therefore, an optimized model is now available: 
         [0000]        y=f (α opt ,β)
 
         [0000]    Note that the optimized model still has inherent uncertainty due to the uncertainty in parameters β. 
         [0019]    A set of solutions to the optimization problem can be plotted in (μ, σ) coordinates, where optimal points corresponding to pre-defined values of λ will form an efficient frontier ( FIG. 7 ). This represents a risk profile of the underlying modeled process. The positive slope of the frontier illustrates the penalty for additional uncertainty (risk). 
         [0020]    From the operational perspective, the goal is to reduce this risk while maintaining the same level of expected performance (represented by μ). In order to reduce the uncertainty, one needs to understand where it is coming from. Therefore, a quantitative link between uncertainties in input parameters (β) and uncertainty in the output is desirable. This link can be quantified using Global Sensitivity Analysis based on variance decomposition. 
         [0021]    Global sensitivity analysis (Saltelli et al., 2004) based on variance decomposition is used to calculate and apportion the contributions to the variance of the measurement signal V(Y) from the uncertain input parameters {X i } of the subsurface model. 
         [0022]    For independent {X i }, the Sobol&#39; variance decomposition (Sobol&#39;, 1993) can be used to represent V(Y) as 
         [0000]        V ( Y )=Σ i=1   N   V   i +Σ 1≦i&lt;j≦N   V   ij   + . . . +V   12 . . . N ,  (1)
 
         [0000]    where V i =V[E(Y|X i )] are the variance in conditional expectations (E) representing first-order contributions to the total variance V(Y) when X i  is fixed i.e., V(X i )=0. Since we do not know the true value of X i  a priori, we have to estimate the expected value of Y when X i  is fixed anywhere within its possible range, while the rest of the input parameters {X ˜i } are varied according to their original probability distributions. Thus, 
         [0000]        S 1 i   =V   i   /V ( Y ) 
         [0000]    is an estimate of relative reduction in total variance of Y if the variance in X i  is reduced to zero. 
         [0023]    Similarly, V ij =V[E(Y|X i , X j )]−V i −V j  is the second-order contribution to the total variance V(Y) due to interaction between X i  and X j . Notice, that the estimate of variance V[E(Y|X i , X j )] when both X i  and X j  are fixed simultaneously should be corrected for individual contributions V i  and V j . 
         [0024]    For additive models Y(X), the sum of all first-order effects S1 i  is equal to 1. This is not applicable for the general case of non-additive models, where second, third and higher-order effects (i.e., interactions between two, three or more input parameters) play an important role. The contribution due to higher-order effects can be estimated via total sensitivity index ST: 
         [0000]        ST   i   ={V ( Y )− V[E ( Y|X   ˜i )]}/ V ( Y ),
 
         [0000]    where V(Y)−V[E(Y|X ˜i )] is the total variance contribution from all terms in Eq. 1 that include X i . Obviously, ST i ≧S1 i , and the difference between the two represents the contribution from the higher-order interaction effects that include X i . 
         [0025]    There are several methods available to estimate S1 i  and ST i  (see (Saltelli et al., 2008) for a comprehensive review). 
         [0026]    In one embodiment, we apply Polynomial Chaos Expansion (PCE) [Wiener, 1938] to approximate the underlying optimized function y=f(α opt ,β). The advantage of applying PCE is that all GSA sensitivity indices can be calculated explicitly once the projection on the orthogonal polynomial basis is computed (Sudret, 2008). 
         [0027]    In another embodiment, GSA sensitivity indices can be calculated using an algorithm developed by Saltelli (2002) that further extends a computational approach proposed by Sobol&#39; (1990) and Homma and Saltelli (1996). The computational cost of calculating both S1 i  and ST i  is N(k+2), where k is a number of input parameters {X i } and N is a large enough number of model calls (typically between 1000 and 10000) to obtain an accurate estimate of conditional means and variances. However, with underlying physical model taking up to several hours to run, this computational cost can be prohibitively high. Therefore, we can use proxy-models that approximate computationally expensive original simulators. Quasi-random sampling strategies such as LPτ sequences (Sobol, 1990) can be employed to improve the statistical estimates of the computed GSA indices. 
         [0028]    Once sensitivity indices are computed, uncertain β-parameters can be ranked according to values of S1. Parameters with the highest values of S1 should be selected for targeted measurement program. Reduction in uncertainty of these parameters will result in largest reduction in uncertainty of predicted model outcome. Parameters with lowest values of ST (typically, below 0.05) can be fixed at their base case value, thus reducing dimensionality of the underlying problem and improving the computational cost of the analysis. 
         [0029]    The summary of the proposed general workflow is given in  FIG. 1 . 
         [0030]    The main steps include:
       1. Identify control variables (α) and uncertain parameters (β). If applicable, define ranges for control variables. Define probability distribution functions (pdfs) for uncertain parameters.   2. Perform optimization under uncertainty (max F (α, β), where F=μ−λσ) and construct relevant points on the efficient frontier for various values of λ.   3. For a given point on the efficient frontier (defined by prescribed value of λ and corresponding values of control parameters α λ ), calculate GSA sensitivity indices and rank uncertain parameters β according to values of S1.   4. Perform additional measurements to reduce uncertainty of parameters β with high values of S1 and redefine pdfs for those parameters.   5. Optional: fix values of parameters β with low (e.g., below 0.05) values of ST to reduce dimensionality of the optimization problem.   6. Perform steps 2-5 until acceptable level of risk is achieved or until the decision is made that the desired level of performance cannot be achieved with the acceptable level of risk.       
 
       Illustrative Example 
     HPP Optimization Under Uncertainty 
       [0037]    We now describe the application to a problem of HPP optimization demonstrating the method. 
         [0038]    The underlying physical model along with the methods and numerical tools developed to simulate it are disclosed in “Method for Predicting Heterogeneous Proppant Placement and Conductivity” (U.S. Provisional Patent Application Ser. No. 61/870,901, filed Aug. 28, 2013 which is incorporated by reference above). Below we provide a short description of the main steps involved in calculating fracture conductivity resulting from HPP. 
         [0039]      FIG. 2  shows a flow diagram highlighting the inputs and outputs utilized in our specific example of fracture conductivity when considering the heterogeneous placement of proppant from a vertical well intersecting a vertical hydraulic fracture as depicted in  FIG. 3 . In this instance, the heterogeneity of proppant in the fracture is achieved through a combination of pulsing of the proppant into the fracture (see  FIG. 4 ) and mixing phenomena (see  FIG. 5 ) that are characterized by a mixing length. 
         [0040]      FIG. 6  shows in more detail how the inputs are broken down into those related to the placement of the proppant and those related to the subsequent deformation and conductivity calculations. The complete list of model inputs utilized by the example application is provided in Table 1 along with descriptions of the inputs, their units and initial ranges used in this example. 
         [0041]    We start by following the steps of the workflow disclosed in  FIG. 1 . 
       Step 1. 
       [0042]    Identify control variables (α) and uncertain parameters (β) and define their ranges and probability distributions. 
         [0000]    Control variables (α) include:
       1. Injector spacing   2. Pumping rate   3. Full cycle length   4. Proppant pulse length   5. Injector phase shift   6. Proppant Young&#39;s modulus   7. Proppant permeability (parameter was fixed in this study since the dominating flow mechanism in successful HPP job should be though the channels formed between the proppant pillars rather than through the proppant itself)       
 
         [0050]    Ranges for control variables are given in Table 1.  FIG. 4  illustrates some of these variables related to heterogeneous placement of proppant and consequently some systems can accommodate a pumping schedule that includes variations in proppant concentration with time. 
         [0000]    Uncertain variables (β) include: 
         [0051]    1. Fracture aperture during placement 
         [0052]    2. Proppant mixing length 
         [0053]    3. Formation Young&#39;s modulus 
         [0054]    4. Formation Poisson ratio 
         [0055]    Ranges for uncertain variables are given in Table 1. All variables were assumed to be uniformly distributed, except for “Proppant mixing length” that was assumed to be uniformly distributed on a log scale. 
         [0000]    
       
         
               
             
               
               
               
               
             
           
               
                 TABLE 1 
               
             
             
               
                   
               
               
                 List of inputs for fracture conductivity calculation applied to injection 
               
               
                 into a vertically oriented fracture from a vertical well. 
               
             
          
           
               
                 Input 
                 Description 
                 Units 
                 Range 
               
               
                   
               
               
                 Injector 
                 Vertical distance between 
                 Length (m) 
                 0.5-3     
               
               
                 spacing 
                 injectors 
               
               
                 Pumping 
                   
                 Volume per 
                 0.1-0.5 
               
               
                 rate 
                   
                 unit time (bpm) 
               
               
                 Full cycle 
                 Length in time of repeated 
                 Time (s) 
                 15-25 
               
               
                 length 
                 cycle of heterogeneous 
               
               
                   
                 injection 
               
               
                 Proppant 
                 Fraction of total injection 
                 Non- 
                 0.25-0.75 
               
               
                 pulse length 
                 period dedicated to 
                 dimensional 
               
               
                   
                 proppant injection 
               
               
                 Injector 
                 The systematic delay be- 
                 Non- 
                 0-1 
               
               
                 phase shift 
                 tween the cycles of the 
                 dimensional 
               
               
                   
                 injectors (as fraction of 
               
               
                   
                 total cycle length) 
               
               
                 Fracture 
                 Fracture assumed to have 
                 Length (mm) 
                 3-7 
               
               
                 aperture 
                 constant aperture during 
               
               
                 during 
                 displacement for this 
               
               
                 placement 
                 demonstration. 
               
               
                 Proppant 
                 The permeability of the 
                 Length*Length 
                 fixed at 
               
               
                 permeability 
                 permeability can be stress 
                 (m 2 ) 
                 10 −10   
               
               
                 under stress 
                 dependent. In this 
               
               
                   
                 demonstration it was 
               
               
                   
                 assumed constant. 
               
               
                 Proppant 
                 Characteristic length scale 
                 Length (m) 
                 0.001-0.25  
               
               
                 mixing 
                 over which proppant and 
               
               
                 length 
                 clean fracturing fluid mix 
               
               
                   
                 during placement 
               
               
                 Proppant 
                 Assumed elastic constant 
                 Stress (MPa) 
                  50-500 
               
               
                 Young&#39;s 
                 characterizing compression 
               
               
                 modulus 
                 of proppant. 
               
               
                 Formation 
                   
                 Stress (GPa) 
                  5-50 
               
               
                 Young&#39;s 
               
               
                 modulus 
               
               
                 Formation 
                   
                 Non- 
                 0.15-0.35 
               
               
                 Poisson 
                   
                 dimensional 
               
               
                 ratio 
               
               
                 Closure 
                   
                 Stress (MPa) 
                 0.1-30  
               
               
                 stress levels 
               
               
                   
               
             
          
         
       
     
       Step 2. 
       [0056]    Perform optimization under uncertainty (max F (α, β), where F=μ−λσ) and construct relevant points on the efficient frontier for various values of λ. 
         [0057]    The underlying quantity to be optimized is fracture conductivity at a predefined closure stress (20 MPa in this example). In general, the objective function can be based on other performance metrics of proppant delivery and placement in the fracture including total hydrocarbon produced through the fracture, hydrocarbon production rate, and a financial indicator characterizing profitability of the fracturing job. Results of the optimization step comprise a risk profile shown in  FIG. 7 . Corresponding values of mean, standard deviation for three λ points along with P10-P50-P90 estimates for facture conductivity corresponding to these three operational scenarios are given in Table 2. 
         [0000]    
       
         
               
             
               
               
               
               
             
               
               
               
               
             
           
               
                 TABLE 2 
               
             
             
               
                   
               
               
                 Results of optimization with initial uncertainty. 
               
             
          
           
               
                   
                 λ = 0 
                 λ = 1 
                 λ = 2 
               
               
                   
                   
               
             
          
           
               
                 Mean fracture conductivity (D · m) 
                 248.15 
                 232.05 
                 193.6 
               
               
                 Mean fracture conductivity (log 10) 
                 −9.605 
                 −9.634 
                 −9.713 
               
               
                 Standard deviation (log10 cycles) 
                 1.21 
                 1.15 
                 1.10 
               
               
                 P90 (D · m) 
                 2.21 
                 2.79 
                 2.83 
               
               
                 P50 (D · m) 
                 562 
                 507 
                 408 
               
               
                 P10 (D · m) 
                 4582 
                 3619 
                 2622 
               
               
                   
               
             
          
         
       
     
       Step 3. 
       [0058]    For a given point on the efficient frontier (defined by prescribed value of λ and corresponding values of control parameters α λ ), calculate GSA sensitivity indices and rank uncertain parameters β according to values of S1. 
         [0059]    We apply Polynomial Chaos Expansion approach to calculate GSA sensitivity indices for optimized models corresponding to values λ=0, 1, 2. The values for first-order sensitivity index (S1) and total sensitivity index (ST) for each uncertain parameter β are given in Table 3. For all three optimal points on the efficient frontier, “Proppant mixing length” is responsible for almost 70% of variance in predicted fracture conductivity. The second largest contributor is “Fracture aperture during placement” (15-20% of variance). 
         [0000]    
       
         
               
             
               
               
               
               
             
               
               
               
               
               
               
               
             
               
               
               
               
               
               
               
             
           
               
                 TABLE 3 
               
             
             
               
                   
               
               
                 GSA sensitivity indices for optimized models (uncertain 
               
               
                 parameters ranked according to S1). 
               
             
          
           
               
                   
                 λ = 0 
                 λ = 1 
                 λ = 2 
               
             
          
           
               
                   
                 S1 
                 ST 
                 S1 
                 ST 
                 S1 
                 ST 
               
               
                   
                   
               
             
          
           
               
                 Proppant mixing length 
                 0.72 
                 0.75 
                 0.69 
                 0.73 
                 0.65 
                 0.70 
               
               
                 Fracture aperture during 
                 0.17 
                 0.18 
                 0.17 
                 0.18 
                 0.19 
                 0.20 
               
               
                 placement 
               
               
                 Formation Young&#39;s modulus 
                 0.08 
                 0.11 
                 0.09 
                 0.13 
                 0.11 
                 0.15 
               
               
                 Formation Poisson ratio 
                 0.00 
                 0.00 
                 0.00 
                 0.00 
                 0.00 
                 0.00 
               
               
                   
               
             
          
         
       
     
       Step 4. 
       [0060]    Perform additional measurements to reduce uncertainty of parameters β with high values of S1 and redefine pdfs for those parameters. 
         [0061]    Based on results of Step 3, “Proppant mixing length” was identified as a single largest contributor to variance of fracture conductivity at 20 MPa. For illustration, we assume that additional measurements were performed to reduce the uncertainty range of this parameter from 0.001 m-0.25 m (slightly more than two log 10 cycles) to 0.005-0.05 (one log 10 cycle) with uniform distribution on log scale. 
       Step 5. 
       [0062]    Optional: fix values of parameters β with low (&lt;0.05) values of ST to reduce dimensionality of the optimization problem. 
         [0063]    Analyzing total-sensitivity values, we notice that “Formation Poisson ratio” has values very close to zero. Therefore, fixing this parameter in the middle of its original uncertainty range (0.15-0.35) will not significantly affect the outcome of the subsequent analysis (Sobol, 2001) while improving its computational cost since the dimensionality of the problem will be reduced. 
       Step 6. 
       [0064]    Perform optimization step 2 with updated ranges of uncertain parameters. 
         [0065]    Results of the optimization step are shown in  FIG. 8 . Three points of the initial efficient frontier are also included for comparison. The updated efficient frontier has moved to the left (desired reduction in uncertainty) and slightly up. We note that the vertical direction of the shift in efficient frontier depends on underlying values in the physical quantity of interest (fracture conductivity) in the updated range of the uncertain parameter (Proppant mixing length). Corresponding values of mean, standard deviation for three λ points along with updated P10-P50-P90 estimates for facture conductivity corresponding to these three operational scenarios are given in Table 4. We observe the significant reduction in standard deviation (on log scale) compared to the initial case. The reduction in P10-P90 range on a linear scale is also noticeable. 
         [0000]    
       
         
               
             
               
               
               
               
             
               
               
               
               
             
           
               
                 TABLE 4 
               
             
             
               
                   
               
               
                 Results of optimization with updated 
               
               
                 uncertainty ranges (based on GSA). 
               
             
          
           
               
                   
                 λ = 0 
                 λ = 1 
                 λ = 2 
               
               
                   
                   
               
             
          
           
               
                 Mean fracture conductivity (D · m) 
                 743.39 
                 698.31 
                 589.72 
               
               
                 Mean fracture conductivity (log 10) 
                 −9.129 
                 −9.156 
                 −9.229 
               
               
                 Standard deviation (log10 cycles) 
                 0.68 
                 0.59 
                 0.53 
               
               
                 P90 (D · m) 
                 81 
                 109 
                 103 
               
               
                 P50 (D · m) 
                 981 
                 943 
                 782 
               
               
                 P10 (D · m) 
                 4494 
                 3010 
                 2219 
               
               
                   
               
             
          
         
       
     
         [0066]    The shift of efficient frontier to the left is expected in most cases. With the rare exception when the local variance underlying of values in the physical quantity of interest in the updated range of the uncertain parameter is higher than that in the initial range. Although even for this exception case, we argue that the disclosed approach provides iterative way to accurately estimate risk-reward profile for a given HPP job and allows one to avoid costly mistakes that would result in an underperforming fracture. 
         [0067]    We disclosed a method for adaptive optimization of heterogeneous proppant placement under uncertainty. A predictive physics-based HPP model is used to estimate fracture conductivity under the desired closure stress. The input parameters of the model are divided into control variables and uncertain variables. The model is first optimized to obtain values of control variables maximizing mean fracture conductivity (at given closure stress) under initial uncertainty of formation properties. An efficient frontier may be obtained at this step to characterize the dependence between the optimized mean value of fracture conductivity and its uncertainty expressed by the standard deviation. Global sensitivity analysis is then applied to quantify and rank contributions from the uncertain input parameters to the standard deviation of the optimized values of fracture conductivity. The uncertain parameters are ranked according to their calculated sensitivity indices and additional measurements can be performed to reduce uncertainty in the high-ranking parameters. Constrained optimization of the model with reduced ranges of uncertain parameters is performed and a new efficiency frontier is obtained. In most cases, the points of the updated efficient frontier will shift to the left indicating reduction in risk associated with achieving the desired fracture conductivity. The disclosed method provides an adaptive GSA-optimization approach that results in iterative improvement of estimated risk-reward profile of an optimized HPP job under uncertainty. 
         [0068]    Some embodiments may use a computer system including a computer processor (e.g., a microprocessor, microcontroller, digital signal processor or general purpose computer) for executing any of the methods and processes described herein. The computer system may further include a memory such as a semiconductor memory device (e.g., a RAM, ROM, PROM, EEPROM, or Flash-Programmable RAM), a magnetic memory device (e.g., a diskette or fixed disk), an optical memory device (e.g., a CD-ROM), a PC card (e.g., PCMCIA card), or other memory device. The memory can be used to store computer instructions (e.g., computer program code) that are interpreted and executed by the processor.