A Method for Projecting Cloud Shadows onto a Central Receiver Field to Predict Receiver Damage

. This work demonstrates methods of mapping high-spatial-resolution direct normal irradiance (DNI) data from satellites, Total Sky Imagers (TSIs), and analogous data sources onto a heliostat field for characterizing the spatial and temporal variation of the incident flux on a central receiver tower during cloud transient events. The mapping methods are incorporated into an optical software module that interfaces with CoPylot–SolarPILOT’s python API– to provide computationally efficient optical simulation of the heliostat field and the solar power tower. Eventually, this optical model will be incorporated into optimization models whereby a plant operator can understand the effects of cloud transient events on overall power production and receiver lifetime due to creep-fatigue damage and therefore make better informed decisions about receiver shutdown events. By more accurately modelling the effects of cloud events on receiver flux maps, this work may determine the magnitude and frequency of thermal cycling on receiver tubes and panels using actual or realistic cloud shapes instead of averaged DNI values–which may undercount the total cycle number. This work may also prevent unnecessary plant shutdowns due to overly precautionary control strategies and characterize the relative impact of various cloud types on receiver life. We plan to eventually integrate this methodology into the System Advisor Model (SAM) to improve performance model accuracy during periods of cloudiness. In this paper, we demonstrate generating DNI maps and mapping them to a solar field in CoPylot using 10 m resolution data from publicly available Sentinel-2 satellite data over the Crescent Dunes plant.


Introduction
In order to maximize the revenue of a central receiver plant and drive down the levelized cost of electricity (LCOE) both during the design and the operation of the plant, the plant operator must consider several factors including the lifetime of the plant and component maintenance, the availability of the plant, and other factors such as the current price of electricity and thermal storage if applicable.One software package used to assist plant operators is the System Advisor Model (SAM) [1] developed at the National Renewable Energy Lab (NREL).SAM contains a module for concentrating solar plants with and without thermal energy storage and integrates a mixed-integer linear program (MIP) to optimally dispatch electricity during times of peak load (and peak price) to the grid while respecting the various physical and economic constraints associated with the plant.Compared to baseline dispatch strategies, the approach implemented in SAM can improve plant profitability by 5-20% [2].A version of SAM has also been developed for real-time operations, including during periods of variability [3].However, SAM relies on historical DNI and weather data points that are either averaged spatially over the entire plant or recorded in one location; there has yet to be a study on how overall plant revenue is affected when spatially resolved DNI data and additional cloud data-such as cloudtype classifications-are used to inform the optimization model during transient events.
Historically, accounting for cloud transient events in both the design and operation of a central receiver plant has been managed by using large safety margins on the receiver tubes and conservative control strategies where the heat transfer fluid (HTF) flow rate would be set to match clear-sky conditions during short-duration large cloud events, which may result in a significant loss in revenue for the plant in both the capital expenditures (CAPEX) from the construction of the plant due to higher material costs and the operating expenses (OPEX) due to lost electricity revenue.The large safety margins on receiver tubes were intended to prevent premature failure due to creep-fatigue damage-assuming a design life of 30 years, and the control strategies were designed to prevent catastrophic failure due to over temperature conditions once a cloud front leaves the field and DNI levels rapidly approach clear-sky conditions.Both design and operation decisions were in large part due to a lack of design standards and test data for CSP plants to account for the creep and fatigue damage due to thermal cycling on the receiver and a lack of detailed high spatial and temporal resolution cloud transient data and predictive modelling.
In a report by Kistler [4] using historical weather data from the Solar One plant, it is clear why such large safety margins on receiver tubes were used in practice.The effects of small cloud transients that only crossed portions of the field were disregarded entirely and only two locations in the field measured DNI; such undercounting necessarily resulted in the use of a large safety factor associated with the allowable number of fatigue cycles, as a single nonhomogenous cloud-front may induce multiple thermal cycles on a given receiver tube.Additionally, in both Kistler's report and a paper by Narayanan et.al [5], they used a modified version of the ASME B&PV Code Case N47 developed for nuclear pressure vessels to calculate the creep-fatigue damage due to thermal cycling.Although appropriate for nuclear applications, this modified approach still likely overpredicts the damage done to a receiver due to cloud events.Other project designs around this period required the receiver to be modular and to anticipate tube replacement as an expected operating cost due to the uncertainty related to calculating the creep-fatigue damage [6].
Likewise, the Solar Two plant employed eight photometers that were pointed towards the receiver panels to detect changes in reflected solar flux due to clouds [7]; with such few datapoints and limited scope, it would be impossible to obtain accurate predictions of the duration and intensity of cloud events in future timesteps, necessitating more conservative control strategies.As such, the control algorithm for the heat transfer fluid (HTF) that was employed in Solar Two included a Cloud Standby (CSB) mode, whereby the plant would increase the HTF flow rate and thereby decrease the outlet temperature to prevent an over temperature condition once the cloud front leaves the field.Even if one part of the field experienced cloud cover due to cumulus clouds, the flow rate would be set to the maximum clear-sky flow rate for both flow paths on either side of the receiver.The use of this control algorithm was also justified with simulations that demonstrated that the reduced temperature variations during cloud events would yield less fatigue damage on the receiver tubes.While optimizing between component cost and power production would be ideal, the lack of high-quality cloud data and modelling necessitated a precautionary control strategy.
Recently, work by Schwager et al. [8] investigated strategies to model the spatial variation of transient cloud events and the resulting receiver flux.They developed complementary aim-point and HTF control strategies, where it was found that not accounting for the spatial variations due to cloud transients can lead to overpredicting the annual solar yield by 2-4% and control strategies that do not account for local spatial DNI variation can result in exceeding the temperature limits of the HTF.Another paper by Crespi et al. examines the effects of different synthetic cloud events on receiver efficiency and finds a 1% increase in receiver efficiency using aiming strategies that account for spatial DNI data [9].Rangel et.al used actual DNI data from CIEMAT's Plataforma Solar de Almeria (PSA) and the Solar Tower Ray Tracing Laboratory (STRAL) tool [10] to find the flux on the receiver and subsequently the creep-fatigue damage using the rainflow-counting algorithm and a simplified elastic stress analysis [11].These studies demonstrate the utility of using spatially resolved DNI data to increase receiver efficiency and more accurately determine annual solar yield.However, these methods have yet to be integrated into annual production simulations that optimize between the increase in electricity revenue due to less conservative control strategies and the subsequent reductions in receiver component lifetime in order to maximize long term revenue.This paper aims to demonstrate methods whereby cloud shadows may be projected onto the heliostat field using the CoPylot Python API [12].While CoPylot is not yet optimized for multi-threading, in the future it may approach the computing time of SolarPILOT, which can simulate the incident flux on a central receiver in <105 ms [13].This would allow for the fast-optical characterization of large sets of cloud data or for real-time operations.Sources of cloud transient data that may be used and explanations of the software methods in this paper are described in the following sections.

