The Impact of Utility-Scale Photovoltaics Plant on Near Surface Turbulence Characteristics in Gobi Areas

: With the rapid deployment of utility-scale photovoltaic (PV) plants, the impact of PV plants on the environment is a new concern of the scientiﬁc and social communities. The exchange of sensible and latent heat energy and mass between land and air in PV plants is crucial to understanding its impact. It is known that the near surface turbulence characteristics rule the exchange. Therefore, it is essential for understanding the impact to study the characteristics of near surface turbulence. However, it is not well recognized. Turbulent ﬂuxes and strength characteristics for the PV plant and the adjacent reference site in the Xinjiang Gobi area were investigated in this study. Various surface layer parameters including friction velocity, stability parameter, momentum ﬂux, and turbulent ﬂux were calculated using eddy correlation system. Results indicate that compared to the reference site, near the surface boundary layer was more unstable during the daytime due to the stronger convection heating, while it was more stable at night in the PV plant. In the PV plant, I u was weakened and I v was strengthened during the daytime, and I u and I v were all weakened at night, while I w was strengthened across the whole day. The signiﬁcant difference between I u and I v in the PV plant indicated that the horizontally turbulence strengths were affected by the plant layout. The turbulent kinetic energy of the PV plant was lower than the reference site and the momentum in the PV plant was higher than the reference site, especially during the daytime. Compared to the reference site, the PV plant had a higher sensible heat ﬂux and less latent heat ﬂux. The turbulent components of wind followed the 1/3 power law in the unstable conditions and stable conditions in the PV plant and the reference site.


Introduction
Solar energy development is considered as one of the efficient countermeasures against environmental crisis. The rapid progress of PV technologies and the dramatic reduction in cost has impelled the development of the PV industry [1,2]. Nowadays, more and more PV plants have been setup worldwide. Consequently, one concern, the environmental impact of PV plants is raising awareness in the public and scientific community. The exchange of sensible and latent energy and mass between land and air in PV plants is crucial to formulate its impact. It is known that the near surface turbulence characteristics rule the exchange.
The transport of momentum, heat, and moisture between the Earth's surface and atmosphere is dominated by turbulence, which, as a key physical process in the atmospheric boundary layer (ABL), plays an important role in the formation and evolution of the planetary boundary layer [3][4][5][6][7]. Topography modifies atmospheric turbulence through various processes including radiative, thermodynamic, and dynamic effects [8,9]. The heterogeneity in surface cover generates a significant impact on the micrometeorological parameters [10,11]. The characteristics of the micrometeorological parameters of surface boundary layer such as turbulent fluxes, turbulent diffusion, and wind profiles have been described in previous studies [12][13][14][15][16][17].
According to the Monin-Obukhov (M-O) similarity theory, the non-dimensional wind speed component can be expressed as a universal function of the non-dimensional stability parameter for the isotropic boundary layer with a steady and horizontally homogeneous underlying surface [18,19]. In recent decades, with the advancement of theoretical analysis and the progress of observational technology, many investigations and verifications of M-O similarity theory over different surface conditions such as urban, plateau, ocean, and mountainous regions have been carried out [20][21][22][23][24][25][26][27][28]. Panofsky et al. found that the relationship between normalized standard deviation of wind speed and the stability parameter on the flat underlying surface was in accordance with the 1/3 power law [20]. M-O similarity theory has failed to be verified for the coastal marine boundary layer in near-neutral and stable conditions [21]. Moreover, the characteristics of turbulence over mountainous terrain demonstrates that M-O similarity theory is also applicable to the anisotropic boundary layer [28].
Although some features of the surface boundary layer over different underlying types have been investigated in the above studies, few studies have been carried out in areas with a utility-scale deployment of PV plants. Generally, darker-colored PV panels change the radiative forcing on the Earth's surface by reducing the surface albedo [29,30] and usually increases surface roughness [31]. In previous studies, some observational studies found that a utility-scale PV plant in barren areas produced a heating effect of about 2-4 K on the near surface air due to low albedo and specific heat capacity of PV panels [32][33][34]. A recent observation found that PV plants increased the sensible heat flux and reduced the storage and release of soil heat flux [35]. Unfortunately, near surface turbulence characteristics, which are crucial to formulate these impacts of PV plants, are not well recognized.
In this study, we investigated the characteristics of the near surface air turbulence, dynamical, and thermal-dynamical properties in a utility-scale PV plant and an undisturbed reference site in the Wujiaqu Gobi area in Xinjiang. The rest of this paper is organized as follows. Section 2 introduces the experimental design and the derived methodology related to the turbulence parameters discussed in this study; Section 3 presents the comparison results and relevant discussion; and Section 4 summarizes the conclusions.

