Technoeconomic Analysis of Changing PV Array Convective Cooling Through Changing Array Spacing

Accuracy in photovoltaic (PV) module temperature modeling is crucial to achieving precision in energy performance yield calculations and subsequent economic evaluations of PV projects. While there have been numerous approaches to PV temperature modeling based on both the steady-state and transient thermal assumptions, there have been few attempts to account for changing convective cooling on PV module surfaces resulting from changes in the PV system layout. Changes in system row spacing, in particular, can have a meaningful impact on module electrical efficiency and subsequent economic performance, even when considering additional costs from the changes in row spacing. Using a heat transfer approach based on the spatial definition of a PV array, technoeconomic analyses of different plant configurations are presented here that show an improved system levelized cost of energy (LCOE) for fixed-tilt PV systems when increasing system row spacing. These LCOE improvements have been found to be as high as 2.15% in climates characterized by low ambient temperatures and higher average annual wind speeds in U.S. climates. While the LCOE improvements are primarily driven by incident irradiance changes for altered row spacing, the waterfall analysis of the different components of changing system LCOE show that modifications in the heat transfer dynamics have a 0.5% contribution to the LCOE reduction for the largest LCOE, compared with a 3.3% reduction from irradiance changes.


I. INTRODUCTION
P HOTOVOLTAIC (PV) module temperature is second only to irradiance in terms of importance to module electrical output. The module max power output is decreased by an average of around 0.3%-0.5% per degree increase in the module temperature [1]. As such, the module temperature must be accurately modeled in PV performance models to account for the thermal effects on PV array electrical output. Many PV performance models rely on steady-state models, such as the Sandia steady-state model and the nominal operating cell temperature (NOCT) model that calculate the temperature based on an assumed thermal equilibrium condition at each time step, and the input parameters of incident irradiance, wind speed, and ambient temperature [2], [3]. These steady-state models have shown that the irradiance is the primary driver of the module temperature, with the wind speed having a secondary effect and the ambient temperature serving as a temperature offset that sets the base environment temperature from which the module temperature increases or decreases. There are also models that use an iterative heat transfer balance of the PV module at each time step to converge on the module temperature that results in thermal equilibrium between the module and the environment. These models take the convection, radiation, and incident irradiance into account when determining the module temperature and typically require numerous input parameters that are often difficult to measure or quantify from module specifications [4].
For models that require finer temporal resolution to account for the transient thermal behavior of the module, transient models such as those in [5] and [6] account for the heat capacity of the module to reduce modeled temperature variability due to abrupt changes in environmental conditions. These models account for changing thermal effects over time, but they are often based on fixed heat transfer assumptions that do not account for the changing heat transfer flow due to changing wind speed and ambient temperature. Both the steady-state and transient modeling approaches have been validated against measured PV system data and are commonly used in PV performance modeling practices.
While there are several approaches to PV temperature modeling that are widely used in research and industry, there have been few attempts to account for the changes in heat transfer This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ flow caused by differences in array configuration parameters, such as row spacing and panel height above the ground. There have been previous attempts to model the changes in heat flow for changing air gaps in rooftop PV systems through the use of temperature offsets based on bins of gap sizes between the panel backside and mounting surface [7], but this approach fails to capture the changing heat transfer balance between the PV modules and the environment in a given time step. Changes to the array configuration can affect the total irradiance that is incident on the module surfaces and, thus, increase or decrease the module temperature. These system modifications also alter the path of heat-removing wind flow, changing convective cooling efficiency for PV systems [8], [9], [10]. This effect is considered via the convective heat transfer coefficient h, an input to cell temperature models. Determining h for a PV power plant is often empirical, using either location-specific quantities dependent on wind conditions as in [11] or flat plate convection correlations based on only panel geometries [12], [13], [14]. Yet, solar module convection is not merely dependent on plate dimensions, and each PV power plant arrangement is uniquely designed for its location-optimized for maximum solar incidence, land availability, and environmental conditions. Thus, a more versatile approach to PV temperature modeling acknowledges module heat removal as a full plant phenomenon and takes array configuration parameters as input. This article presents a new perspective on performance and cost modeling for PV power plants, where the convective heat transfer coefficient h is derived considering 3-D geometric features. Here, we discuss a lacunarity-based length scale, inspired by forest and urban canopy research, describing the spatial heterogeneity of 3-D PV power plant arrangements. The lacunarity scale describes the spatial arrangement of the PV systems by taking parameters, such as row spacing, panel height, and tilt angle into account. This parameter is calculated for a number of fixed-tilt PV power plants and embedded within a modified convection correlation, resulting in a value for h, which takes all the geometric plant features into account. The result is then applied to the open-source performance and financial models in the NREL System Advisor Model (SAM) Version 2021.12.02 [7], finding a dependence of both levelized cost of energy (LCOE) and annual energy per module on lacunarity for a given PV power plant.

