Calculate the storm’s forward speed

First, the code determines the forward speed of the storm at each observation, in meters per second. The forward speed of the storm needs to be removed from the observed maximum wind speed to get an estimate of the maximum wind speed associated just with the rotational movement of the storm, which is what needs to go into the Willoughby model (Willoughby, Darling, and Rahn 2006). This asymmetry needs to be removed before we convert winds to gradient level. The code adds back in a forward motion component to the surface winds at grid points near the end of the the modeling process.

The code uses the Haversine method with great circle distance to calculate the distance, in kilometers, between the storm’s latitude and longitude coordinates for the current observation and the following observation. Then the time difference is determined between those two observations. From these, the storm’s forward speed is calculated. This speed is converted to meters per second (from kilometers per hour), so it’s in the same units as \(V_{max}\) from the imputed storm tracks.

Here is the equation used in the code for the Haversine method with great circle distance to calculate the distance between two locations based on their latitudes and longitude:

\[ hav(\gamma) = hav(\phi_1 - \phi_2) + cos(\phi_1)*cos(\phi_2)*hav(L_1 - L_2) \] \[ D = R_{earth} * \gamma \]

where:

• \(\phi_{1,rad}\): Latitude of first location, in radians
• \(L_{1,rad}\): Longitude of first location, in radians
• \(\phi_{2,rad}\): Latitude of second location, in radians
• \(L_{2,rad}\): Longitude of second location, in radians
• \(\gamma\): Intermediary result
• \(hav(\gamma)\): The haversine function, \(hav(\gamma) = sin^2 \left(\frac{\gamma}{2}\right)\)
• \(R_{earth}\): Radius of the earth, here assumed to be 6378.14 kilometers
• \(D\): Distance between the two locations, in kilometers

In the package, there is a helper function to convert between degrees and radians (degrees_to_radians), as well as a function called latlon_to_km that uses the Haversine method to calculate the distance between two storm center locations (tclat_1 and tclat_2 are the two storm center latitudes; tclon_1 and tclon_2 are the two storm center longitudes):

stormwindmodel:::degrees_to_radians

## function (degrees) 
## {
##     degrees * pi/180
## }
## <bytecode: 0x7fe8d4a79638>
## <environment: namespace:stormwindmodel>

stormwindmodel:::latlon_to_km

## function (tclat_1, tclon_1, tclat_2, tclon_2, Rearth = 6378.14) 
## {
##     tclat_1 <- degrees_to_radians(tclat_1)
##     tclon_1 <- degrees_to_radians(tclon_1)
##     tclat_2 <- degrees_to_radians(tclat_2)
##     tclon_2 <- degrees_to_radians(tclon_2)
##     delta_L <- tclon_1 - tclon_2
##     delta_tclat <- tclat_1 - tclat_2
##     hav_L <- sin(delta_L/2)^2
##     hav_tclat <- sin(delta_tclat/2)^2
##     hav_gamma <- hav_tclat + cos(tclat_1) * cos(tclat_2) * hav_L
##     gamma <- 2 * asin(sqrt(hav_gamma))
##     dist <- Rearth * gamma
##     return(dist)
## }
## <bytecode: 0x7fe8d49d2608>
## <environment: namespace:stormwindmodel>

The package uses these functions in a function called calc_forward_speed that calculates the distance between two storm locations, calculates the time difference in their observation times, and from that determines the forward speed in kilometers per hour and converts it to meters per second.

stormwindmodel:::calc_forward_speed

## function (tclat_1, tclon_1, time_1, tclat_2, tclon_2, time_2) 
## {
##     dist <- latlon_to_km(tclat_1, tclon_1, tclat_2, tclon_2) * 
##         1000
##     time <- as.numeric(difftime(time_2, time_1, units = "secs"))
##     tcspd <- dist/time
##     return(tcspd)
## }
## <bytecode: 0x7fe8d4275580>
## <environment: namespace:stormwindmodel>

Calculate direction of storm movement

Next, the code calculates the direction of the motion of the storm (storm bearing). Later code will use this to add back in a component for the forward motion of the storm into the wind estimates. The following equation is used to calculate the bearing of one point from another point based on latitude and longitude. So, this is calculating the bearing of a later storm observation, as seen from an earlier storm observation. The function restricts the output to be between 0 and 360 degrees using modular arithmetic (%% 360).

