Wind Forecasts Using Large Eddy Simulations for Stratospheric Balloon Applications

(1)

Wind Forecasts Using Large Eddy Simulations for Stratospheric Balloon

Applications

Ludvig Sjöberg

Space Engineering, master's level 2019

Luleå University of Technology

Department of Computer Science, Electrical and Space Engineering

(2)

(3)

i

LULEÅ UNIVERSITY OF TECHNOLOGY

Abstract

Department of Computer Science, Electrical and Space Engineering

Master of Science in Engineering in Space Engineering

Wind Forecasts Using Large Eddy Simulations for Stratospheric Balloon Applications by Ludvig SJÖBERG

The launch of large stratospheric balloons is highly dependant on the meteorological conditions at ground level, including wind speed. The balloon launch base Esrange Space Center in northern Sweden currently uses forecasts delivered through the Swedish Mete- orological and Hydrological Institute to predict opportunities for balloon launches. How- ever the staff at Esrange Space Center experience that the current forecasts are not accurate enough. For that reason the Weather Research and Forecasting model is used to improve the forecast. The model performs a Large Eddy Simulation over the area closest to Esrange Space Center to predict wind speed and turbulence. During twelve hypothetical launch days the improved forecast have an overall accuracy of 93% compared to the old forecast accuracy of 69%. With some improvements and the right computational power the system is thought to be operationally viable.

(4)

(5)

iii

Acknowledgements

Thank you to Ricardo Fonseca for all the help with introducing me to the Weather Re- search and Forecasting Model, as well as answering any questions along the way.

I would also give large thanks to Roberto Mantas Nakhai for all the discussions and for helping me out in a pinch.

Thank you to my supervisor Kent Andersson for supporting me throughout this thesis and for the feedback provided.

Last but not least thank you to my examiner Mathias Milz for providing feedback as well as taking the time to help with secure computational resources.

The simulations were performed on resources provided by the Swedish National Infras- tructure for Computing (SNIC) at High Performance Computing Center North (HPC2N).

(6)

(7)

v

List of Figures

1.1 The different programs used in constructing the WRF forecast. . . 2 2.1 Comparison between the best resolution default dataset and the topography

adapted for WRF. . . 7 2.2 Comparison between MODIS lake category and the lake category adapted

for WRF. Yellow color represents land and blue color represents water.. . . . 8 3.1 Comparison of actual and predicted 10-meter wind speed during the launch

window on 2016-01-19 . . . 13 3.2 Comparison of actual and predicted 10-meter wind speed during the launch

window on 2017-02-15. . . 14 3.3 Comparison of actual and predicted 10-meter wind speed during the launch

window on 2017-12-07. . . 24 3.13 The angle represents variance similarity and the radial distance represents

absolute value of the correlation. The optimal model performance is marked by the black star. . . 26 3.14 Plot representing the bias and error variance connection to the normalized

root mean square error. The contours represent normalized RMSE in increments of 0.2. The optimal model performance is marked with a black star. . 27 4.1 Comparison of actual and predicted temperature in the Troposphere on 2017-

05-31. Sounding balloon released 07:43 UTC and WRF forecast with target time 07:43 UTC. . . 30

(10)

viii

4.2 Comparison of actual and predicted water vapor content on 2017-05-31.

Sounding balloon released 07:43 UTC and WRF forecast with target time 07:43 UTC. . . 31 4.3 Comparison between actual and predicted 100-meter wind speed before,

during and after the launch window on 2016-09-01.. . . 32 4.4 As figure 4.3 but with a simulation using a longer spin up time added. . . 33 4.5 Comparison of actual and predicted temperature of the lower atmosphere

on 2017-12-07. Sounding balloon released 09:57 UTC, WRF forecast with target time 09:57 UTC and ERA5 Reanalysis with target time 10:00 UTC.. . . 34 4.6 A buildup of the strong inversion after the initialization of the highest reso-

lution domain on 2017-12-07. Profiles taken every 30 minutes from start to end of domain. . . 35 4.7 Measured temperature at BPW compared with the finest resolution forecast

from domain initialization on 2017-12-07. . . 36 4.8 Average wind in the temperature inversion during the start and end of the

microscale forecast on 2017-12-07.. . . 37

(11)

ix

List of Tables

2.1 Domain setup used in the study. . . 5

2.2 Specifications of the WT instrument. . . 9

2.3 Specifications of the BPW instrument.. . . 10

2.4 A summary of the case days. All times are in UTC. . . 10

3.1 Evaluation of balloon launch criterion from the different sources on 2016- 01-19. A green color indicates a pass and a red color indicates a fail for the given hour. . . 14

(12)

x

3.13 Number of correct calls compared to measurements for the different balloon launch criteria. . . 25 3.14 Color key for figures 3.13 and 3.14 . . . 27

(13)

xi

List of Abbreviations

ACM2 Asymmetric Convective Model 2 AGL Above Ground Level

ARW Advanced Research WRF BP Balloon Pad

BPW Balloon Pad West

ECMWF European Center for Medium-Range Weather Forecasts ERA5 ECMWF ReAnalysis 5

ESC Esrange Space Center GFS Global Forecasting System LES Large Eddy Simulation

NCAR National Center for Atmospheric Research NWP Numerical Weather Prediction

PBL Planetary Boundry Layer RH Radar Hill

RMSE Root Mean Squared Error

SMHI Swedish Meteorological and Hydrological Institute TI Terra Incognita

TKE Turbulent Kinetic Eenergy

WRF Weather Research and Forecasting Model

WT Wind Tower

(14)

(15)

1

1 Introduction

This thesis is written as a final step to obtain a Master of Science in Engineering from Luleå University of Technology. The thesis is done in cooperation with SSC, a worldwide provider of advance space services. SSC owns and operates Esrange Space Center(ESC) which is a launching base for sounding rockets and balloons.

1.1 Background

The following sections provide information regarding the context of the thesis.

1.1.1 Esrange Space Center