Site Overview
This experiment was conducted in a utility-scale PV plant in the Wujiaqu Gobi area of Xinjiang. Wujiaqu is located at the intersection of the desert and oasis in the southeast of the Jungar Basin. It is cold in winter and hot in summer, and it is dry all year round, a typical continental temperate arid climate. The PV plant had a length of 1120 m from north to south and a width of 1030 m from east to west, covering an area of about 1.15 km 2 and a capacity of 70 MW. The PV panel lies with the inclination angle of 37 • , the row spacing of the PV array was 7.5 m, and the top edge of the PV panel was 2.5 m above the ground. The panel material was polysilicon and the maximum power was 250 W.
We set up two observation systems in and out of the plant, respectively. One (Figure 1c) was located out of the plant and taken as the reference site as it was free of PV panels, compared to the one inside the plant (Figure 1a) tude of 437 m. The terrain is relatively flat. The observation systems contained two meteorological observation towers of 11 m high, respectively. The underlying surface is Gobi with sparse vegetation. The soil type is sandy soil, and the soil texture is relatively coarse. The various turbulence parameters were derived from the ultrasonic open-circuit vortex covariance system, which were installed at the height of the 5.5 m level in the PV plant and at a 3 m level in the reference site, respectively. The equipment of the ultrasonic open-circuit vortex covariance system used in our experiment was IRGASON-IC-BB by Campbell, USA. It integrates an open-circuit infrared gas analyzer and a three-dimensional ultrasonic anemometer. The drift error was less than ±8cm/s for Ux, Uy, and Uz. IRGASON-IC-BB can measure the vertical flux of heat and water vapor in the atmosphere and various pulsations. The data collection and sampling frequency was 10 Hz, and a 30 min-averaged wind speed was used for calculating the related surface layer parameters. The data period was from 8 June to 31 October 2019, and the time system used in this paper was Beijing time.

Flux Data Processing and Calculation Methods
The high frequency ultra-sonic anemometer observational data were processed by Easyflux, which is an online processing software used for post-data processing and thee quality control of ultrasonic observations. The main procedures include outlier elimination, double rotation and plane fitting method coordinate rotation [36], frequency delay time correction [37], ultrasonic virtual temperature correction [38], and WPL correction [39]. Finally, 30 min-averaged turbulent fluxes were calculated, and the data quality evaluation such as the turbulence stability test and development adequacy ITC (Integral Turbulence Characteristics) test were used to check the accuracy of the fluxes.
The surface layer parameters analyzed in this study mainly include friction velocity * , temperature scale * , stability parameter , roughness length 0 , turbulence strength , turbulent kinetic energy , momentum flux , sensible heat flux , and latent heat flux . The eddy correlation method was employed to derive these parameters. Normalized standard deviation of 3D wind were analyzed as a function of .
In the boundary layer, after introducing dimensionless variables through scale parameters, the contours of many variables can be expressed in a general form [40]. Among them, the scaling speed, as a basic surface speed scaling parameter, also called friction speed, can be obtained by the following formula: The various turbulence parameters were derived from the ultrasonic open-circuit vortex covariance system, which were installed at the height of the 5.5 m level in the PV plant and at a 3 m level in the reference site, respectively. The equipment of the ultrasonic open-circuit vortex covariance system used in our experiment was IRGASON-IC-BB by Campbell, USA. It integrates an open-circuit infrared gas analyzer and a threedimensional ultrasonic anemometer. The drift error was less than ±8 cm/s for Ux, Uy, and Uz. IRGASON-IC-BB can measure the vertical flux of heat and water vapor in the atmosphere and various pulsations. The data collection and sampling frequency was 10 Hz, and a 30 min-averaged wind speed was used for calculating the related surface layer parameters. The data period was from 8 June to 31 October 2019, and the time system used in this paper was Beijing time.

