- Email: [email protected]
Multivariate simulation for assessing the joint wind and ice hazard in the United States
Journal of Wind Engineering & Industrial Aerodynamics 184 (2019) 436–444
Contents lists available at ScienceDirect
Journal of Wind Engineering & Industrial Aerodynamics journal homepage: www.elsevier.com/locate/jweia
Multivariate simulation for assessing the joint wind and ice hazard in the United States Hung Nguyen Sinh a, *, Franklin T. Lombardo b, Chris Letchford c a b c
BMT Fluid Mechanics Ltd, New York City, NY, USA Civil & Environmental Engineering Department, University of Illinois at Urbana-Champaign, Urbana, IL, USA Civil & Environmental Engineering Department, Rensselaer Polytechnic Institute, Troy, NY, USA
A R T I C L E I N F O
A B S T R A C T
Keywords: Joint hazard Multivariate simulation Performance-based design Wind/ice hazard Freezing rain Ice load
Multi-hazard assessment is receiving increasing attention as often commonly occurring hazards may contain a second or even a third associated hazard. Generally, multi-hazards could be joint hazards (i.e. these hazards occur at the same time), for example, a hurricane usually comes with strong wind, storm surge, and heavy rainfall. They can also be one hazard preceding another, for example, an earthquake generating a tsunami. The jointly occurring hazards in some cases can be more devastating than a single hazard. This paper will describe the assessment of two jointly occurring hazards, wind and ice formation by using an advance numerical simulation. First, multivariate (hourly wind, temperature and precipitation) simulation for many years in the Midwestern US is undertaken. By setting a condition for freezing rain, these parameters were input into a simple ice-accretion model to estimate ice-thickness. A joint hazard curve for wind and ice was then constructed and compared to the existing prescribed treatment of these two hazards in design loading guidelines, such as ASCE 7 in the United States, ISO 12494 (International Standard); and the results from our previous model.
1. Introduction For locations in the US, Canada, and other higher latitude countries, wind and ice (e.g., freezing rain) are hazards whose properties are of interest for design of transmission lines and other energy infrastructure (e.g., wind turbines). Freezing rain is one of the major weather hazards in the US causing annual losses of $314 million (2000 US$) during 1949–2000 (Changnon et al., 2003). As pointed out in our previous paper (Nguyen Sinh et al., 2016), although there have been numerous scientific and engineering studies of freezing rain, few considered the combined effect of the wind hazard and the ice hazard occurring simultaneously. Nguyen Sinh et al. (2016) used a simple approach to deal with the joint wind and ice hazard during freezing rain events in which the wind speed and ice data based on actual observations over a long period (~300 years using the superstation approach (Peterka, 1992)) were simulated separately to assess the joint wind and ice hazard. The details for this approach is shown in Fig. 1. Firstly, the meteorological data from measurement stations in the Midwestern US were collected from the ISH (2001) database. Quality control measures to remove any spurious (i.e. not physically plausible) data points were undertaken (e.g. freezing rain that occurred at high temperatures (>0 Celsius degree), unrealistic wind
speeds and precipitations (wind speed values > 30 m/s, precipitations > 15 mm/h for winter and >40 mm/h for summer, …). Then the simple ice-accretion model (Goodwin et al., 1983; Jones, 1998) was used to estimate the ice thickness for each freezing rain event. The superstation approach (Peterka, 1992) was then used to group these stations into a “superstation”. The wind speeds (occurring during freezing rain) were fitted by a Weibull distribution while the Generalized Pareto (GP) distribution was used to fit the ice thickness using the peaks-over-threshold method. The wind speeds and ice thickness were then simulated separately and their joint histogram was established to create the joint wind and ice hazard. Details and results of this approach are described in Nguyen Sinh et al. (2016). Even though this simplified approach was conservative, it still confirmed that the ASCE 7 (2016) design values for 50 year Mean Recurrence Interval (MRI) for joint wind-ice hazard belong to much higher MRI (~150 years). In other words, the 50 year ice thickness values from ASCE 7 (2016) must be combined with much lower wind speed than the concurrent 3-s gust wind speed in the “Equivalent radial ice thicknesses due to freezing rain with concurrent 3-s gust speeds, for a 50-year Mean Recurrence Interval” map from ASCE 7 (2016) and vice versa or both these values must be reduced to match the 50 year joint hazard
* Corresponding author. E-mail addresses: [email protected] (H. Nguyen Sinh), [email protected] (F.T. Lombardo), [email protected] (C. Letchford). https://doi.org/10.1016/j.jweia.2018.12.012 Received 8 October 2018; Received in revised form 12 December 2018; Accepted 13 December 2018 Available online 20 December 2018 0167-6105/© 2018 Elsevier Ltd. All rights reserved.
H. Nguyen Sinh et al.
Journal of Wind Engineering & Industrial Aerodynamics 184 (2019) 436–444
associated with each freezing rain event. Finally, the joint wind and ice hazard results are established and compared with the results from the simple model (Nguyen Sinh et al., 2016) as well as ASCE 7 (2016)) and ISO 12494 (2017) design values. Since the freezing rain data was only used to set conditions for freezing rain, actual missing data during freezing rain is no longer problematic. Although the simulation procedure in this approach is very different to the one in the simple model (Nguyen Sinh et al., 2016), the joint hazard curves were established using the joint probability distribution of ice and wind which is similar to the simple model (Nguyen Sinh et al., 2016). The difference between this joint hazard approach and ASCE 7 was discussed in detail in Nguyen Sinh et al. (2016). 2. Data and methodology The meteorological data used for this study is collected from the National Climate Data Center (NCDC) within NOAA's Integrated Surface Global Hourly (ISH) database, which can be found at ftp://ftp.ncdc.noaa. gov/pub/data (ISH, 2001). Table 1 shows the station numbers and names, locations and the associated time periods used for this study. Fig. 3 shows that these 8 stations have very similar ice thickness for the 50 year MRI according to ASCE 7 (2016). They are the same data set used in our previous paper (Nguyen Sinh et al., 2016). Table 1 also shows the annual frequency (λ) of freezing rain for each station determined from the total number of freezing rain events divided by the total number of years of data for that station. It is seen that all have approximately 2 events per year. Noting that, as indicated in Nguyen Sinh et al. (2016), the ‘hourly’ wind speed data and the ‘hourly’ temperature in this paper are actually the 2-min mean prior to every top-of-hour. The detailed procedure employed for this simulation is shown in Fig. 4. Firstly, all the data (wind speeds, precipitations and the weather observations) from the 8 Midwestern US stations (Table 1) were collected from the ISH (2001) database. These historical data are used to estimate parameters for the simulation discussed in more detail below. The simulation process is time dependent within each variable and the three variables (wind speed, temperature and precipitation) are inter-dependent. For instance, one might expect that the ‘hourly’ temperature during a wet hour in summer is more likely to be lower than during a dry hour (Richardson, 1981). Also the correlation between wind speed and precipitation, or wind speed and temperature, may not be obvious but still present. The main idea is to consider precipitation as the primary variable and then conditioning the other two variables (wind speed and temperature) for a given hour on whether that hour is wet or dry. A Markov chain is used to simulate the precipitation (P) (Richardson, 1981). Then multivariate simulation (wind, temperature) for that given hourly precipitation is achieved by using a multi-variable Markov generating process (Matalas, 1967). By establishing the condition for freezing rain, these parameters were input into a simple ice-accretion model (Goodwin et al., 1983; Jones, 1998) to estimate the ice
Fig. 1. The procedure for the simple model (Nguyen Sinh et al., 2016). (GP: Generalized Pareto distribution).
contour line developed in Nguyen Sinh et al. (2016) and shown in Fig. 2. However, one of the significant problems for this approach was missing observations that commonly occur during freezing rains, leading to an incomplete dataset to drive simulations. In this paper, a more sophisticated model is introduced in which all the meteorological variables needed for an ice-accretion model (Goodwin et al., 1983; Jones, 1998); total hourly precipitation, hourly temperature and mean hourly wind speed, are simulated simultaneously such that the variables are time dependent and their inter-dependence established. Conditions for freezing rain are then set and freezing rain events estimated from the simulated meteorological data. These freezing rain event meteorological variables are input into a simple ice-accretion model (Goodwin et al., 1983; Jones, 1998) to estimate the ice thickness
Table 1 Information on stations used for Midwestern, US. Station Number
Station Name
725450 725460 725510
Cedar Rapids Des Moines Lincoln Municipal Grand Island County Sioux City Muni Madison Sioux Fall Minneapolis/ St. Paul
725520 725570 726410 726510 726580
Fig. 2. Joint gust wind and ice hazard for Midwestern US superstation considering wind effects after freezing rain events (Nguyen Sinh et al., 2016): a) Joint hazard curves; b) Marginal distribution for ice thickness during freezing rain events; c) Marginal distribution for wind speed during freezing rain events. 437
Length of recorded data (years)
FR frequency (λ)
2010 2010 2010
33 38 38
2.03 2.16 1.63
1978
2010
33
2.18
IA
1973
2010
38
1.53
WI SD MN
1973 1973 1973
2010 2010 2010
38 38 38
1.68 1.97 1.92
State
Time period Start
End
IA IA NE
1978 1973 1973
NE
H. Nguyen Sinh et al.
Journal of Wind Engineering & Industrial Aerodynamics 184 (2019) 436–444
Fig. 3. ASCE 50 year MRI ice map (adapted from ASCE 7 (2016)).
Markov Chain (Richardson, 1981) is employed. Basically, the occurrence of wet or dry hours was modeled by a first-order Markov Chain (Bailey, 1964) and then on the rainy hour, the best-fit distribution for precipitation intensity is used to estimate the amount of precipitation in that hour. In the past, Gabriel and Neuman (1962), Caskey (1963), Weiss (1964), Hopkins and Robillard (1964) and Haan (1976) used Markov chain models for simulating daily precipitation. “Smith and Schreiber (1973) found that a first-order Markov chain was superior to a Bernoulli model for describing the occurrence of wet or dry days in southeastern Arizona, US” (Richardson, 1981). This study will use a similar first-order Markov chain approach. Fig. 5 shows the detail of applying the Markov chain to generating hourly total precipitation. Let W (wet) and D (dry) denote for the state of an hour, P (W/D) is the probability of a W hour given a condition D at the previous hour. Similarly for the definition of P (D/D), P(D/W) and P(W/ W) (Richardson, 1981). Then: P(D/W)¼1-P(W/W)
(1)
P(D/D)¼1-P(W/D)
(2)
Therefore the time series of precipitation is fully defined given P (W/ W) and P (W/D) and the wet or dry state at the first hour. In other words, if these conditional probabilities can be estimated from the historical data, by assuming the initial state for the first hour, condition W or D of all hours in a time series can be generated. Based on the historical events, these probabilities are indicated in Table 2. Another probability distribution is used to characterize the amount of precipitation at each wet hour. Inverse cumulative density function (CDF) (Benjamin and Cornell, 1970) is used to determine the best-fit distribution for the precipitation intensity. This approach assumes that the CDF values (u) of a variable (X) are initially uniformly distributed over [0 1], then different types of probability distributions are selected. For each type of probability distribution, the variable X'i is calculated by taking the inverse CDF of u (i.e. X'i ¼ CDF1 i (u)). All of the points (X, X'i)
Fig. 4. The procedure for the sophisticated model.
thickness. The joint histogram and joint exceedance probability of wind and ice thickness are then constructed to establish the joint hazard curves for wind and ice. 2.1. Precipitation simulation Many models have been proposed to simulate precipitation: Buishand (1978), Chin (1957) and Gabriel and Neuman (1962). In this research, a 438
H. Nguyen Sinh et al.
Journal of Wind Engineering & Industrial Aerodynamics 184 (2019) 436–444
Fig. 6. Best-fit distribution for hourly total precipitation at Des Moines, IA.
temperature based on the precipitation that has been simulated above. First the residual variables were calculated by removing the mean and standard deviation of each variable as expressed in Equation (3). xi ðpÞ ¼ Fig. 5. The procedure for the Markov chain.
Cedar Rapids Des Moines Lincoln Municipal Grand Island County Sioux City Muni Madison Sioux Fall Minneapolis/St. Paul
State
P
Best fit distribution for P
P(W/ W)
P(W/ D)
P(W)
IA IA NE NE
0.838 0.78 0.794 0.809
0.025 0.025 0.021 0.019
0.135 0.102 0.091 0.089
Weibull Log normal Log normal Log normal
IA WI SD MN
0.789 0.767 0.795 0.768
0.021 0.029 0.02 0.024
0.092 0.11 0.09 0.095
Weibull Weibull Weibull Log normal
(3)
where XðpÞ and σ ðpÞ are the mean and standard deviation and xi(p) is the residual component for variable p, at time ‘i’. Then the weakly-stationary generating process suggested by Matalas (1967) was used to generate a residual series of all variables.
Table 2 Estimated probabilities for Precipitation. Station Name
Xi ðpÞ XðpÞ σ ðpÞ
xi ðpÞ ¼ Axi1 ðpÞ þ Bεi ðpÞ
(4)
where xi ðpÞ, xi1 ðpÞ are residual variables at time ‘i’ and ‘i-1’ of variable p; εi ðpÞ are independent random components that are normally distributed with a zero mean and a unit variance associated with variable p. The A and B matrices are defined such that the new sequences have the desired serial correlation and cross-correlation coefficients. Equation (4) implies that the residuals of variables are normally distributed and that the serial correlation of each variable may be described by a first-order linear autoregressive model (Matalas, 1967). The A and B matrices can be determined from the following equations.
for each distribution are then plotted and the distribution closest to the line of slope 1 is the best-fit distribution. Fig. 6 shows that Weibull is the best-fit distribution for the precipitation in Des Moines, IA within six types of selected distributions (Normal, Lognormal, Extreme Type 1 (ET1), Extreme Type 2 (ET2), Weibull (WB) and Exponential (EXP)). Among methods for finding best-fit probability distributions, the inverse CDF method (Benjamin and Cornell, 1970) was selected because it provides a simple, graphical technique that can be generalized for any distribution form without the need for specific types of plotting papers and furthermore can be easily programmed.
A ¼ M1 M 1 0
(5)
T BBT ¼ M0 M1 M 1 0 M1
(6)
where the superscripts 1 and T denote the inverse and transpose of a matrix respectively. M0 and M1 are the lag-0 and lag-1 covariance matrices. For the case of 2 variables (p ¼ 2), M0 and M1 are expressed by following equations M0 ¼
2.2. Lag-one Markov generating process
M1 ¼
A lag-one Markov process - a multivariate weakly-stationary generating process (Matalas, 1967) - was used to simulate the wind speed and
1 ρ0 ð2; 1Þ
ρ0 ð1; 2Þ
(7)
1
ρ1 ð1; 1Þ ρ1 ð1; 2Þ ρ1 ð2; 1Þ ρ1 ð2; 2Þ
(8)
where ρ0 ði; jÞ and ρ1 ði; jÞ are the lag-0 and lag-1 cross-correlation coeffi439
H. Nguyen Sinh et al.
Journal of Wind Engineering & Industrial Aerodynamics 184 (2019) 436–444
cient between variables ‘i’ and ‘j’, respectively. In this application, ‘i’ will represent wind speed and ‘j’ will represent temperature.
3. Analysis Using Des Moines, IA data this section illustrates the simulation procedure in more detail. The probabilities P (W/W) and P (W/D) were estimated from historical data and are shown in Table 2. The initial state (wet or dry) was determined based on a Bernoulli trial with P (W) ¼ 0.102 and P (D) ¼ 1- P (W) ¼ 0.898. All other states (at each following hour) were determined by a Markov chain with P (W/ W) ¼ 0.78 and P (W/D) ¼ 0.025. At all wet hours, the amount of precipitation was determined using the best-fit distribution for the precipitation in Des Moines, IA. Although, the best-fit distribution varies for each location (Table 2) and both Log normal and Weibull distributions are best fit for data, the Log normal distribution usually predicted slightly higher values than historical data (i.e. observed data). For consistency and conservatism, a Log normal distribution was used for all locations.
2.3. Multivariate hourly wind speed and hourly temperature conditioned on hourly precipitation series In this section, application of the multi-variable lag-one Markov generating process for simulating ‘hourly’ wind speed and ‘hourly’ temperature will be discussed. Noting that, Matalas' model (Matalas, 1967) (i.e. Equation (4)) works best when variables (and hence residual variables) are normally distributed. For the ‘hourly’ temperature, Lombardo et al. (2014) showed that temperature did indeed follow a normal distribution for each season. For the ‘hourly’ wind speed, the Weibull distribution is most often the best-fit distribution. However, using the Inverse CDF method (Benjamin and Cornell, 1970), Fig. 7 shows that a log-normal distribution is as good a fit in the region of 0–15 m/s which is the typical range of wind speeds during freezing rains (Roberts and Stewart, 2008) for this and other locations. Hence, it is reasonable to apply Matalas' model (Matalas, 1967) with the first parameter being the natural log of wind speed (i.e. X (1) ¼ ln (V), where V is the ‘hourly’ wind speed) and the second parameter being temperature (i.e. (X (2) ¼ T). There are hours with calm (i.e. no wind) that lead to undefined natural log, hence it is necessary to add a value of 0.01 (m/s) to these values. Since the accuracy resolution for anemometers is 1 knot (~0.514 m/s), this adding has no effect on the statistical parameter estimating. In addition, monthly analysis (i.e. these statistical parameters are estimated for each month) is undertaken to ensure these variables (wind speed and temperature) follow the normal distribution. Fig. 8 outlines the detailed simulation procedure. In summary, first, a time series for precipitation is constructed using a Markov chain. The initial values (i.e. at first hour) were selected as the values at the beginning of the first available year from historical data. Then at each following hour, by using the Matalas' model (Matalas, 1967) (or using Equation (4)) the residual vector of two parameters are calculated (the first parameter (i.e. X(1)) is the natural log of ‘hourly’ wind speed and the second parameter (i.e. X(2)) is ‘hourly’ temperature). Finally, the residual vector is converted to the original parameter vector using Equation (3) and the mean and standard deviation of each parameter for the condition (wet or dry) of this hour.
The input parameters for the simulation: X w , σW, X D , σD were estimated from historical data. Noting that the first parameter in X is natural log of wind speed (Ln (V)) and the second parameter is temperature (T). The mean and standard deviation of X (during wet hours (denoted as subscript W) or dry hours (denoted as subscript D)) are monthly values (i.e. were estimated for each month) to better ensure X (i.e. X (1) and X (2)) are Gaussian distributions as noted above. The matrices A and B were determined from Equations (5)–(8). For Des Moines, IA, they become: M0 ¼ M1 ¼
1 ρ0 ð2; 1Þ
ρ0 ð1; 2Þ
¼
1
1 0:0096 0:0096 1
ρ1 ð1; 1Þ ρ1 ð1; 2Þ 0:6779 0:0169 ¼ ρ1 ð2; 1Þ ρ1 ð2; 2Þ 0:0086 0:9955
(9) (10)
Hence A¼ B¼
0:6778 0:0104 0:0001 0:9955 0:7350 0:0012 0:0091 0:0942
(11)
(12)
Noting that, to obtain values for ρ0 ð1; 2Þ, ρ1 ð1; 1Þ, ρ1 ð1; 2Þ and ρ1 ð2; 1Þ, the temperature and wind speed data must be continuous (i.e. no missing data). Hence, the best period (with least missing data and as long as possible) for each station was chosen. Also, for any missing data of up to 4 h, linear interpolation was used to fill this gap. Since the biggest missing gap in these periods was selected to be less than 4 h, linear interpolation for the missing mean hourly data is acceptable. Also, as noted in (Matalas, 1967), there are more than one matrix B satisfying Equation (9) since if B is a root of Equation (8) and Γ is a matrix satisfying ΓΓT ¼ I, where I is the identity matrix, B* ¼ B*Γ is also a root of Equation (9). However, as pointed out in (Matalas, 1967), “Insofar as synthetic hydrology is concerned, the elements of B carry no physical significance” and any matrix B that satisfies Equation (8) could be used. The initial value for wind speed and temperature were selected as the values on January 1st of the first available year of data set. Following the procedure shown in Fig. 8 and 1000 years of data was simulated for each station. Table 3 compares statistical parameters between historical data (1973–2010) and simulated data for Des Moines, IA over the same 38year long block. The averages and standard deviations of these parameters are determined from 26 38-year long blocks (~1000 years). In general, there is very good agreement across nearly all parameters.
3.1. Setting conditions for freezing rain events The presence of freezing rain is ultimately determined by the vertical profile of temperature. Since there is only one measurement height in the data, here freezing rain is conditioned on the surface temperature alone.
Fig. 7. Best-fit distribution for ‘mean hourly’ wind speed at Des Moines, IA. 440
H. Nguyen Sinh et al.
Journal of Wind Engineering & Industrial Aerodynamics 184 (2019) 436–444
Fig. 8. Detail simulation procedure (*i ¼ 8784 for leap year).
ice-accretion model (Goodwin et al., 1983; Jones, 1998) to estimate ice thickness. It is crucial to check how the freezing rain data from simulated data compare with the ASCE 7 (2016) values. Fig. 9 shows the freezing rain data from 8000 years of simulation for the whole Midwestern superstation (8 stations). From this figure, the ice thickness for 50 year and 500 year MRI are 17.7 and 39.6 mm (~0.7 and 1.6 inches), respectively. These are quite consistent with the 50 year ice thickness in ASCE 7 (2016) (0.75 inches) for this region. Also as indicated in ASCE 7 (2016), the ratio between 500 year ice thickness and 50 year values is 2 which is quite close to this results (39.6/17.7–2.2). Another interesting result is that the ratio of 3-year ice thickness and 50-year ice thickness from this simulation procedure is 0.26 (4.6/17.7) which is very close with the recommended value of 0.3 from ISO 12494 (2017).
At any wet hour, the probability that freezing rain occurs will depend on the actual temperature being in ranges (preferred range -5 C < T < 0 C or adverse range 15 C < T < -5 C). If the temperature is outside of these ranges (T > 0 C or T < 15 C), then freezing rain is assumed not to occur. Probabilities were estimated from historical data and depended on location. After the freezing rain commences, the probability that freezing rain will continue in the next hour is denoted as P (FR/FR) and is indicated in Table 4. Comparing with Table 2, these probabilities (P (FR/ FR)) are quite similar to those for subsequent wet hours P (W/W). Table 4 also compares simulated data and historical data probabilities. In general, freezing rain rates in the simulated data are very close to freezing rain rates in the historical data. After setting a condition for freezing rain, the data (wind speed and precipitation) during freezing rain are extracted and placed in the simple 441
H. Nguyen Sinh et al.
Journal of Wind Engineering & Industrial Aerodynamics 184 (2019) 436–444
are 0.4 and 0.45 for glaze class G1 and G2, respectively. From Fig. 9, the ice thickness for 3-year and 50-year event are 4.6 mm and 17.7 mm, respectively. Hence the 3-year event will fall into G1 class and the 50year event will fall into G2 class. Typically the wind load varies the with square of wind speed, hence the k factor for wind load could be converted to the reduced wind speed factor of sqrt(0.4) and sqrt(0.45). Using the gust wind speed data from these stations and applying the superstation method (Peterka, 1992) provides 3-year and 50-year wind speed of 27.0 and 37.5 m/s, respectively. So the design values following ISO 12494 (2017) are approximately (23.7 m/s and 4.6 mm) or (18.1 m/s and 17.7 mm). Fig. 11 shows the 50 year MRI joint gust wind and ice thickness for the different locations, the earlier result from the simple approach (Nguyen Sinh et al., 2016), the design value pair from ASCE 7 (2016) and ISO 12494 (2017). In general all the joint hazard contours confirm the conservatism of the 50 year design values in ASCE 7 (2016) and ISO 12494 (2017). It is also shown that these stations have quite similar characterizations of joint wind and ice hazard. Interestingly, the simple model (Nguyen Sinh et al., 2016) represents a reasonable upper bound for much of the simulated data. Fig. 12 compares the results for a superstation (8 stations) between the earlier simplified approach (10,000 years simulation) (Nguyen Sinh et al., 2016) and the new sophisticated approach (8000 years of simulation – 1000 years for each station). Overall for MRI greater than 50 years the simplified approach is quite a lot more conservative than the new sophisticated approach especially for the ice thickness. The reason for this is in the simple approach an ill-fitting extreme-value distribution was used for ice thickness (GP distribution) while for the new approach, ice thickness was calculated based on jointly simulated meteorological parameters.
Table 3 Comparison of historical and simulated data for Des Moines, IA. Statistical parameters
Mean wind speed (m/s) Maximum wind speed (m/s) Mean precipitation (during wet hour) (mm/h) Maximum precipitation (mm/h) Mean temperature (oC) Maximum temperature (oC) Minimum temperature (oC)
Historical data
Simulated data Average of 26 38-year blocks
Standard deviation of 26 38-year blocks
4.5 25.2
4.5 20.4
0.01 1.81
1.3
1.3
0.01
152.4
162.2
30.27
10.1 41.7
10.3 42.4
0.27 1.93
31.7
32.6
2.25
Setting conditions for freezing rain events.
3.2. Joint hazard analysis The total ice thickness associated with each freezing rain event and the maximum mean wind speed during each freezing rain event were then extracted to run the joint hazard analysis. In some cases, high wind speeds can occur after the freezing rain event ends while the ice remains un-melted. To account for this, as discussed in (Nguyen Sinh et al., 2016), the wind speed associated with each freezing rain event will be the maximum mean wind speed from the beginning of the freezing rain event until the temperature rose above 1 C. Also, since gust wind speeds are typically required for strength design the maximum ‘hourly’ wind speeds were converted to 10 m height 3s-gust wind speeds utilizing a gust factor, G ¼ 1.33, assuming open country exposure (ASCE 7, 2016; Nguyen Sinh et al., 2016). Again noting that this ‘hourly’ wind speeds are actually the 2-min mean prior to every top-of-hour resulting in the gust factor of 1.33. Joint annual histograms for gust wind speed and ice thickness (i.e. the joint histogram of gust wind speed and ice thickness multiply with the annual frequency of freezing rain) were then constructed by counting the number of joint events with wind speed and ice thickness within specific ranges, dividing that number by the total number of data pairs and multiplying the outcome by the annual frequency of freezing rain (λ). The joint annual exceedance probability of gust wind and ice thickness can be determined from the histogram and is shown in Fig. 10 for Des Moines, IA. Once the joint annual exceedance probability distribution for wind and ice thickness is established, contours corresponding to specific hazard levels can be drawn (i.e. each contour representing different annual exceedance probabilities). 4. Results Noting that ISO 12494 suggests two combinations (for wind and ice) should be checked for design: a) a “reduced” 50-year wind action combined with the 3-year ice load and b) a “reduced” 3-year wind action combined with the 50-year ice load. The “reduced” factor k for wind load
Fig. 9. Freezing rain simulated data from superstation Midwestern US.
Table 4 Probability of FR occurrences. Station Name
Cedar Rapids Des Moines Lincoln Municipal Grand Island County Sioux City Muni Madison Sioux Fall Minneapolis/St. Paul
P(FR/-5
P(FR/-15 < T < 5)
P(FR/FR)
Historical data
Simulated data
Historical data
Simulated data
Historical data
Simulated data
0.015 0.019 0.011 0.016 0.012 0.0092 0.016 0.013
0.015 0.019 0.012 0.016 0.012 0.0094 0.016 0.013
0.0014 0.0021 0.0012 0.0021 0.0022 0.00044 0.0024 0.0018
0.0015 0.0021 0.0012 0.0021 0.0021 0.00044 0.0025 0.002
0.774 0.749 0.783 0.753 0.722 0.731 0.774 0.743
0.776 0.744 0.788 0.754 0.720 0.713 0.770 0.737
442
H. Nguyen Sinh et al.
Journal of Wind Engineering & Industrial Aerodynamics 184 (2019) 436–444
5. Conclusions In this paper, the joint wind-ice hazard for freezing rain events was established using a new approach. This study used data from 8 weather stations in US Midwest to simulate 1000 years of wind speed, temperature and precipitation data simultaneously for each station to establish the joint wind and ice thickness hazard. The results are compared with our earlier simplified approach (Nguyen Sinh et al., 2016), ASCE 7 (2016) values and ISO 12494 (2017) values. Firstly, both the simplified and new sophisticated approaches confirm that the design values in ASCE 7 (2016) and ISO 12494 (2017) are conservative for this set of stations despite ASCE 7 (2016) and ISO 12494 (2017) considering the joint effect of wind and ice hazard. The recommended 50 year coupled design values in ASCE 7 (2016) and ISO 12494 (2017) for gust wind speed and ice thickness actually belong to much higher MRI. In other words, the 50 year ice thickness values from ASCE 7 (2016) must be combined with much lower wind speed than the concurrent 3-s gust wind speed in the “Equivalent radial ice thicknesses due to freezing rain with concurrent 3-s gust speeds, for a 50-year Mean Recurrence Interval” map from ASCE 7 (2016) and vice versa, or both these values should be reduced following the 50 year joint hazard contour line in Fig. 11 or 12. The reduction factor k in ISO 12494 could be reduced to lower values (~0.35 for these stations). Secondly, one of the significant problems for the simple approach was missing observations that commonly occur during freezing rains, leading to an incomplete dataset to drive simulations. To overcome this problem, the current sophisticated approach described here was developed to examine freezing rain events and better match historical data as well as ASCE 7 (2016) values. Moreover this methodology could be used to for other types of precipitations (e.g. snow, rainfall) by setting different meteorological conditions for them, and to examine the impact of climate change. Also, it was shown that the simple approach represents a reasonable upper bound for many of the results obtained from the sophisticated model. Finally, this procedure used the simple ice-accretion model (Goodwin et al., 1983; Jones, 1998) in which the important role of icicles is entirely neglected. Hence the ice thickness values indicated in this paper may not predict realistic values due to neglecting icicle growth. However the simple ice-accretion model for horizontal circular cylinders is based on three assumptions (Jones, 1998): “a) the collision efficiency of the raindrops with the cylinder is 100%; b) all rain water impinging on the cylinder sticks to the cylinder and freezes; and c) the ice accretes uniformly around the circumference of the cylinder”. These assumptions are conservative, hence the simple ice-accretion model is conservative in this sense. Another important note when designing the structures under icing load is that for high elevations region as well as for tall towers in general, not only freezing precipitation is needed to be considered but also the in-cloud icings as indicated in ISO 12494 (2017).
Fig. 10. Joint annual exceedance rate distribution at Des Moines, IA.
Fig. 11. 50 year MRI Joint Hazard Curves for different Midwestern US stations.
Appendix A. Supplementary data Supplementary data to this article can be found online at https://doi. org/10.1016/j.jweia.2018.12.012. References ASCE, 2016. Minimum Design Loads for Buildings and Other Structures. American Society of Civil Engineers, Reston, VA. Bailey, N.T.J., 1964. The Elements of Stochastic Processes. John Wiley, New York. Benjamin, J.R., Cornell, C.A., 1970. Probability, Statistics and Decision for Civil Engineers. McGraw-Hill, New York. Buishand, T.A., 1978. Some remarks on the use of daily rainfall models. J. Hydrol. 36, 295–308. Caskey Jr., J.E., 1963. A Markov chain model for the probability of precipitation occurrence in intervals various lengths. Monthly Weather Review-American Meteorological Society 91, 298–301. Changnon, S.A., Creech, T.G., 2003. NOTES and correspondence, sources of data on freezing rain and resulting damages. J. Appl. Meteorol. 42, 1514–1518.
Fig. 12. Comparing sophisticated and simplified approaches for the Midwestern US superstation. 443
H. Nguyen Sinh et al.
Journal of Wind Engineering & Industrial Aerodynamics 184 (2019) 436–444 Lombardo, F.T., Nguyen Sinh, H., Letchford, C.W., Rosowsky, D.V., 2014. “Assessing the Joint Wind and Temperature Hazard for the United States”. Structure Congress 2014, Boston, MA. Matalas, N.C., 1967. Mathematical assessment of synthetic hydrology. J. Water Resour. Res. 3 (4), 937–945. Nguyen Sinh, H., Lombardo, F.T., Letchford, C.W., Rosowsky, D.V., 2016. Characterization of the joint wind and ice hazard in the midwestern United States. Nat. Hazards Rev. 17 (3) ( J. ASCE). Peterka, J.A., 1992. Improved extreme wind prediction for the United States. J. Wind Eng. Ind. Aerod. 41–44, 533–541. Richardson, C.W., 1981. Stochastic simulation of daily precipitation, temperature, and solar radiation. J. Water Resour. Res. 17 (1), 182–190. Roberts, E., Stewart, R.E., 2008. On the occurrence of freezing rain and ice pellets over the Eastern Canadian Arctic. Atmos. Res. 89 (1–2), 93–109. Smith, R.E., Schreiber, H.A., 1973. Point processes of seasonal thunderstorm rainfall, 1, distribution of rainfall event. J. Water Resour. Res. 9 (4), 871–884. Weiss, L.L., 1964. Sequence of wet and dry days described by a Markov chain probability model. Monthly Weather Review- American Meteorological Society 92, 169–176.
Chin, E.H., 1957. Modeling daily precipitation occurrence process with Markov chain. J. Water Resour. Res. 13 (6), 949–956. Gabriel, R., Neuman, J., 1962. A Markov chain model for daily rainfall occurrence at Tel Aviv, Israel. Q. J. R. Meteorol. Soc. 88, 90–95. Goodwin III, E.J., Mozer, J.D., Di Gioia Jr., A.M., Power, B.A., 1983. Predicting ice and snow loads for transmission lines. In: The Proceedings of the First IWAIS, pp. 267–273. Haan, C.T., Allen, D.M., Street, J.O., 1976. A Markov chain model of daily rainfall. J. Water Resour. Res. 12 (3), 443–449. Hopkins, J.W., Robillard, P., 1964. Some statistics of daily rainfall occurrence for the Canadian Prairi provinces. Journal of Applied Meteorology and Climatology 3, 600–602. Integrated surface global hourly data (ISH), 2001. National Oceanic and Atmospheric Administration (NOAA). Retrieved from. ftp://ftp.ncdc.noaa.gov/pub/data. (Accessed 1 January 2001). ISO 12494, 2017. International Standard. Atmospheric Icing of Structures. Jones, K.F., 1998. A simple model for freezing rain ice loads. Atmos. Res. 46, 87–97.
444
Contents lists available at ScienceDirect
Journal of Wind Engineering & Industrial Aerodynamics journal homepage: www.elsevier.com/locate/jweia
Multivariate simulation for assessing the joint wind and ice hazard in the United States Hung Nguyen Sinh a, *, Franklin T. Lombardo b, Chris Letchford c a b c
BMT Fluid Mechanics Ltd, New York City, NY, USA Civil & Environmental Engineering Department, University of Illinois at Urbana-Champaign, Urbana, IL, USA Civil & Environmental Engineering Department, Rensselaer Polytechnic Institute, Troy, NY, USA
A R T I C L E I N F O
A B S T R A C T
Keywords: Joint hazard Multivariate simulation Performance-based design Wind/ice hazard Freezing rain Ice load
Multi-hazard assessment is receiving increasing attention as often commonly occurring hazards may contain a second or even a third associated hazard. Generally, multi-hazards could be joint hazards (i.e. these hazards occur at the same time), for example, a hurricane usually comes with strong wind, storm surge, and heavy rainfall. They can also be one hazard preceding another, for example, an earthquake generating a tsunami. The jointly occurring hazards in some cases can be more devastating than a single hazard. This paper will describe the assessment of two jointly occurring hazards, wind and ice formation by using an advance numerical simulation. First, multivariate (hourly wind, temperature and precipitation) simulation for many years in the Midwestern US is undertaken. By setting a condition for freezing rain, these parameters were input into a simple ice-accretion model to estimate ice-thickness. A joint hazard curve for wind and ice was then constructed and compared to the existing prescribed treatment of these two hazards in design loading guidelines, such as ASCE 7 in the United States, ISO 12494 (International Standard); and the results from our previous model.
1. Introduction For locations in the US, Canada, and other higher latitude countries, wind and ice (e.g., freezing rain) are hazards whose properties are of interest for design of transmission lines and other energy infrastructure (e.g., wind turbines). Freezing rain is one of the major weather hazards in the US causing annual losses of $314 million (2000 US$) during 1949–2000 (Changnon et al., 2003). As pointed out in our previous paper (Nguyen Sinh et al., 2016), although there have been numerous scientific and engineering studies of freezing rain, few considered the combined effect of the wind hazard and the ice hazard occurring simultaneously. Nguyen Sinh et al. (2016) used a simple approach to deal with the joint wind and ice hazard during freezing rain events in which the wind speed and ice data based on actual observations over a long period (~300 years using the superstation approach (Peterka, 1992)) were simulated separately to assess the joint wind and ice hazard. The details for this approach is shown in Fig. 1. Firstly, the meteorological data from measurement stations in the Midwestern US were collected from the ISH (2001) database. Quality control measures to remove any spurious (i.e. not physically plausible) data points were undertaken (e.g. freezing rain that occurred at high temperatures (>0 Celsius degree), unrealistic wind
speeds and precipitations (wind speed values > 30 m/s, precipitations > 15 mm/h for winter and >40 mm/h for summer, …). Then the simple ice-accretion model (Goodwin et al., 1983; Jones, 1998) was used to estimate the ice thickness for each freezing rain event. The superstation approach (Peterka, 1992) was then used to group these stations into a “superstation”. The wind speeds (occurring during freezing rain) were fitted by a Weibull distribution while the Generalized Pareto (GP) distribution was used to fit the ice thickness using the peaks-over-threshold method. The wind speeds and ice thickness were then simulated separately and their joint histogram was established to create the joint wind and ice hazard. Details and results of this approach are described in Nguyen Sinh et al. (2016). Even though this simplified approach was conservative, it still confirmed that the ASCE 7 (2016) design values for 50 year Mean Recurrence Interval (MRI) for joint wind-ice hazard belong to much higher MRI (~150 years). In other words, the 50 year ice thickness values from ASCE 7 (2016) must be combined with much lower wind speed than the concurrent 3-s gust wind speed in the “Equivalent radial ice thicknesses due to freezing rain with concurrent 3-s gust speeds, for a 50-year Mean Recurrence Interval” map from ASCE 7 (2016) and vice versa or both these values must be reduced to match the 50 year joint hazard
* Corresponding author. E-mail addresses: [email protected] (H. Nguyen Sinh), [email protected] (F.T. Lombardo), [email protected] (C. Letchford). https://doi.org/10.1016/j.jweia.2018.12.012 Received 8 October 2018; Received in revised form 12 December 2018; Accepted 13 December 2018 Available online 20 December 2018 0167-6105/© 2018 Elsevier Ltd. All rights reserved.
H. Nguyen Sinh et al.
Journal of Wind Engineering & Industrial Aerodynamics 184 (2019) 436–444
associated with each freezing rain event. Finally, the joint wind and ice hazard results are established and compared with the results from the simple model (Nguyen Sinh et al., 2016) as well as ASCE 7 (2016)) and ISO 12494 (2017) design values. Since the freezing rain data was only used to set conditions for freezing rain, actual missing data during freezing rain is no longer problematic. Although the simulation procedure in this approach is very different to the one in the simple model (Nguyen Sinh et al., 2016), the joint hazard curves were established using the joint probability distribution of ice and wind which is similar to the simple model (Nguyen Sinh et al., 2016). The difference between this joint hazard approach and ASCE 7 was discussed in detail in Nguyen Sinh et al. (2016). 2. Data and methodology The meteorological data used for this study is collected from the National Climate Data Center (NCDC) within NOAA's Integrated Surface Global Hourly (ISH) database, which can be found at ftp://ftp.ncdc.noaa. gov/pub/data (ISH, 2001). Table 1 shows the station numbers and names, locations and the associated time periods used for this study. Fig. 3 shows that these 8 stations have very similar ice thickness for the 50 year MRI according to ASCE 7 (2016). They are the same data set used in our previous paper (Nguyen Sinh et al., 2016). Table 1 also shows the annual frequency (λ) of freezing rain for each station determined from the total number of freezing rain events divided by the total number of years of data for that station. It is seen that all have approximately 2 events per year. Noting that, as indicated in Nguyen Sinh et al. (2016), the ‘hourly’ wind speed data and the ‘hourly’ temperature in this paper are actually the 2-min mean prior to every top-of-hour. The detailed procedure employed for this simulation is shown in Fig. 4. Firstly, all the data (wind speeds, precipitations and the weather observations) from the 8 Midwestern US stations (Table 1) were collected from the ISH (2001) database. These historical data are used to estimate parameters for the simulation discussed in more detail below. The simulation process is time dependent within each variable and the three variables (wind speed, temperature and precipitation) are inter-dependent. For instance, one might expect that the ‘hourly’ temperature during a wet hour in summer is more likely to be lower than during a dry hour (Richardson, 1981). Also the correlation between wind speed and precipitation, or wind speed and temperature, may not be obvious but still present. The main idea is to consider precipitation as the primary variable and then conditioning the other two variables (wind speed and temperature) for a given hour on whether that hour is wet or dry. A Markov chain is used to simulate the precipitation (P) (Richardson, 1981). Then multivariate simulation (wind, temperature) for that given hourly precipitation is achieved by using a multi-variable Markov generating process (Matalas, 1967). By establishing the condition for freezing rain, these parameters were input into a simple ice-accretion model (Goodwin et al., 1983; Jones, 1998) to estimate the ice
Fig. 1. The procedure for the simple model (Nguyen Sinh et al., 2016). (GP: Generalized Pareto distribution).
contour line developed in Nguyen Sinh et al. (2016) and shown in Fig. 2. However, one of the significant problems for this approach was missing observations that commonly occur during freezing rains, leading to an incomplete dataset to drive simulations. In this paper, a more sophisticated model is introduced in which all the meteorological variables needed for an ice-accretion model (Goodwin et al., 1983; Jones, 1998); total hourly precipitation, hourly temperature and mean hourly wind speed, are simulated simultaneously such that the variables are time dependent and their inter-dependence established. Conditions for freezing rain are then set and freezing rain events estimated from the simulated meteorological data. These freezing rain event meteorological variables are input into a simple ice-accretion model (Goodwin et al., 1983; Jones, 1998) to estimate the ice thickness
Table 1 Information on stations used for Midwestern, US. Station Number
Station Name
725450 725460 725510
Cedar Rapids Des Moines Lincoln Municipal Grand Island County Sioux City Muni Madison Sioux Fall Minneapolis/ St. Paul
725520 725570 726410 726510 726580
Fig. 2. Joint gust wind and ice hazard for Midwestern US superstation considering wind effects after freezing rain events (Nguyen Sinh et al., 2016): a) Joint hazard curves; b) Marginal distribution for ice thickness during freezing rain events; c) Marginal distribution for wind speed during freezing rain events. 437
Length of recorded data (years)
FR frequency (λ)
2010 2010 2010
33 38 38
2.03 2.16 1.63
1978
2010
33
2.18
IA
1973
2010
38
1.53
WI SD MN
1973 1973 1973
2010 2010 2010
38 38 38
1.68 1.97 1.92
State
Time period Start
End
IA IA NE
1978 1973 1973
NE
H. Nguyen Sinh et al.
Journal of Wind Engineering & Industrial Aerodynamics 184 (2019) 436–444
Fig. 3. ASCE 50 year MRI ice map (adapted from ASCE 7 (2016)).
Markov Chain (Richardson, 1981) is employed. Basically, the occurrence of wet or dry hours was modeled by a first-order Markov Chain (Bailey, 1964) and then on the rainy hour, the best-fit distribution for precipitation intensity is used to estimate the amount of precipitation in that hour. In the past, Gabriel and Neuman (1962), Caskey (1963), Weiss (1964), Hopkins and Robillard (1964) and Haan (1976) used Markov chain models for simulating daily precipitation. “Smith and Schreiber (1973) found that a first-order Markov chain was superior to a Bernoulli model for describing the occurrence of wet or dry days in southeastern Arizona, US” (Richardson, 1981). This study will use a similar first-order Markov chain approach. Fig. 5 shows the detail of applying the Markov chain to generating hourly total precipitation. Let W (wet) and D (dry) denote for the state of an hour, P (W/D) is the probability of a W hour given a condition D at the previous hour. Similarly for the definition of P (D/D), P(D/W) and P(W/ W) (Richardson, 1981). Then: P(D/W)¼1-P(W/W)
(1)
P(D/D)¼1-P(W/D)
(2)
Therefore the time series of precipitation is fully defined given P (W/ W) and P (W/D) and the wet or dry state at the first hour. In other words, if these conditional probabilities can be estimated from the historical data, by assuming the initial state for the first hour, condition W or D of all hours in a time series can be generated. Based on the historical events, these probabilities are indicated in Table 2. Another probability distribution is used to characterize the amount of precipitation at each wet hour. Inverse cumulative density function (CDF) (Benjamin and Cornell, 1970) is used to determine the best-fit distribution for the precipitation intensity. This approach assumes that the CDF values (u) of a variable (X) are initially uniformly distributed over [0 1], then different types of probability distributions are selected. For each type of probability distribution, the variable X'i is calculated by taking the inverse CDF of u (i.e. X'i ¼ CDF1 i (u)). All of the points (X, X'i)
Fig. 4. The procedure for the sophisticated model.
thickness. The joint histogram and joint exceedance probability of wind and ice thickness are then constructed to establish the joint hazard curves for wind and ice. 2.1. Precipitation simulation Many models have been proposed to simulate precipitation: Buishand (1978), Chin (1957) and Gabriel and Neuman (1962). In this research, a 438
H. Nguyen Sinh et al.
Journal of Wind Engineering & Industrial Aerodynamics 184 (2019) 436–444
Fig. 6. Best-fit distribution for hourly total precipitation at Des Moines, IA.
temperature based on the precipitation that has been simulated above. First the residual variables were calculated by removing the mean and standard deviation of each variable as expressed in Equation (3). xi ðpÞ ¼ Fig. 5. The procedure for the Markov chain.
Cedar Rapids Des Moines Lincoln Municipal Grand Island County Sioux City Muni Madison Sioux Fall Minneapolis/St. Paul
State
P
Best fit distribution for P
P(W/ W)
P(W/ D)
P(W)
IA IA NE NE
0.838 0.78 0.794 0.809
0.025 0.025 0.021 0.019
0.135 0.102 0.091 0.089
Weibull Log normal Log normal Log normal
IA WI SD MN
0.789 0.767 0.795 0.768
0.021 0.029 0.02 0.024
0.092 0.11 0.09 0.095
Weibull Weibull Weibull Log normal
(3)
where XðpÞ and σ ðpÞ are the mean and standard deviation and xi(p) is the residual component for variable p, at time ‘i’. Then the weakly-stationary generating process suggested by Matalas (1967) was used to generate a residual series of all variables.
Table 2 Estimated probabilities for Precipitation. Station Name
Xi ðpÞ XðpÞ σ ðpÞ
xi ðpÞ ¼ Axi1 ðpÞ þ Bεi ðpÞ
(4)
where xi ðpÞ, xi1 ðpÞ are residual variables at time ‘i’ and ‘i-1’ of variable p; εi ðpÞ are independent random components that are normally distributed with a zero mean and a unit variance associated with variable p. The A and B matrices are defined such that the new sequences have the desired serial correlation and cross-correlation coefficients. Equation (4) implies that the residuals of variables are normally distributed and that the serial correlation of each variable may be described by a first-order linear autoregressive model (Matalas, 1967). The A and B matrices can be determined from the following equations.
for each distribution are then plotted and the distribution closest to the line of slope 1 is the best-fit distribution. Fig. 6 shows that Weibull is the best-fit distribution for the precipitation in Des Moines, IA within six types of selected distributions (Normal, Lognormal, Extreme Type 1 (ET1), Extreme Type 2 (ET2), Weibull (WB) and Exponential (EXP)). Among methods for finding best-fit probability distributions, the inverse CDF method (Benjamin and Cornell, 1970) was selected because it provides a simple, graphical technique that can be generalized for any distribution form without the need for specific types of plotting papers and furthermore can be easily programmed.
A ¼ M1 M 1 0
(5)
T BBT ¼ M0 M1 M 1 0 M1
(6)
where the superscripts 1 and T denote the inverse and transpose of a matrix respectively. M0 and M1 are the lag-0 and lag-1 covariance matrices. For the case of 2 variables (p ¼ 2), M0 and M1 are expressed by following equations M0 ¼
2.2. Lag-one Markov generating process
M1 ¼
A lag-one Markov process - a multivariate weakly-stationary generating process (Matalas, 1967) - was used to simulate the wind speed and
1 ρ0 ð2; 1Þ
ρ0 ð1; 2Þ
(7)
1
ρ1 ð1; 1Þ ρ1 ð1; 2Þ ρ1 ð2; 1Þ ρ1 ð2; 2Þ
(8)
where ρ0 ði; jÞ and ρ1 ði; jÞ are the lag-0 and lag-1 cross-correlation coeffi439
H. Nguyen Sinh et al.
Journal of Wind Engineering & Industrial Aerodynamics 184 (2019) 436–444
cient between variables ‘i’ and ‘j’, respectively. In this application, ‘i’ will represent wind speed and ‘j’ will represent temperature.
3. Analysis Using Des Moines, IA data this section illustrates the simulation procedure in more detail. The probabilities P (W/W) and P (W/D) were estimated from historical data and are shown in Table 2. The initial state (wet or dry) was determined based on a Bernoulli trial with P (W) ¼ 0.102 and P (D) ¼ 1- P (W) ¼ 0.898. All other states (at each following hour) were determined by a Markov chain with P (W/ W) ¼ 0.78 and P (W/D) ¼ 0.025. At all wet hours, the amount of precipitation was determined using the best-fit distribution for the precipitation in Des Moines, IA. Although, the best-fit distribution varies for each location (Table 2) and both Log normal and Weibull distributions are best fit for data, the Log normal distribution usually predicted slightly higher values than historical data (i.e. observed data). For consistency and conservatism, a Log normal distribution was used for all locations.
2.3. Multivariate hourly wind speed and hourly temperature conditioned on hourly precipitation series In this section, application of the multi-variable lag-one Markov generating process for simulating ‘hourly’ wind speed and ‘hourly’ temperature will be discussed. Noting that, Matalas' model (Matalas, 1967) (i.e. Equation (4)) works best when variables (and hence residual variables) are normally distributed. For the ‘hourly’ temperature, Lombardo et al. (2014) showed that temperature did indeed follow a normal distribution for each season. For the ‘hourly’ wind speed, the Weibull distribution is most often the best-fit distribution. However, using the Inverse CDF method (Benjamin and Cornell, 1970), Fig. 7 shows that a log-normal distribution is as good a fit in the region of 0–15 m/s which is the typical range of wind speeds during freezing rains (Roberts and Stewart, 2008) for this and other locations. Hence, it is reasonable to apply Matalas' model (Matalas, 1967) with the first parameter being the natural log of wind speed (i.e. X (1) ¼ ln (V), where V is the ‘hourly’ wind speed) and the second parameter being temperature (i.e. (X (2) ¼ T). There are hours with calm (i.e. no wind) that lead to undefined natural log, hence it is necessary to add a value of 0.01 (m/s) to these values. Since the accuracy resolution for anemometers is 1 knot (~0.514 m/s), this adding has no effect on the statistical parameter estimating. In addition, monthly analysis (i.e. these statistical parameters are estimated for each month) is undertaken to ensure these variables (wind speed and temperature) follow the normal distribution. Fig. 8 outlines the detailed simulation procedure. In summary, first, a time series for precipitation is constructed using a Markov chain. The initial values (i.e. at first hour) were selected as the values at the beginning of the first available year from historical data. Then at each following hour, by using the Matalas' model (Matalas, 1967) (or using Equation (4)) the residual vector of two parameters are calculated (the first parameter (i.e. X(1)) is the natural log of ‘hourly’ wind speed and the second parameter (i.e. X(2)) is ‘hourly’ temperature). Finally, the residual vector is converted to the original parameter vector using Equation (3) and the mean and standard deviation of each parameter for the condition (wet or dry) of this hour.
The input parameters for the simulation: X w , σW, X D , σD were estimated from historical data. Noting that the first parameter in X is natural log of wind speed (Ln (V)) and the second parameter is temperature (T). The mean and standard deviation of X (during wet hours (denoted as subscript W) or dry hours (denoted as subscript D)) are monthly values (i.e. were estimated for each month) to better ensure X (i.e. X (1) and X (2)) are Gaussian distributions as noted above. The matrices A and B were determined from Equations (5)–(8). For Des Moines, IA, they become: M0 ¼ M1 ¼
1 ρ0 ð2; 1Þ
ρ0 ð1; 2Þ
¼
1
1 0:0096 0:0096 1
ρ1 ð1; 1Þ ρ1 ð1; 2Þ 0:6779 0:0169 ¼ ρ1 ð2; 1Þ ρ1 ð2; 2Þ 0:0086 0:9955
(9) (10)
Hence A¼ B¼
0:6778 0:0104 0:0001 0:9955 0:7350 0:0012 0:0091 0:0942
(11)
(12)
Noting that, to obtain values for ρ0 ð1; 2Þ, ρ1 ð1; 1Þ, ρ1 ð1; 2Þ and ρ1 ð2; 1Þ, the temperature and wind speed data must be continuous (i.e. no missing data). Hence, the best period (with least missing data and as long as possible) for each station was chosen. Also, for any missing data of up to 4 h, linear interpolation was used to fill this gap. Since the biggest missing gap in these periods was selected to be less than 4 h, linear interpolation for the missing mean hourly data is acceptable. Also, as noted in (Matalas, 1967), there are more than one matrix B satisfying Equation (9) since if B is a root of Equation (8) and Γ is a matrix satisfying ΓΓT ¼ I, where I is the identity matrix, B* ¼ B*Γ is also a root of Equation (9). However, as pointed out in (Matalas, 1967), “Insofar as synthetic hydrology is concerned, the elements of B carry no physical significance” and any matrix B that satisfies Equation (8) could be used. The initial value for wind speed and temperature were selected as the values on January 1st of the first available year of data set. Following the procedure shown in Fig. 8 and 1000 years of data was simulated for each station. Table 3 compares statistical parameters between historical data (1973–2010) and simulated data for Des Moines, IA over the same 38year long block. The averages and standard deviations of these parameters are determined from 26 38-year long blocks (~1000 years). In general, there is very good agreement across nearly all parameters.
3.1. Setting conditions for freezing rain events The presence of freezing rain is ultimately determined by the vertical profile of temperature. Since there is only one measurement height in the data, here freezing rain is conditioned on the surface temperature alone.
Fig. 7. Best-fit distribution for ‘mean hourly’ wind speed at Des Moines, IA. 440
H. Nguyen Sinh et al.
Journal of Wind Engineering & Industrial Aerodynamics 184 (2019) 436–444
Fig. 8. Detail simulation procedure (*i ¼ 8784 for leap year).
ice-accretion model (Goodwin et al., 1983; Jones, 1998) to estimate ice thickness. It is crucial to check how the freezing rain data from simulated data compare with the ASCE 7 (2016) values. Fig. 9 shows the freezing rain data from 8000 years of simulation for the whole Midwestern superstation (8 stations). From this figure, the ice thickness for 50 year and 500 year MRI are 17.7 and 39.6 mm (~0.7 and 1.6 inches), respectively. These are quite consistent with the 50 year ice thickness in ASCE 7 (2016) (0.75 inches) for this region. Also as indicated in ASCE 7 (2016), the ratio between 500 year ice thickness and 50 year values is 2 which is quite close to this results (39.6/17.7–2.2). Another interesting result is that the ratio of 3-year ice thickness and 50-year ice thickness from this simulation procedure is 0.26 (4.6/17.7) which is very close with the recommended value of 0.3 from ISO 12494 (2017).
At any wet hour, the probability that freezing rain occurs will depend on the actual temperature being in ranges (preferred range -5 C < T < 0 C or adverse range 15 C < T < -5 C). If the temperature is outside of these ranges (T > 0 C or T < 15 C), then freezing rain is assumed not to occur. Probabilities were estimated from historical data and depended on location. After the freezing rain commences, the probability that freezing rain will continue in the next hour is denoted as P (FR/FR) and is indicated in Table 4. Comparing with Table 2, these probabilities (P (FR/ FR)) are quite similar to those for subsequent wet hours P (W/W). Table 4 also compares simulated data and historical data probabilities. In general, freezing rain rates in the simulated data are very close to freezing rain rates in the historical data. After setting a condition for freezing rain, the data (wind speed and precipitation) during freezing rain are extracted and placed in the simple 441
H. Nguyen Sinh et al.
Journal of Wind Engineering & Industrial Aerodynamics 184 (2019) 436–444
are 0.4 and 0.45 for glaze class G1 and G2, respectively. From Fig. 9, the ice thickness for 3-year and 50-year event are 4.6 mm and 17.7 mm, respectively. Hence the 3-year event will fall into G1 class and the 50year event will fall into G2 class. Typically the wind load varies the with square of wind speed, hence the k factor for wind load could be converted to the reduced wind speed factor of sqrt(0.4) and sqrt(0.45). Using the gust wind speed data from these stations and applying the superstation method (Peterka, 1992) provides 3-year and 50-year wind speed of 27.0 and 37.5 m/s, respectively. So the design values following ISO 12494 (2017) are approximately (23.7 m/s and 4.6 mm) or (18.1 m/s and 17.7 mm). Fig. 11 shows the 50 year MRI joint gust wind and ice thickness for the different locations, the earlier result from the simple approach (Nguyen Sinh et al., 2016), the design value pair from ASCE 7 (2016) and ISO 12494 (2017). In general all the joint hazard contours confirm the conservatism of the 50 year design values in ASCE 7 (2016) and ISO 12494 (2017). It is also shown that these stations have quite similar characterizations of joint wind and ice hazard. Interestingly, the simple model (Nguyen Sinh et al., 2016) represents a reasonable upper bound for much of the simulated data. Fig. 12 compares the results for a superstation (8 stations) between the earlier simplified approach (10,000 years simulation) (Nguyen Sinh et al., 2016) and the new sophisticated approach (8000 years of simulation – 1000 years for each station). Overall for MRI greater than 50 years the simplified approach is quite a lot more conservative than the new sophisticated approach especially for the ice thickness. The reason for this is in the simple approach an ill-fitting extreme-value distribution was used for ice thickness (GP distribution) while for the new approach, ice thickness was calculated based on jointly simulated meteorological parameters.
Table 3 Comparison of historical and simulated data for Des Moines, IA. Statistical parameters
Mean wind speed (m/s) Maximum wind speed (m/s) Mean precipitation (during wet hour) (mm/h) Maximum precipitation (mm/h) Mean temperature (oC) Maximum temperature (oC) Minimum temperature (oC)
Historical data
Simulated data Average of 26 38-year blocks
Standard deviation of 26 38-year blocks
4.5 25.2
4.5 20.4
0.01 1.81
1.3
1.3
0.01
152.4
162.2
30.27
10.1 41.7
10.3 42.4
0.27 1.93
31.7
32.6
2.25
Setting conditions for freezing rain events.
3.2. Joint hazard analysis The total ice thickness associated with each freezing rain event and the maximum mean wind speed during each freezing rain event were then extracted to run the joint hazard analysis. In some cases, high wind speeds can occur after the freezing rain event ends while the ice remains un-melted. To account for this, as discussed in (Nguyen Sinh et al., 2016), the wind speed associated with each freezing rain event will be the maximum mean wind speed from the beginning of the freezing rain event until the temperature rose above 1 C. Also, since gust wind speeds are typically required for strength design the maximum ‘hourly’ wind speeds were converted to 10 m height 3s-gust wind speeds utilizing a gust factor, G ¼ 1.33, assuming open country exposure (ASCE 7, 2016; Nguyen Sinh et al., 2016). Again noting that this ‘hourly’ wind speeds are actually the 2-min mean prior to every top-of-hour resulting in the gust factor of 1.33. Joint annual histograms for gust wind speed and ice thickness (i.e. the joint histogram of gust wind speed and ice thickness multiply with the annual frequency of freezing rain) were then constructed by counting the number of joint events with wind speed and ice thickness within specific ranges, dividing that number by the total number of data pairs and multiplying the outcome by the annual frequency of freezing rain (λ). The joint annual exceedance probability of gust wind and ice thickness can be determined from the histogram and is shown in Fig. 10 for Des Moines, IA. Once the joint annual exceedance probability distribution for wind and ice thickness is established, contours corresponding to specific hazard levels can be drawn (i.e. each contour representing different annual exceedance probabilities). 4. Results Noting that ISO 12494 suggests two combinations (for wind and ice) should be checked for design: a) a “reduced” 50-year wind action combined with the 3-year ice load and b) a “reduced” 3-year wind action combined with the 50-year ice load. The “reduced” factor k for wind load
Fig. 9. Freezing rain simulated data from superstation Midwestern US.
Table 4 Probability of FR occurrences. Station Name
Cedar Rapids Des Moines Lincoln Municipal Grand Island County Sioux City Muni Madison Sioux Fall Minneapolis/St. Paul
P(FR/-5
P(FR/-15 < T < 5)
P(FR/FR)
Historical data
Simulated data
Historical data
Simulated data
Historical data
Simulated data
0.015 0.019 0.011 0.016 0.012 0.0092 0.016 0.013
0.015 0.019 0.012 0.016 0.012 0.0094 0.016 0.013
0.0014 0.0021 0.0012 0.0021 0.0022 0.00044 0.0024 0.0018
0.0015 0.0021 0.0012 0.0021 0.0021 0.00044 0.0025 0.002
0.774 0.749 0.783 0.753 0.722 0.731 0.774 0.743
0.776 0.744 0.788 0.754 0.720 0.713 0.770 0.737
442
H. Nguyen Sinh et al.
Journal of Wind Engineering & Industrial Aerodynamics 184 (2019) 436–444
5. Conclusions In this paper, the joint wind-ice hazard for freezing rain events was established using a new approach. This study used data from 8 weather stations in US Midwest to simulate 1000 years of wind speed, temperature and precipitation data simultaneously for each station to establish the joint wind and ice thickness hazard. The results are compared with our earlier simplified approach (Nguyen Sinh et al., 2016), ASCE 7 (2016) values and ISO 12494 (2017) values. Firstly, both the simplified and new sophisticated approaches confirm that the design values in ASCE 7 (2016) and ISO 12494 (2017) are conservative for this set of stations despite ASCE 7 (2016) and ISO 12494 (2017) considering the joint effect of wind and ice hazard. The recommended 50 year coupled design values in ASCE 7 (2016) and ISO 12494 (2017) for gust wind speed and ice thickness actually belong to much higher MRI. In other words, the 50 year ice thickness values from ASCE 7 (2016) must be combined with much lower wind speed than the concurrent 3-s gust wind speed in the “Equivalent radial ice thicknesses due to freezing rain with concurrent 3-s gust speeds, for a 50-year Mean Recurrence Interval” map from ASCE 7 (2016) and vice versa, or both these values should be reduced following the 50 year joint hazard contour line in Fig. 11 or 12. The reduction factor k in ISO 12494 could be reduced to lower values (~0.35 for these stations). Secondly, one of the significant problems for the simple approach was missing observations that commonly occur during freezing rains, leading to an incomplete dataset to drive simulations. To overcome this problem, the current sophisticated approach described here was developed to examine freezing rain events and better match historical data as well as ASCE 7 (2016) values. Moreover this methodology could be used to for other types of precipitations (e.g. snow, rainfall) by setting different meteorological conditions for them, and to examine the impact of climate change. Also, it was shown that the simple approach represents a reasonable upper bound for many of the results obtained from the sophisticated model. Finally, this procedure used the simple ice-accretion model (Goodwin et al., 1983; Jones, 1998) in which the important role of icicles is entirely neglected. Hence the ice thickness values indicated in this paper may not predict realistic values due to neglecting icicle growth. However the simple ice-accretion model for horizontal circular cylinders is based on three assumptions (Jones, 1998): “a) the collision efficiency of the raindrops with the cylinder is 100%; b) all rain water impinging on the cylinder sticks to the cylinder and freezes; and c) the ice accretes uniformly around the circumference of the cylinder”. These assumptions are conservative, hence the simple ice-accretion model is conservative in this sense. Another important note when designing the structures under icing load is that for high elevations region as well as for tall towers in general, not only freezing precipitation is needed to be considered but also the in-cloud icings as indicated in ISO 12494 (2017).
Fig. 10. Joint annual exceedance rate distribution at Des Moines, IA.
Fig. 11. 50 year MRI Joint Hazard Curves for different Midwestern US stations.
Appendix A. Supplementary data Supplementary data to this article can be found online at https://doi. org/10.1016/j.jweia.2018.12.012. References ASCE, 2016. Minimum Design Loads for Buildings and Other Structures. American Society of Civil Engineers, Reston, VA. Bailey, N.T.J., 1964. The Elements of Stochastic Processes. John Wiley, New York. Benjamin, J.R., Cornell, C.A., 1970. Probability, Statistics and Decision for Civil Engineers. McGraw-Hill, New York. Buishand, T.A., 1978. Some remarks on the use of daily rainfall models. J. Hydrol. 36, 295–308. Caskey Jr., J.E., 1963. A Markov chain model for the probability of precipitation occurrence in intervals various lengths. Monthly Weather Review-American Meteorological Society 91, 298–301. Changnon, S.A., Creech, T.G., 2003. NOTES and correspondence, sources of data on freezing rain and resulting damages. J. Appl. Meteorol. 42, 1514–1518.
Fig. 12. Comparing sophisticated and simplified approaches for the Midwestern US superstation. 443
H. Nguyen Sinh et al.
Journal of Wind Engineering & Industrial Aerodynamics 184 (2019) 436–444 Lombardo, F.T., Nguyen Sinh, H., Letchford, C.W., Rosowsky, D.V., 2014. “Assessing the Joint Wind and Temperature Hazard for the United States”. Structure Congress 2014, Boston, MA. Matalas, N.C., 1967. Mathematical assessment of synthetic hydrology. J. Water Resour. Res. 3 (4), 937–945. Nguyen Sinh, H., Lombardo, F.T., Letchford, C.W., Rosowsky, D.V., 2016. Characterization of the joint wind and ice hazard in the midwestern United States. Nat. Hazards Rev. 17 (3) ( J. ASCE). Peterka, J.A., 1992. Improved extreme wind prediction for the United States. J. Wind Eng. Ind. Aerod. 41–44, 533–541. Richardson, C.W., 1981. Stochastic simulation of daily precipitation, temperature, and solar radiation. J. Water Resour. Res. 17 (1), 182–190. Roberts, E., Stewart, R.E., 2008. On the occurrence of freezing rain and ice pellets over the Eastern Canadian Arctic. Atmos. Res. 89 (1–2), 93–109. Smith, R.E., Schreiber, H.A., 1973. Point processes of seasonal thunderstorm rainfall, 1, distribution of rainfall event. J. Water Resour. Res. 9 (4), 871–884. Weiss, L.L., 1964. Sequence of wet and dry days described by a Markov chain probability model. Monthly Weather Review- American Meteorological Society 92, 169–176.
Chin, E.H., 1957. Modeling daily precipitation occurrence process with Markov chain. J. Water Resour. Res. 13 (6), 949–956. Gabriel, R., Neuman, J., 1962. A Markov chain model for daily rainfall occurrence at Tel Aviv, Israel. Q. J. R. Meteorol. Soc. 88, 90–95. Goodwin III, E.J., Mozer, J.D., Di Gioia Jr., A.M., Power, B.A., 1983. Predicting ice and snow loads for transmission lines. In: The Proceedings of the First IWAIS, pp. 267–273. Haan, C.T., Allen, D.M., Street, J.O., 1976. A Markov chain model of daily rainfall. J. Water Resour. Res. 12 (3), 443–449. Hopkins, J.W., Robillard, P., 1964. Some statistics of daily rainfall occurrence for the Canadian Prairi provinces. Journal of Applied Meteorology and Climatology 3, 600–602. Integrated surface global hourly data (ISH), 2001. National Oceanic and Atmospheric Administration (NOAA). Retrieved from. ftp://ftp.ncdc.noaa.gov/pub/data. (Accessed 1 January 2001). ISO 12494, 2017. International Standard. Atmospheric Icing of Structures. Jones, K.F., 1998. A simple model for freezing rain ice loads. Atmos. Res. 46, 87–97.
444