Esrange Space Center is located in the north of Sweden, approximately 200 km north of the arctic circle. The first rocket was launched from ESC in 1966 and the first stratospheric balloon was launched in 1974. Since the beginning more than 500 rockets and more than 600 balloons have been launched from ESC with the number increasing every year.

1.1.2 Stratospheric Balloons

Stratospheric balloons can be used for a wide range of purposes, lifting payloads that weighs only a couple of kilos up to several tons. After launch the balloons rise to a float- ing altitude that is generally between 25-40 km and stay there drifting with the wind for the duration of the mission. A mission can be a short test of a drop body that only have a minimal flight time or a long circumpolar flight taking multiple weeks.

The launch of the balloons at ESC takes place from the Balloon Pad(BP). It is a large gravel field with an area of approximately 250 000 m². The balloon is laid out on the ground and inflated with helium. After the inflation process is complete the balloon and payload are released. Launching a stratospheric balloon requires good meteorological conditions. A launch with precipitation is generally avoided and the winds have to be both light and stable in direction. The exact meteorological requirements for a balloon launch differ depending on the type of balloon, but generally larger balloons requires better conditions for a successful launch.

1.1.3 Weather Prediction

The weather of the Earth is the result of many different forces and interactions. Different parameters govern the interactions depending on what length scale is considered. There- fore the dynamics are generally divided into three scales: the synoptic, meso and micro scales.

(16)

2 Chapter 1. Introduction

The synoptic scale meteorology focuses on large scale circulations driven by high and low pressure systems in the atmosphere. These systems are of global scale with wavelengths of 1000-4000 kilometers (American Meteorological Society,2019).

The mesoscale meteorology focuses on regional scale features such as warm or cold fronts as well as thunderstorms and weather generated by topography. The length scales of these systems are at least a few kilometers and up to what can be classified as synoptic (Ameri- can Meteorological Society,2019).

The microscale meteorology focuses on small scale features with length scales of 2 kilo- meters or less (American Meteorological Society,2019). On this scale the local atmospheric circulations are resolved. Local features such as topography or land use are generally very impactful on the microscale.

1.1.4 Weather Research and Forecasting Model

The Weather Research and Forecasting Model(WRF) is a Numerical Weather Prediction(NWP) system capable of many different types of meteorology. It contains the solver core Ad- vanced Research WRF(ARW) that is maintained by the American National Center for At- mospheric Research(NCAR). It consists of five different programs arranged in the flow shown in figure1.1.

Ungrib unpacks input meteorological data and transforms it into a format that WRF can use.

Geogrid construct a geographical domain to conduct the simulation in, the data is taken from libraries of static data.

Metgrid takes the output from "Ungrib" and "Geogrid" and constructs a combined earth- atmosphere model. It also interpolates data points horizontally.

Real interpolates the data points vertically and creates the initial and boundary conditions required.

WRF, which is the main program, performs the actual numerical integration and con- structs the forecast.

FIGURE 1.1: The different programs used in constructing the WRF forecast.

WRF-ARW is very flexible and can be configured for domain sizes ranging from global down to a few hundred meters. It also features nesting which allows for higher resolution domains to be driven by the output of lower resolution domains (Skamarock et al.,2008).

(17)

1.2. Problem Statement 3

1.2 Problem Statement

Currently the light wind required for the launch of a stratospheric balloon are predicted using synoptic and mesoscale forecasts delivered through the Swedish Meteorological and Hydrological Institute(SMHI). These forecasts are generally not accurate enough and it is common that a launch attempt has to be aborted due to higher than forecasted winds.

A theory is that the forecasts delivered by SMHI are too coarse in spatial resolution and do not capture the local features of ESC well enough to accurately forecast the wind. Can the forecast be improved using a model with finer resolution?

1.3 Objectives

The objective of the thesis is to develop a model using WRF that can predict conditions for balloon launches at ESC. The model output shall be verified against on-site measurements and compared to already available forecasts.

1.4 Research in the field

The technique of using nested NWP-models to produce an accurate local forecast has been used in many different situations. Many studies of microscale meteorology using WRF used a technique called Large Eddy Simulation(LES). This technique resolves some of the turbulent eddies in the flow, in contrast to a Planetary Boundary Layer(PBL) parameterization which does not resolve any turbulence (Wyngaard,2004). But as Wyngaard continued to point out, there exists a regime between the PBL parameterization and LES for which no method is designed. He coined this region Terra Incognita(TI) and a clear strategy on how to handle it does not exist.

In Mazzaro et al.,2017an attempt was made to investigate simulations through the TI.

The paper focused on convective structures and how they develop when going from a mesoscale simulation with a PBL parameterization to a LES simulation. The paper concluded that using the finest possible mesoscale resolution reduced the impact of the TI.

However only one LES domain was considered and the grid scale factor ^∆_∆^Meso

LE S was relatively large. As the paper pointed out there does not yet exist a clear best practise when scaling down from mesoscale simulations to LES, however a ratio of 3 is commonly used.

Liu et al.,2011used WRF with LES to downscale the synoptic scale Global Forecasting System(GFS) forecast through six domains ranging in size from 30 km down to 0.123 km, making the last domain a clear example of microscale meteorology. They ran the finest two nests in LES mode, while keeping other model variables the same for all nests. This was a similar approach to many other microscale WRF simulations. The practice of using at least two LES nests was recommended by J. Mirocha, Branko Kosovi´c, and Kirkil,2014.

By doing this the model was able to better predict the wind speed, which was verified using in-situ measurements. The paper also emphasized the importance of using a spin up time to let the model settle before using the output.

Chu et al.,2014used WRF with LES nested down to a resolution of 100 m to simulate several properties of a storm over complex terrain. The model was verified using measurements from both ground based and air based systems. The model was able to capture the turbulent boundary layer as well as the structure of the storm. However it was seen that

(18)

4 Chapter 1. Introduction

the model was too dissipating, something that has been consistently seen in other studies of microscale WRF-LES simulations (e.g. Kirkil et al.,2012).