Data Sources
The three primary sources of readily available DNI data that are spatially and temporally resolved enough to be applicable to a solar power tower field are: Total Sky Imagers (TSIs), correlated DNI data from photovoltaic (PV) panels or sensors on heliostats, and satellite imagery.Additional sources of data may come from radar, lidar, or drones, but have yet to be validated in the literature on heliostat fields.
TSIs or cloud cameras use a fisheye lens and a shadow arm to image the entire sky at ~10 m spatial and 30 s temporal resolution and have successfully been used for DNI nowcasting [14].Nouri et al. used 4 TSIs to calculate cloud height and estimate cloud transmittance by generating 3D cloud voxels [15], [16].Other methods that use multiple TSIs to generate shadow masks from 3D cloud shapes and project them onto heliostat or PV fields are documented in [17] and [18].Examples of publicly available datasets include those from the Atmospheric Radiation Measurement (ARM) Climate Research Facility [19] and a benchmark dataset released by Coimbra et al. [20].To our knowledge there are no publicly available TSI datasets that consist of two or more cameras for implementing 3D stereoscopic cloud shape algorithms, which means that cloud height must be estimated using some other method.Some CSP plants are opting to use battery-powered heliostats powered by PV panels in-stead of using an underground wired power system as this reduces cost and may potentially avoid issues during electrical storms [21].PV panels are powered by global horizontal irradiance (GHI), which can be correlated to DNI using a variety of models if the diffuse and ground reflected radiation are known.The Megalim plant operated by BrightSource Energy utilizes this method.Given the high spatial and temporal resolution that results from this approach, this method is likely the most accurate data source outside of an array of pyrheliometers for determining the DNI over a field [22].Unfortunately, the range is limited to where the hardware can be placed, meaning other sources are necessary for accurate nowcasting and new plant site determination.
Satellite DNI data is mostly available at the mesoscale (~1-4 km) from geostationary satellites at short time resolutions (5-30 mins) [23].For example, the National Solar Radiation Database (NSRDB) has 5 min time resolution and 2 km spatial resolution (as of 2018) for anywhere in the United States using sensors from a collection of satellites such as GOES and MODIS to calculate DNI, GHI, and a variety of other variables [24].Similar systems to evaluate DNI exist for other countries such as EUMETSAT's MSG satellite; in [25], Sirch et al. used the SEVIRI instrument aboard MSG to determine and forecast the DNI reduction from both water and cirrus ice clouds with a 15 min temporal resolution and 3 km spatial resolution.
To achieve higher spatial resolutions for satellite data, constellations of lower orbiting sun synchronous polar satellites are often used at the expense of temporal resolution (re-imaging the same location every 1-12 days).The two publicly available satellite datasets with the highest spatial resolution are the LANDSAT and Sentinel-2 constellations with resolutions of 30 and 10 meters, respectively.Private satellite datasets, often used for surveillance or land monitoring purposes, offer higher spatial (0.3-1 m) and temporal resolution (<1 day) at the expense of less available spectral bands.While algorithms such as those described in [26] are used to identify clouds and their shadows in Sentinel-2 imagery, for example, there are no physical models to calculate DNI from these datasets.Nevertheless, as this is the only publicly available dataset with comparable resolution to heliostat PV-DNI correlations, this resource is used by the authors in the following sections.

