Specs:Cloud model

Cloud Model Specification -- The purpose of this document is to describe the specifications of the Cloud model in the climate module.

Overview
The cloud model for GridLAB-D is intended to replicate the effects that clouds produce on a feeder, allowing for the simulation of broader weather-related phenomena. The initial intended application is the cloud-induced power transients from PV panels. Other applications (such as heating of residences) may be possible but will likely have to be explicitly supported through modification of the existing implementation. It is expected that the transient nature of these cloud-based effects will manifest themselves as system transients and that a typical GridLAB-D study will likely simulate a single day or less. The model will support high temporal and spatial resolutions to accurately replicate these transitory effects.

Up until this point, the outputs provided by the climate module have been global in nature: one temperature, one humidity level, one wind speed, one irradiance. Clouds are definitionally not a global phenomena and will force two important changes to the existing GridLAB-D architecture: any object wishing to use the cloud model MUST have a defined spatial location, and irradiance can no longer be considered a global value. The former is already available in GridLAB-D but is rarely used; typically the only geographical information included in feeder models is line length. The latter, though, is a completely new concept in the weather model and will require new code development in the climate module but also modification to the code of objects currently using the solar information from the climate module.

Qualitative Specification
The cloud module will be responsible for:
 * Cloud Pattern Generation : Generating a realistic cloud pattern of an appropriate size to cover all points in the feeder that need irradiance data (at a specified minimum spatial resolution).
 * Cloud Pattern Density Variation : Varying the density of the cloud pattern to follow the measured TMY data currently in use.
 * Cloud Pattern Location Variation : Moving the cloud pattern across the feeder according to the wind speed and direction of the TMY data currently in use.
 * Irradiance Value Calculation : Modifying the sun irradiance value provided to solar PV installations and/or residences based on the coverage determined by the current state of the cloud pattern.

Cloud Pattern Generation
The cloud pattern will be generated using the square-diamond algorithm outlined in [6] which produces cumulus-like cloud patterns. The algorithm is a well-known and well-used technique for the fractal generation of clouds and is not difficult to implement. In addition to the square-diamond algorithm, a technique also outlined in [6] will be used to add varying density to the individual clouds in the pattern with the center of the clouds having a higher opacity than the edges.

The size of the pattern will be determined based on the coordinates of the objects that will be using the pattern. To this end, as part of its initialization function, the cloud model will collect the locations of all solar PV objects.

Cloud Size
Studies such as [11], [8], [4], [9] and [5] have examined the size of the cumulus clouds and as may be expected, the range of typical cloud sizes is quite large, ranging from 0.1 km2 to 13 km2. If we assume a square cloud, this translates into approximately 300 m to 3600 m per side. Choosing a cloud pattern pixel to represent 20 m to a side (400 m2), a typical cloud would take up 15 to 180 pixels per side. The limitations of the square-diamond algorithm for generating the cloud pattern require the pattern size to conform to $$2^n +1$$ pixels per side (where $$n$$ is an integer). To ensure the pattern can entirely contain a single cloud, the minimum cloud pattern tile generated must be 257 pixels to a side ($$n$$ = 8).

To determine if a larger cloud pattern is needed the maximum difference in latitude and longitude between all given locations in the feeder model must be found. Once these values are determined, they must be translated into metric values. Based on [1], latitude distances are constant throughout the globe where

$$ 1^\circ$$ latitude $$ = 111.32 $$ km

Longitude values, on the other hand, have variable spacing based on the latitude of the point of interest.

$$ 1^\circ $$ longitude $$ = 111.32 \cos ( latitude ) $$ km

Cloud Pattern Density Variation
Naturally occurring cloud coverage density changes in time and the cloud model must be able to accomplish this. The technique outlined in [6] offers a feasible method for doing so. The fractal pattern generated by the diamond-square technique produces a surface ($$xy$$ location with a magnitude). To transform this into a cloud pattern, this surface is sliced at an arbitrary elevation; all points on the surface below this cut are considered cloud and points above it are considered clear sky. To cover more of the sky with cloud the cutting elevation is increased and to create more blue sky the elevation is decreased. To alter the cloud coverage density, the cutting elevation will be adjusted up or down;it will be necessary to search for the correct cutting elevation to produce the target level of coverage. The cloud coverage pattern generation is heavily randomized and there does not exist a general relationship between cutting elevation and cloud coverage.

In the expected use case, the cloud coverage density signal will come from the "Opaque Sky Cover" signal found in TMY datasets. The cumulus clouds are typically opaque and the use of the TMY measurement will allow some degree of correlation between the cloud pattern and other weather-related phenomena being simulated. Since TMY data is only updated hourly and the expected use-case calls for much higher time resolution, the cloud coverage density signal will be interpolated using GridLAB-D's built-in interpolation functionality for TMY datasets. A default static value will be assumed if the specified field in the dataset is not available.

Cloud Pattern Location Variation
Clouds change in both shape and location over time. As in item 2, the wind speed and direction from the user-provided TMY dataset will be used to move the clouds over the feeder model. Cumulus clouds are low elevation (less than 1000 m) and though this is much higher than the 10 m at which the wind TMY data was collected, the wind speed at the higher elevation can be estimated using an terrain roughness measurement/estimate and the measured wind speed at 10 m as outlined in the European Wind Atlas [10].

$$ u_z = u_{ref} \frac { \ln \left( \frac {z} {z_0} \right) } { \ln \left( \frac { z_{ref} } { z_0} \right) } $$