Fonseca, Martín-Torres, and Andersson,2018have investigated wind forecasts at ESC. The paper focused on rocket applications and used WRF to forecast sudden wind shifts that can cause delays just before the launch of a sounding rocket. The paper both tests which PBL parameterization produces the best forecasts at the ESC as well as test the feasibility of using WRF in an operational context at ESC. This study does not employ LES and had a finest resolution grid size of 900 m, a grid size which do not completely capture the terrain around ESC. The study found that out of the tested PBL schemes, Asymmetric Convective Model 2(ACM2) produced the best forecasts. It also found that forecasts can be produced with a speed and accuracy that makes them feasible for operational use.

Another case of using WRF-LES in both ideal and real conditions can be found in Talbot, Bou-Zeid, and Smith,2012. The study first evaluates different LES setups in an ideal case.

It then used a six-domain WRF simulation with three LES-domains for a real case. The finest grid had a resolution of 50 m. The simulation was validated using both a sounding balloon and lidar measurements. The study concluded that the synoptic scale forecast used to force the WRF simulation has a significant impact on the result even at a very fine resolution grid. The authors suggested that a possible solution to eliminate some of the synoptic forecast bias could be to run WRF in a mesoscale simulation setup before initiating the LES.

(19)

5

2 Setup

In this thesis WRF of version 3.8.0 is used. A total of five domains is run, of those the innermost two conducts LES. This follows the recommendations of J. Mirocha, Branko Kosovi´c, and Kirkil,2014.

TABLE2.1: Domain setup used in the study.

Domain PBL Domain Size Grid Size Time Step Sound Steps per Time Step

01 ACM2 220x160 8100 m 27 s 4

02 ACM2 220x160 2700 m 9 s 4

03 ACM2 220x160 900 m 3 s 8

04 LES 220x160 300 m 1 s 16

05 LES 220x160 100 m 0.33 s 8

2.1 Physics

The physics used are similar to what Fonseca, Martín-Torres, and Andersson,2018used for simulation over the ESC area. The differences are mostly due to the different resolution. The microphysics scheme is the Goddard six-class (Tao, Simpson, and McCumber, 1989). It contains mixed-phase cloud processes which makes it suitable for high resolution simulations (Skamarock et al.,2008). The radiation scheme for both shortwave and longwave radiation is the Rapid Radiative Transfer Model for Global Circulation Models (J. Iacono et al.,2008) which is called every minute. Slope radiation and topology shad- ing are used for the finest three nests. The surface treatment is the unified Noah land surface model (Chen and Dudhia,2001). The PBL parameterization in the coarsest three nests is chosen to be ACM2 (Pleim,2007) after the result of Fonseca, Martín-Torres, and Andersson,2018. In the finest two nest no PBL parameterization is used. No cumulus parameterization scheme is employed in any of the nests.

Since the simulation in the two innermost domains are done using LES the dynamic for these domains are of special interest. The eddy viscosity is calculated from the Turbulent Kinetic Energy(TKE) using the equation that is described in Skamarock et al.,2008and follows

K = Cklp

e (2.1)

where e is the TKE, l is the length scale (∆x∆y∆z)^1/3and Ckis a constant that in this setup is kept at the default value of 0.15. It is important to not confuse the eddy viscosity with the molecular viscosity of air, the latter is smaller by several orders of magnitude. Instead, the eddy viscosity is a function of the flow and is not constant in time or space. The TKE eddy diffusivity approach also includes a way to model Sub Filter Stress(SFS) (Lilly,1967). SFS is a parameterization of the turbulent structures that are too small to be resolved (Kirkil et al.,2012). However, this model has a number of serious flaws. For example it assumes that production and dissipation of turbulent energy is balanced locally (J. D. Mirocha,

(20)

6 Chapter 2. Setup

Lundquist, and B. Kosovi´c,2010). This is obviously not true in a domain with any kind of heterogeneity.

For this reason a more advanced SFS model is used, the NBA(Nonlinear Backscatter Anisotropic) SFS model (Kosovic,1997). This SFS model is more complex and performs better than the linear models in an ideal test case (J. D. Mirocha, Lundquist, and B. Kosovi´c,2010).

2.2 Static Data

The static data used are mostly the default with a few exceptions. For the topography a dataset with a resolution of 50 m has been adapted for use with WRF. The data is sourced from Lantmäteriet (Lantmäteriet,2019) and consists of laser measurements from an air- craft with a mean error of about 1 m. The resolution can be compared to the resolution of WRF’s most higher resolution topography which is 30⁰⁰or about 900 m in the latitudinal axis. A comparison between WRF default topography and the added data can be seen in figure2.1. The difference between the two datasets can be as high as ±100 m. The largest differences can be found around the more complex terrain. As expected of a poor resolution dataset the default data smooths out the terrain and makes hills lower and valleys shallower. The high resolution topography is applied to the two finest nests. The land use model is primarily the IGPD-Modified MODIS (A. Friedl et al.,2010). In the innermost two domains the water category from MODIS is replaced by a dataset sourced from SMHI’s registry of lakes in Sweden (SMHI,2019). The dataset have been modified and adapted for use with WRF. This new dataset is far more accurate when it comes to the size and position of the lakes in the ESC area. A comparison over the innermost nest can be seen in figure2.2. The total lake area in the innermost domain is 1.715 km²in the MODIS dataset and 2.140 km²is the new lake category.

(21)

2.2. Static Data 7

(A) Best resolution WRF topography in the ESC area.

(B) Improved topography in the ESC area.

(C) Difference between default and improved topography in the ESC area.

FIGURE2.1: Comparison between the best resolution default dataset and the topography adapted for WRF.

(22)

8 Chapter 2. Setup

(A) MODIS lake category, plotted in the ESC area.

(B) New lake category implemented in WRF, plotted in the ESC area.

FIGURE2.2: Comparison between MODIS lake category and the lake category adapted for WRF. Yellow color represents land and blue color repre-

sents water.

(23)

2.3. Meteorological Data 9

2.3 Meteorological Data