II. HEAT TRANSFER BACKGROUND
Convective cooling for PV systems is generally a function of inflow conditions and is related to the Reynolds number (Re = U ∞ L/ν), describing the ratio of inertial to viscous forces in a flow. The term U ∞ is the inflow wind velocity, ν is the kinematic viscosity of the air, and L is a length scale, which characterizes the object or space which the air moves around or through. The Reynolds number is then typically used to derive the Nusselt number that describes the ratio of convective to conductive cooling on a surface through empirical relations. The knowledge of the Nusselt number and thermal conductive properties allows for the calculation of the heat transfer coefficient h to determine the rate of heat transfer between the surface and the surrounding environment. For individual flat plates, the characteristic length L often represents a reference dimension for a single surface length. For example, PV cooling models employ flat plate convection theory, where the characteristic length often used is L c = 4 A/P , with A being the surface area enclosed by perimeter P [15]. However, convection varies significantly for PV plants with the identical module dimensions but varied plant configurations-suggesting that the full array geometry is responsible for flow behavior and relative cooling [8], [16]. To account for these effects, we employ a length scale L Λ based on the 3-D lacunarity to quantify the heterogeneity of solar arrays as canopy flows [17], [18]. As further described in Section VI, L Λ characterizes the space through which flow is able to move in a solar array as a measure of "gappiness" based on the full plant geometry [19]. Embedding this parameter as the length scale in Re allows for a Nusselt number correlation of potential PV cooling Nu = f (Re, P r), which encompasses all the unique physical characteristics present in solar arrays [18]. This Nusselt number correlation was derived from operating field data, wind tunnel measurements, and numerical simulations and is shown in the following equation [9], [18]: where L H is the height of the solar array canopy. Fit coefficients found in this study of a = 0.09, m = 1/5, n = 1/12, and b = 1.91 resulted in a logarithmic relation with the coefficient of determination over R 2 = 90%. The works [9], [16], and [18] include more of the detailed heat transfer analysis used to determine the Nusselt number relationship used in this article. With this relation, a given PV plant canopy and flow configuration could estimate convective cooling through solving for h and use this value for h in existing thermal models to more accurately calculate power production.