Flux Data Processing and Calculation Methods
The high frequency ultra-sonic anemometer observational data were processed by Easyflux, which is an online processing software used for post-data processing and thee quality control of ultrasonic observations. The main procedures include outlier elimination, double rotation and plane fitting method coordinate rotation [36], frequency delay time correction [37], ultrasonic virtual temperature correction [38], and WPL correction [39]. Finally, 30 min-averaged turbulent fluxes were calculated, and the data quality evaluation such as the turbulence stability test and development adequacy ITC (Integral Turbulence Characteristics) test were used to check the accuracy of the fluxes.
The surface layer parameters analyzed in this study mainly include friction velocity u * , temperature scale T * , stability parameter ζ, roughness length Z 0 , turbulence strength I, turbulent kinetic energy e, momentum flux τ, sensible heat flux H, and latent heat flux LE. The eddy correlation method was employed to derive these parameters. Normalized standard deviation of 3D wind were analyzed as a function of ζ.
In the boundary layer, after introducing dimensionless variables through scale parameters, the contours of many variables can be expressed in a general form [40]. Among them, the scaling speed, as a basic surface speed scaling parameter, also called friction speed, can be obtained by the following formula: The temperature scale and length scale parameters are as follows: Among them, the fluctuations of the three-dimensional wind speed component and temperature are obtained by subtracting the average value from the instantaneous value of each parameter, that is, Karman constant, and g is the acceleration of gravity.
Atmosphere stability parameter ζ is a measure of the buoyancy convection ability and is obtained by: where z is the measurement height (5.5 m for the PV plant and 3 m for the reference site), d is the zero-plane displacement. The roughness length Z 0 depends not only on the surface type, but also on the surface roughness texture, and can be derived from a wind profile under neutral conditions [41]: Turbulence strength is defined as the ratio of the standard deviation of the three wind speed components to the average wind speed: where σ u , σ v , σ w are the standard deviation of wind speed fluctuations: The turbulent parameters are computed using the following formulae: where ρ and C p are the density and the constant specific heat capacity of air, respectively; q is the fluctuation of water vapor density; and λ is the latent heat of evaporation. According to the M-O similarity theory, the normalized standard deviation of the turbulent velocity component can be expressed by the universal function of the stability parameter ζ: where a and b are the site-specific modeling coefficients.

Atmospheric Stability Parameter and Aerodynamic Roughness
Atmospheric stability ζ and aerodynamic roughness Z 0 are important factors that control the energy exchange between the Earth and atmosphere. The diurnal variation of ζ was analyzed for the PV plant and the reference site in Figure 2. The positive ζ during early morning and at night indicated stable conditions (ζ > 0) of near-surface atmosphere due to the stratification of the atmosphere caused by surface radiational cooling, while the negative ζ during the daytime from 9:00 to 17:00 indicated unstable conditions (ζ < 0), which resulted from the accumulation of surface solar radiation. Although similar characteristics of ζ were found at both sites, there were a few differences. During the daytime, the near surface air has higher instability with the average minimum ζ of smaller than −1.5 at 12:00 in the PV plant compared with the reference site (Figure 2a,b). The highly unstable conditions of near surface air may be caused by the strong convection exchange in the PV plant reported in the previous study [35]. For nocturnal hours, larger ζ in the PV plant indicated more stable near surface atmospheric stratification than that of the reference site, which may have resulted from the thermal insulation effect of PV modules on the ground surface at night. The magnitude of daily variation of ζ in the PV plant was larger than that of the reference site, as evident from the length of the whiskers. unstable conditions of near surface air may be caused by the strong convection exchange in the PV plant reported in the previous study [35]. For nocturnal hours, larger in the PV plant indicated more stable near surface atmospheric stratification than that of the reference site, which may have resulted from the thermal insulation effect of PV modules on the ground surface at night. The magnitude of daily variation of in the PV plant was larger than that of the reference site, as evident from the length of the whiskers.
Aerodynamic roughness 0 refers to the height with a wind speed of zero, and changes dynamically because its magnitude not only depends on the structure and shape of the surface roughness element, but also depends on the near-surface airflow conditions to a certain extent [42]. 0 in the PV plant with the value of 0.08~0.1 m was larger than the reference site with a value slightly larger than 0.04 m (Figure 2c), indicating that the PV underlying surface was rougher than the flat Gobi surface. More obvious diurnal variation of 0 in the PV plant also revealed a more complicated influence of panels on the airflow compared with the reference site. The averaged 0 was 0.089 m in the PV plant and 0.041 m in the reference site, respectively, and the order of magnitude was similar to the previous studies in grassland and shrub sandbags in New Mexico [43], but an order of magnitude larger than studies in the Heihe Desert (0.0045 m), Heihe Gobi (0.0045 m), and Dunhuang Gobi (0.0019 m).