Meteorological data is sourced from the European Center for Medium-Range Weather Forecasts(ECMWF) dataset named ECMWF ReAnalysis 5(ERA5) (Copernicus Climate Change Service Climate Data Store (CDS),2019). The data is a reanalysis, which means that measurement data have been incorporated in the simulation. Some amount of post-processing may also have been done. The goal of a reanalysis is to produce the most accurate representation of the atmosphere at that time. The use of reanalysis data instead of forecast data have some implications when it comes to the interpretation of the results.

Since the boundary conditions that comes from the ERA5 reanalysis can be considered as better in quality compared to a standard forecast, the overall quality of the WRF forecast is expected to rise slightly. However since the coarsest domain is so large in size and the forecast run time is so short the effect of the boundary conditions will not be that great.

The input data that has the most impact on the forecast is the data used to initialize the simulation at the first time step. This initial data will always consist of analysis data, even when producing real forecasts in an operational environment.

However the ERA5 data used for the comparison in the result section is also from the reanalysis, where in an operational environment it will consist of forecast data. This will lead to an increase in quality of that data as well. It is also reasonable to believe that the increase in quality between a reanalysis dataset and a forecast dataset is higher than the difference in quality between a WRF forecast forced with either reanalysis or forecast data.

This leads to the difference in quality between ERA5 reanalysis and WRF forecast might be underestimated, but to what margin is unknown.

The choice to use reanalysis data instead of forecast data is made for a few reasons.

1. The data is of higher quality, possibly making is easier to evaluate WRF.

2. The data is a single, consistent dataset over every point in time.

3. From a license standpoint, the data is easier to use and publish.

2.4 Measurement Data

There are a number of measurement sites around the ESC area. Those used in this thesis are detailed in table2.2and2.3. They will be used as a verification of the modeled results. The specifications for the instruments are taken from Fonseca, Martín-Torres, and Andersson,2018.

TABLE2.2: Specifications of the WT instrument.

Name: Wind Tower(WT)

Latitude: 67.893° N

Longitude: 21.106° E

Altitude: 10,25,45,65,85 and 100 m Above Ground Level(AGL)

Sensors: Wind

Resolution: 0.01 m s⁻¹, 1°

Range: 0-60 m s⁻¹, 0-359°

Accuracy: 2%, ±2°

(24)

10 Chapter 2. Setup

TABLE2.3: Specifications of the BPW instrument.

Name: Balloon Pad West(BPW)

Latitude: 67.890° N

Longitude: 21.080° E

Altitude: 3.5 m AGL

Sensors: Wind, temperature, pressure and humidity

2.5 Case Description

The study is conducted over 12 days, one in each month of the year. The dates are chosen to coincide with sounding balloon launches at ESC. Since these takes place quite irreg- ularly the study days are spread out over several years. To verify the simulation results comparisons are made to meteorological instruments placed on ESC.

TABLE2.4: A summary of the case days. All times are in UTC.

Simulation Started Simulation Ended Launch Window 2016-01-17 21:00 2016-01-19 09:00 06:00-09:00 2017-02-13 18:00 2017-02-15 06:00 03:00-06:00 2015-03-08 18:00 2015-03-10 06:00 03:00-06:00 2017-04-23 21:00 2017-04-25 09:00 06:00-09:00 2017-05-29 21:00 2017-05-31 09:00 06:00-09:00 2018-06-07 00:00 2018-06-08 12:00 09:00-12:00 2015-07-22 06:00 2015-07-23 18:00 15:00-18:00 2017-08-10 15:00 2017-08-12 03:00 00:00-03:00 2016-08-30 21:00 2016-09-01 09:00 06:00-09:00 2018-10-15 18:00 2018-10-17 06:00 03:00-06:00 2017-10-31 18:00 2017-11-02 06:00 03:00-06:00 2017-12-06 00:00 2017-12-07 12:00 09:00-12:00

2.5.1 Balloon Launch Criterion

As mentioned in section1.1a launch of a large stratospheric balloon can only take place under specific meteorological conditions. These can be complex and significantly varied depending on payload, balloon type and other factors. For the purpose of this study a simplified balloon launch criterion shall be defined. Due the the difference in models and measurement positions, completely equivalent parameters cannot be obtained. However they should still remain comparable.

The last three hours of a simulation is termed the "Launch Window". It is in this time the hypothetical balloon launch shall take place. For each hour all the three different criterias are evaluated. They are the low altitude wind speed, high altitude wind speed and the maximum gust.

The low altitude wind is defined as the average horizontal component of the wind speed at an altitude of 10 m AGL. This shall be below a 1-hour average of 4 m s⁻¹ for the entire launch window. The variable in WRF is the two horizontal wind components, named

"U10" and "V10". ERA5 also have the same variable, so for this the two models can be

(25)

2.6. Statistical Analysis 11

compared completely equivalently. The reference measurement is the 10 m level of the wind tower at ESC.

The high altitude wind is defined as the average horizontal component of the wind speed at an altitude of 100 m AGL. This shall be below a 1-hour average of 6 m s⁻¹for the entire launch window. The variable in ERA5 data is two horizontal wind components named

"U100" and "V100". WRF does not have this variable so instead the variables "U" and "V"

are used at the 9th model level (about 102 m AGL). This is deemed to be comparable. The reference measurement is the 100 m level of the wind tower at ESC.

The maximum gust is defined as the maximum horizontal component of the wind speed at an altitude of 10 m AGL. This shall be below 6 m s⁻¹for the entire launch window. The variable in ERA5 data is called "Instantaneous 10m wind gust" and is defined as the maximum 3-second average wind at 10 m AGL. In WRF the simple maximum value of the 10- meter horizontal wind speed is taken. The reason for this is that in the treatment of the wind field, features with a spatial wavelength of less than twice the grid size are implicitly filtered out. This makes the simple maximum approximately equal to the 3-second average. The primary reference measurement is the 10 m level of the WT where the maximum gust is taken as the maximum 3-second average wind.

2.6 Statistical Analysis

