A Four Dimensional Variational Data Assimilation Framework for Wind Energy Potential Estimation

In this paper, we propose a Four-Dimensional Variational (4D-Var) data assimilation framework for wind energy potential estimation. The framework is defined as follows: we choose a numerical model which can provide forecasts of wind speeds then, an ensemble of model realizations is employed to build control spaces at observation steps via a modified Cholesky decomposition. These control spaces are utilized to estimate initial analysis increments and to avoid the intrinsic use of adjoint models in the 4D-Var context. The initial analysis increments are mapped back onto the model domain from which we obtain an estimate of the initial analysis ensemble. This ensemble is propagated in time to approximate the optimal analysis trajectory. Wind components are post-processed to get wind speeds and to estimate wind energy capacities. A matrix-free analysis step is derived from avoiding the direct inversion of covariance matrices during assimilation cycles. Numerical simulations are employed to illustrate how our proposed framework can be employed in operational scenarios. A catalogue of twelve Wind Turbine Generators (WTGs) is utilized during the experiments. The results reveal that our proposed framework can properly estimate wind energy potential capacities for all wind turbines within reasonable accuracies (in terms of Root-Mean-Square-Error) and even more, these estimations are better than those of traditional 4D-Var ensemble-based methods. Moreover, large variability (variance of standard deviations) of errors are evidenced in forecasts of wind turbines with the largest rate-capacity while homogeneous variability can be seen in wind turbines with the lowest rate-capacity.


Introduction
Data Assimilation (DA) is the process by which an imperfect numerical forecast {x b k } G k=0 is adjusted according to real noisy observations {y k } G k=0 [1,2], where x b k ∈ R n×1 and y k ∈ R m×1 are the background state and the observations at step k, for 0 ≤ k ≤ G, respectively, n is the model size (model resolution), m denotes the number of observations per assimilation step, and G is the size of the assimilation window (the number of times wherein observations are available). In strong constraint Four-Dimensional Variational (4D-Var) methods, cost functions of the form [3,4]: are employed to perform the assimilation process, where B 0 ∈ R n×n and R k ∈ R m×m are the covariance matrix of the initial background errors and the estimated data error covariance matrix at step k, respectively. Likewise, where M : R n×1 → R n×1 is a numerical model which, for instance, mimics the behavior of the ocean and/or the atmosphere and even more, we assume that the model can predict wind components (or wind speeds). We then seek the initial condition which best fit the data: J (x 0 ) , subject to (2) .
To solve the optimization problem (3), for instance, we can make use of adjoint models (i.e., transpose of linearization of the numerical model) or an ensemble of model realizations. Regardless of which one is chosen, the initial condition (3) can provide a forecast of relevant physical variables (depending on the numerical model) such as wind components, temperature, and humidity. Forecasts of wind components can be exploited in the context of Wind Turbine Generators (WTGs) to estimate the potential energy capacities of wind turbines and then to employ green sources of energy in cities, countries, and even more, remote places wherein these are the unique option. Hence, we can make use of numerical models and 4D-Var optimization problems to estimate wind components (in particular, wind speeds since they are key inputs to size WTG), then by using parameters from wind turbines such as cut-in speed, cut-out speed, and rated wind speed, we can estimate their potential wind power capacity for a specific place. WTG parameters allow choosing the best WTG for a specific place according to its wind-speed values. This can be exploited in places such as the Latin American and the Caribbean (LAC) countries since these are widely known for their large power generation capacity using renewable energies. This makes them highly attractive for clean energy investment [5]. Recent studies indicate that the full deployment of this capacity can be almost seven times larger than the current world installed one, and even more, it can constitute a near-zero carbon emissions option for developing countries [6]. This could provide substantial societal benefits, including energy security, local and global environmental benefits, domestic job creation, and improved balance of payments, amongst others [7]. Given LAC geography, wind turbines (wind farms) can be exploited as clean energy sources. The economic benefits of wind farms are better than those of traditional sources such as solar farms, which make the first option desirable for long-term plans at national levels. The proper planning and scheduling of wind power systems can lead to almost no impact on LAC ecosystems, neither visual or audible. Besides, this can serve as a complement of the hydro-dominated electricity grids [8], as the winds are stronger during the dry season when hydroelectric generation is most limited. Moreover, wind farms can provide a full or complementary source of energy in some areas of difficult access; the application of wind turbines is primarily in windmills that are used to generate electricity [9]. These wind turbines can be used to avail off-grid electricity in remote regions (i.e., some islands).
The structure of this paper is as follows: Section 2 discusses DA formulations and wind turbine generators (WTG). In Section 3, we propose a novel framework for electrical power estimation via WTGs and 4D-Var ensemble DA. In Section 4, numerical simulations are employed to show how our framework can be employed; we generate synthetic scenarios by using an Atmospheric General Circulation Model. Lastly, Section 5 states the conclusions of this research.