\[ S = cos(\phi_{2,rad}) * sin(L_{1,rad} - L_{2,rad}) \]

\[ C = cos(\phi_{1,rad}) * sin(\phi_{2,rad}) - sin(\phi_{1,rad}) * cos(\phi_{2,rad}) * cos(L_{1,rad} - L_{2,rad}) \]

\[ \theta_{F} = atan2(S, C) * \frac{180}{\pi} + 90 \]

where:

• \(\phi_{1,rad}\): Latitude of first location, in radians
• \(L_{1,rad}\): Longitude of first location, in radians
• \(\phi_{2,rad}\): Latitude of second location, in radians
• \(L_{2,rad}\): Longitude of second location, in radians
• \(S\), \(C\): Intermediary results
• \(\theta_F\) is the direction of the storm movement, in degrees

This calculation is implemented in our package through the calc_bearing function:

stormwindmodel::calc_bearing

## function (tclat_1, tclon_1, tclat_2, tclon_2) 
## {
##     tclat_1 <- degrees_to_radians(tclat_1)
##     tclon_1 <- degrees_to_radians(-tclon_1)
##     tclat_2 <- degrees_to_radians(tclat_2)
##     tclon_2 <- degrees_to_radians(-tclon_2)
##     S <- cos(tclat_2) * sin(tclon_1 - tclon_2)
##     C <- cos(tclat_1) * sin(tclat_2) - sin(tclat_1) * cos(tclat_2) * 
##         cos(tclon_1 - tclon_2)
##     theta_rad <- atan2(S, C)
##     theta <- radians_to_degrees(theta_rad) + 90
##     theta <- theta%%360
##     return(theta)
## }
## <bytecode: 0x7fe8d4069778>
## <environment: namespace:stormwindmodel>

Calculate u- and v-components of forward speed

Next, the code uses the estimated magnitude and direction of the storm’s forward speed to calculate u- and v-components of this forward speed. Later, these two components are added back to the modeled surface wind speed at grid points, after adjusting with the Phadke correction factor for forward motion (Phadke et al. 2003).

To calculate the u- and v-components of forward motion, \(F_{u}\) and \(F_{v}\), the function uses:

\[ F_{u} = \left|\overrightarrow{F}\right| cos(\theta_F) \]

\[ F_{v} = \left|\overrightarrow{F}\right| sin(\theta_F) \]

where \(\theta_F\) is the direction of the storm movement (0 for due east, 90 for due north, etc.).

tcspd_u = tcspd * cos(degrees_to_radians(tcdir))
tcspd_v = tcspd * sin(degrees_to_radians(tcdir))

Adjust wind speed to remove forward motion of storm

Next, there’s a function called remove_forward_speed that uses the estimated forward speed to adjust wind speed to remove the component from forward motion, based on equation 12 (and accompanying text) from Phadke et al. (2003).

Because here the code is adjusting the maximum sustained wind to remove the component from forward speed, we can assume that we are adjusting wind speed at the radius of maximum winds and for winds blowing in the same direction as the direction of the forward motion of the storm. Therefore, the correction term for forward motion is \(U = 0.5F\) (Phadke et al. 2003), or half the total forward speed of the storm. This correction factor is subtracted from the maximum sustained wind speed to remove the forward component, and because we assume that the maximum winds are blowing in the same direction as the direction of the storm’s forward motion, this correction term is directly subtracted from the magnitude of the wind speed (rather than breaking into u- and v-components for a vector addition).

If \(V_{max}\) after removing the forward storm motion is ever negative, the remove_forward_speed function resets the value to 0 m / s.

stormwindmodel:::remove_forward_speed

## function (vmax, tcspd) 
## {
##     vmax_sfc_sym <- vmax - 0.5 * tcspd
##     vmax_sfc_sym[vmax_sfc_sym < 0] <- 0
##     return(vmax_sfc_sym)
## }
## <bytecode: 0x7fe8d3aaa670>
## <environment: namespace:stormwindmodel>

Convert symmetric surface wind to gradient wind