A short statistical analysis will also be performed on the simulation results. This is to have an evaluation metric that is independent from the operationally important, but somewhat arbitrary, balloon launch criterion. The analysis will be focused on the scalar wind speed and will be done according the methods proposed by Koh, Wang, and Bhatt,2012.

The input to the analysis is the wind speed at the 10-meter level. A sample average over 20 seconds is taken every minute from the 10-meter level of the WT. This averaging is done to smooth out any variations that is not resolvable in WRF to produce two somewhat comparable datasets. This is then compared to samples taken every minute from the WRF forecast. The length of each sampling vector is 60 samples, or an hour in time.

This approach gives a total of 2160 samples for the analysis. The error will be divided as a mean error and a pattern error. This report will only contain a brief description of the error parameters, for a full derivation the reader is referred to Koh, Wang, and Bhatt,2012.

2.6.1 Pattern Error

The pattern error is the error that remains after the mean error is removed. It can be divided into two main parts: the correlation and the variance similarity.

The correlation is denotedρ and can be viewed as a measure of the phase agreement between the observation and the forecast. The optimum score is 1 and indicates that the forecast is completely in phase with the observations. With the vector of forecast samples being F and the vector of the observational samples being O the equation describing the correlation is defined as

ρ = 1

σOσF〈(F − 〈F〉)· (O − 〈O〉)〉 , −1 ≤ ρ ≤ 1 (2.2) As is common,σOandσF denotes sample standard deviation of the corresponding vec- tors and the angular brackets denotes a mean operation.

(26)

12 Chapter 2. Setup

The variance similarity is denotedη and can be viewed as a measure of the amplitude agreement between the observation and the forecast. The optimum score is 1 and indicates that the forecast completely predicts the observed variability. The variance similarity can be described as

η = 2σOσF

σ²_O+ σ²_F , 0 ≤ η ≤ 1 (2.3)

Finally the normalized error variance can be defined. It is denotedα and has an optimum score of 0 which indicates perfect correlation and variance similarity. The normalized error variance is defined as

α = 1 − ηρ , 0 ≤ α ≤ 2 (2.4)

2.6.2 Mean Error

The mean error indicates if there is any systematic bias in the simulation. The main com- ponent is called normalized bias and is denotedµ. It is defined with the equation

µ =〈D〉

σD

(2.5)

where D is the discrepancy vector D = F − O. Together with the normalized error variance α the common metric of normalized root mean squared error(RMSE), here denoted with δ can be defined as

δ =q

α(1 + |µ|²) (2.6)

(27)

13

3 Results

3.1 Launch Analysis

In this section the results relating to the balloon launch criterion are presented for the different launch days. Only data from the launch window, i.e. the last three hours of a simulation are presented. For each day the table shows the predictions made both by the input data ERA5 and the WRF model together with the measured data. A green value indicates that the balloon launch criterion is met and a red value indicates a violation. An optimal result would be that the correct prediction is made for every criteria at every point in time. This includes to correctly predict a violation of the balloon launch criteria when that is the case.

The 10-meter wind speed as a function of time is also included to serve as a context for the formal results.

3.1.1 2016-01-19

FIGURE 3.1: Comparison of actual and predicted 10-meter wind speed during the launch window on 2016-01-19

(28)

14 Chapter 3. Results

TABLE 3.1: Evaluation of balloon launch criterion from the different sources on 2016-01-19. A green color indicates a pass and a red color in-

dicates a fail for the given hour.

Time 09-10 10-11 11-12

Mean Speed 10 meter WT 1.6 1.2 1.6 Mean Speed 10 meter WRF 0.4 0.4 1.0 Mean Speed 10 meter ERA5 0.9 1.5 1.1 Mean Speed 100 meter WT 4.4 3.4 3.9 Mean Speed 100 meter WRF 0.3 0.3 0.6 Mean Speed 100 meter ERA5 0.5 1.5 0.9 Max Gust 10 meter WT 3.2 2.9 3.5 Max Gust 10 meter WRF 1.0 0.6 1.9 Max Gust 10 meter ERA5 1.5 1.3 1.3

3.1.2 2017-02-15

FIGURE 3.2: Comparison of actual and predicted 10-meter wind speed during the launch window on 2017-02-15.

(29)

3.1. Launch Analysis 15

Time 03-04 04-05 05-06

Mean Speed 10 meter WT 0.8 0.5 0.4

Mean Speed 10 meter WRF 1 1 1.1

Mean Speed 10 meter ERA5 3.0 2.9 3.0 Mean Speed 100 meter WT 1.3 1.8 2.1 Mean Speed 100 meter WRF 0.8 1.0 1.6 Mean Speed 100 meter ERA5 5.8 5.6 5.9 Max Gust 10 meter WT 2.0 1.4 1.4 Max Gust 10 meter WRF 1.7 1.5 1.6 Max Gust 10 meter ERA5 5.0 4.9 5.1

3.1.3 2015-03-10

(30)

Time 03-04 04-05 05-06

3.1.4 2017-04-25

(31)

Time 06-07 07-08 08-09

Mean Speed 10 meter WT 0.9 0.9 1.5 Mean Speed 10 meter WRF 0.5 0.7 0.8 Mean Speed 10 meter ERA5 2.1 2.1 2.2 Mean Speed 100 meter WT 0.8 1.3 1.6 Mean Speed 100 meter WRF 0.4 0.9 1.2

Mean Speed 100 meter ERA5 3 3 3.1

Max Gust 10 meter WT 2.9 3.2 4.1 Max Gust 10 meter WRF 0.9 1.1 1.2 Max Gust 10 meter ERA5 4.2 4.9 5.3

3.1.5 2017-05-31

(32)

Time 06-07 07-08 08-09

3.1.6 2018-06-08

(33)

Time 09-10 10-11 11-12

3.1.7 2015-07-23

(34)

Time 15-16 16-17 17-18

3.1.8 2017-08-12

(35)

Time 00-01 01-02 02-03

3.1.9 2016-09-01

(36)

Time 06-07 07-08 08-09

3.1.10 2018-10-17

(37)