where
 * $$u_z$$ is the wind speed at target height $$z$$
 * $$u_{ref}$$ is the wind speed at the measured height $$z_{ref}$$
 * $$z$$ is the target height
 * $$z_0$$ is the terrain roughness length
 * $$z_{ref}$$ height where wind speed was measured

Terrain roughness class is used to determine the roughness length and can be will be qualitatively determined based on descriptions in [?] but in many applications will likely land between 2.5 and 3.5. The equation relating roughness class and roughness length is [?]

$$ class = 3.91249 + \frac { \ln z_0 } { \ln 3.3333 } $$

(for roughness lengths over 0.03) where Assuming a roughness class of 3, the roughness length is 0.4 meters. Applying this to equation \ref{eq: wind height translation} seeking a target height of 1000m based on the TMY2 measurements at 10m, said equation becomes: $$ \begin{aligned} u_z & = u_{ref} \frac { \ln \left( \frac { 1000 } { 0.4 } \right) } { \ln \left( \frac { 10 } { 0.4 } \right) } \\ & = u_{ref} \frac { \ln 2500 } { \ln 25 } \\ & = 2.43 ~ u_{ref} \end{aligned} $$
 * $$class$$ is the terrain roughness class
 * $$z_0$$ is the terrain roughness length

A quick investigation of a few TMY2 datasets in Fargo ND, Goodland KS, Miami FL, and Springfield MO revealed a peak recorded wind speed of approximately 25 m/s. Translating this up to 1000 m we have a wind speed of almost 61 m/s. It is likely that this is not the highest wind speed recorded in all TMY2 data sets; a maximum wind speed at 1000 m of 80 m/s will be assumed.

The movement of the clouds will require that new cloud patterns be generated once the edge of the pattern reaches the edge of the feeder model. Based on the upcoming wind speed and direction, the cloud module will generate additional patterns that will be stitched onto the existing pattern. The diamond square algorithm easily accommodates this by defining the leading edge of the new pattern being added with values from the trailing edge of the existing pattern. The algorithm uses these as seed values as it generates the new pattern addition, ensuring that it seamlessly integrates with the former pattern.

Just as in the the case of the changing cloud coverage density, if the wind speed and direction signal is updated at a slower rate than the simulation, values will be interpolated for all times between the recorded data points using GridLAB-D's built-in interpolation functionality.

Irradiance Value Calculation
Once the cloud pattern has been created, cut at the appropriate elevation to define the coverage density correctly, and translated to the correct position over the feeder based on the wind speed and direction, each location in the feeder needing irradiance information must be evaluated against the current cloud pattern to determine if it is shaded or not. For simplicity, the cloud pattern is assumed to be located directly and immediately over the location of the PV panel such that the location of the sun in the sky does not impact the shading of the panel. If the cloud pattern shows the PV installation is covered, the shadow of the cloud will fall on the the installation, regardless of the time of day.

As part of this calculation, both the clear sky and shaded irradiance values must be defined. There are several papers that provide relatively simple algebraic models for "clear " sky solar irradiance such as \cite{Carroll:1985hg}, \cite{Atwater:1981hl}, \cite{Atwater:1978dp}, and \cite{Troen:1989wc}. Using judicious assumptions and TMY2 extraterrestrial solar radiation values, the non-clouded irradiance values can be calculated (see below). The shading value used for locations shaded by clouds will be a constant fraction of the calculated non-clouded value, optionally defined by the user. If the user does not define this value, a default value will be used.

Using references \cite{Atwater:1978dp} and \cite{Carroll:1985hg} as an example, the ``clear sky" solar radiation can be calculated as follows:

$$ F_d = I_0 \cos Z ((P_R P_A)_d - a_w) P_a (P_c)_d $$

where
 * $$F_d$$ is the direct horizontal radiation reaching the ground
 * $$I_0$$ is the extraterrestrial solar radiation, assumed constant at 1353 W/m2
 * $$Z$$ is the zenith angle, a function of hour of day and day of year and calculated using the existing solar angles library based on date, time and location on the globe.
 * $$(P_R P_A )_d$$ is the transmissivity after attenuation due to Rayleigh scattering and gaseous attenuation, defined as $$1.041 - 0.15 \sqrt[]{\frac{p ~ 949 \times 10^{-5} +0.051}{M}}$$ where
 * $$p$$ is the atmospheric pressure, likely assumed constant.
 * $$M$$ is the optical path length defined as $$\frac{p}{35 \times 101.3} \sqrt[]{1224 cos Z^2 +1}$$
 * $$a_w$$ is the absorptivity of water vapor defined as $$0.077\left(\frac{u}{M}\right)^{0.3}$$ where
 * $$u$$ is the precipitable water in cm, assumed to be a moderate value of 3 to 4.5 cm
 * $$P_a$$ is the transmissivity after aerosol attenuation, defined as $$e^{-\alpha M}$$ assumed to be 1 as in [3]. ($$\alpha$$ is the aerosol volume absorption coefficient.)
 * $$(P_c)_d$$ is the transmissivity after cloud attenuation. For this model we will assume a constant factor of 30 % (based on estimates in [1], figure 12) which can be over-ridden by the user. This factor will only be used when a given location is calculated to be covered by clouds.

The above irradiance value must be calculated at every solar installation location and will be impacted by the clouds moving over the simulated feeder. The above explanation also shows that very few assumptions need to be made that those that do ($$a_w$$, $$P_a$$, etc) have justifiable values readily available. It is not anticipated that implementing the required functions (including the above) will be particularly difficult.