III. TECHNOECONOMIC ANALYSIS
Changes in the PV array configuration that alter the convective heat transfer flow over the PV module surfaces must be evaluated for their impact on the technoeconomic analysis (TEA) performance of PV systems, as increased row spacing or ground coverage ratio (GCR) results in an increased balance of system costs (such as wiring costs) and increased land costs in addition to the changing energy performance. The GCR is defined as the ratio of the diagonal length of the side of a row of PV modules over the row-to-row horizontal distance between identical points on the module in adjacent rows. The GCR in this analysis is determined from the following relation [7]: where GCR is the unitless ground coverage ratio, l is the length of the side of the array in meters, and d is the distance between rows in meters. Additional land lease costs for increased acreage are needed to install a system with the same rated energy output for the increasing row spacing associated with decreasing GCR values. Changing the array layout can also affect the amount of incident irradiance on the surfaces of the module and the amount of energy lost due to wiring losses. The TEA allows for the comparison of the performance gains from changing irradiance and convection conditions against the additional project costs in order to evaluate the cost-effectiveness of the thermal cooling considerations for PV power plant projects. The analysis is performed using the SAM, an open-source NREL tool with detailed PV performance models and financial models [7], [20]. The thermal modeling with the proposed convection handling is evaluated against existing convective heat transfer coefficient calculations used in the SAM. The existing convective heat transfer coefficient is calculated based on the following equations: where Nu is the Nusselt number, Re is the Reynolds number, P r is the Prandtl number, k air is the thermal conductivity of air at the ambient conditions in W/mK, and L c is the characteristic length in meters of the PV module, which is defined as where A is the module surface area in square meters and P is the module perimeter in meters. The proposed convection based on array lacunarity length scale presents a modified Nusselt number fit and an h calculation based on the lacunarity length scale L sc rather than the characteristic length L c . The performance analysis considers a 1000-kWDC South-facing PV system with a fixed tilt angle of 30 • over different GCR conditions, with weather data from different climates to evaluate the seasonal effects of the changing heat transfer balance for different GCR values. The additional balance of system costs scale with the GCR based on previous sensitivity analysis on PV system costs for array configuration parameters [21]. An annual land lease cost of 0.054 $/m 2 was also applied to the analysis to evaluate the increasing land cost for increased acreage needed for the system.

A. Convection Cooling Gains
Before analyzing the impact of changing array configuration on overall plant performance in different climates, it is helpful to show how the different modeling approach presented in this article impacts energy output from PV systems as compared with the existing convection modeling practices used in the SAM. For both the modeling approaches, increasing the distance between rows increases the amount of incident irradiance on a PV module surface through increased ground reflections and decreased row-to-row shading, in addition to the increased convective cooling that comes with increased spacing. The newly defined convection model takes the amount of row spacing into account when calculating the amount of convection cooling due to the wind flow, which can result in lower module temperatures and, thus, higher energy output than the existing convection model. This is shown in Fig. 1, which shows per module annual energy outputs for the reference system in Phoenix, AZ, USA, for both the convection modeling approaches. The annual energy in kilowatt for each GCR value is normalized to the value at When comparing results with those found at a GCR value of 0.46 for lacunarity convection calculations, the annual per module output at a GCR of 0.58 decreases by 3.4% due to increased row-to-row shading and less convective cooling flow on the module surface. Conventional convection modeling results in an energy decrease of 2.64%, as the changing convective cooling flow is not considered. When increasing row spacing to a GCR value of 0.35, the energy output increases by 1.7% due to increased incident irradiance and convective cooling. This result indicates increased predicted energy output over the existing convection heat transfer modeling, which results in a 0.98% annual energy increase for the same GCR change. Comparing these two modeling approaches across a wide range of GCR values gives a better understanding of the impact the convective cooling dynamics in the lacunarity approach can have on the energy output of a system.