Description of software
The two methods used to map cloud shadows onto the heliostat field will be referred to as the Raster and Vector methods, respectively.The Raster method converts the DNI data into a numerical array using the numpy library in Python and scales/masks the array onto the heliostat field by adjusting each heliostat's soiling or optical efficiency parameter.The array can either be a binary array (corresponding to turning the DNI on or off) or an array of floats normalized between 0 and 1 (corresponding to reducing the DNI).Eq. ( 1) gives the values in the array, where ̅ is the position vector, usually provided in Cartesian coordinates with respect to the tower: In [25], Sirch et al. uses the Lambert-Beer Law and the "strict definition" of DNI that is often used in radiative transfer models to derive an approximate equation for the DNI over a given point described in Eq. ( 2) where  0 is the extra-terrestrial solar irradiance integrated over the entire spectrum (~1367 W/m 2 ) [27], τ  is the optical thickness of the atmosphere, which depends on factors such as water vapor content (note that Sirch et al. does not account for aerosols, but this term may be added for better accuracy τ  ), τ  represents cloud optical thickness (sometimes referred to as cloud optical depth), and  0 is the solar zenith angle.Sirch et al. also notes that cloud optical thickness is often split between high-altitude thin ice clouds ( , ) and low-altitude water clouds ( , ), where the low-altitude water clouds often have such high optical thickness values that they reduce DNI to below levels usable for CSP (<200 W/m 2 ).As such, using a binary cloud shadow mask is often appropriate and simple to implement in the presence of low-altitude water clouds.
The Vector method draws a set of polygons around cloud shadows on the field.The vertices of the polygons are then passed to CoPylot and used to adjust the soiling parameter of the heliostats within them.The polygons are generated using the Sci-kit image library [27], which uses a special 2D case of the "marching cubes" algorithm [28] called the "marching squares" algorithm, the details of which can be found in [29].To generate the polygons from DNI data, a threshold value must be passed; in the case of binary cloud shadow masks for low-altitude water clouds, the selection of a threshold value is straightforward.However, in the presence of thin cirrus clouds, which reduce DNI gradually based on optical thickness, multiple threshold values need to be selected depending on the bit depth of the DNI values, which offsets any computation time savings compared to the Raster method.
Figure 1 demonstrates a synthetic cloud shadow mask composed of a 500 x 500 m cloud moving across the SAM default field and its effect on incident receiver flux using the Raster method.The cloud event duration is 15-minutes simulated using 30 s time steps at a 1 m spatial resolution.Figure 2 shows histograms of the computation time for projecting the cloud shadow in Figure 1 onto the field as it moves from West to East for the two methods, respectively.Using one core on an i7-7700HQ 2.8 GHZ CPU, the total time for calculating the receiver flux for the cloud event was 190 s and 175 s for the Raster and Vector methods, respectively.

Satellite Data Use Case
To demonstrate how satellite data or an analogous data source may be used to project shadows onto the heliostat field, data from the Sentinel-2 satellite was retrieved over the Crescent Dunes plant using the Google Earth Engine (GEE) and GEEMap libraries [30].This demonstration only shows how the software may be used and is not meant to accurately estimate DNI.Sentinel-2 contains 13 spectral bands and a variety of data products such as aerosol optical thickness (AOT) and water vapor content, which may be used in a physical model to more accurately calculate DNI [31].However, for simplicity, precomputed cloud and shadow masks are taken from [32]. Figure 3 demonstrates how the precomputed cloud shadow mask is projected onto the heliostat field in Cartesian coordinates using a UTM projection library.As the Sentinel-2 time resolution is not sufficient for capturing transient cloud events, one method from Crespi et al. [9] can be employed where the cloud shadows may be translated in time both before and after the local image time using the local wind velocity vector and the Hellman equation for altitude correction.Local wind velocity vectors may be taken from a nearby weather station or from GOES satellite data products such as the Wind Integration National Dataset (WIND) Toolkit [33].

Conclusion
With optimization for multithreading, the methods described in this paper may be used to compute the incident flux on the receiver during cloud transient events at computational speeds similar to SolarPILOT (<105 ms), which would be suitable for characterizing the effects of large datasets of clouds or to examine real-time plant operations.While data from Sentinel-2 was demonstrated on an external cylindrical receiver, there is no limitation to the data sources that may be used to compute DNI maps or the types of receivers that may be modeled, so long as they are within SolarPILOT.This paper also highlights the need for physics-based DNI models for polar orbiting satellites such as Sentinel-2 and public repositories for stereoscopic TSI imagery, as these resources would allow for more accurate assessment of the local solar resource when selecting sites for future plants.

Figure 1 .
Figure 1.(a-b) synthetic 2500 m 2 square cloud moving over heliostat field from west to east using Raster method (c-d) receiver incident solar flux map from CoPylot.

Figure 2 .
Figure 2. (a) Raster Method histogram of computation time for each time step for square cloud moving across field (b) Vector Method histogram (Note: to account for multiple small cloud shapes, longer computing time is required).