Scaling Parameters
Friction velocity * , as an important dynamic characteristic parameter, characterizes the interaction between near-surface airflow and surface roughness element, and reflects the diffusion and transport capacity in all directions in space [44]. The diurnal variation of * was analyzed for the two sites (Figure 3a). The obvious diurnal characteristics of * indicated strong wind speed during the daytime, and greater * in the PV plant suggests that the airflow is subjected to stronger friction due to the presence of PV panels, which is Aerodynamic roughness Z 0 refers to the height with a wind speed of zero, and changes dynamically because its magnitude not only depends on the structure and shape of the surface roughness element, but also depends on the near-surface airflow conditions to a certain extent [42]. Z 0 in the PV plant with the value of 0.08~0.1 m was larger than the reference site with a value slightly larger than 0.04 m (Figure 2c), indicating that the PV underlying surface was rougher than the flat Gobi surface. More obvious diurnal variation of Z 0 in the PV plant also revealed a more complicated influence of panels on the airflow compared with the reference site. The averaged Z 0 was 0.089 m in the PV plant and 0.041 m in the reference site, respectively, and the order of magnitude was similar to the previous studies in grassland and shrub sandbags in New Mexico [43], but an order of magnitude larger than studies in the Heihe Desert (0.0045 m), Heihe Gobi (0.0045 m), and Dunhuang Gobi (0.0019 m).

Scaling Parameters
Friction velocity u * , as an important dynamic characteristic parameter, characterizes the interaction between near-surface airflow and surface roughness element, and reflects the diffusion and transport capacity in all directions in space [44]. The diurnal variation of u * was analyzed for the two sites (Figure 3a). The obvious diurnal characteristics of u * indicated strong wind speed during the daytime, and greater u * in the PV plant suggests that the airflow is subjected to stronger friction due to the presence of PV panels, which is not only affected by the surface roughness element, but also related to the wind speed and the strength of the mechanical turbulence.

Scaling Parameters
Friction velocity * , as an important dynamic characteristic parameter, characterizes the interaction between near-surface airflow and surface roughness element, and reflects the diffusion and transport capacity in all directions in space [44]. The diurnal variation of * was analyzed for the two sites (Figure 3a). The obvious diurnal characteristics of * indicated strong wind speed during the daytime, and greater * in the PV plant suggests that the airflow is subjected to stronger friction due to the presence of PV panels, which is not only affected by the surface roughness element, but also related to the wind speed and the strength of the mechanical turbulence.
The characteristic temperature * (Figure 3b) has the opposite trend related to * . At both the PV plant and reference site, negative * decreased after 8:00 and reduced to a minimum value of −0.65 in the PV plant and −0.47 in the reference site in the afternoon, respectively. Positive * gradually increased after 20:00. Across the whole day, greater * and * in the PV plant indicated that the PV panel strengthened the ability of the vertical transport of turbulence in contrast to the reference site.  The characteristic temperature T * (Figure 3b) has the opposite trend related to u * . At both the PV plant and reference site, negative T * decreased after 8:00 and reduced to a minimum value of −0.65 in the PV plant and −0.47 in the reference site in the afternoon, respectively. Positive T * gradually increased after 20:00. Across the whole day, greater u * and T * in the PV plant indicated that the PV panel strengthened the ability of the vertical transport of turbulence in contrast to the reference site.