B. LCOE Evaluations
After quantifying the potential performance gains from the convective cooling provided by increased row spacing, the convection heat transfer coefficient approach outlined in previous sections must be evaluated for its impact on technoeconomic performance for a variety of system types and climates. This analysis includes the effects of irradiance, shading, and convective cooling for changing GCR, as well as cost consideration stemming from changing balance of system costs and land lease costs. The annual energy from the system archetypes described in the previous subsection was evaluated along with the LCOE for a variety of climate conditions. The LCOE is a metric that encapsulates the life-cycle costs of the PV project in the numerator and the life-cycle energy generation of the system in the denominator for a resulting single value in units of cents/kWh, as shown in the following equation [7], [20]: where C 0 is the initial project investment in $, C n is the annual project cost in year n, Q n is the electricity delivered by the system in kWh in year n, and d is the discount rate of the project. Evaluating the shifts in the LCOE for changing GCR allows one to determine how the performance gains from irradiance gain and convective cooling change compared with the cost increases.
Decreasing LCOE would indicate that the performance gains outweigh the additional costs, whereas increasing LCOE would indicate that the additional costs outweigh the performance gains.
The lacunarity length scale of a PV plant was varied from 2 to 11 m, which corresponds to GCR values of 0.73 to 0.08. Each system had linear increases in the balance of system costs ranging from $0.16/kWDC for the high end of the GCR scale to $0.32/kWDC for the lower end of the GCR scale. The annual land lease cost was held at a constant $ 0.054/m 2 as the system acreage increases with increasing row spacing based on estimates motivated by cost numbers used in [21]. Typical meteorological years (TMY) for the climates of Portland, OR, and Phoenix, AZ, were simulated in the SAM to get annual energy and LCOE outputs. Each climate was modeled using the heat transfer balance approach that is used in the SAM [7]. The SAM heat transfer model used rack mounting configuration and module dimensions assumptions to best model the effects of the changing convection flow on the module for large-scale fixed-racking systems. The convective heat transfer is calculated using the approach defined in the previous section, with the lacunarity length scale input changing based on the changing system GCR. The annual energy output of the system is normalized to the output of a single module on an annual basis. The graphical results of the analysis for the different climate conditions are shown in Fig. 2. Results in this figure are normalized for both the LCOE and module annual energy output values found at a GCR value of 0.46, as was done in Fig. 1 Table I to show the relevant cost increases resulting from the changing GCR and lacunarity length scale. Results reveal that the system GCR value where the performance gains no longer outperform the increase in costs changes depending on the climate conditions for the system.
Each plot in the figure shows increasing energy production per module for decrease in the system GCR due to the increased incident irradiance and increased convective cooling. The increased spacing also requires increased system costs due to greater land usage and higher balance of system costs for system components such as wiring. The LCOE values decrease to a minimum point along the GCR axis before eventually increasing and showing decreased system economic performance. The point of LCOE minimum represents the point at which the increased energy output for decreasing GCR no longer outweighs the increase in costs associated with increased array row spacing, as increased energy output decreases the system LCOE, while increased costs increases the system LCOE. The point along with GCR axis where the LCOE reaches its optimum (minimum) value varies depending on the climate of the PV array site. Optimizing LCOE is not the only consideration in the PV system design, but should be taken into consideration along with available land area, annual energy requirements, and capital cost constraints. Accounting for the convective cooling changing for changing system layout allows for more accuracy in system energy yield modeling when taking these investment considerations into account.
To further evaluate the impact of the changing convective cooling for increased row spacing, the LCOE change from a system with a GCR of 0.46 to a GCR of 0.35 was evaluated for the 50 U.S. states using the latest TMY weather data for each state's capital. These GCR values were chosen to represent  changing system GCR from a value more indicative of fixed tilt systems to one more representative of the row spacing used for single-axis tracking systems. The analysis was performed for monofacial and bifacial systems with the same system parameters as those used to generate Fig. 2 [22]. For the bifacial systems, a 0.05 $/Wdc cost increase was applied to simulate the increased module costs for the bifacial panels. The 0.05 $/Wdc was a conservative estimate based on literature review. Some literature, such as [23], indicates that cost increases closer to 0.02 $/Wdc may be acceptable. The resulting heatmaps of the percentage improvement (decrease) in the LCOE for the change in the GCR are shown in Fig. 3 for the monofacial system and Fig. 4 for the bifacial system. These results show that for both the bifacial and monofacial systems, the greatest LCOE changes occur in climates characterized by low average annual ambient temperatures and moderate-to-high average annual wind speeds in U.S. climates. Warmer climates with lower average annual wind speeds see little to no improvement in the system LCOE as the cost increases negate the energy performance gains in the LCOE calculation. Wind direction was not found to have a meaningful impact on results. For the monofacial systems, it can be seen that in certain locations, the LCOE actually increases or worsens for the decrease in the system GCR, indicating that the improvements in convective cooling and incident irradiance are outweighed by the cost increases and increased wiring losses caused by the GCR change. The energy gains are greater for the bifacial systems due to the increased ground reflected irradiance on both the module front and rear surfaces for increased module row spacing. A preliminary analysis of tilt angle effects on the system LCOE and energy production show that increasing the tilt angle increases the convective heat transfer coefficient due to increased turbulent wind flow incident on the module surface. This increased convective cooling, along with increased energy production for tilt angles closer to the site latitude for most states, leads to improvements in the monofacial system LCOE across all but three states when going from 30 • tilt assumption to 41 • tilt assumption. More detailed model sensitivities will be analyzed in future efforts.
Waterfall model analysis was also performed on a subset of systems to isolate the effects of the changing convective cooling on the system LCOE from the effects of other factors, such as changing incident irradiance, balance of system cost increases, land lease cost increases, and dc wiring loss increases. This was done by first modeling the GCR 0.46 and 0.35 systems normally as was done in the previous figures. To isolate the effects of convective cooling on the temperature, the module temperature outputs from the GCR 0.35 results were applied to the GCR 0.46 model with all the other inputs remaining the same. Similar modeling approaches were used to isolate the remaining changing model factors. The incident irradiance of the 0.35 system was applied in a model with the module temperatures from the GCR 0.46 system. Balance of system costs, land area costs, and additional dc wiring losses were applied individually to the GCR 0.46 system to isolate those effects that occurred when changing the GCR value to 0.35. The dc wiring losses associated with changing cable length were assumed to increase from a 2% dc loss for the GCR 0.46 system to a 2.2% loss for the GCR 0.35 system. The resulting LCOE changes from each of these isolated theoretical simulations were then weighted to determine the individual effect each factor had on the LCOE change. The waterfall plot for this modeling approach is shown in Fig. 5 for North Dakota, the state with the highest LCOE change in Figs. 3 and 4. The blue in the plots indicates reduction or improvement in the system LCOE, while red indicates the increasing LCOE. While the incident irradiance is the highest weighted factor on the changing LCOE as the primary variable in PV performance modeling, it can be seen that the temperature changes do result in an approximate 0.5% reduction in the system LCOE compared with a 3.3% improvement from irradiance and 1.4% increase in the LCOE from increased wiring costs. Similarly for the Arizona system waterfall plot in Fig. 6, it can be seen that while the LCOE is relatively constant for the GCR change, the temperature changes from changing convective cooling have a much greater weight on the LCOE changes than in the colder climate of North Dakota. The Arizona system has a 0.67% improvement in the LCOE from both the temperature changes and irradiance changes, and the same 1.4% increase in the LCOE due to increased wiring costs. While the wiring cost, wiring loss, and land cost changes are identical for the North Dakota and Arizona systems, the difference in annual energy between the systems results in different LCOE behavior. These results show that the convective cooling considerations of the PV array must be a part of the consideration when designing PV systems.