The extended tracks database gives maximum winds as 1-minute sustained wind speeds 10 meters above the ground. To make calculations easier (by not having to deal with friction), the code converts the symmetric 1-minute sustained wind speeds at 10 meters (\(V_{max,sym}\)) to gradient wind speed (\(V_{max,G}\)), using the following equation:

\[ V_{max,G} = \frac{V_{max,sym}}{f_r} \]

where:

• \(V_{max,G}\): Maximum gradient-level 1-minute sustained wind (m/s)
• \(V_{max,sym}\): Maximum 10-m 1-minute sustained wind with motion asymmetry removed (m/s)
• \(f_r\): Reduction factor for \(r \le\) 100 km. At this point, the winds are calculated for \(r = R_{max}\).

The code uses the reduction factor from Figure 3 of Knaff et al. (2011). We assume that \(R_{max}\) is always 100 km or less. Then, \(f_r\) is 0.9 if the storm’s center is over water and 80% of that, 0.72, if the storm’s center is over land (Knaff et al. 2011).

stormwindmodel::calc_gradient_speed

## function (vmax_sfc_sym, over_land) 
## {
##     reduction_factor <- 0.9
##     if (over_land) {
##         reduction_factor <- reduction_factor * 0.8
##     }
##     vmax_gl <- vmax_sfc_sym/reduction_factor
##     return(vmax_gl)
## }
## <bytecode: 0x7fe8d381ca20>
## <environment: namespace:stormwindmodel>

To determine if the storm center is over land or over water, the package includes a dataset called landmask, with locations of grid points in the eastern US and offshore, with a factor saying whether each point is over land or water:

data(landmask)
head(landmask)

## # A tibble: 6 x 3
##   longitude latitude land 
##       <dbl>    <dbl> <fct>
## 1       260     15   water
## 2       260     15.2 water
## 3       260     15.4 water
## 4       260     15.6 water
## 5       260     15.8 water
## 6       260     16   water

The check_over_land function finds the closest grid point for a storm location and determines whether the storm is over land or over water at its observed location.

stormwindmodel:::check_over_land

## function (tclat, tclon) 
## {
##     lat_diffs <- abs(tclat - stormwindmodel::landmask$latitude)
##     closest_grid_lat <- stormwindmodel::landmask$latitude[which(lat_diffs == 
##         min(lat_diffs))][1]
##     lon_diffs <- abs(tclon - (360 - stormwindmodel::landmask$longitude))
##     closest_grid_lon <- stormwindmodel::landmask$longitude[which(lon_diffs == 
##         min(lon_diffs))][1]
##     over_land <- stormwindmodel::landmask %>% dplyr::filter(.data$latitude == 
##         closest_grid_lat & .data$longitude == closest_grid_lon) %>% 
##         dplyr::mutate(land = .data$land == "land") %>% dplyr::select(.data$land)
##     over_land <- as.vector(over_land$land[1])
##     return(over_land)
## }
## <bytecode: 0x7fe8d5de2730>
## <environment: namespace:stormwindmodel>

Here is an example of applying this function to the tracks for Hurricane Floyd:

floyd_tracks$land <- mapply(stormwindmodel:::check_over_land,
                            tclat = floyd_tracks$latitude,
                            tclon = -floyd_tracks$longitude)
ggplot(landmask, aes(x = longitude - 360, y = latitude, color = land)) +
  geom_point() + 
  geom_point(data = floyd_tracks, aes(x = longitude, y = latitude, 
                                     color = NULL, shape = land)) + 
  scale_color_discrete("Land mask") + 
  scale_shape_discrete("Track over land")

In the add_wind_radii function, there is code to check whether the point is over land or water and apply the proper reduction factor:

over_land = mapply(check_over_land, tclat, tclon),
vmax_gl = mapply(calc_gradient_speed,
                vmax_sfc_sym = vmax_sfc_sym,
                over_land = over_land)

Calculate radius of maximum wind speed

Next, the package code calculates \(R_{max}\), the radius of maximum wind speed, using Eq. 7a from Willoughby et al. (2006):

\[ R_{max} = 46.4 e^{- 0.0155 V_{max,G} + 0.0169\phi} \]

where:

