VISITOR ARRIVAL FORECASTS AMID COVID-19: A perspective from the Asia and Pacific team

Keywords:

It is important to provide scientific assessments concerning the future of tourism under the uncertainty surrounding COVID-19. To this purpose, this paper presents a two-stage threescenario forecast framework for inbound-tourism demand across 20 countries. The main findings are as follows: in the first-stage ex-post forecasts, the stacking-based algorithms are more accurate and robust, especially when combining five single models. The second-stage ex-ante forecasts are based on three recovery scenarios: a mild case assuming a V-shape recovery, a medium one with a V/U-shape, and a severe one with an L-shape. The forecast results show a wide range of recovery (10%-70%) in 2021 compared to 2019. This two-stage three-scenario framework contributes to improvement in the accuracy and robustness of tourism demand forecasting.

Keyword: COVID-19,

INTRODUCTION

Covid-19 has had devastating effects on many facets of society and the economy of every country in the world. The impact has been particularly felt in the tourism and hospitality sector. To reduce the spread of the virus, nations closed their borders and limit the mobility of its residents. This resulted in a drastic economic and social crisis. Hotels and restaurants closed, airlines cancelled flights and grounded planes, travel agencies and tour operators ceased operation, and tourist attractions shut their doors. This all happened within a few weeks, if not days when countries realized the gravity of the pandemic. The downturn in the global economy has been far more reaching and far deeper than other shocks in the recent past, such as the September 11 terrorist attacks of 2001, the SARS outbreak in 2003, and the global financial crisis of 2008/9. Several studies have sort to estimate the economic impact of COVID-19 (see [START_REF] Farzanegan | International Tourism and Outbreak of Coronavirus (COVID-19): A Cross-Country Analysis[END_REF][START_REF] Mariolis | The COVID-19 multiplier effects of tourism on the Greek economy[END_REF][START_REF] Yang | Coronavirus pandemic and tourism: Dynamic stochastic general equilibrium modeling of infectious disease outbreak[END_REF] yet few studies have sort to model a geographically comprehensive, methodologically rigorous post-COVID-19 recovery. The need for this research is important, given that tourism destinations and tourism and hospitality businesses need to plan when and how to re-open. This has significant implications for the wellbeing of residents as businesses decide to reemploy staff. At the time of writing, the pandemic is still ongoing. Some countries have opened again, some remain fully closed and others lie somewhere in between. How and when the tourist and hospitality market will recover remains unknown. However, accurate forecasting of the impact of COVID-19 on the tourism industry and its recovery is critical for strategic planning of tourist destinations and tourism-related businesses.

The evidence of the importance of forecasting tourism demand can be seen by the way tourism demand forecasting has received considerable attention in the literature. Numerous review papers have been published summarizing research published to date. These reviews include [START_REF] Crouch | The study of international tourism demand: A survey of practice[END_REF], [START_REF] Witt | Forecasting tourism demand: A review of empirical research[END_REF], [START_REF] Lim | Review of international tourism demand models[END_REF][START_REF] Lim | A meta-analytic review of international tourism demand[END_REF], [START_REF] Li | Recent developments in econometric modeling and forecasting[END_REF], [START_REF] Song | Tourism demand modelling and forecasting-A review of recent research[END_REF], and more recently [START_REF] Wu | New developments in tourism and hotel demand modeling and forecasting[END_REF] and [START_REF] Song | A review of research on tourism demand forecasting: Launching the Annals of Tourism Research Curated Collection on tourism demand forecasting[END_REF]. As noted by Song & Li (2008, p. 217) "It is

To achieve the research objectives and as will be elaborated on below, we undertake a tourism forecasting competition, along the lines of [START_REF] Athanasopoulos | The tourism forecasting competition[END_REF]. This 'competition' is a contest among various methods for the 'best' forecast of post-COVID-19 tourism demand. There are two stages to the forecasting. The first stage involves ex-post forecasting of international visitor arrivals before COVID-19, to identify the most accurate forecasting method(s) during this period. In the second stage, ex-ante judgemental-adjusted scenario forecasting of visitor arrivals during and after COVID-19 are made up until the end of 2021. The purpose of this two-stage forecasting procedure is to identify effective forecasting framework and evaluation procedures in a crisis.

There are numerous measures to assess the accuracy of forecasts. [START_REF] Hyndman | Another look at measures of forecast accuracy[END_REF] advocate for using the Mean Absolute Scaled Error (MASE) above others. They argue that this measure is less sensitive to outliers and it is independent of the scale of the data. An additional benefit is the ease of interpretation: if MASE is greater than one, it is a poorer forecast than the average one-step naïve forecast computed in-sample. Conversely, if MASE is less than one it signifies an improved forecast than the average one-step naïve forecast computed in-sample.

The remainder of the paper is structured as follows: the next section describes the modelling strategies of both stages while the subsequent section presents and discusses the results. The final section concludes the paper, and highlights the methodological and theoretical contributions of the research and the policy implications of the work.

MODELLING STRATEGIES

The data used in the present competition is collected from 20 tourist destinations across the world, which covers all UNWTO regions. They are Australia, Bulgaria, Canada, Chile, Czech, Finland, Indonesia, Japan, South Korea, Malaysia, Mauritius, Mexico, New Zealand, Singapore, South Africa, Sweden, Thailand, Tunisia, UK, and USA. For each destination, the total volume of international visitor arrivals and visitor arrivals from five key source markets are collected as measures of tourism demand [START_REF] Song | A review of research on tourism demand forecasting: Launching the Annals of Tourism Research Curated Collection on tourism demand forecasting[END_REF]. A total of 120 time-series are adopted for the forecasting practice. Quarterly data is used to the end of 2019 (2019Q4), representing the most frequently modelled time interval [START_REF] Song | Tourism demand modelling and forecasting-A review of recent research[END_REF]. Quarterly data has the advantage of being long enough to help policymakers assess trends and not overreact to random fluctuations but short enough to help policymakers see the effects of their decisions. The starting points of the tourism demand series vary from 1991Q1 to 2010Q1, resulting in the longest series having observations and the shortest series with 40 observations. To incorporate the influence of external factors, economic variables, such as GDP, CPI, and exchange rates of relevant countries/regions (markets hereafter), are collected from the Global Economic Monitor (GEM) database of the World Bank and the International Financial Statistics (IFS) database of IMF. To facilitate forecast adjustments of the COVID-19 pandemic, the cumulative number of cases and deaths of each relevant markets are collected from Coronavirus Disease COVID-19 Dashboard (WHO, 2020), and the travel ban (border control) information is collected from the Coronavirus Government Response Tracker (University of Oxford, 2020) and the Policy Responses to the Coronavirus Pandemic (Our World in Data, 2020a).

The main framework of the forecasts is a two-stage-process. In the first stage, ex-post forecasting of international visitor arrivals before COVID-19 is examined and evaluated. We use data up until the end of 2019, with the data until the end of 2019 as training and data of 2019 as testing. single models and 26 stacking models are examined and compared. In the second stage, judgmental adjusted forecasting under three recovery scenarios and one set of baseline forecasting which assumes there is no COVID-19 are generated separately to estimate the impact of the COVID-19 on tourism demand. Figure 1 illustrates this framework.

[Figure 1 about here] First Stage -ex post forecasting before COVID-19 Individual forecasting models The difficulty of the first stage forecast is in trying to find one algorithm that not only provides accurate forecasts of the tourism demand series but also in finding an algorithm that performs stably across various origin-destination pairs and different time horizons. Ensemble methods have proved effective in improving the forecast accuracy of constituent methods [START_REF] Jaganathan | A combination-based forecasting method for the M4competition[END_REF]. In particular, bagging, boosting, and stacking techniques have emerged to provide outstanding forecasts in tourism demand analysis (e.g. [START_REF] Athanasopoulos | Bagging in tourism demand modeling and forecasting[END_REF][START_REF] Cankurt | Tourism demand forecasting using ensembles of regression trees[END_REF]. For the present competition, stacking-based algorithms are examined and compared, and the best algorithm is proposed for forecasting the 120 time-series of international visitor arrivals. 11 preliminary single models are adopted to generate forecasts of the tourism demand series, namely seasonal Naïve (SN), seasonal autoregressive integrated moving average (SARIMA), exponential smoothing (ETS), seasonal and trend decomposition using Loess (STL), exponential smoothing state space model with Box-Cox transformation, ARMA errors, trend, and seasonal components (TBATS), autoregressive distributed lag (ARDL), static regression with time-varying parameter (SR-TVP), univariate multilayer perceptron (MLP), multivariate multilayer perceptron (MLPX), univariate extreme learning machine (ELM), and multivariate extreme learning machine (ELMX). These models cover the three major families of quantitative forecasting approaches commonly adopted in tourism demand forecasts: time series methods; econometric models; and AI-based techniques [START_REF] Song | A review of research on tourism demand forecasting: Launching the Annals of Tourism Research Curated Collection on tourism demand forecasting[END_REF]. These models include both univariate and multivariate methods, which ensure both historical patterns of tourism demand series and the influence of economic factors are covered.

Both the tourism demand series and the economic factor series are pre-processed before the model estimation. Missing values in the middle of the time series are imputed using Kalman Smoothing on structural time series models [START_REF] Harvey | Forecasting, Structural Time Series Models and the Kalman Filter[END_REF] with the "na_kalman" function of "imputeTS" package in R [START_REF] Moritz | imputeTS: Time series missing value imputation in R[END_REF]. Missing values at the end of the time series (i.e. series end earlier than 2019Q4) are extended by exponential smoothing using the "ets" function of "forecast" package in R [START_REF] Hyndman | Automatic time series forecasting: The forecast package for R[END_REF][START_REF] Hyndman | forecast: Forecasting functions for time series and linear models[END_REF]. The imputation of missing values provides a complete time series for model estimation while preserving the characteristics of the original time series. In statistical analysis, extreme values and outliers may have a significant influence on model estimation. In the tourism demand context, these extreme values and outliers are usually caused by one-off events or sudden fluctuations in the model variables. These disturbances are usually treated by incorporating dummy variables into the statistical model, but their treatment can be subjective and inflexible. Furthermore, in a large-scale forecast project like the present competition, it is extremely costly to implement such treatments. In contrast, we adopt an outlier smoothing process. The extreme values in both tourism demand series and economic variables are firstly identified by "tsoutlier" function of "forecast" package in R. The identified outliers remain in the time series in multivariate model estimations if the fluctuations of tourism demand series and economic variables happen simultaneously. In the single variate model estimations, or in the case that identified outliers in the tourism demand series and the economic variables do not occur simultaneously, extreme values are replaced by STL decomposition smoother. This process removes extreme values in the data while preserving the most out of original information.

The estimations of the 11 preliminary single models are performed in R. The time series models (SN, SARIMA, ETS, STL, and TBATS) are estimated by functions of "forecast" package ("snaive", "auto.arima", "ets", "stlm", and "tbats", respectively). All the time series methods, except for SN, are optimized within the family. For example, SARIMA models with different AR lags and MA lags are considered and compared using the Akaike information criterion (AIC). Estimations of economic models are conducted on tourism demand series whenever the corresponding economic variables available. The ARDL model is estimated using "auto_ardl" function of "ARDL" package in R [START_REF] Natsiopoulos | ARDL: ARDL, ECM and Bounds-Test for cointegration[END_REF]. The lag of each variable in the ARDL model is derived automatically through an AIC comparison with the maximum of lags set to eight (two years in the current context). The optimal ARDL models are also selected so that the correct signs correspond with the coefficients of economic variables (i.e. positive for income (GDP) and negative for prices (exchange rate adjusted CPIs)). The SR-TVP model is estimated using "tvLM" function of "tvReg" in R [START_REF] Casas | tvReg: Time-varying coefficients linear regression for single and multi-equations in R[END_REF]2020). Both of the AIbased techniques are trained with "nnfor" package in R [START_REF] Kourentzes | nnfor: Time series torecasting with neural networks[END_REF], using "mlp" function for MLP and MLPX, and "elm" function for ELM and ELMX.

Stacking models

To train the eleven preliminary single models as well as the stacking-based algorithms. Each tourism demand series is divided into two training sets and one testing set. The testing set of every tourism demand series has identical timeframes, from 2019Q1 to 2019Q4. The division of two training sets varies according to the starting point of the particular tourism demand series. The estimation of each preliminary single model uses the data of the first fold of the training set as the inputs and the width of the second fold of training data as the rolling window. Multiple forecasts generated through rolling are used as the basis of stacking in the following steps.

The 11 preliminary single models are stacked in 26 different ways to explore the best stacking algorithm in the current context. In particular, one regression-based stacking, 24 "best n" stacking, and a "seasonal naïve" stacking are considered. The regression-based stacking takes the forecasts of the 11 preliminary single models as explanatory variables and fits them to the actual tourism demand. The preliminary single models with statistically significant coefficients are chosen and weighted by their standardized coefficients. The preliminary single models with insignificant coefficients are ignored in the weighing process. If no preliminary single model has significant coefficients, the regression-based stacking does not generate a forecast. The "best n" type stacking selects several preliminary single models according to certain criteria and combines the selected models using different averaging methods. In the present competition, we consider the selection of three and five algorithms, using MAPE, RMSE, or MASE as criteria, and combine these algorithms using simple average, error weighted average, square-error weighted average, and accuracy measure weighted average. Therefore, a total of 24 (two by three by four) "best n" type of stacking are generated. The "seasonal naïve" examines all 25 stacking algorithms above and allows different stacking algorithms to forecast different quarters.

For each origin-destination pair, 10 absolute scaled errors (ASE) are calculated as

ASE = |𝑒 𝑑,𝑗 | 1 𝑇 -1 ∑ |𝑌 𝑑,𝑡 -𝑌 𝑑,𝑡-𝑚 | 𝑇 𝑡=2
, where 𝑒 𝑑,𝑗 is the difference between the forecast value and the actual value of international visitor arrivals of origin-destination pair 𝑑 at time 𝑗; 𝑌 𝑑,𝑡 represents the actual value of international visitor arrivals of origin-destination pair 𝑑 at time 𝑡, and 𝑚 is the seasonal period. The 10 ASEs correspond to 1-step-ahead forecasts of 2019Q1, 2019Q2, 2019Q3, and 2019Q4, 2-step-ahead forecasts of 2019Q2, 2019Q3, and 2019Q4, 3-step-ahead forecasts of 2019Q3 and 2019Q4, and 4-step-ahead forecast of 2019Q4. The arithmetic average of these 10 ASEs comprises the mean absolute scaled error (MASE) of the specific origin-destination pair. The arithmetic average and the standard deviation of origin-destination specific MASE are used as the standard to select the best algorithm among all 37 algorithms (11 preliminary single models and 26 stacking-based algorithms).

Second Stageex ante scenario forecasting under COVID-19

The best performing model from the first stage is used to generate a baseline ex-ante forecast for the period 2020Q1-2021Q4.

The COVID-19 pandemic generated a sudden sharp shock to world travel and tourism in 2020.

The most recent data available, as of September 15 th (UNWTO, 2020a,b), show that international visitor arrivals decreased by 65% in the first quarter of 2020 over to the same period last year and that a bottom occurred during the second quarter of 2020 with a 95% decline relative the same quarter in the previous year.

By comparison, world visitor arrivals had declined by only -0.4% in 2003 because of the SARS epidemic, and by -4.0% in 2009 because of the global economic crisis (UNWTO, 2020c). As a result, international arrivals took 11 months after SARS and 19 months after the global economic crisis to recover to their pre-crisis levels (UNWTO, 2020d).

There is considerable uncertainty about the magnitude of the COVID-19 impact on world tourism and the speed of its recovery. In the three scenarios presented by the UNWTO (2020d), the decrease in international visitor arrivals in 2020 ranges from -58% to -78%, depending on when (July, September, or December) borders gradually reopen and travel restrictions are lifted. A recovery of international tourism arrivals to their 2019 level is estimated to take from 2.5 to 4 years, according to the UNWTO (2020a), and up to five years, according to Tourism Economics (2020).

A strong rebound in 2021 is "based on the assumption of a reversal of the evolution of the pandemic, significant improvement in traveller confidence and major lifting of travel restrictions by the middle of the year" UNWTO (2020a).

We define three scenarios to reflect mild, medium, and severe impacts based on this recovery while taking into account national specificities.

-Scenario 1: Mild (a V-shaped pattern)

The pandemic fades out at the end of 2020, the international travel of each specific original-destination pair starts to recover as soon as the bilateral travel bans are lifted, and the recovery to the 2019 average level takes 12 quarters (three years). -Scenario 2: Medium (a deep V/U-shaped pattern)

The pandemic fades out in the first half of 2021, the recovery of international travel for each specific original-destination pair starts slowly after travel restrictions are removed, and the recovery to the 2019 average level takes 16 quarters (four years). -Scenario 3: Severe (an L shape pattern)

The influence of the pandemic does not fade away until the vaccine is available to the general public (approaching the end of 2021), international travel for each specific original-destination pair remains at its bottom for two more quarters relative to scenario 2, and the recovery to the 2019 average level takes 16 quarters (four years).

Our forecasting methodology is based on the following steps for each origin-destination pair:

Step Step 2. Determine the number of tourism arrivals in 2020Q1 and 2020 Q2, using annual percentage changes (over the same quarter in 2019) for each destination, from the "UNWTO World Tourism Barometer and Statistical Annex, August/September" (2020a), or national sources, if not already available in the original database.

Step 3. Determine the number of arrivals during the period of severe travel restrictions, by replicating the bottom value of bilateral flows for 2020Q2.

Step 4. Determine the number of arrivals during the first two quarters of recovery, by using the following annual percentage changes over the average value for 2019:

• -70.0% for the first quarter, then -40.0% for the second quarter for scenario 1 • -90.0% for the first quarter, then -67.0% for the second quarter for scenario 2 • -95.2% for the first quarter, then -94.0% for the second quarter for scenario 3 These parameters were calibrated so that the 2019 to 2020 annual percentage changes would be equal to -58% for scenario 1, -70% for scenario 2, -78% for scenario 3 (similar to UNWTO's three scenarios), using the actual annual percentage changes for 2020Q1 (-27,8%), for 2020Q2 (-95,2%), and the values above for 2020Q3 and 2020Q4. For consistency, the results are also compared with two benchmarks (an annual forecast based on a gravity model, and a forecast from Tourism Economics, 2020).

Step 5. Determine the number of arrivals at the date (quarter) of recovery (= start date for recovery + 11 quarters for scenario 1 or 15 quarters for scenarios 2 & 3), to be equal to its average value for 2019.

Step 6. Determine the number of arrivals during the period of recovery, by linearly linking the value of the second quarter of recovery to the value of the date of recovery, to form a trend line of adjustment;

Step 7. Seasonally adjust the trend line with observed values prior to the events.

The following figure illustrates the adjustment process.

[Figure 2 about here]

RESULTS AND DISCUSSION Empirical Results from the First Stage

Table 1 indicates the performance of all 37 models including 11 preliminary single models and 26 stacking algorithms. The MASE values reflect the forecasting accuracy of each model, and the standard deviations of the MASEs reflect the robustness of the forecasting accuracy for each model. Ideally, the models with smaller MASEs and smaller standard deviations of MASEs are more favourable. Amongst these 37 models, regression stacking based on all 11 single models performs extremely poor with the largest MASE and the standard deviation (4.1632 and 23.8915 respectively). After inspections of origin-destination specific forecasts, it is found out that regression stacking is extremely sensitive to the quality and volumes of the data. The algorithm generates unreliable forecasts when the tourism demand series is short (i.e. less data). Given its very poor performance, this model is dropped from the subsequent discussion.

[Table 1 about here]

The results of the rest 36 models are further shown in Figure 3, in which the blue line refers to the MASEs and the orange dashed line refers to the standard deviations of these MASEs. From the figure, it is seen that overall stacking algorithms outperform the single models on average. The average MAPEs for the 11 single models and the 25 stacking algorithms are 1.4287 and 1.0122 respectively, and the average MASE standard deviations for the 11 single models and the 25 stacking algorithms are 1.3462 and 0.9406 respectively (see Table 1). Therefore, stacking algorithms are more effective in generating more accurate and more robust forecasts of tourism demand. This conclusion is consistent with the existing studies that combined forecasts, based on a number of single models, often outperform the corresponding single models (see, for example, [START_REF] Li | The combination of interval forecasts in tourism[END_REF][START_REF] Song | A review of research on tourism demand forecasting: Launching the Annals of Tourism Research Curated Collection on tourism demand forecasting[END_REF][START_REF] Wu | New developments in tourism and hotel demand modeling and forecasting[END_REF]. In Figure 3 it is also observed that stacking models exhibit less variation whereas the performance of the 11 single models varies more widely.

[Figure 3 about here] Table 1 also shows that among the 11 single models, five single models of SARIMA, ETS, STL, ELM, and ELMX outperform the benchmark model of SNAIVE, with lower MASEs. In considering robustness among these five models, three models have smaller standard deviations than the benchmark of SNAIVE. These models are SARIMA, ETS, and STL (Figure 4). The favourable performance of these three models provides empirical evidence for the adoption of these models in the tourism demand forecasting using single models.

[Figure 4 about here]

Where stacking algorithms are concerned, Figure 5 shows that the MASEs do not vary much across models. Their standard deviations are also quite similar. The MASEs vary between 0.997 and 1.027, and the standard deviations vary between 0.906 and 0.985 for 25 stacking models (see Table 1). Compared to the benchmark model of SNAIVE, all 24 "best n" stacking algorithms outperform the benchmark for both MASEs and their standard deviations. Amongst these algorithms, four stacking algorithms, namely bn5rs, bn5re, bn5r2, and bn5rc perform best, having both the smallest MASEs and the smallest MASE standard deviations. We therefore conclude that the optimal number of single models to be combined in this study is five.

[Figure 5 about here] Table 2 reports the forecasting accuracy of the stacking algorithm "bn5re", which combines the five best single models by minimising the RMSE, for each destination. The stacking algorithm "bn5re" is the best performed one among 11 single models and 26 stacking algorithms, with the lowest MASE of 0.9969 (see Table 1). It is observed that on average the forecasts are more accurate for lower forecasting horizons. In addition, the forecasting accuracy varies over origindestination pairs. For example, when the one quarter ahead forecasting is concerned, the MASEs vary from 0.2198 (the case of South Africa) to 2.6232 (the case of Chile) (The results for destinations are available upon request). This observation is consistent with the existing studies [START_REF] Song | A review of research on tourism demand forecasting: Launching the Annals of Tourism Research Curated Collection on tourism demand forecasting[END_REF][START_REF] Wu | New developments in tourism and hotel demand modeling and forecasting[END_REF].

[Table 2 about here]

In conclusion, this stage focuses on the forecasting competition before COVID-19 based on 11 single models and 25 stacking-based algorithms. Some interesting findings are discovered. Firstly, the forecasting performance of the 11 single models varies. Further, SARIMA, ETS, and STL perform best taking into consideration both forecasting accuracy and robustness. Secondly, compared to single models, the stacking-based algorithms provide more accurate and more robust forecasts on average. Thirdly, the optimal number of stacked single models is five. Fourthly, RMSE is the best criteria in selecting the single models to be included in the stacking-based algorithms. These findings provide helpful information for future tourism demand forecasting practices.

Empirical Results from the Second Stage

The COVID-19 pandemic spreads first in Asia Pacific (China, South Korea, Japan, Thailand, Australia, Indonesia, Malaysia, New Zealand, Singapore), then to North America (USA, Canada, Mexico, Chile), then to Europe (UK, Bulgaria, Sweden, Finland, Czech Republic), and finally Africa (Mauritius, South Africa, Tunisia), in the order of first confirmed COVID-19 deaths all occurring in the first quarter of 2020 (WHO).

Table 3 provides a summary of second stage forecast indicators (deepest impact, bottom date, 2021 recovery rate) for the three scenarios across the 20 nations.

[Table 3 about here] Almost all destinations closed their borders (or banned travellers from high-risk regions, for Japan, Mexico, Singapore, South Korea, Sweden, the UK, and the USA) by the end of March 2020 (Our World in Data, 2020b). However, these policy responses to COVID-19 (IMF, 2020) differed across the countries in the sample. Most European nations (and Mauritius) opened their borders to travel in July 2020 (i.e. bottom date = 2020Q2) to generate some tourism revenues from Northern Hemisphere summer vacationers, whereas South Africa and Malaysia plan to relax their travel restrictions in 2020Q4, and the remainder of the destinations are waiting for 2021Q1 to do so. This bottom date is often the same for all origin countries for a given destination, but in some cases, tourists from high-risk regions are not permitted to enter a destination (e.g. visitors from the USA to Europe) while other source markets are permitted.

The overwhelming majority of destinations experienced an almost total collapse of their international tourists during 2020Q2 (their deepest impact ranging from -95% to -100%); the exception being Bulgaria, Indonesia, Mexico, and Sweden. The latter is a special case in Europe since the Swedish government decided to not to resort to a lockdown in the first place and have relied on voluntary social distancing and a few other light constraints, in response to COVID-19.

The 2021 recovery rate ranges from 53% to 70% for the mild scenario 1 (after a 58% fall in 2020), from 29% to 45% for the medium scenario 2 (after a 70% decline in 2020), and only from 9% to 23% for the severe scenario 3 (after a 78% collapse in 2020). The latest UNWTO World Tourism Barometer August/September 2020 issue, which was published on September 15 th , comments: "current trends point to a decline in international arrivals closer to 70% for 2020" (UNWTO, 2020b). Yet, a great deal of uncertainty remains about the evolution of the pandemic and thus the most likely scenario.

We select the following five destinations, Australia, Sweden, Thailand, Tunisia, and the USA, to demonstrate our results.

Australia

Australia is one of the destinations that has implemented a continuous ban of international tourism arrivals since March 2020 to significantly reduce the likelihood of importing COVID-19 cases. The Australian Government (2020) is nevertheless establishing a safe travel zone with New Zealand, a very low-risk country, starting from October 16 th , 2020. Australia also indicated they would gradually lift other travel restrictions in 2021, notably for international students. The following six graphs clearly show the absence of international arrivals for the last three quarters of 2020, and the beginning of a soft recovery in 2021, where rates are equal to 55% for scenario 1 (red line), 31% for scenario 2 (blue line), and 11% for scenario 3 (green line), respectively.

[Figure 6a about here] Sweden Sweden is also a special case due to its unique strategy of pursuing herd immunity without any lockdown or strict sanitary measures. Contrary to most of its European neighbours, Sweden experienced an early wave of infection in June and July 2020, in which COVID-19 almost disappeared, in contrast to a second wave of infections appearing in the rest of the European Union. As a result, forecasted tourism arrivals in Sweden show a rapid and strong rebound, where recovery rates are some of the highest, i.e. 70% for scenario 1 (red line), 45% for scenario 2 (blue line), and 23% for scenario 3 (green line).

[Figure 6b about here] Thailand Thailand is one of the hardest-hit economies due to its dependence on tourism revenues. The IMF (2020b) concludes that "in Thailand, a decrease in tourism due to COVID-19 could bring the country's overall exports down by eight percentage points of GDP and have a direct net impact of about six percentage points of GDP on its current account balance in 2020. That could erode part of the seven percent overall current account surplus the country had in 2019". Since October 1 st , Phuket has reopened its tourism, but under strict rules such as a 14-day hotel quarantine, which will discourage many tourists from coming to Thailand. The graphs related to forecasted tourism flows to Thailand show recovery starting in 2021, whose rates are equal to a meager 58% for scenario 1 (red line), 33% for scenario 2 (blue line), and 11% for scenario 3 (green line).

[Figure 6c about here] Tunisia Tunisia also reopened its borders to international tourists at the end of June 2020, but a color-coded system based on COVID-19 risk (green for low risk; orange for medium risk; red for high risk). Low risk involves completing a form. Medium risk requires a negative PCR test to enter Tunisia while high risk necessitates a mandatory 14-day quarantine at a public center. As a result, tourism is likely to recover slowly in 2021, with a moderate recovery of 58% for scenario 1 (red line), 33% for scenario 2 (blue line), and 11% for scenario 3 (green line), and with strong seasonality, as shown on the six graphs related to this destination.

[Figure 6d about here] The United States The United States closed its borders to tourists from Canada, China, Europe, and Mexico in 2020. A gradual recovery is expected in 2021, as shown on the six corresponding US graphs, with rates of 56% for scenario 1 (red line), 32% for scenario 2 (blue line), and 11% for scenario 3 (green line). The US elections will probably play a role in the policies to be implemented against COVID-19 in 2021, and this will also affect recovery.

[Figure 6e about here]

In conclusion, these second-stage forecasts are highly dependent on the policies implemented by countries regarding the pandemic and its evolution. There is a great deal of uncertainty regarding these policies including the timing of a vaccine and the extent to which tourists will regain confidence in international travel.

CONCLUSIONS

Being totally unanticipated, COVID-19 brought international travel to an immediate halt, severely impacting the global economy and in particular, the tourism and hospitality sector. The impact is far-reaching. Both private and public sectors have instituted various countermeasures, both travelrelated as well as health-and safety-related, to combat the spread of the contagious virus. At the time of writing (October 2020), the full effects of these measures remain to be unseen. Thus, it is important for tourism economists to estimate evidence-based forecasts in how the tourism and hospitality industry will recover. In this paper, we forecast inbound tourism demand for 20 countries using a two-stage forecast framework: ex post forecasting before COVID-19 and ex ante forecasting post-COVID-19.

The main findings are as follows: regarding the first stage ex post forecast models, among the single models SARIMA, ETS and STL performed well in terms of accuracy and robustness. Furthermore, the stacking-based algorithms provide more accurate and robust forecasts, based on MASEs. The optimal number of single models that could be used for the combination model is five. These methodological findings contribute to constructing more accurate and robust forecast framework.

Regarding the second stage ex ante forecasting, three scenarios of different patterns of inbound tourism recovery were modelled: a mild scenario displaying a V-shaped rebound, a medium scenario which is a deep V/U-shape, and a severe scenario showing an L-shaped recovery. These scenarios are adapted for each country-wise situation. The second stage results demonstrate wide variations in the 2021 recovery rates depending on the three scenarios, allowing for country-wise adjustments. Compared to the same quarter in 2019, the mild scenario shows international arrivals ranging from 53% to 70% of the same level over the comparative period. The medium scenario has international arrivals being between 29% and 45% of the previous level while the recovery in the severe scenario forecasts tourism arrivals ranging from 9% to 23% compared to the same quarter in 2019. Although there is considerable uncertainty regarding the likely scenario, a country's travel policies will influence recovery speeds. For instance, Sweden has opted for a unique open-door no-lockdown herd-immunity policy, so may experience faster recovery than those countries with strict travel-control policies.

Finally, the COVID-19 pandemic is unprecedented so generating accurate and robust forecasts is a most challenging task for economists. Nevertheless, this paper has used state-of-the-art forecasting techniques to provide policymakers and tourism businesses with the tourism demand forecasts under different scenarios, with the available information at the time and forecasting competition results.. This framework is also particularly useful for impact evaluation of crises on tourism demand and its recovery estimation.

Once a vaccine is found and tourism has converged to its new normal, further research could analyse international arrivals and improve their forecasts by integrating models based on its fundamental long-term determinants (such as gravity models), which are of little use for the shortterm forecasting exercise in this paper, but could become again valuable tools for future investigations. Another direction for research could be an analysis of the impact of greater uncertainty on confidence intervals in forecasts, due to exogenous shocks such as the pandemic and its evolution as a result of various policies. Indeed, in a matter of weeks, the mild scenario has now become much less probable, after a second wave of infections which appears to be suddenly more difficult to manage than expected, as the time of writing (October 2020). Note: The first two digits, "bn", stand for "best n" type of stacking methods, with the third digit representing the number of preliminary single models used in stacking, the fourth digit representing the criteria of selecting the best algorithms (p -MAPE, r -RMSE, a -MASE), and the fifth digit representing the method of combining the selected algorithms (ssimple average, eerror weighting, 2squared error weighting, cselection criteria based weighting).

  1. Determine important dates (which quarter and number of quarters) for the start of COVID-19, duration of travel restrictions, the start of recovery, and length of recovery:• The start date for COVID-19 (2020Q1 for all countries): cf. WHO (2020a) • Period of severe travel restrictions, defined by the "international travel controls" index (a component of the Oxford Coronavirus Government Response Tracker: cf. Our World in Data, 2020b) equal to its maximum value on a range from 0 to 4, or from the following sources: France Diplomatie; International Monetary Fund; Travel Off Path; U.S. Department of State; World Nomads. • The bottom date for COVID-19 (last quarter of severe travel restrictions = start date for COVID-19 + number of quarters in the period of severe travel restrictions) • The start date for recovery (= bottom quarter + 1) • Date of recovery to the 2019 average level (= start date for recovery + 11 quarters for scenario 1 or 15 quarters for scenarios 2 & 3)

  Figure 3. Forecasting performance of all 36 models

Table 1 . MASEs and standard deviation of MASEs for all models MASE Std. Deviation of MASEs MASE

 1

	Std. Deviation
	of MASEs

Table 2 Forecasting performance (MASE) for all destinations

 2

		1Q-ahead	2Q-ahead 3Q-ahead	4Q-ahead
	Australia	0.6456	0.4715	0.5475	0.6310
	Bulgaria	0.5535	0.6797	0.6734	3.0527
	Canada	1.1054	1.2251	2.1464	1.8193
	Chile	2.6232	1.7403	1.4582	1.8180
	Czech	0.7521	0.8005	0.7496	0.6241
	Finland	0.9505	0.7558	1.5307	1.1113
	Indonesia	0.8908	1.1560	0.8358	1.5825
	Japan	1.2069	2.0178	1.6670	1.4687
	South Korea	0.5486	0.6263	1.4568	0.9181
	Malaysia	1.1541	0.9210	0.7220	0.8479
	Mauritius	1.1529	1.1363	1.3825	1.4843
	Mexico	0.8510	0.6409	0.6064	0.7628
	New Zealand	0.4812	0.7479	0.9676	0.9961
	Singapore	0.7303	0.5787	0.8325	0.9167
	South Africa	0.2198	1.1358	0.9472	0.7696
	Sweden	0.5929	0.7497	1.5836	1.3647
	Thailand	0.2597	0.5225	0.4272	0.8387
	Tunisia	0.4068	0.8598	1.4806	0.6688
	UK	1.2188	1.1081	0.9035	0.9434
	USA	0.2944	0.5090	0.5238	0.6737
	Overall	0.8319	0.9191	1.0721	1.1646

Note: "Deepest impact" refers to the deepest decline of visitor arrivals in percentage against the same quarter in 2019; "Bottom time" is the quarter when the deepest impact takes place; "2021 recovery rate" refers to the ratio of the predicted annual total visitor arrivals in 2021 against the annual visitor arrivals in 2019.


Lauren T W Fey

Marie Z Anne

Tan

Thomas D Swinburne

Danny Perez

Dallas R Trinkle

Accelerated molecular dynamics simulations of dislocation climb in nickel

The mechanical behavior of materials operating under high temperatures is strongly influenced by creep mechanisms such as dislocation climb, which is controlled by the diffusion of vacancies. However, atomistic simulations of these mechanisms have traditionally been impractical due to the long time scales required. To overcome these time scale challenges, we use Parallel Trajectory Splicing (ParSplice), an accelerated molecular dynamics method, to simulate dislocation climb in nickel. We focus on modeling the activity of a vacancy near a jog on an edge dislocation in order to observe vacancy pipe diffusion and vacancy absorption at the jog. From rigorously constructed trajectories encompassing more than 2000 vacancy absorption events over a simulation time of more than 4 μs at 900 K, a comprehensive sampling of available atomistic mechanisms is collated and analyzed further with molecular statics calculations. We estimate average rates for pipe diffusion and vacancy absorption into the jog using data from the dynamic and static calculations, finding very good agreement. Our results strongly suggest that the dominant mechanism for vacancy absorption by jogs is via biased diffusion to the dislocation core followed by fast pipe diffusion to the jog.

I. INTRODUCTION

At high temperatures, crystalline materials deform through a time-dependent process known as creep [START_REF] Meyers | Dynamic Behavior of Materials[END_REF][START_REF] Haasen | Physical Metallurgy[END_REF]. Understanding creep is critical for high-temperature applications, such as in the nickel superalloy turbine blades in jet engines [START_REF] Caron | Improvement of creep strength in a nickel-base single-crystal superalloy by heat treatment[END_REF][START_REF] Pollock | Creep resistance of CMSX-3 nickel base superalloy single crystals[END_REF][START_REF] Murakumo | Creep behaviour of Ni-base single-crystal superalloys with various γ volume fraction[END_REF][START_REF] Rae | Primary creep in single crystal superalloys: Origins, mechanisms and effects[END_REF][START_REF] Smith | Diffusion processes during creep at intermediate temperatures in a Ni-based superalloy[END_REF]. These engines are more fuel efficient at higher temperatures, yet the operating temperatures are limited by the creep behavior of the component materials. Therefore, understanding creep mechanisms and rates is critical for advancing these applications. One mechanism of creep is dislocation climb, which involves the absorption of vacancies into dislocation cores, allowing dislocations to move out of their slip plane and glide around obstacles.

Simulations of dislocation climb are challenging due to the long time scales required. Dislocation climb requires the flux of vacancies through the bulk toward dislocations, which is a slow process compared to the nanosecond time scales for atomistic simulations. Dislocation jogs, which can be created thermally or by the intersection of dislocations, then absorb the vacancies in order to climb. These atomistic mechanisms are difficult to observe experimentally, so simulation techniques are necessary to advance our understanding. Most models of dislocation climb ignore the atomistic details of the dislocation climb mechanism but instead describe climb phenomenologically. Several groups have studied dislocation climb using dislocation dynamics methods [START_REF] Mordehai | Introducing dislocation climb by bulk diffusion in discrete dislocation dynamics[END_REF][START_REF] Gao | 3D discrete dislocation dynamics study of creep behavior in Ni-base single crystal superalloys by a combined dislocation climb and vacancy diffusion model[END_REF][START_REF] Danas | Plane-strain discrete dislocation plasticity with climb-assisted glide motion of dislocations[END_REF]. Mordehai et al. was the first to couple dislocation dynamics with * dtrinkle@illinois.edu vacancy diffusion [START_REF] Mordehai | Introducing dislocation climb by bulk diffusion in discrete dislocation dynamics[END_REF], assuming that the rate-limiting step in dislocation climb is the bulk diffusion of vacancies to the dislocation lines. However, these models assume that vacancies can be absorbed anywhere along the dislocation line, ignoring the jog density and the rate of pipe diffusion, i.e., the potential acceleration of vacancy diffusion along the dislocation line. Pipe diffusion has been incorporated into some discrete dislocation dynamics models [START_REF] Gao | Investigations of pipe-diffusion-based dislocation climb by discrete dislocation dynamics[END_REF], but the interactions of vacancies with the jogs is not accounted for. Other modeling efforts for dislocation climb include phase field models [START_REF] Ke | Phase field microelasticity model of dislocation climb: Methodology and applications[END_REF][START_REF] Geslin | A phase field model for dislocation climb[END_REF]. These are also mesoscale models which do not explicitly include the atomistic mechanisms involved in dislocation climb, such as pipe diffusion and vacancy absorption at jogs. To parametrize and improve these models, the rates and mechanisms of the individual steps of dislocation climb must be understood. Kabir et al. used atomistic methods to calculate the vacancy migration barriers near a dislocation and used the results to run kinetic Monte Carlo simulations [START_REF] Kabir | Predicting Dislocation Climb and Creep from Explicit Atomistic Details[END_REF]. Baker and Curtin coupled a discrete diffusion model to an accelerated molecular dynamics simulation to study the diffusion of vacancies to a dislocation and the nucleation of jogs [START_REF] Baker | Multiscale diffusion method for simulations of long-time defect evolution with application to dislocation climb[END_REF]. However, this model did not include pipe diffusion due to limited length scales. Density functional theory has also been used to calculate the vacancy migration barriers for pipe diffusion [START_REF] Wirth | Density functional study of self-diffusion along an isolated screw dislocation in fcc Ni[END_REF], which were then used to parametrize a kinetic Monte Carlo model. These results show that pipe diffusion is orders of magnitude faster than bulk diffusion and can contribute to the overall creep rate when the dislocation density is high.

In this paper, we use ParSplice [START_REF] Perez | Long-time dynamics through parallel trajectory splicing[END_REF], an accelerated molecular dynamics technique, to study dislocation climb in nickel. ParSplice manages the production of many short, independently generated, molecular dynamics trajectories in parallel, which can then be rigorously "spliced" to form a single long trajectory between metastable states, thereby allowing for timewise parallelization and hence for long time simulations. The method is especially efficient to accelerate the production of escape trajectories from large superbasins of energy minima. This is ideally suited to the problem we study here, namely, the thermally activated motion of a vacancy near an edge dislocation jog, a prototypical process in dislocation climb. From these simulations, we extract rates for vacancy pipe diffusion near and far from the jog site and further analyze found mechanisms for vacancy absorption and emission in the context of harmonic transition state theory via the use of static energy barrier calculations. It is important to emphasize that ParSplice generates long-time trajectories in a totally unbiased fashion. Indeed, every local basin of attraction on the potential energy surface is an allowed state and no restrictions are placed on possible transitions. The only factor that controls the accuracy of the method is the so-called correlation time τ c used in the construction of the segments [START_REF] Perez | Long-time dynamics through parallel trajectory splicing[END_REF][START_REF] Bris | A mathematical formalization of the parallel replica dynamics[END_REF]; when states are defined as basins of attractions on the potential energy surface, τ c is a characteristic vibrational timescale. It can be shown that the accuracy of the spliced trajectories increases exponentially with τ c [START_REF] Bris | A mathematical formalization of the parallel replica dynamics[END_REF][START_REF] Perez | The parallel replica dynamics method-Coming of age[END_REF], so that the resulting trajectories can be made essentially MD-accurate if τ c is chosen conservatively.

Our main result is that the vacancy flux to jog sites is dominated by biased vacancy diffusion to the dislocation line, followed by rapid pipe diffusion to the jog site. This is to be distinguished from vacancy diffusion to the jog site directly from the bulk. The implications of our findings for the multiscale modeling of climb and therefore creep are discussed.

II. METHODS

In order to simulate the process of dislocation climb, an edge dislocation with a jog was created. All calculations are made using LAMMPS [START_REF] Plimpton | Fast parallel algorithms for short-range molecular dynamics[END_REF] with the Mishin EAM potential for nickel [START_REF] Mishin | Interatomic potentials for monoatomic metals from experimental data and ab initio calculations[END_REF]. This potential was chosen for its known ability to reproduce experimental values of the equilibrium lattice parameter, cohesive energy, elastic constants, intrinsic stacking fault energy, and bulk vacancy formation energy of face-centered cubic (FCC) nickel. An initially perfect edge dislocation with a total line vector of N t + j was created in a cylinder, where t is the line direction of the pure edge segment, j is the line direction of the jog segment, and N is an integer Fig. 1(a). The dislocation has a Burgers vector of a 2 [110] and t in the a 2 [1 12] direction, so it lies in the ( 111) slip plane. Since dislocation jogs can be formed by the intersection of dislocations, the jog has a line direction of a 2 [1 10], a Burgers vector on a different slip plane in FCC. On relaxation, the initially straight dislocation line splits into pure edge segments and a jog segment, with the edge components dissociating into partial dislocations separated by a stacking fault [Fig. 1(b)]. This partial decomposition lowers the line energy of the dislocation, but it leaves a stacking fault region between the two partial dislocations. At the jog, the two edge dislocations are closer together due to the higher stacking fault energy on the (001) plane. In the stacking fault region, the stacking of atomic planes resembles hexagonal close packing instead of the equilibrium FCC packing of nickel. For this reason, it is useful to use the centrosymmetry parameter as a visualization tool. The centrosymmetry parameter measures the inversion symmetry around each atom, giving zero for perfect FCC crystals and higher values near defects such as dislocations or vacancies; removing atoms with low centrosymmetry allows us to isolate the defective regions. Interestingly, depending on the plane in which the jog exists, two different jog configurations can be formed. These will be referred to as Jog A and B [Fig. 1(c)]. The two jog configurations arise due to the two distinct (001) planes on which the jog, with a line direction of a 2 [1 10], can intersect the dislocation line. Both jog configurations have the same partial splitting distance of 19.3 Å far from the jog, but Jog A is slightly more compact than Jog B, as can be seen by counting the number of atoms shown in the jogged region. When a jog absorbs a vacancy, the jog is displaced by one half a dislocation threading vector. Two vacancy absorption events are required to return to the original configuration and displace the jog a full threading vector.

The total energy of the system includes contributions from the cohesive bulk energy, the dislocation line energy, and the jog energy. The jog energy was estimated for each jog configuration by calculating the total energy of the system for different N, which controls the jog-jog distance. Extrapolating to an isolated jog with a distance of 0 gives the energy of the jog only without the edge dislocation segments (Fig. 2). The energy for both Jog A and Jog B is 3.52 eV. Therefore, which configuration a jog forms is determined by geometric considerations only.

Two different system sizes were created: A large system with approximately 120,000 atoms and a smaller system with approximately 16,000 atoms. The smaller system was used for the accelerated molecular dynamics simulations to reach longer timescales, while the larger system was used to verify the results through molecular statics simulations. During all simulations, periodic boundary conditions are applied along the dislocation line, while fixed boundary conditions are used perpendicular to the dislocation line. The fixed boundary conditions are implemented by fixing an outer ring of atoms at positions prescribed by the anisotropic elastic solution for a perfect edge dislocation [START_REF] Bacon | Anisotropic continuum theory of lattice defects[END_REF]. In the large system, this corresponds to a 60 Å radius cylinder being simulated, while the smaller system has a radius of 30 Å. The lengths of the dislocations lines were also varied by adjusting N, the number of pure edge segment threading vectors. For the large system, N is 15, and for the small system, N is 8. To ensure that the small system size is still representative of a dislocation in bulk, the energy barrier for pipe diffusion was calculated for both system sizes using the nudged elastic band method (NEB) [START_REF] Henkelman | Improved tangent estimate in the nudged elastic band method for finding minimum energy paths and saddle points[END_REF][START_REF] Henkelman | Climbing image nudged elastic band method for saddle points and minimum energy paths[END_REF]. The initial and final states of the NEB were taken directly from the ParSplice simulations. The average of the forward and reverse barriers is 0.86 eV for the large system and 0.87 eV for the small system, suggesting that pipe diffusion is not significantly impacted by the change in system size. Additionally, the jog splitting distance, defined as the distance between the two partial dislocations at the jog, is similar for the two system sizes. Using Ovito's dislocation analysis tool [START_REF] Stukowski | Visualization and analysis of atomistic simulation data with OVITO-the Open Visualization Tool[END_REF], the two atomic sites where the dislocation lies closest to the jog were identified, and the distance between these two sites is defined as the jog splitting distance. For Jog A, this distance is 11.7 and 11.6 Å for the small and large systems, respectively. For Jog B, the respective values are 14.0 and 14.1 Å. Therefore, we feel confident that performing the accelerated MD in the small system will still yield results that are representative of behavior in the bulk.

We used ParSplice [START_REF] Perez | Long-time dynamics through parallel trajectory splicing[END_REF], an accelerated molecular dynamics method, to simulate a vacancy near the jogged dislocation and calculate rates of vacancy pipe diffusion and absorption rates. The simulations were run at 900 K using a Langevin thermostat with a timestep of 2 fs. Each ParSplice simulation used 200 replicas with four cores assigned to each replica. The correlation time used was 2 ps. In each simulation, a vacancy was added to the system at an atomic site near or inside the partial dislocations. Each vacancy site was chosen manually using the centrosymmetry parameter to identify the dislocation line. Figure 1(c) shows these partial dislocation sites in dark blue. Each state identified by ParSplice was automatically analyzed to determine the position of the vacancy. This is done by identifying a cluster of atoms that are all missing a nearest neighbor in approximately the same location. An atom is considered to be missing a nearest neighbor if the sum of its vectors to nearest neighbors has a magnitude greater than 0.75 Å. The position of the vacancy is estimated as the mean position of these missing neighbors. If the vacancy position moves less than 0.1 Å between consecutive ParSplice events, those states are considered the same for the purposes of vacancy diffusion. This automatic vacancy identification method allows for many simulations to be statistically analyzed together to calculate vacancy diffusion rates. In addition to the 900 K simulations, selected simulations were also run at 800 and 1000 K to evaluate entropic effects.

To support the ParSplice simulation results, we also perform molecular statics calculations on the larger jogged dislocation system. Vacancy-dislocation interaction energies near the dislocation, which measure the energy change of a vacancy due to the nearby dislocation, were calculated. The interaction energy can be predicted using elasticity theory for vacancy positions far from the dislocation; however, near the dislocation, elasticity theory breaks down and atomic simulations are required to accurately determine these interaction energies [START_REF] Clouet | The vacancy-edge dislocation interaction in fcc metals: A comparison between atomic simulations and elasticity theory[END_REF]. The interaction energies were calculated for a variety of lattice sites by measuring the total energy of the system with vacancy at that site. The interaction energy is given by

E int D+V = E D+V -E D + E 0 -E F V , (1)
where E int D+V is the interaction energy, E D+V is the total energy of the dislocation-vacancy system, E D is the energy of system with no vacancy, E 0 is the cohesive energy, and E F V is the bulk vacancy formation energy.

Additionally, the energy barriers for vacancy movement within the dislocation and jog were calculated with the NEB method [START_REF] Henkelman | Improved tangent estimate in the nudged elastic band method for finding minimum energy paths and saddle points[END_REF][START_REF] Henkelman | Climbing image nudged elastic band method for saddle points and minimum energy paths[END_REF]. Starting and ending configurations were chosen for specific vacancy events: pipe diffusion along the dislocation line, vacancy absorption into the jog, vacancy movement out of the slip plane, etc. An initial path is created by linearly interpolating between the starting and ending configurations.

III. RESULTS AND DISCUSSION

Using the relaxed dislocation system, we performed Par-Splice simulations and observed the vacancy diffusing along the partial dislocations and being absorbed into the jogs. By using this accelerated MD technique, we are able to achieve total simulation times of over 2.7 and 1.6 μs for vacancies diffusing and absorbing into Jog A and Jog B, respectively. In each simulation, the vacancy diffuses randomly along the partial dislocation cores until it reaches the jog and is absorbed. The vacancy spends the vast majority of the simulation time on the compressive side of the dislocation, and any hops across the slip plane are quickly followed by hops back to the compressive side. After the jog absorbs the vacancy, it will occasionally emit vacancies, which are always immediately reabsorbed. In total, 1843 vacancy hops along the partials and 18 initial vacancy absorptions were observed during the combined simulation time. Including vacancy reabsorption, there were a total of 510 vacancy absorptions.

We are first interested in how the vacancy diffuses within the partial dislocations. The vacancy escape time is defined as the amount of time between vacancy hops within the partial dislocations. The escape times for every simulation are calculated, not including hops into or out of the jog, and are plotted in Fig. 3. Occasionally, the observed vacancy diffusion took place in the bulk or outside of the slip plane instead of within the dislocation because the vacancy jumped out of the dislocation line. In these cases, the vacancy always returns to the dislocation line as opposed to diffusing away into the bulk since vacancy has a lower energy within the dislocation. These outliers were removed in order to study only the pipe diffusion mechanism. The overall mean escape time was 2.34 ns. For Jog A, the mean escape time was 2.31 ns ± 0.10 ns, and the mean escape time for Jog B was 2.39 ns ± 0.19 ns. Therefore, there is no discernible difference between pipe diffusion behavior in dislocations with Jog A versus Jog B.

If we assume that vacancy diffusion within the partial cores is a constant rate Poisson process where the energy barriers for pipe diffusion are constant and not significantly impacted by the presence of the jog, then the escape times should follow an exponential distribution,

F (t ) = 1 -e -t/τ , ( 2
)
where t is the escape time, F (t ) is the cumulative probability distribution, and τ is the average time between events. The vast majority of the 1843 vacancy escape times recorded fall below 15 ns. In this region, the plot is linear with a slope corresponding to τ = 3.14 ns, supporting the assumption that vacancy diffusion is a Poisson process. The long tail of escape times greater than 15 ns can be attributed to more rare events such as vacancy hopping in and out of the partial dislocation cores. These events will have a different average escape time and therefore do not fall on the linear fit. Additionally, there is a cluster of escape times less than 0.025 ns that also do not fall on the linear fit, which is attributed to the increased hopping rate for vacancies near the jog. Therefore, the linear portion with τ = 3.14 ns corresponds to pipe diffusion along the dislocation line independent of the presence of the jog.

To better quantify the impact of the jog on vacancy diffusion, we study how the distance from the jog affects the direction of vacancy flow. At each point along the dislocation line, the number of vacancies passing by in the positive and negative z directions is calculated. Figure 4 shows the fraction of these vacancies that are observed traveling in the positive z direction. For random vacancy diffusion, this fraction should be near 0.5, which can be seen for positions far from the jog. However, within about 10 Å of the jog, this value deviates from 0.5 as the jog biases the direction of diffusion. Because the jog attracts vacancies, vacancies are more likely to flow toward the jog. While vacancy pipe diffusion is random and vacancy hops in both directions occur on either side of the jog, the jog creates a diffusion flux for vacancies near the jog.

To estimate the free-energy difference between the jog and the displaced jog + vacancy state, we ran simulations starting with a jog only to observe vacancy emission, reaching a total of 38.5 μs of simulation time (Table I). The free-energy difference can be estimated by comparing the amount of time spent in each state:

G = k B T ln t jog t jog+vac , ( 3
)
where t jog is the simulation time spent in the jog only state and t jog+vac is the time spent in the displaced jog + vacancy state. To investigate possible entropy effects, the same simulations with a jog only were repeated at 800 and 1000 K. Table I shows that Jog B undergoes a larger free-energy change on emitting a vacancy at all three temperatures studied. The results do not show a strong temperature dependence for either jog configuration, suggesting that entropic effects are minimal.

To complement the ParSplice results, we also performed molecular statics calculations. First, vacancy-dislocation interaction energies were calculated for a variety of atomic sites in a straight edge dislocation and jogged dislocation. The interaction energy describes the difference in energy for a vacancy at a specific site relative to a vacancy in the bulk, so the vacancies will be attracted to areas with more negative interaction energies. For both dislocation types, the interaction energies outside of the dislocation correspond to the stress state around the dislocation, with tensile regions having positive interaction energies and compressive regions having negative interaction energies. Inside the dislocation cores, the interaction energies are much larger. In a straight edge dislocation, the interaction energy of a site within the core is -0.34 eV, while in the jogged dislocation, the interaction energy at a site within the core but far from the jog is -0.37 eV. This suggests that the jog has a small effect on vacancies at sites far from the jog. As was observed in the ParSplice simulations, the vacancy is attracted to the dislocation cores due to their negative interaction energies, and the vacancies FIG. 4. The fraction of vacancies hopping in the positive z direction as a function of distance along the dislocation line. The shaded region shows the standard error. Far from the jog, this fraction is about 0.5, so there is no directional preference for vacancy hopping. Closer to the jog, the vacancies are more likely to hop toward the jog than away, showing that the jog biases the vacancy diffusion direction.

were unlikely to hop out of the dislocation core due to the associated energy penalty. The highest magnitude interaction energies are found at sites directly adjacent to the jog. For Jog A, there are two equivalent sites on either side of the jog, each with an interaction energy of -0.58 eV. Between these two sites, there are four sites where, if the vacancy is placed there, it is absorbed and the jog displaces. These four sites correspond to the four atoms behind the turquoise atoms in Fig. 1(c). For Jog B, there are also two sites on either side of the jog with an interaction energy of -0.66 eV. However, due to symmetry, there are only three sites where the vacancy is absorbed, and these atoms are also found behind the row of turquoise atoms in Fig. 1(c). These strongly negative interaction energies suggest that the vacancy will be attracted to the jog. Overall, the calculated interaction energies show that vacancies are naturally attracted to the partial dislocations and jog.

The energy barriers for vacancy hops in and around the dislocation were calculated for a variety of positions in the system. Four distinct types of hops are identified and their energy barriers are shown in Fig. 5. The energy barrier for vacancy migration in the bulk is 1.29 eV, and all four of the vacancy events shown have much lower barriers so will therefore occur much more quickly than bulk diffusion. The first hop shows a vacancy diffusing along the length of the partial dislocation, referred to as pipe diffusion. The energy barrier for pipe diffusion is nearly symmetric with a barrier of approximately 0.86 eV. We find that the jog has little effect on the barrier for pipe diffusion, as the barrier in a jog-free edge dislocation was also calculated to be 0.86 eV. For a hop across the slip plane, the energy increases since the vacancy moves from a compressive region to a tensile region, which is consistent with the interaction energy calculations performed previously. The final two barriers show a vacancy being absorbed from either the bulk or the stacking fault. Both show that the total energy is lowered when the vacancy is in dislocation cores. Together, these results support the ParSplice observation of vacancy attraction to the partial dislocation core followed by pipe diffusion.

These calculated barriers can be compared to the ParSplice results by transition state theory, which states that the rate is related to the energy barrier by

= ν 0 e -E kT , ( 4
)
where is the rate, ν 0 is a frequency prefactor, and E is the energy barrier. Assuming a barrier of 0.86 eV for pipe diffusion and an average time between events of 3.14 ns, transition state theory predicts ν 0 to be 2.1 × 10 13 s -1 at 900 K. This FIG. 5. NEB-calculated energy barriers for vacancy hops in and around the partial dislocations. The arrow shows the direction of motion for the atom that replaces the vacancy. Compared to the bulk vacancy diffusion barrier of 1.29 eV, all of these vacancy hopping events have much lower barriers. The total energy is lowered when the vacancy is absorbed into the partial dislocation from either the stacking fault or the bulk, which is consistent with the interaction energy calculation. The vacancy also prefers to be on the compressive side of the slip plane, as shown by the barrier for movement across the slip plane. The barrier for pipe diffusion is approximately 0.86 eV, which suggests that pipe diffusion along the partial dislocations will be much faster than bulk diffusion. value of ν 0 is in good agreement with the typical estimation of ν 0 = 10 13 s -1 , showing that the assumption of a Poisson process is valid and consistent with the NEB-computed energy barrier of 0.86 eV.

The NEB method was also used to calculate the energy barrier for a vacancy being emitted from the jog. The energy profile associated with emitting a vacancy from each jog configuration is slightly different, as shown in Fig. 6. Since the two jogs have the same energy, both curves begin at 0 eV, but the ending energies are different due to the difference in vacancy interaction energies for sites near the jog as discussed previously. The barriers for vacancy emission are comparable to pipe diffusion, 0.779 eV for a vacancy emitted from Jog A and 0.916 eV for a vacancy emitted from Jog B. Based on these calculations, Jog A would be expected to emit vacancies more often than Jog B, which is consistent with the ParSplice emission rates shown in Table I.

The barriers for the reverse process, vacancy absorption, are relatively small. The absorption barrier for a vacancy being absorbed by Jog A is 0.120 and 0.053 eV for a vacancy being absorbed by Jog B. We can compare the average absorption time predicted by transition state theory to that calculated with ParSplice. Assuming ν 0 is 10 13 s -1 , the average absorption time is 0.470 ps for Jog A and 0.198 ps for Jog B at 900 K. This is consistent with ParSplice calculations as the vast majority of vacancies emitted are reabsorbed within one ParSplice segment (1 ps). Since the absorption time is on the scale of picoseconds while the vacancy hopping time is on the nanosecond scale, the rate-limiting step for dislocation climb will be vacancy diffusion to the jog, not absorption into the jog.

The NEB-calculated energy difference at 0 K can be compared to the free-energy difference calculated by ParSplice (Table I). The 0 K energy differences are 0.726 eV for Jog A and 0.796 eV for Jog B. At 900 K the ParSplice-calculated free-energy difference is only 0.092 eV higher for Jog A, and 0.088 eV higher for Jog B. These small differences can be attributed to entropic effects. The agreement between these two methods of calculation again supports the ParSplice results.

IV. CONCLUSION

Accelerated molecular dynamics and molecular statics were used to study a vacancy near a jogged edge dislocation in nickel in order to observe dislocation climb. Two different jog configurations with identical energies are formed due to the geometry of the lattice, and absorption of a vacancy causes the dislocation to climb and change jog configurations. We used ParSplice, an accelerated molecular dynamics method, to simulate a vacancy diffusing along the dislocation line to the jog. At 900 K, the mean time for pipe diffusion was 3.14 ns, which is consistent with transition state theory estimates using the NEB-calculated barriers of 0.86 eV for pipe diffusion. By comparison, the TST-estimated bulk diffusion is 1.67 μs.

Once the vacancy has diffused to the jog, it is absorbed very quickly, on the picosecond scale, due to the low-energy barriers for vacancy absorption. Therefore, assuming vacancies are present in the dislocation line, pipe diffusion will be the ratelimiting step. The jog can also emit vacancies, but these were always reabsorbed in our simulations due to the high attraction between the jog and the vacancy. The free-energy differences between the jog state and the vacancy plus displaced jog state were calculated at 800, 900, and 1000 K. The free-energy difference for Jog B emitting a vacancy is greater than Jog A, and the free-energy differences do not show a strong temperature dependence. The energy differences calculated with ParSplice are consistent with those calculated via molecular statics and strongly suggest that vacancy absorption to jogs is dominated by biased diffusion to the dislocation line followed by fast pipe diffusion to the jog. These results can be used to parametrize future mesoscale models to better understand dislocation climb in complex systems.

FIG. 1 .

 1 FIG. 1. (a) Schematic of the initial, unrelaxed dislocation system. A perfect edge dislocation with b = a 2 [110] and line direction N t + j (purple line) is introduced in a cylindrical simulation cell. The atoms in the gray region are then allowed to relax while all other atoms are held fixed. The simulation cell is aligned such that the periodic z direction is parallel to N t + j and the y direction is parallel to b = a 2 [110]. (b) Schematic of the relaxed dislocation structure in the small system. The dislocation relaxes into an edge dislocation with a jog. The edge components dissociate into partial dislocations (purple) with b 1 = a 6 [12 1] and b 2 = a 6 [211]. The stacking fault region between the partials is shown in red. The green dislocation lines represent the jog, which has a line direction of a 2 [1 10]. The extent of dissociation at the jog is much less than along the edge components, which has the effect of constricting the stacking fault region. (c) Comparison of the atomic configurations for the two different jog types. Atoms are colored by the centrosymmetry parameter with very low centrosymmetry atoms removed. The dark blue and green atoms surround the partial dislocations, and the yellow and orange atoms show the stacking fault regions. The two partial dislocations pinch together at the jog, which is shown by the red and turquoise atoms. Jog A is slightly more compact than Jog B, as can be seen by the turquoise atoms in the jog.

FIG. 2 .

 2 FIG.2. The total energy of the jogged edge dislocation as a function of the total dislocation length. Extrapolating the linear fit to a dislocation of length zero gives the energy of the jog only. The energy for both Jog A and Jog B is 3.52 eV.

FIG. 3 .

 3 FIG. 3. Distribution of vacancy hopping times. F(t) is the cumulative probability distribution. The majority of the data show a linear relation, showing that the main vacancy diffusion process follows a Poisson distribution. The slope of the plotted black line is -τ -1 , where τ is 3.14 ns.

FIG. 6 .

 6 FIG. 6. NEB-calculated energy profile of a vacancy being emitted from each jog configuration. The barriers for vacancy emission are 0.916 eV for Jog B and 0.779 eV for Jog A. Conversely, the barriers for vacancy absorption are relatively low at 0.120 eV for a vacancy absorbed into Jog A and 0.053 eV for a vacancy absorbed into Jog B.

TABLE I .

 I Vacancy emission from the jog. At all temperatures, Jog A emitting a vacancy has a lower free-energy difference than Jog B emitting a vacancy. This causes Jog A to have a significantly higher vacancy emission rate than Jog B.

	Process	Jog A → Jog B + Vacancy	Jog B → Jog A + Vacancy
	Temperature (K)	800	900	1000	800	900	1000
	Total simulation time (μs)	11.5	16.1	5.76	18.3	22.4	11.1
	Number of vacancies emitted	81	402	162	27	162	257
	Vacancy emission rate (μs -1 )	7.04	25.0	28.1	1.48	7.23	23.2
	Free-energy difference (eV)	0.814 ± 0.007	0.818 ± 0.004	0.828 ± 0.004	0.918 ± 0.011	0.884 ± 0.023	0.910 ± 0.007

ACKNOWLEDGMENTS

This work was partially supported by the Science Undergraduate Laboratory Internships program from the Department of Energy Office of Science. This work was also supported by the Illinois Scholars Undergraduate Research Program, and the Department of Energy National Nuclear Security Administration Stewardship Science Graduate Fellowship, which is provided under cooperative agreement number DE-NA0003960. Work at Los Alamos National Laboratory was supported by the Laboratory Directed Research and Development program of Los Alamos National Laboratory under Project No. 20150557ER. Los Alamos National Laboratory is operated by Triad National Security, LLC, for the National Nuclear Security Administration of U.S. Department of Energy (Contract No. 89233218CNA000001).

This research is part of the Blue Waters sustained-petascale computing project, which is supported by the National Science Foundation (awards OCI-0725070 and ACI-1238993) and the state of Illinois. Blue Waters is a joint effort of the University of Illinois at Urbana-Champaign and its National Center for Supercomputing Applications.