IV. CONCLUSION
This article described a PV module temperature modeling approach that accounts for the changes in convective cooling flow on the module surfaces that come with changes in array spacing and configuration. The convective heat transfer curve was generated through computational flow simulations and wind tunnel experiments that allowed for convective heat transfer to be described for a lacunarity length scale value that describes the spacing of the entire PV array through a single length unit. Implementing this unique approach to convective heat transfer with the PV models in the SAM revealed that accounting for the changing convective cooling flow could lead to around 1.5% difference in annual energy values when compared with conventional convection modeling approaches. TEA performed in the SAM reveals that the changing GCR can improve the system LCOE, with the primary energy increases from increased incident irradiance and secondary convective cooling effects from increased row spacing outweighing the increased land and wiring costs in certain climates for monofacial and bifacial PV systems.
Future work for this convective modeling approach should include the sensitivity analysis and reduced -order modeling of the lacunarity length scale definitions of PV arrays to allow for an easier implementation of the heat transfer modeling presented here in full PV model workflows. This would also allow for less computational expense in the lacunarity calculations, which would aid in applying these heat transfer considerations to tracking axis systems that would require a recalculation of the heat transfer paradigms at each time step. Module performance degradation was also not modified in these analyses, but improvements in module thermal performance would in theory lead to reduced degradation rates that could increase the reduction in the system LCOE over system lifetime. Initial estimates for the temperature changes associated with the GCR changes in Figs. 3-6 indicate estimated service life increases of 10%-52% for North Dakota following the methodologies presented in [24]. Finally, a more detailed sensitivity analysis of the effects of system parameters, such as tilt angle, would offer more insights into the balance of heat transfer, system energy output, and system costs that users should consider in system planning.