Time 03-04 04-05 05-06

3.1.11 2017-11-02

(38)

Time 03-04 04-05 05-06

3.1.12 2017-12-07

(39)

Time 09-10 10-11 11-12

3.1.13 Summary

The summarized results is presented in table3.13. The maximum score in each category is 36, making a total maximum score of 108. This means that in this trial the WRF model made the correct call 93% of the time while ERA5 made the correct call 69% of the time.

TABLE3.13: Number of correct calls compared to measurements for the different balloon launch criteria.

Criterion WRF ERA5

Mean Speed 10 meter 33 29 Mean Speed 100 meter 33 22 Max Gust 10 meter 34 23

Total 100 74

3.2 Statistical Analysis

In figures3.13and3.14the result of the statistical analysis can be seen. As described in section2.6the analysis is done for each hour of every launch window simulated. The result for the first hour in each launch window is marked by a square, the second hour is marked by a diamond and the third hour is marked by a star. The color key for the symbols can be found in table3.14.

(40)

FIGURE3.13: The angle represents variance similarity and the radial distance represents absolute value of the correlation. The optimal model per-

formance is marked by the black star.

(41)

FIGURE 3.14: Plot representing the bias and error variance connection to the normalized root mean square error. The contours represent normalized RMSE in increments of 0.2. The optimal model performance is

marked with a black star.

TABLE3.14: Color key for figures3.13and3.14

Date Color

2016-01-19 2017-02-15 2015-03-10 2017-04-25 2017-05-31 2018-06-08 2015-07-23 2017-08-12 2016-09-01 2018-10-17 2017-11-02 2017-12-07

(42)

(43)

29

4 Discussion

This section will provide some discussion of the overall results of the study. It will also provide comments and additional information on some of the more interesting simulation results. An effort is made to try and explain the model behavior with a focus on identi- fying the strengths and weaknesses of using microscale LES to predict winds in a specific location.

4.1 General Discussion of Results

The results in section3.1gives positive indications of the effect of using WRF to improve forecasts. The prediction of all three launch criteria were improved compared to ERA5.

The WRF forecast had a very similar success rate on all three criteria, while ERA5 had a much larger success rate in forecasting the 10-meter wind. This indicates that the improvement gained using WRF is lower for 10-meter wind than for 100-meter wind. This could be seen as a surprising result as the ERA5 reanalysis is expected to be better at a higher altitude above ground.

Looking at the 10-meter wind as a function of time in figures3.1to3.12it can be seen that the forecast often produces far less turbulence than observed. There also seem to be more turbulence generated during the simulation between May-July. The lack of turbulence is most likely at least in part caused by the relatively large grid size of 100 m. As the assump- tion is made that the energy containing eddies are resolved a large grid size can lead to a lack of turbulence. This was known already in the beginning of the study, however the resolution could not be increased due to a lack of computational power. Doing another nest down to a grid size of 33 m produced a large increase in resolved TKE using NBA in an ideal study case by J. Mirocha, Branko Kosovi´c, and Kirkil,2014. This makes it likely that a smaller grid size would increase the resolved turbulence also in this case, likely producing a more accurate solution.

The statistical analysis in section3.2gives further indications of the forecast skill. In figure 3.13it can be seen the forecast often have very small or even negative correlation to the observations. This is to a large degree expected results. Since the statistical analysis is done using 1-minute samples the resulting signal frequencies are high, as it should be in a microscale forecast. If the analysis was done on a mesoscale forecast a sampling time of 10 or 30 minutes might have been used. With a longer sampling time the more low frequency shifts in the atmosphere would have been analyzed and a higher correlation would have been expected.

Instead the variance similarity is expected to be close to one for an accurate forecast. This will indicate that the forecast produced a correct amount of variation within the hour, variation that often can be attributed to turbulence. The variance similarity for many forecast points are very close to one, showing that the forecast is somewhat successful in simulation variance.

(44)

30 Chapter 4. Discussion

When looking at3.14it can be seen that the forecast seems to have a positive bias at the 10 meter level. This means that the forecast often overestimates the wind speed compared to the observations. There can be multiple causes for this including too coarse grid size, incorrect surface parameterization or an incorrect input from the mesoscale domains.

It can also be seen that in most cases the bias in the primary contributor to the normalized root mean square error. The error variance is so high due to low correlation that any reasonable change would only have a marginal effect on the normalized RMSE. Instead the normalized bias can be said to influence most of the error.

A comparison between ERA5 and WRF for the statistical skill scores can not be made due to the low temporal resolution of ERA5.

No clear correlation between forecast hour and forecast skill can be observed. This is in- dicative of that the spin up period is sufficient.

4.2 2017-05-31 - Large overestimation of low altitude wind speed

The only occasion where the WRF forecast failed to make the correct prediction for all three hours was on 2017-05-31. As seen in figure3.5the turbulence was well developed but the forecasted wind speed was consistently higher than the measured. However the mean 100-meter wind speed and the maximum gust were well predicted. Data from a sounding balloon launched 07:43 UTC are available to help investigate how the simulation might have went wrong.

FIGURE4.1: Comparison of actual and predicted temperature in the Tro- posphere on 2017-05-31. Sounding balloon released 07:43 UTC and WRF

forecast with target time 07:43 UTC.

(45)

4.2. 2017-05-31 - Large overestimation of low altitude wind speed 31

In figure4.1the forecasted temperature is compared with the measured. The important thing that stands out is that the small inversion that is measured at approximately 2000 m altitude is not present in the forecast. The forecast does not contain any temperature inversions at all in the troposphere, which is a strange behavior. Another indication on what might be the problem can be seen in figure4.2. There the water vapor content of the lower troposphere can be seen. It can be seen that a portion of the atmosphere contains a rather high moisture gradient.

This can lead to problems first investigated by Yamaguchi and Feingold,2012. It is related to how WRF treats potential temperature, which determines if air is rising or falling. The problem occurs since the potential temperature used in the small advection time steps assumes that there is no moisture perturbation. In situations with a high moisture gradient this can lead to artificially high mixing, especially on the top of clouds. It is also at the cloud tops that a temperature inversion is usually found.