• \(R_{max}\): Radius from the storm center to the point at which the maximum wind occurs (km)
• \(V_{max,G}\): Maximum gradient-level 1-min sustained wind (m / s)
• \(\phi\): Latitude, in decimal degrees

stormwindmodel::will7a

## function (vmax_gl, tclat) 
## {
##     Rmax <- 46.4 * exp(-0.0155 * vmax_gl + 0.0169 * tclat)
##     return(Rmax)
## }
## <bytecode: 0x7fe8d5d62a78>
## <environment: namespace:stormwindmodel>

Calculate parameters for Willoughby model

Next, the code calculates \(X_1\), which is a parameter needed for the Willoughby model. This is done using equation 10a from Willoughby et al. (2006):

\[ X_1 = 317.1 - 2.026V_{max,G} + 1.915 \phi \]

where:

• \(X_1\): Parameter for Willoughby wind model
• \(V_{max,G}\): Maximum gradient-level 1-min sustained wind (m / s)
• \(\phi\): Latitude, in decimal degrees

stormwindmodel::will10a

## function (vmax_gl, tclat) 
## {
##     X1 <- 317.1 - 2.026 * vmax_gl + 1.915 * tclat
##     return(X1)
## }
## <bytecode: 0x7fe8d5eb19c8>
## <environment: namespace:stormwindmodel>

Next, the code calculates another Willoughby parameter, \(n\). This is done with equation 10b from Willoughby et al. (2006):

\[ n = 0.4067 + 0.0144 V_{max,G} - 0.0038 \phi \]

where:

• \(n\): the exponential for the power law inside the eye
• \(V_{max,G}\): Maximum gradient-level 1-min sustained wind (m / s)
• \(\phi\): Latitude, in decimal degrees

stormwindmodel::will10b

## function (vmax_gl, tclat) 
## {
##     n <- 0.4067 + 0.0144 * vmax_gl - 0.0038 * tclat
##     return(n)
## }
## <bytecode: 0x7fe8d5ff7d80>
## <environment: namespace:stormwindmodel>

Next, the code calculates another Willoughby parameter, \(A\), with equation 10c from Willoughby et al. (2006):

\[ A = 0.0696 + 0.0049 V_{max,G} - 0.0064 \phi \] \[ A = \begin{cases} 0 & \text{ if } A < 0\\ A & \text{ otherwise } \end{cases} \]

where:

• \(A\): A parameter for the Willoughby model
• \(V_{max,G}\): Maximum gradient-level 1-min sustained wind (m / s)
• \(\phi\): Absolute value of latitude (degrees)

stormwindmodel::will10c

## function (vmax_gl, tclat) 
## {
##     A <- 0.0696 + 0.0049 * vmax_gl - 0.0064 * tclat
##     A[A < 0 & !is.na(A)] <- 0
##     return(A)
## }
## <bytecode: 0x7fe8d89d4758>
## <environment: namespace:stormwindmodel>

Determine radius of start of transition region

Now, the code uses a numerical method to determine the value of \(R_1\), the radius to the start of the transition region, for the storm for a given track observation. The requires finding the root of this equation (Willoughby, Darling, and Rahn 2006):

\[ w - \frac{n ((1 - A) X_1 + 25 A)}{n ((1 - A) X_1 + 25 A) + R_{max}} = 0 \]

where:

• \(w\): The weighting function (a function of \(R_{max}\), \(R_1\), and \(R_2 - R_1\); see below)
• \(n\): A parameter for the Willoughby model
• \(R_{max}\): Radius at which the maximum wind occurs (km)
• \(A\): A parameter for the Willoughby model
• \(X_1\): A parameter for the Willoughby model
• \(X_2\): A parameter for the Willoughby model, set to 25 (Willoughby, Darling, and Rahn 2006)

The weighting function, \(w\), is:

\[ w = \begin{cases} 0 & \text{if } \xi < 0 \\ 126 \xi^5 - 420 \xi^6 + 540 \xi^7- 315 \xi^8 + 70 \xi^9 & \text{if } 0 \le \xi \le 1\\ 1 & \text{if } \xi > 1 \end{cases} \]

where:

\[ \xi = \frac{R_{max} - R_1}{R_2 - R_1} \]

and where:

• \(w\): The weighting function
• \(\xi\): A nondimensional argument
• \(R_{max}\): radius at which the maximum wind occurs (km)
• \(R_1\): lower boundary of the transition zone (in km from the storm center)
• \(R_2\): upper boundary of the transition zone (in km from the storm center)

To find the root of this equation, the package has a function that uses the Newton-Raphson method (Press et al. 2002; Jones, Maillardet, and Robinson 2009). This function starts with an initial guess of the root, \(x_0\), then calculates new values of \(x_{n+1}\) using:

\[ x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)} \]

This process iterates until either the absolute value of \(f(x_{n+1})\) is smaller than some threshold (\(\epsilon\); default of 0.001 in our function) or the maximum allowed number of iterations is reached (default of 100 iterations in our function).

In this case, we’re using it to try to find the value of \(\xi\) that is a root of the following function:

\[ f(\xi) = 126 \xi^5 - 420 \xi^6 + 540 \xi^7- 315 \xi^8 + 70 \xi^9 - \frac{n ((1 - A) X_1 + 25 A)}{n ((1 - A) X_1 + 25 A) + R_{max}} \]

The derivative of this function, which is also needed for the Newton-Raphson method, is:

\[ f'(\xi) = 5 * 126 \xi^4 - 6 * 420 \xi^5 + 7 * 540 \xi^6 - 8 * 315 \xi^7 + 9 * 70 \xi^8 \]

stormwindmodel::will3_right

## function (n, A, X1, Rmax) 
## {
##     eq3_right <- (n * ((1 - A) * X1 + 25 * A))/(n * ((1 - A) * 
##         X1 + 25 * A) + Rmax)
##     return(eq3_right)
## }
## <bytecode: 0x7fe8d8a98308>
## <environment: namespace:stormwindmodel>

stormwindmodel::will3_deriv_func

## function (xi, eq3_right) 
## {
##     deriv <- 70 * 9 * xi^8 - 315 * 8 * xi^7 + 540 * 7 * xi^6 - 
##         420 * 6 * xi^5 + 126 * 5 * xi^4
##     func <- 70 * xi^9 - 315 * xi^8 + 540 * xi^7 - 420 * xi^6 + 
##         126 * xi^5 - eq3_right
##     deriv_func <- c(deriv, func)
##     return(deriv_func)
## }
## <bytecode: 0x7fe8d8aeae78>
## <environment: namespace:stormwindmodel>

stormwindmodel::solve_for_xi

## function (xi0 = 0.5, eq3_right, eps = 0.001, itmax = 100) 
## {
##     if (is.na(eq3_right)) {
##         return(NA)
##     }
##     else {
##         i <- 1
##         xi <- xi0
##         while (i <= itmax) {
##             deriv_func <- will3_deriv_func(xi, eq3_right)
##             if (abs(deriv_func[2]) <= eps) {
##                 break
##             }
##             xi <- xi - deriv_func[2]/deriv_func[1]
##         }
##         if (i < itmax) {
##             return(xi)
##         }
##         else {
##             warning("Newton-Raphson did not converge.")
##             return(NA)
##         }
##     }
## }
## <bytecode: 0x7fe8d8b43348>
## <environment: namespace:stormwindmodel>

While the Newton-Raphson method can sometimes perform poorly in finding global maxima, in this case the function for which we are trying to find the root is well-behaved. Across tropical storms from 1988 to 2015, the method never failed to converge, and identified roots were consistent across storms (typically roots are for \(\xi\) of 0.6–0.65). We also tested using the optim function in the stats R package and found similar estimated roots but slower convergence times than when using the Newton-Raphson method.

Now an equation from the Willoughby et al. 2006 paper can be used to calculate \(R_1\) (Willoughby, Darling, and Rahn 2006):

\[ R_1 = R_{max} - \xi(R_2 - R_1) \]

For this function, the package code assumes that \(R_2 - R_1\) (the width of the transition region) is 25 kilometers when \(R_{max}\) is larger than 20 kilometers and 15 kilometers otherwise.

stormwindmodel::calc_R1

## function (Rmax, xi) 
## {
##     R2_minus_R1 <- ifelse(Rmax > 20, 25, 15)
##     R1 <- Rmax - xi * R2_minus_R1
##     return(R1)
## }
## <bytecode: 0x7fe8d8cba168>
## <environment: namespace:stormwindmodel>

Determine radius for end of transition region

Next, the code estimates the radius for the end of the transition period, \(R_2\). We assume that smaller storms (\(R_{max} \le 20\mbox{ km}\)) have a transition region of 15 km while larger storms (\(R_{max} > 20\mbox{ km}\)) have a transition region of 25 km:

\[ R_2 = \begin{cases} R_1 + 25 & \text{ if } R_{max} > 20\mbox{ km}\\ R_1 + 15 & \text{ if } R_{max} \le 20\mbox{ km} \end{cases} \]

where:

• \(R_1\): Radius to the start of the transition region (km)
• \(R_2\): Radius to the end of the transition region (km)

R2 = ifelse(Rmax > 20, R1 + 25, R1 + 15)

Get the radius from the storm center to the grid point

The first step is to determine the distance between the grid point and the center of the storm (\(C_{dist}\)). This is done using the Haversine method of great circle distance to calculate this distance. For this, the package uses the latlon_to_km function given in the code shown earlier.

Determine gradient wind speed at each location

Next, the package calculates \(V_G(r)\), the gradient level 1-minute sustained wind at the grid point, which is at radius \(r\) from the tropical cyclone center (\(C_{dist}\) for the grid point). Note there are different equations for \(V_G(r)\) for (1) the eye to the start of the transition region; (2) outside the transition region; and (3) within the transition region.

First, if \(r \le R_1\) (inside the transition region), \(V_G(r)\) is calculated using (equation 1a, Willoughby, Darling, and Rahn 2006):

\[ V_G(r) = V_i = V_{max,G} \left( \frac{r}{R_{max}} \right)^n, (0 \le r \le R_1) \]

Next, if \(r \ge R_2\) (outside the transition region), \(V_G(r)\) is calculated using (equation 4, dual-exponential replacement of equation 1b, Willoughby, Darling, and Rahn 2006):

\[ V_G(r) = V_o = V_{max,G}\left[(1 - A) e^\frac{R_{max} - r}{X_1} + A e^\frac{R_{max} - r}{X_2}\right], R_2 < r \]

Finally, if \(r\) is between \(R_1\) and \(R_2\), first \(\xi\) must be calculated for the grid point (Willoughby, Darling, and Rahn 2006):

\[ \xi = \frac{r - R_1}{R_2 - R_1} \]

and then:

\[ w = \begin{cases} 0 & \text{if } \xi < 0 \\ 126 \xi^5 - 420 \xi^6 + 540 \xi^7- 315 \xi^8 + 70 \xi^9 & \text{if } 0 \le \xi \le 1\\ 1 & \text{if } \xi > 1 \end{cases} \]

where:

• \(w\): Weighting variable
• \(\xi\): A nondimensional argument
• \(r\): radius from the storm center to the grid point (in km)
• \(R_{max}\): radius at which the maximum wind occurs (km)
• \(R_1\): lower boundary of the transition zone (in km from the storm center)
• \(R_2\): upper boundary of the transition zone (in km from the storm center)

The, \(V(r)\) is calculated using (equation 1c, Willoughby, Darling, and Rahn 2006):

\[ V_G(r) = V_i (1 - w) + V_o w, (R_1 \le r \le R_2) \]

All of this is wrapped in a function:

stormwindmodel:::will1

## function (cdist, Rmax, R1, R2, vmax_gl, n, A, X1, X2 = 25) 
## {
##     if (is.na(Rmax) || is.na(vmax_gl) || is.na(n) || is.na(A) || 
##         is.na(X1)) {
##         return(NA)
##     }
##     else {
##         Vi <- vmax_gl * (cdist/Rmax)^n
##         Vo <- vmax_gl * ((1 - A) * exp((Rmax - cdist)/X1) + A * 
##             exp((Rmax - cdist)/X2))
##         if (cdist < R1) {
##             wind_gl_aa <- Vi
##         }
##         else if (cdist > R2) {
##             wind_gl_aa <- Vo
##         }
##         else {
##             eps <- (cdist - R1)/(R2 - R1)
##             w <- 126 * eps^5 - 420 * eps^6 + 540 * eps^7 - 315 * 
##                 eps^8 + 70 * eps^9
##             wind_gl_aa <- Vi * (1 - w) + Vo * w
##         }
##         wind_gl_aa[wind_gl_aa < 0 & !is.na(wind_gl_aa)] <- 0
##         return(wind_gl_aa)
##     }
## }
## <bytecode: 0x7fe8d7db3f38>
## <environment: namespace:stormwindmodel>

Estimate surface level winds from gradient winds

Now, to get estimated symmetric surface-level tangential winds from gradient-levels winds, the code uses a method from Knaff et al. (2011).

In this case, we can always assume the point is over land, not water, so the code always uses a reduction factor that is reduced by 20%. This method uses a reduction factor of 0.90 up to a radius of 100 km, a reduction factor of 0.75 for any radius 700 km or greater, and a linear decreasing reduction factor for any radius between those two radius values (Knaff et al. 2011). If the point is over land (true for any county), this reduction factor is further reduced by 20%.

Here is the function the code uses for this part:

stormwindmodel:::gradient_to_surface

## function (wind_gl_aa, cdist) 
## {
##     if (cdist <= 100) {
##         reduction_factor <- 0.9
##     }
##     else if (cdist >= 700) {
##         reduction_factor <- 0.75
##     }
##     else {
##         reduction_factor <- 0.9 - (cdist - 100) * (0.15/600)
##     }
##     reduction_factor <- reduction_factor * 0.8
##     wind_sfc_sym <- wind_gl_aa * reduction_factor
##     return(wind_sfc_sym)
## }
## <bytecode: 0x7fe8d4e1a998>
## <environment: namespace:stormwindmodel>

Here is a reproduction of part of Figure 3 from Knaff et al. (2011) using this function:

rf_example <- data.frame(r = 0:800,
                         rf = mapply(
                           stormwindmodel:::gradient_to_surface, 
                           wind_gl_aa = 1, cdist = 0:800))
ggplot(rf_example, aes(x = r, y = rf)) + 
  geom_line() + 
  theme_classic() + 
  xlab("Radius (km)") + 
  ylab("Reduction factor") + 
  ylim(c(0.5, 0.9))

Calculate the direction of gradient winds at each location

So far, the model has just calculated the the rotational component of the wind speed at each grid location (tangential wind). Now the code adds back asymmetry from the forward motion component of the wind speed (translational wind) at each location. The total storm winds will be increased by forward motion of the storm on the right side of the storm. Conversely, on the left side of the storm, the forward motion of the storm will offset some of the storm-related wind.

We’ve already determined the direction of the forward direction of the storm. Now, we need to get the direction of the storm winds at the grid point location.

First, the function determines, for each storm location, the direction from the storm center to the location. For this, the package uses the calc_bearing function described earlier to determine the bearing of the storm at each storm observation. Now, it is used to calculate the angle from the storm center to the grid point:

stormwindmodel:::calc_bearing

## function (tclat_1, tclon_1, tclat_2, tclon_2) 
## {
##     tclat_1 <- degrees_to_radians(tclat_1)
##     tclon_1 <- degrees_to_radians(-tclon_1)
##     tclat_2 <- degrees_to_radians(tclat_2)
##     tclon_2 <- degrees_to_radians(-tclon_2)
##     S <- cos(tclat_2) * sin(tclon_1 - tclon_2)
##     C <- cos(tclat_1) * sin(tclat_2) - sin(tclat_1) * cos(tclat_2) * 
##         cos(tclon_1 - tclon_2)
##     theta_rad <- atan2(S, C)
##     theta <- radians_to_degrees(theta_rad) + 90
##     theta <- theta%%360
##     return(theta)
## }
## <bytecode: 0x7fe8d4069778>
## <environment: namespace:stormwindmodel>

Next, the function calculates the gradient wind direction based on the bearing of a location from the storm. This gradient wind direction is calculated by adding 90 degrees to the bearing of the grid point from the storm center.

gwd = (90 + chead) %% 360

Calculate the surface wind direction

The next step is to change from the gradient wind direction to the surface wind direction. To do this, the function adds an inflow angle to the gradient wind direction (making sure the final answer is between 0 and 360 degrees). This step is necessary because surface friction changes the wind direction near the surface compared to higher above the surface.

The inflow angle is calculated as a function of the distance from the storm center to the location and the storm’s \(R_{max}\) at that observation point (eq. 11a–c, Phadke et al. 2003):

\[ \beta = \begin{cases} & 10 + \left(1 + \frac{r}{R_{max}}\right) \text{ if } r < R_{max}\\ & 20 + 25\left(\frac{r}{R_{max}} - 1 \right ) \text{ if } R_{max} \le r \le 1.2R_{max}\\ & 25 \text{ if } r \ge 1.2R_{max} \end{cases} \]

Over land, the inflow angle should be about 20 degrees more than it is over water. Therefore, after calculating the inflow angle from the equation above, the function adds 20 degrees to the value since all modeled locations are over land. The final calculation, then, is:

\[ \theta_S = \theta_G + \beta + 20 \]

Here is the function for this step:

add_inflow

## function (gwd, cdist, Rmax) 
## {
##     if (is.na(gwd) | is.na(cdist) | is.na(Rmax)) {
##         return(NA)
##     }
##     if (cdist < Rmax) {
##         inflow_angle <- 10 + (1 + (cdist/Rmax))
##     }
##     else if (Rmax <= cdist & cdist < 1.2 * Rmax) {
##         inflow_angle <- 20 + 25 * ((cdist/Rmax) - 1)
##     }
##     else {
##         inflow_angle <- 25
##     }
##     overland_inflow_angle <- inflow_angle + 20
##     gwd_with_inflow <- (gwd + overland_inflow_angle)%%360
##     return(gwd_with_inflow)
## }
## <bytecode: 0x7fe8d9816fd0>
## <environment: namespace:stormwindmodel>

Add back in wind component from forward speed of storm

Next, to add back in the storm’s forward motion at each grid point, the code reverses the earlier step that used the Phadke correction factor (equation 12, Phadke et al. 2003). The package calculates a constant correction factor (correction_factor), as a function of r, radius from the storm center to the grid point, and \(R_{max}\), radius from storm center to maximum winds.

\[ U(r) = \frac{R_{max}r}{R_{max}^2 + r^2}F \] where:

• \(U(r)\): Phadke correction factor at a radius \(r\) from the storm’s center
• \(r\): Radius from the storm’s center (km)
• \(R_{max}\): Radius of maximum winds (km)
• \(F\): Forward speed of the storm (m / s)

Both the u- and v-components of forward speed are corrected with this correction factor, then these components are added to the u- and v-components of tangential surface wind. These u- and v-components are then used to calculate the magnitude of total wind associated with the storm at the grid point:

stormwindmodel:::add_forward_speed

## function (wind_sfc_sym, tcspd_u, tcspd_v, swd, cdist, Rmax) 
## {
##     wind_sfc_sym_u <- wind_sfc_sym * cos(degrees_to_radians(swd))
##     wind_sfc_sym_v <- wind_sfc_sym * sin(degrees_to_radians(swd))
##     correction_factor <- (Rmax * cdist)/(Rmax^2 + cdist^2)
##     wind_sfc_u <- wind_sfc_sym_u + correction_factor * tcspd_u
##     wind_sfc_v <- wind_sfc_sym_v + correction_factor * tcspd_v
##     wind_sfc <- sqrt(wind_sfc_u^2 + wind_sfc_v^2)
##     wind_sfc <- ifelse(wind_sfc > 0, wind_sfc, 0)
##     return(wind_sfc)
## }
## <bytecode: 0x7fe8d5123ef8>
## <environment: namespace:stormwindmodel>

Calculate 3-second gust wind speed from sustained wind speed

The last step in the code calculates gust wind speed from sustained wind speed by applying a gust factor:

\[ V_{sfc,gust} = G_{3,60}*V_{sfc} \]

where \(V_{sfc}\) is the asymmetric surface wind speed, \(V_{sfc,gust}\) is the surface gust wind speed, and \(G_{3,60}\) is the gust factor.

Here is a table with gust factors based on location (Harper, Kepert, and Ginger 2010):

The stormwindmodel package uses the “in-land” gust factor value throughout.

gustspeed = windspeed * 1.49