Preliminaries
In this section, we briefly discuss concepts related to 4D-Var ensemble DA and wind turbine generators. We primarily focus on the necessary topics for the derivation of our 4D-Var ensemble-based method.

Data Assimilation
To solve the optimization problem (3), we can employ adjoint models to approximate gradients and to condunct optimization steps via, for instance, line-search or trust region methods. However, these models can be labor-intensive to develop and computationally expensive to run. For instance, the adjoint model of the High Resolution Limited Area Modelling (HIRLAM) 4D-Var [10,11] was developed in 10 years in which most of the time was spent to detect and to fix errors in the tangent and the adjoint models [12]. To avoid the use of adjoint models, we can employ an ensemble of model realizations [13] as follows [14,15]: ∈ R n×1 stands for the e-th ensemble member, for 1 ≤ e ≤ N, at time k, for 0 ≤ k ≤ G. Then, the ensemble mean: and the ensemble covariance matrix: act as estimates of the forecast state x b k and the forecast error covariance matrix B k , respectively, where the matrix of member deviations reads: The model trajectory in (1) can be constrained to the space spanned by the background ensemble members (4a), this is: where w ∈ R N×1 is a vector in redundant coordinates to be determined later. This is equivalent to: therefore, the estimation of analysis increments is performed onto the sub-space given by a low-rank square root approximation of the background error covariance matrices (4c) at observation times (times where observations are available). By replacing (5) into the Equation (1) one obtains [16]: where d k = y k − H k · x b k ∈ R m×1 is the innovation vector and Q k = H k · ∆X b k ∈ R m×N . Note that, this cost function does not rely in the numerical model (2) anymore. The optimal value of the control variable w is then seek: The gradient of (6) equals: and from here, the optimal weight (7a) can be approximated as follows: from which the initial analysis state can be estimated: Since in (5), x a k represents an approximation rather than an exact analysis trajectory, the initial analysis is recovered and then, it is evolved in time by using the numerical model (2) from which we obtain an estimate of the optimal trajectory of (3): Notice, in the 4D-EnKF, all computations are performed onto the ensemble space (5) and therefore, the computational cost of estimating (7d) is linearly bounded regadring n and m [17]: Readily, posterior ensemble members at the initial time can be estimated via the implicit covariance matrix in (7c): where: In practice, model dimensions range in the order of millions while ensemble sizes are constrained by the hundreds and as a direct consequence, undersampling degrades the quality of analysis corrections onto the space spanned by (4d). To counteract the effects of sampling noise, localizations methods are commonly employed [18,19], in practice. For instance, methods such as covariance matrix localization (B-localization) [20], domain localization, and observation localization (R-localization) [21][22][23] are employed in operational DA scenarios. Yet another possible choice is to make use of precision matrix estimation. In this context, for instance, the use of the spatial-predecessors concept can be employed to obtain sparse estimators of precision matrices [24]. The predecessors of model component i, from now on Π(i, δ), for 1 ≤ i ≤ n and a radius of influence δ ∈ Z + , are given by the set of components whose labels are lesser than that of the i-th one. Of course, this will depend on the format employed to label components on a numerical grid. In practice, column major and row major format are commonly employed. This idea is exploited in the EnKF formulation proposed in [25,26] wherein the following estimator is employed to approximate precision matrices [27]: where the Cholesky factor L k ∈ R n×n is a lower triangular matrix, whose (sub-diagonal) elements β i,g,k are estimated by fitting linear models: where the variance σ 2 k is unknown, and the diagonal matrix Γ k ∈ R n×n holds the variance of residuals: where the empirical and the actual variances are denoted by var(•) and var(•), respectively.