FIGURE4.2: Comparison of actual and predicted water vapor content on 2017-05-31. Sounding balloon released 07:43 UTC and WRF forecast with

target time 07:43 UTC.

However, a better treatment of potential temperature in the presence of sharp moisture gradients require a rewrite of the ARW core. In the version of WRF used for this study (3.8.0) this is not the default treatment. Further study is suggested with an updated version of the solver. It is also important to ensure that the simulation is using the proper vertical resolution needed to resolve these small inversions.

(46)

4.3 2016-09-01 - Time shift in forecast

When looking at table3.9it can be seen that in the first hour the 100-meter wind is significantly underestimated. The comparison between actual and predicted wind for a long period of 12 hours can be seen in figure4.3.

FIGURE 4.3: Comparison between actual and predicted 100-meter wind speed before, during and after the launch window on 2016-09-01.

Here it can be seen that a very clear decrease in measured wind speed occurs approximately between the times 03:30 to 06:00 UTC. A similar decrease in wind speed can be seen in the WRF forecast approximately between 05:30 UTC to 07:30 UTC while it can not be observed at all in the ERA5 reanalysis.

This observation makes it easy to believe that the WRF forecast is time shifted. One reason for this could be that the large dip in wind speed happens during the first three hours of the LES simulation. This time is considered as a spin up time and the result is not expected to be accurate.

To investigate the difference between different spin up times the simulation is done again.

The setup is identical with the exception the the domains four and five is initiated 24 hours into the simulation instead of 30 hours. The 100 m wind obtained from the simulation can be seen in figure4.4.

(47)

4.4. 2017-12-07 - Winter case with strong surface inversion 33

FIGURE4.4: As figure4.3but with a simulation using a longer spin up time added.

The figure shows that in this case the difference between the regular simulation and the simulation with the longer spin up time is small. The time shift previously observed is still present. The most likely cause for the time shift is the input from the mesoscale domains.

An investigation into why that happened is however outside the scope of this discussion.

4.4 2017-12-07 - Winter case with strong surface inversion

A temperature inversion is when the temperature in the atmosphere increase with increasing altitude. They form due to an imbalance in radiation as well as advection of cold air (Zhang et al.,2011). They are common in the Arctic area. One of the hypothesis in the beginning of the study was that microscale forecasts would significantly increase the ability to simulate the surface based inversions due to a higher resolution of the terrain allowing pools of cold air to form. The explicit vertical advection used in LES could possible also be a factor contributing to a better representation of inversions. Since a temperature inversion works as an insulation from the free atmosphere above it will have an effect on the wind inside the inversion (Bradley, T. Keimig, and Diaz,1992).

During the launch window of the simulation conducted 2017-12-07 a sounding balloon was released at ESC. Using that measurement an accurate temperature profile can be con- structed for lower atmosphere. In4.5this measurement can be seen together with simulation outputs from both WRF and ERA5.

(48)

FIGURE4.5: Comparison of actual and predicted temperature of the lower atmosphere on 2017-12-07. Sounding balloon released 09:57 UTC, WRF forecast with target time 09:57 UTC and ERA5 Reanalysis with target time

10:00 UTC.

It can clearly be seen that in this case WRF treats the inversion in a much better way than ERA5. Further indication of the effect that the microscale forecast have on wind can be seen in figure4.6. There the temperature profile predicted by WRF is plotted every half hour from the initialization or the finest domain. There it can be seen that simulating with LES results in lower temperatures in the inversion. The problem here is that the sounding is only done at one point in time, so it is hard to judge the development of the forecast with that as an only source.

(49)

FIGURE4.6: A buildup of the strong inversion after the initialization of the highest resolution domain on 2017-12-07. Profiles taken every 30 minutes

from start to end of domain.

For that reason a comparison is also made with the time series temperature measurement of the BPW instrument described in table2.3. The comparison can be seen in figure4.7.

Here it can clearly be seen that the input value fed into the domain at start-up is overes- timating the temperature. It can also be seen that the model starts correcting the overes- timated temperature and seems to reach an equilibrium after 250-300 simulated minutes after a short undershoot. Figure4.7definitely illustrates the need for a spin up period when initializing a new domain and also raises questions if the spin up period of three hours employed in the study is too short in cases with a strong inversion.

(50)

FIGURE4.7: Measured temperature at BPW compared with the finest resolution forecast from domain initialization on 2017-12-07.

The question remains on how this affects the wind in the simulation. In figure3.12it can be seen that the 10-meter wind is more correct towards the end of the simulation where the ground temperature is also more correct. However this may just be a coincidence. In 4.8a more detailed look at the forecasted wind inside the surface inversion can be obtained. The lines represent the average wind speed profile obtained by both WRF and ERA5 at the start and end of the microscale forecast. It can be noted that the wind shear near the surface of the WRF simulation increased when the inversion was stronger. This is in line with Monin-Obukhov similarity theory (Monin and Obukhov,1954) which predicts that the wind shear should increase with increasing stability. No similar increase in wind shear can be seen the the ERA5 reanalysis.

(51)

FIGURE 4.8: Average wind in the temperature inversion during the start and end of the microscale forecast on 2017-12-07.

While this study does not focus on the treatment of surface based temperature inversion it is still an important topic for all weather simulations conducted in the Arctic area. An improved inversion model might improve the wind forecast at ESC but it is a subject for further investigations.

(52)

(53)

39

5 Conclusions

5.1 Summary

The goal of this thesis was to develop a method to use the Weather Research and Fore- casting model to improve the quality of wind forecasts at Esrange Space Center. More specifically the study evaluates the prediction of favorable or unfavorable conditions to launch a stratospheric balloon.

The WRF setup used includes an improved terrain model, using adapted data from Lant- mäteriet. It also included an improved lake model, using adapted data from SMHI.

The simulations were run using a simultaneous nesting approach. The three coarsest nests where run using a planetary boundary layer parameterization while the finest two were run using a large eddy simulation approach.