APPENDIX LACUNARITY THEORY AND APPLICATION
This section contains detailed information about lacunarity calculations for solar arrays. Section IV-A expands on background and theory for fractal lacunarity. Links to MATLAB code and descriptive instructions for user application are available in Section VI-B.

A. Theory
Originally coined by Mandelbrot in 1983, lacunarity describes fractal geometries with respect to their encompassed empty space, where a larger lacunarity represents increasingly heterogeneous arrangements with more "gapping" regions in a given pattern [19]. This analysis has since been adapted to characterize a variety of self-similar data for a very diverse set of applications [17]. For example, lacunarity has been used to assess surface roughness for material manufacturing, quantify damage on regions in the brain affected by neurological disease, and even characterize large-scale canopies such as heat islands in urban environments and effects of deforestation [25], [26], [27], [28]. We can also consider PV arrays as a type of canopy, where flow is able to pass through and around the pattern-like rows of heated objects in space [18]. The method used here to ascribe lacunarity values to a theoretical PV array is adapted from forest canopy research as described in [17] as an expansion on the box-counting technique introduced by Plotnick et al. [29].
While applicable in 1-, 2-, or 3-D space, here, we consider the solar array as a 3-D pattern within a finite volume, so changes in module height and inclination can also be accounted for. To calculate lacunarity, a small finite "box" of some finite size r and some shape is moved throughout the plant volume, collecting a sum ni at each relative location based on instances of occupied space. The lacunarity value for this size r is then calculated as Λ(r) = s 2 (r) s 2 (r) + 1 where s 2 (r) = 1 N N i=1 is the mean occupation and s 2 (r) is the relative variance. This process is repeated for subsequently larger r values until a desired spatial limit is reached, resulting in a vector of nondimensional lacunarity Λ(r) respective of each box size r. At a large enough r, Λ(r) will asymptote, signifying the maximum length scale (r max ) that contributes to the heterogeneity of a certain pattern. Using all box size values up to r max (R = [r 1 , r 2 , . . ., r max ]), we can define a lacunarity-based length scale for each unique PV array as where Λ R = [Λ(r 1 ), Λ(r 2 ), . . ., Λ(r max )] and N R is the number of values considered in R and Λ R . This length scale L Λ is then applied as the length scale in Re Λ = U ∞ L Λ /ν. This modification on Re is the key feature in quantifying convection for PV arrays, as it characterizes the heat-removing flow in terms of the space through which it can move within the array. Note that in the cases used for this study, r max generally coincided with row-to-row spacing for a given plant arrangement.

B. Application
The following GitHub repository contains MATLAB code for calculating PV array lacunarity given geometric array parameters as written and employed for [17] and [18]. 1 The associated README file describes the necessary function calls and dependencies and contains all the relevant/necessary geometry values and author/contact information.