Wind Energy Potential
The effects of climate change have triggered alarms to employ alternatives and to reduce Carbon Dioxide (CO2) emissions around the world. In many countries, regulation and CO2 reduction goals promote the substitution of fossil energy sources with Renewable Energy Sources (RES) [28]. For instance, China, the largest energy consumer worldwide, has an economic motivation to execute such substitution [29]: traditional power systems (mainly composed of nuclear, hydro, and thermal generators) are drastically decreasing, and now, they are trying to integrate RES as a shock absorber of this situation. However, RES integration is not straightforward since it brings new issues and challenges that need to be analyzed and addressed. One of the main challenges comes from the intermittency of RES [30]. Intermittency combines variability and uncertainty. The former is produced by the movement of large cloud systems owing to high and low-pressure areas. Uncertainty, also known as unpredictability, comes from the forecast error, which in turn depends on the numerical model (2). Thus, uncertainty amplification relies on model errors (i.e., physics simplifications to make numerical models computationally feasible to run). For instance, if the accuracy of the numerical model is poor, and no Data Assimilation is performed, the bias on the resulting estimate will be large concerning the actual wind speed. Thus, wind speeds can be poorly estimated, and as a direct consequence, wind energy potentials can be underestimated. Hence, Data Assimilation can be employed in this context to mitigate the impact of poor potential energy estimations via real noisy observations of wind speeds.
The potential energy p(v) in MegaWatts (MW) of a wind turbine given a wind speed v (km/h) can be estimated as follows: where v c , v r , v f , and v p are the cut-in wind speed, the rated wind speed, the cut-out wind speed, and the rated power of wind turbine, respectively. Table 1 shows the 12 wind turbine generators types assumed and utilized in many case studies [31]. The outage rate of each wind turbine reads 0.04. Commonly, the useful life of a wind turbine is about 25 years, this does not depend on its size. We also report the capital cost, and the maintenance and operating cost for each turbine, these are taken from [32,33]. Based on the Table 1, places with wind speeds below 10 km/h do not have the chance to generate electrical power from wind speeds since these are lower than the minimum cut-in wind speed across all Wind Turbine Generators (WTGs). Similarly, places with wind speeds greater than 90 km/h cannot produce electrical power because wind speeds exceed the maximum cut-out wind speed (i.e., WTG 6). Although, we do not consider economic impacts of wind-farm placements, it is important to note that wind-speed constraints have economic implications. For instance, for a place with bimodal wind speeds of 40 km/h and 11 km/h, WTGs 1 and 2 can be employed while the rest of them must be discarded in spite of the last are cheaper.

Proposed Framework
In this section, we develop an adjoint-free 4D-Var framework for potential energy estimation. The framework is divided into four stages. First, we build an ensemble of snapshots at observation times by employing a numerical model which can forecast wind components. Second, these snapshots are employed to build control spaces via a modified Cholesky decomposition. Third, the control spaces are utilized to obtain initial conditions whose wind forecasts fit a set of time spaced observations. Lastly, forecasts of wind components are employed to estimate forecasts of wind speeds, which in turn allow us to forecast potential energies of Wind Turbine Generators (WTGs). Since, in practice, model resolutions range in the order of the millions, we develop a matrix-free analysis formulation to avoid the direct inversion of linear systems during assimilation steps. All these stages are clearly detailed next.