The input data used to force the model was the ERA5 reanalysis. The resulting forecast was also compared to the ERA5 dataset using measurement data from Esrange Space Center as verification.

Using this method one test case were run during each of the 12 months in a year. Each test case was evaluated for a "launch window" of three hours where three metrics were evaluated each hour, the 10-meter wind, the 100-meter wind and the 10-meter gust.

The results were that the percentage of correct forecasts rose from 69% for the input ERA5 data to 93% from the resulting WRF forecast. An increase in skill compared to ERA5 was seen in all three categories, but least significant for 10-meter wind. Some turbulence could be produced by the LES approach with in some cases generated realistic gusts. The statistical analysis shows that the variance in time was able to be better simulated than correlation, and that the model generally produced a stronger 10-meter wind speed than observed. Overall the results are viewed as a significant improvement of the forecasting ability.

5.2 Operational viability

The simulations were done on 80 cores of AMD Opteron 6238 using Intel MPI and no OpenMP. Computing was done on the HPC2N Abisko cluster. Using this setup a simulation took about 20 hours, making it too slow for operational use. However the Abisko cluster is an old computing facility and it is estimated that using a similar number of cores on a modern cluster would reduce computing time significantly, making it operationally viable.

(54)

40 Chapter 5. Conclusions

5.3 Possible improvements and further research

During the course of the project some potential improvements to the model were found that were not investigated fully due to time and computing restrictions.

The first and most obvious improvement is running the model at a higher resolution. As explained in section4.1there are strong indications that a notable increase in model performance could be achieved using an additional nest.

An increase in vertical levels would probably also be beneficial for the model. 60 vertical levels were used in this study, and more research is required to find the optimal amount.

A switch to using a hybrid sigma-pressure vertical coordinate system could also be investigated.

As mentioned in section4.2the current treatment of potential temperature is not accurate. Any further study or operational use should use a updated version of WRF which includes the effect of moisture perturbations on potential temperature.

Finally the forecasting ability of other meteorological data affecting balloon launches could be investigated, for example wind direction, precipitation or temperature. The model could also be extended to cover all phases of the balloon flight. In that case an investigation into the forecasting ability of high altitude winds using LES is very interesting.

(55)

41

Bibliography

A. Friedl, Mark et al. (2010). “MODIS Collection 5 Global Land Cover: Algorithm Refine- ments and Characterization of new Datasets”. In: Remote Sensing of Environment 114, pp. 168–182.DOI:10.1016/j.rse.2009.08.016.

American Meteorological Society, ed. (2019). American Meteorological Society.URL:http:

//glossary.ametsoc.org/wiki/.

Bradley, Raymond, Frank T. Keimig, and Henry Diaz (1992). “Climatology of Surface-Based Inversions in the North American Arctic”. In: Journal of Geophysical Research 97, pp. 15699–.

DOI:10.1029/92JD01451.

Chen, Fei and Jimy Dudhia (2001). “Coupling an Advanced Land Surface–Hydrology Model with the Penn State–NCAR MM5 Modeling System. Part I: Model Implementation and Sensitivity”. In: Monthly Weather Review 129.4, pp. 569–585. DOI: 10 . 1175 / 1520 - 0493(2001)129<0569:CAALSH>2.0.CO;2.

Chu, Xia et al. (2014). “A Case Study of Radar Observations and WRF LES Simulations of the Impact of Ground-Based Glaciogenic Seeding on Orographic Clouds and Precipi- tation. Part I: Observations and Model Validations”. In: Journal of Applied Meteorology and Climatology 53.10, pp. 2264–2286.DOI:10.1175/JAMC-D-14-0017.1.

Copernicus Climate Change Service Climate Data Store (CDS) (2019). Copernicus Climate Change Service (C3S) (2017): ERA5: Fifth generation of ECMWF atmospheric reanalyses of the global climate. ECMWF.URL:https://cds.climate.copernicus.eu/cdsapp#

!/home.

Fonseca, Ricardo, Javier Martín-Torres, and Kent Andersson (2018). “Wind Forecasts for Rocket and Balloon Launches at the Esrange Space Center Using the WRF Model”. In:

Weather and Forecasting 33.3, pp. 813–833.DOI:10.1175/waf-d-18-0031.1.

J. Iacono, Michael et al. (2008). “Radiative Forcing by Long-Lived Greenhouse Gases: Cal- culations with the AER Radiative Transfer Models”. In: Journal of Geophysical Research 113.DOI:10.1029/2008JD009944.

Kirkil, Gokhan et al. (2012). “Implementation and Evaluation of Dynamic Subfilter-Scale Stress Models for Large-Eddy Simulation Using WRF”. In: Monthly Weather Review 140.1, pp. 266–284.DOI:10.1175/MWR-D-11-00037.1.

Koh, T.-Y, S Wang, and Bhuwan Bhatt (2012). “A diagnostic suite to assess NWP perfor- mance”. In: Journal of Geophysical Research (Atmospheres) 117, pp. 13109–.DOI:10 . 1029/2011JD017103.

Kosovic, Branko (1997). “Subgrid-scale modelling for the large-eddy simulation of high- Reynolds-number boundary layers”. In: Journal of Fluid Mechanics 336, pp. 151–182.

DOI:10.1017/S0022112096004697.

Lantmäteriet (2019). Lantmäteriet.URL:https://www.lantmateriet.se/sv/Kartor- och-geografisk-information/Hojddata/GSD-Hojddata-grid-50-/.

Lilly, D. K. (1967). “The Represetation of Small-Scale Turbulence in Numerical Simula- tion Experiments”. In: Proc. IBM Scientific Computing Symp. on Environmental Sciences.

DOI:10.5065/D62R3PMM.

Liu, Yubao et al. (2011). “Simultaneous nested modeling from the synoptic scale to the LES scale for wind energy applications”. In: Journal of Wind Engineering and Industrial

Wind Forecasts Using Large Eddy Simulations for Stratospheric Balloon Applications