Turbulence Strength and Turbulent Kinetic Energy
Averaged diurnal variations of the turbulence strength components I u (longitudinal wind), I v (lateral wind), and I w (vertical wind) show a similar behavior for the two sites ( Figure 4). Strong turbulence occurs under unstable conditions with strong wind speed and solar heating during the daytime, while turbulence at night is significantly weakened. For the horizontal turbulence strength I u and I v in the PV plant, the former was lower while the latter was higher than that of the reference site. The significant difference between I u and I v in the PV plant indicated that the horizontally turbulence strengths were affected by the plant layout. Southward inclination of the PV panel diminished the homogeneity of the north-south direction while it increased homogeneity in the east-west direction, resulting in weaker I u and stronger I v . For nocturnal hours with stable stratification conditions, convective exchange near the ground was suppressed, and turbulence was dominated by mechanical dynamics, and I u and I v in the PV plant was smaller than the reference site, mainly due to the lower wind speed compared to the reference site. In the PV plant, I w was significantly higher, which was a result of the increased vertical turbulence mixing caused by increased surface area and roughness compared to the reference site. The averages of I u , I v , and I w in the PV plant (reference site) were 0.29 (0.32), 0.34 (0.34), and 0.12 (0.10), indicating that the turbulence strength was dominated by horizontal flow. The turbulence strength of wind components was higher than that of the Badain Jaran Desert in previous studies [45].
was dominated by mechanical dynamics, and u and v in the PV plant was smaller than the reference site, mainly due to the lower wind speed compared to the reference site. In the PV plant, w was significantly higher, which was a result of the increased vertical turbulence mixing caused by increased surface area and roughness compared to the reference site. The averages of u , v , and w in the PV plant (reference site) were 0.29 (0.32), 0.34 (0.34), and 0.12 (0.10), indicating that the turbulence strength was dominated by horizontal flow. The turbulence strength of wind components was higher than that of the Badain Jaran Desert in previous studies [45]. Figure 5 presents the development of turbulence is subjected to wind speed ( Figure  5). Scatters of turbulence strengths and wind speed showed a similar trend for all three directions. The higher turbulence strengths corresponded with the lower wind speed. In the state of free convection with wind speeds lower than 2 m/s, the turbulence strengths of three directions decreased sharply with the increase in wind speed. The turbulence developed vigorously with a wind speed lower than 4 m/s while it tended to be steady with a wind speed greater than 6 m/s. The variation amplitude of the turbulence strength with the wind speed of the reference site was smaller than that of the PV plant, especially in the vertical direction.   Figure 5 presents the development of turbulence is subjected to wind speed ( Figure 5). Scatters of turbulence strengths and wind speed showed a similar trend for all three directions. The higher turbulence strengths corresponded with the lower wind speed. In the state of free convection with wind speeds lower than 2 m/s, the turbulence strengths of three directions decreased sharply with the increase in wind speed. The turbulence developed vigorously with a wind speed lower than 4 m/s while it tended to be steady with a wind speed greater than 6 m/s. The variation amplitude of the turbulence strength with the wind speed of the reference site was smaller than that of the PV plant, especially in the vertical direction. Turbulent kinetic energy (TKE, e) is an important indicator to measure the development and decline of turbulence, involving the transport of momentum, heat, and water vapor in the boundary layer. During the daytime, the TKE increased rapidly after 8:00 and reached a peak value around 16:00 at the two sites (Figure 6a,b), the lower TKE in the PV plant may be caused by a deficit of mean TKE throughout the panel array, which is related to the inclination of panels. At night, the TKE of the two sites ranged from 0.2 to 0.4, and the difference was very slight. With regard to the relationship between the TKE and wind speed (Figure 6c,d), the TKE increased rapidly with the increase in wind speed with a value lower than 4 m/s and tended to be stable when the wind speed was greater than 4 m/s in the two sites, and TKE at the PV plant was more discrete.