Building an Ensemble of Snapshots
Initially, we choose a numerical model which mimics the dynamics of wind components in places of interest. For this purpose, numerical models such as the Atmospheric General Circulation Model (AT-GCM Speedy) [34] and the Weather Research Forecast (WRF) Model [35,36] can be employed. Once the numerical model is chosen, snapshots of an ensemble of model realizations (4a) are taken at G + 1 observation times. At step k, for 0 ≤ k ≤ G, the background ensemble X b k (4a) is employed to estimate a full-rank square-root approximation of the precision matrix of background errors B −1 k via a modified Cholesky decomposition (8a): At this step, we choose a radius of influence (localization radius) δ to compute the factor V T k . Beyond the scope of this radius (and the predecessors of model components) all components of V T k are assumed zero. We exploit the fact that, when the error correlations of two model components are conditionally independent (given a radius of influence δ), their corresponding entry in the precision matrix of background errors is zero. This results in a sparse Cholesky factor V T k and even more, a localized square-root precision matrix. In this manner, the impact of sampling errors can be mitigated in the square-root approximations (10). Some structures of V k are shown in Figure 1 for a one dimensional grid and different values of δ, cyclic boundary conditions are assumed for physics/dynamics. The square-root approximations (10) serve as control spaces onto which analysis increments can be estimated, therefore, the analysis increment at observation time k can be written as follows: or equivalently: where ∈ R n×n , and α ∈ R n×1 is a vector in redundant coordinates to be determined later. We assume that: Atmosphere 2020, 11, 167 8 of 24 Note that, since the square root approximations (10) are full-rank, the dimension of the spaces (11) equal those of the range of B 1/2 . This differs from what is usually employed in the literature: a control space whose dimension equals the ensemble size (5) and therefore, analysis increments can be highly impacted by sampling noise. We then expect to capture all error dynamics onto the spaces (10).
Since the initial background error covariance matrix B 0 onto the control space (11) is nothing but the identity matrix, the following error statistics hold for the prior weights α b(e) : Due to this, the 4D-Var cost function (1) onto the control space (11) can be written as follows: where Again, this cost function does not rely on the numerical model (2).

Adjoint-Free 4D-Var Optimization
Once the control spaces are estimated across observation times, the adjoint-free optimization problem to solve reads: The gradient of this cost function can be written as follows: whose root reads: and therefore an estimate of the initial analysis state (3) can be computed as follows: whose model trajectory provides a forecast which accounts for the given data into the assimilation window. Note that, the closed form expression (13b) for the optimal weights (13a) is possible since we consider linear observation operators in our formulation. The posterior ensemble onto the control space can then be built by using a square root approximation of the information matrix in (13b), it can be easily shown that the posterior error statistics read: with corresponding analysis members in the model space: x a[e] Then, the analysis members of the initial ensemble are propagated in time , for 1 ≤ f ≤ G and 1 ≤ eleN , from which an estimate of the optimal trajectory and his uncertainty (i.e., by employing a modified Cholesky decomposition on the ensemble members at time k) can be obtained. Note that, the posterior mode (13c) can be written as follows: which is nothing but a linear transformation of the prior increment to the posterior one. In this sense, the analysis step is similar to that of square root filter formulations. However, we compute the analysis increments of the initial ensemble members by using synthetic data, which is statistically consistent with the posterior error distribution: Readily: , and therefore, in spite of the posterior mode of the analysis distribution can be estimated via a linear transformation of the initial background increments, the analysis increments of the initial ensemble are actually computed by employing synthetic data. This places our proposed filter formulation into the family of stochastic formulations of data assimilation methods. Notice, given the special structure of our estimator B −1/2 k , the Woodbury matrix identity can be exploited to avoid direct inversions [37]. We denote this filter implementation Four Dimensional Variational Data Assimilation via a Modified Cholesky Decomposition (4D-Var-MC).

Post-Processing of Data, Potential Energy Estimation
Once the model trajectory is computed for each ensemble member, we proceed to map wind fields to wind energy potentials in two steps: whenever is necessary, the wind components of ensemble members are mapped to wind-speeds, 2.
this subset of information is exploited to estimate the wind energy potential of each analysis ensemble member: where w : R n×1 → R h×1 is a function that maps model states to potential energy states (this is, for each ensemble member, its wind-speed components are mapped to wind energy potentials), where h is the number of wind-speed components (with h ≤ n), and x a[e] k ∈ R h×1 is the k-th transformed member. The mapping process depends on the wind turbine employed, for instance, one can consider the wind turbines discussed in Section 2.2.
Note that, the empirical moments of the samples (16) can be exploited to estimate mean and standard deviations of wind energy potential capacities. Besides, covariances of such samples can be estimated via a modified Cholesky decomposition to understand better (and to estimate) their uncertainties.

Further Comments: Matrix-Free Formulation of the 4D-Var-MC
In practice, the number of model components n range in the order of the millions and therefore, matrix computations can be constrained by computational resources. For instance, the direct inversion of (13b) is prohibitive. Thus, it is mandatory to count with a matrix-free implementation of any data assimilation process. Following the ideas discussed in [38], we can develop a matrix-free equation for the analysis step of the 4D-Var-MC implementation. We can proceed as follows, consider: where Ω k = Q k · R −1/2 k ∈ R n×m , and O = m · G, the precision matrix in (14) can be written as follows, where ω [o] ∈ R n×1 is the o-th column of matrix Ω, for 0 ≤ o ≤ O, and I = T T · C · T is the Cholesky decomposition of I (all factors equal the identity matrix). Consider the sequence of matrices, A (2) = A (1) + ω [2] · ω [ . . .
where V (0) ∈ R n×n and Γ (0) ∈ R n×n are the factors of the Cholesky decomposition of the identity matrix I. Among iteration steps o, for 0 ≤ o ≤ O, one can see that: where Via the Cholesky factors of, and, by considering Equation (18), the matrix A (j) can be decomposed as follows, and therefore, the next relations hold: and

8: end function
Now, we are ready to test our proposed framework.

Numerical Results
In this section, we employ our proposed framework by using the Atmospheric General Circulation Model (AT-GCM) Speedy [39]. This model is a general circulation model that mimics the behavior of the atmosphere across different pressure levels [40]. The number of numerical layers in this model is 7, and we employ a T-30 spectral model resolution (96 × 48 grid components) for the space discretization of each model layer [41,42]. The number of physical variables is 5. These are detailed in the Table 2 with their corresponding units and number of numerical layers. Note that the total number of model components to be estimated reads n = 133, 632. We let the number of model realizations (ensemble size) as N = 30 for all experimental scenarios. In this case, the model resolution is approximately 4454 times larger than the sample size (n N), which is very common in operational DA scenarios. Additional details of the experimental settings are described below, some of them are similar to those detailed in [43]: • Starting with a system in equilibrium, the model is integrated over a long time period to obtain an initial condition whose dynamics are consistent with those of the SPEEDY model.

•
The initial condition is perturbed N times and propagated over a long-time period from which the initial background ensemble is obtained.

•
We employ the trajectory of the initial condition as the reference one. This reference trajectory serves to build synthetic observations. Besides, we will consider that the actual potential capacities of WTGs are based on this solution.

•
The experiments are performed under perfect model assumptions.

•
The number of assimilation steps reads G = 15. Thus, the simulation times is 7.5 days.

•
We use the wind turbines discussed in Section 2.2 for computing wind potential energies.

•
To estimate wind speeds, the wind fields (zonal and meridional components) are taken from the numerical grid at the pressure level 100 mb.

•
Our results are compared with those obtained by the 4D-EnKF formulation.

•
We employ the L − 2 error norm as a measure of accuracy for the estimation of wind energy potential: where p * (v) k is the reference wind energy potential, and p k (v) is the estimated one by a filter implementation. Likewise, k stands for observation time and v for wind speed.

•
The Root-Mean-Square-Error (RMSE) provides an estimate of the performance of a filter for a given assimilation window: • We estimate the potential energy capacities of Wind Turbines Generators (WTGs) discussed in Section 2.2.

•
Our numerical results are compared with those of the 4D-EnKF formulation discussed in Section 2.

Results with p = 50% of Observations from the Model State
The L − 2 error norms (21) of wind-energy-potential estimations for an ensemble size of N = 20 is shown in Figure 3. We employ a log scale in the y axis to render the text easier to read. As can be seen, for all WTGs, the compared filter implementations provide better estimates of potential generations than those obtained by pure forecasts, as should be expected. Thus, regardless of the employed DA method, the accuracy of forecasts can be improved by injecting real observations of the dynamical system. This can be beneficial for taking actions on whether to employ or not green sources of energy during, for instance, industrial operations. In all cases, on average, the estimated analysis trajectories in the 4D-Var-MC context outperform those computed by the 4D-EnKF formulation. In the 4D-Var-MC, the dimension of control spaces equals those of model one; therefore, we have enough degrees of freedom to capture most of the directions where errors grow faster. This allows our proposed implementation to properly correct initial background states with the information brought by observations in time. Besides, the initial analysis state (initial condition of the initial value problem) relies on the quality of the estimated background error correlations: As is proven in [26] [Theorem 1], the precision matrix estimator (8a) converges to the actual precision matrix B −1 0 as long as log(n)/N goes to zero. This value, under the current experimental settings, reads ∼0.170, which can explain as well why the accuracy of the 4D-EnKF-MC method is better than that of the 4D-EnKF. On the other hand, the control space in the 4D-EnKF formulation relies on the ensemble size, whose dimension is much lesser than that of the model one. Consequently, this sub-space can be highly sensitive to sampling noise, which can create spurious correlations among distant model components. Besides, there is no guaranty that such sub-space can capture the leading directions where errors grow faster. This results in the poor estimation of the analysis increments of the initial ensemble mean and, as a direct consequence, the analysis members of the initial ensemble. For this reason, the benefits of increasing the number of model realizations are just evident in the 4D-EnKF context; for instance, the accuracy of this formulation improves drastically as the ensemble size increases. The Table 3 provides an overview of the compared filter implementations in terms of performance (RMSE values) and all parameter configurations, RMSE values are computed based on the analysis trajectory (estimated initial condition). It is clear that, on average, our proposed filter implementation outperforms the traditional 4D-EnkF one in terms of accuracy. In general, both filters formulations can improve their performance as the ensemble size is increased. In Figure 4, snapshots of the estimated initial wind-energy-potential are shown for the proposed 4D-Var-MC method. Their corresponding standard deviations of errors (based on analysis ensembles) are shown in Figure 5. Recall that this initial state is our estimate of the initial condition in the optimization problem (3). As can be expected, most of the wind-energy-potential is produced on the ocean where wind speeds get the largest rise for all wind turbines. This serves as a validation test since no wind farms (turbines) can be placed under such a place. However, countries well-known for their potential capacities are just evident in these results, for instance, countries such as those from Latin American and the Caribbean and Africa. Moreover, by taking a close look at standard deviations, one can see that forecasts are obtained with low uncertainties for all filter implementations. This, together with the RMSE results, show that the proposed framework can be employed to estimate wind energy potentials with high accuracy and low variations.
Notice, green sources of energy such as those based on wind speeds are impacted by three observable conditions (which can be implicitly evidenced in our numerical results): variability, unpredictability, and placement. Variability obeys to the fact that as time moves forward, wind speeds can drastically vary, this impacts the potential energy that can be generated by WTGs. Regardless the DA method employed to estimate wind potential capacities, unpredictability is always present: numerical forecasts are imperfect and even more, uncertain. Placement of WTGs is crucial, as can be seen in Figure 4, no all WTGs can properly work in different zones of the globe, this is, electrical power via WTGs can drastically vary from one place to other. We can stood out the importance of employing WTGs as sources of energy based on wind speeds in different regions of the globe but, we cannot argue which turbine is better than others. To do this analysis, we should consider other relevant factors such as wind speed variability, WTGs constraints, and economic considerations, the last two are out of the scope of our analysis.
Consider, again, the WTG parameters reported in the Table 1 and the initial snapshots reported in Figures 4 and 5. Note that, in most places in the globe, WTGs 1 and 2 can be installed to guaranty electrical power from wind speeds; these WTGs are the ones with rate-capacity 0.5. Note that, as the cut-in wind speed and the rate-capacity increase, the electrical energy generation of WTGs can be impacted. For instance, near the poles, the power generation is almost null for WTGs with the largest cut-in wind speed parameters. WTGs with rated-capacity of 1 are the ones that have large variability across different places in the world. Lastly, the largest amount of energy across different places in the domain can be obtained for WTGs with the largest rate-capacity values. However, WTGs with low rate-capacity values are the ones whose numerical forecasts are obtained with lesser variability (i.e., in Figure 5 the standard deviation of errors has a homogeneous behavior across different regions of the domain). As the rate-capacity increases, the variance of the standard deviation of errors increases as well. This means, more variability of errors can be evidenced across different parts of the world. Thus, WTGs with large rate capacities provide forecasts with a large amount of clean energy, but these come with large uncertainties in certain regions of the world, which can difficult decision making. Time Step (ℓ) 4D-Var-MC 4D-EnKF Background (l) Wind Turbine 12

Single Observations across Observation Times
In this section, we briefly discuss the performance of our proposed 4D-EnKF-MC method by using a single observation test. We hold the same experimental settings as those in Section 4.1 and report the estimation errors in the initial conditions via L − 2 norms (21). The single observation, across all observation times, is placed as is shown in Figure 6. In the Table 4, we can clearly seen the advantages of employing the control space (11) instead of the traditional approach based on the ensemble sub-space (5). For instance, for WTG with low rate-capacity, error differences are of one order of magnitude. In Figures 7 and 8, we report the potential energy estimation and the uncertainty of each component (as the standard deviation of errors from the initial members of the analysis ensemble). As can be seen, low-rate capacity WTGs such as the WTG 1 and the WTG2 provide estimates whose error dispersion is small. Again, as the rate-capacity increases, the spread of ensemble members grow. The accuracy of the proposed method obeys to the fact that the precision matrix is full-rank, well-conditioned, and even more localized. Thus, the impact of spurious correlations is mitigated in the analysis increments of the initial ensemble.

Conclusions
We propose a 4D-Var ensemble-based data assimilation framework for wind energy potential estimation. In this formulation, in the 4D-Var context, the intrinsic need of adjoint models is avoided via the use of an ensemble of model realizations. These ensembles are employed to build control spaces onto which analysis increments are estimated. Control spaces are built via a modified Cholesky decomposition. The particular structure of this estimator allows for a matrix-free implementation of the proposed filter formulation. Experimental tests are performed, making use of wind turbines catalogs and the Atmospheric General Circulation Model Speedy. The results reveal that our proposed framework can properly estimate wind energy potential capacities within reasonable accuracies in terms of Root-Mean-Square-Error, and even more, these estimations are better than those of traditional 4D-Var ensemble-based methods. Besides, Wind Turbine Generators (WTGs) with low rate-capacity are the ones which provide homogeneous behavior of error estimations around the globe. As the rate-capacity increases, the potential energy increases as well, but the error dispersion of ensemble members grow, which can difficult decision-making processes. Of course, rate-capacity is just a single parameter of many in the WTG context, and we do not consider, for instance, economic aspects in our study, which can be crucial for deciding whether or not to employ green sources of energy.