Momentum Flux and Heat Flux
The diurnal variation of momentum flux, as shown in Figure 7, as used to investigate the surface layer dynamics difference between the PV plant and the reference site. Momentum was positive across the whole day and showed obvious diurnal variation in the two sites. The momentum fluxes during the daytime were higher than that at nighttime, and high values of momentum fluxes corresponded to the strong convective turbulence and unstable conditions. PV panels enhanced the vertical interaction with the air flow above the panels, hence the wind shear as the main contributor for mechanical turbulence was stronger in the PV plant than the reference site. The peak daytime momentum fluxes in the PV plant and the reference site were 0.17 and 0.14 kg·m −1 ·s −2 , respectively. The lower momentum flux at night was associated with the weak development of turbulence, and the momentum flux in the PV plant was slightly higher than that of the reference site in the late night.

The Relationship between Turbulence Variance and Stability
The normalized standard deviations of the wind speed components, which describe the relationships between the near-surface variance and the fluxes, and represent the overall characteristics of turbulence, were applied to a variety of atmospheric diffusion modes. These relationships depend on the spatial and temporal in-homogeneities of the atmospheric motions and the atmospheric stability [46]. Additional, turbulent strength relationships may be affected by geographical location, dynamical, and the thermal characteristics of the atmosphere [47].
Here, we fitted the relationship between σ u /u * , σ v /u * , σ w /u * , and ζ, respectively, for the PV plant and the reference site according to Equation (11), and the best similarity coefficients are shown in Table 1. In the near neutral range (−0.1 < ζ < 0.1), turbulence is assumed to be isotropic and turbulent motion is decoupled from the surface, so the statistical characteristics of turbulence are independent of height and tend to be a constant. Under unstable conditions, the relationship between the horizontal and vertical wind speed normalized deviations and ζ can conform to the 1/3 law in the PV plant and the reference site ( Figure 9). Compared with the horizontal velocity, the vertical velocity variance had better similarity, and the dispersion of scatters was small, the main reason being that the energy of vertical movement is buoyancy; and the wind velocity variance similarity of the three directions of the PV plant was better than that of the reference site, which may be caused by the difference in air flow between the PV plant and reference site. Under stable conditions, the wind speed variance of the PV plant and the reference site can also basically meet the similarity theory, but the dispersion of scatters was greater than that of the unstable condition ( Figure 10). The similarity of the wind speed variance of the PV plant was better than that of the reference site, especially in the vertical direction.

Conclusions
To explore the impact of PV plants on the turbulence field, the observation of turbulence characteristics in a utility-scale PV plant and adjacent reference site over the Gobi area from 8 June to 31 October 2019 in Xinjiang, China were compared. As far as we know, this is the first time to explore the atmospheric turbulence characteristics of a utility-scale PV plant. Results indicate that in the study region, near surface air stratification is more unstable during the daytime while more stable at night in the PV plant than the reference site, ζ varies from −1.5 to 1.25 in the PV plant while it varied from −0.8 to 0.45 in the reference site in the averaged diurnal variation. In the PV plant, I u was lower, I v and I w were higher than the reference site during the daytime, I u and I v were lower, and I w was higher than the reference site at night. The averages of I u , I v , and I w in the PV plant (reference site) were 0.29 (0.32), 0.34 (0.34), and 0.12 (0.10), respectively. The turbulent kinetic energy in the PV plant during the daytime was lower than the reference site.
With regard to the turbulent fluxes, the momentum in the PV plant was higher than the reference site, especially during the daytime. Compared to the reference site, the lower heat capacity and increased surface area of the PV panels generated stronger convective conditions, resulting in higher sensible heat flux while the latent heat flux was reduced due to the wind resistance and shielding effect of PV panels in the PV plant.
Site specific turbulent strength relationships were developed for the PV plant and the reference site by analyzing wind components as a function of ζ. The turbulent components of wind followed the 1/3 power law under the unstable conditions and stable conditions in the two sites. Comparison of the derived relationships of the PV plant with the reference site indicates that the turbulent strength for vertical winds in the PV plant is relatively strong, especially under unstable conditions.