Multi-decadal trends in the tropical Pacific western boundary currents retrieved from historical hydrological observations

As large-scale ocean circulation is a key regulator in the redistribution of oceanic energy, evaluating the multi-decadal trends in the western Pacific Ocean circulation under global warming is essential for not only understanding the basic physical processes but also predicting future climate change in the western Pacific. Employing the hydrological observations of World Ocean Atlas 2018 (WOA18) from 1955 to 2017, this study calculated the geostrophic currents, volume transport and multi-decadal trends for the North Equatorial Current (NEC), the North Equatorial Countercurrent (NECC), the Mindanao Current (MC), the Kuroshio Current (KC) in the origin and the New Guinea Coastal Undercurrent (NGCUC) within tropical western Pacific Ocean over multi-decades. Furthermore, this study examined the contributions of temperature and salinity variations. The results showed significant strengthening trends in NEC, MC and NGCUC over the past six decades, which is mainly contributed by temperature variations and consistent with the tendency in the dynamic height pattern. Zonal wind stress averaged over the western Pacific Ocean in the same latitude of each current represents the decadal variation and multi-decadal trends in corresponding ocean currents, indicating that the trade wind forcing plays an important role in the decadal trend in the tropical western Pacific circulation. Uncertainties in the observed hydrological data and trends in the currents over the tropical western Pacific are also discussed. Given that the WOA18 dataset covers most of the historical hydrological sampling data for the tropical western Pacific, this paper provides important observational information on the multi-decadal trend of the large-scale ocean circulation in the western Pacific.


Introduction
The ocean is the largest heat storage in the earth's climate system (Chen and Tung, 2014;Cheng et al., 2017). Under warmer climate influenced by continuously emitted greenhouse gases, the energy budget of the climate system is out of balance due to additional heat obtained from solar radiation that more than 90% of which enters the ocean, leading to sustained ocean warming (Trenberth and Fasullo, 2013). Redistributing the heat into the ocean which the ocean dynamic process formulates the basic pattern of the ocean environment, mainly resulting from the large-scale ocean circulation that is a basic physical process of ocean energy redistribution (Ganachaud and Wunsch, 2000;Wu et al., 2012;Forget and Ferreira, 2019;Wu, 2020). Therefore, it is crucial to investigate how the large-scale ocean circulation varies over the past decades, i.e., a long-term trend, for understanding and predicting future climate variations. Given long-term climate warming, the global mean sea surface wind speed accelerates (Reguero et al., 2019;Young and Ribal, 2019), resulting in an accelerating trend in the global mean ocean circulation over the last 20 years (Hu et al., 2020). However, the response of ocean circulation to global climate variation is not evenly distributed, i.e., considerably varying with regional characteristics (Hu et al., 2015).
Intense atmospheric convection processes take place in the tropical western Pacific, the world's largest warm pool, which is the main source of typhoons and locale of El Niño-Southern Oscillation (ENSO) events. Furthermore, the tropical western Pacific exerts significant impacts on the climate of global tropics and China. The upper tropical Pacific Western Boundary Current (WBC) system consists of the Mindanao Current (MC), the Kuroshio Current (KC), the New Guinea Coastal Current/Undercurrent (NGCC/ NGCUC), the South Equatorial Current (SEC), the North Equatorial Current (NEC) and the North Equatorial Countercurrent (NECC) near the western boundary. The tropical Pacific WBCs, shuttling through the Indo-Pacific Warm Pool, play a role in the climate system (e.g., Guan et al., 2019), the variations of sea surface temperature (SST) in the region of the Western Pacific Warm Pool (Qu et al., 1997), the occurrence of the splitting of the western Pacific Warm Pool Wirasatriya et al., 2020), typhoon processes (Zhou et al., 2019), Asian monsoon (Huang and Li, 1988), ENSO (e.g., Jin, 1997) and other tropical climate processes (Hu et al., 2015).
Under the context of global warming, the tropical western Pacific has undergone notable changes and the Indo-Pacific Warm Pool has experienced expansion (Weller et al., 2016;, consequently, alterations in the western tropical Pacific WBCs and serious consequence on the Indo-Pacific climate are expected (Hu et al,. 2015). According to the assimilation data analysis over the past nearly 60 years from 1950 to 2008, the bifurcation of NEC in the tropical Pacific Ocean (10°-20°N) has a tendency of moving southward under the wind stress curl regulation, leading to the strengthened Kuroshio in the origin (Chen and Wu, 2012). Based on multiple numerical simulations under the background of global warming, the volume transports of NEC, MC and KC might decrease while the bifurcation of NEC might move northward (Duan et al., 2017). However, the multi-decadal trend of MC was not significant comparing with the reanalysis data due to not only a higher variability and uncertainties in the wind data but also the effects of multi-decadal variations of the WBCs in the tropical Pacific and the influence of Pacific Decadal oscillation (PDO) (Duan et al., 2019a). Synthesis of reanalysis data and model simulation indicates that the MC has been intensified for 1974-1982 and 2002-2008 while weakened for 1983-1990 and 1994-2001, resulting from tropical north-west Pacific winds forcing off the equator (Duan et al., 2019b). According to hydrological, assimilation and reanalysis data on the 137°E section, the volume transport of the SEC reached the maximum during 1980/1981, 1994/1995 and 2004/2005 while to the minimum during 1989/1990 and 1999/2000. Thereinto, the multi-decadal variation of the NEC is induced by local Ekman dynamic processes and propagation of the baroclinic Rossby wave caused by multi-decadal wind forcing in the tropical western Pacific (Zhai et al., 2013). The low-frequency variations and trends of the western Pacific sea surface height (SSH) and WBCs such as MC are also resulted from the propagation of baroclinic Rossby waves (Cheng et al., 2016;. The massive accumulation of satellite altimeter data and hydrological data over the 137°E section indicate that the rapid rising rate of SSH in the western tropical Pacific since the 1990s is 3 times faster than the global mean, along with a significant moving southward and accelerating trend of NEC and NECC, and the trend of accelerating Pacific Walker circulation is considered to be the main reason for the strengthening of NEC (Qiu and Chen, 2012). The accelerating NEC over the past 20 years has been corroborated by direct measurement. Shipboard acoustic Doppler current profiler (ADCP) measurements at the 137°E section show that the NEC increased by about 7 Sv (1 Sv≡10 6 m 3 s −1 ) from 1993 to 2008, accounting for 17% of the mean transport 71 Sv of the NEC during the period . Moreover, observations and reanalysis data indicate that the Indonesian Throughflow (ITF) and Makassar Throughflow (MTF), located in the downstream of the MC and NGCUC, have a significant strengthening trend, mainly arising from wind forcing and the effect of salinity associated with PDO (Liu et al., 2015;Hu and Sprintall, 2017;Gordon et al., 2019). Table 1 briefly summarizes key findings from the literature relevant to the trend in the western Pacific circulation. It is evident that there are contradictrory key findings and differences between the previous studies. It is questionable whether the tropical Pacific western boundary current system is strengthening or not, under the background of global warming. To answer this question accurately, detecting evidences from observational data is a crucial prerequisite. As shown in the literature, model simulation and reanalysis data have been used extensively, however, it is difficult to estimate the reliability of these data from models due to relatively high uncertainty in model-based data. Therefore, it is essential to investigate the historical observational data for better understanding the multi-decadal trend of the tropical Pacific WBC system. This study investigates the statistical characteristics and multidecadal trends of tropical Pacific WBCs using a hydrological dataset accumulated since the 1950s.

Data
This study employed the decadal mean historical temperature and salinity data (Garcia et al., 2019) released by the World Ocean Atlas 2018 (WOA18). The first decade with enough observational data to calculate the decadal mean climatological field is from 1955 to 1964. Therefore, WOA18 calculated the decadal mean climatological field every 10 years from 1955 to 2017 (e.g., 1955-1964, 1965-1974, 1975-1984, 1985-1994, 1995-2004, and 2005-2017) and accordingly there are 6 decadal climatological mean temperature and salinity fields (Garcia et al., 2019). Comparing with the WOA 2013 available until 2012, we found that the longer period of 2005-2017 as a decade in WOA2018 had minor influence on the results. Each mean field contains a large amount of observational data at a relatively finer spatial resolution, a horizontal resolution of 0.25°×0.25°and 102 levels in vertical 0-5500 m with 47 levels in the upper 1000 m and a mean depth of 54 m. Hence, WOA18 can commendably distinguish the temperature and salinity fields in the region of WBCs. Figure 1a and 1b present a spatial distribution of sampling stations that measure surface temperature and salinity in the tropical western Pacific and corresponding error estimates during the first decade from 1955 to 1964. It is prominent that there are still quite a few observation stations close to the western boundary and the number of salinity stations is apparently less than that of temperature. In the tropical western Pacific (120°E-180°, 10°S-20°N), the total number of temperature stations is 6090 while just 671 for salinity. Since the first decade, however, the number of observational stations for both temperature and salinity has been rapidly increased, exceeding more than 5000 stations. Currently, 12601 temperature observations and 8500 salinity observations are available, which are sufficient to estimate the decadal mean fields of temperature and salinity ( Figure 1c).
Systematic quality control was implemented for all gridded data (Garcia et al., 2019), including the correction of expendable bathythermograph (XBT) data (Cheng et al., 2014). The error distribution of salinity and temperature ( Figure 1) is relatively uniform with a significant error distribution for salinity less than 0.1 psu and temperature between 0.15 and 0.25°C. The error in the temperature and salinity measurements during each decade was generally distributed within this range, indicating that the data quality is acceptable to proceed with further analysis.

Methods
As previous studies have demonstrated that the main component of the WBC system in the tropical Pacific is geostrophic currents, it is an effective and widely adopted approach to study the structure and variations of the WBCs by retrieving geostrophic currents (e.g., Lindstrom et al., 1987;Hu and Cui, 1989;Qu et al., 2012). In this study, decadal mean temperature and salinity fields from WOA18 were used to calculate the geostrophic currents over the sections of 130°E, 18°N, 8°N and 147°E. Given two hydrological stations A and B, the velocity V between the isobaric surfaces P 1 and P 2 , which is perpendicular to the section formed by station A and station B, is where f is the Coriolis parameter, L is the horizontal distance between station A and B, the direction of V is perpendicular to the section AB, and α A and α B are the specific heat capacity, a function of ocean temperature (t), salinity (S) and pressure (p): S t p = ( , , ). (

2)
A A Assuming that the bottom layer of the temperature and salinity stations is a surface without flow, the velocity for each layer can be obtained by eq. (1). Whereas the different vertical distribution of temperature and salinity in the western boundary of the tropical Pacific, the deepest layers with no flow were chosen as a reference surface, i.e., 3245, 1819, 2937 and 1008 dbar (1 dbar=0.1 MPa) for sections 130°E, 18°N, 8°N and 147°E, respectively. Multi-decadal trend of MC is not significant Duan et al., 2019bIn 1974-1982and 2002, the MC intensified relatively In 1983-1990and 1994-2001, the MC weakened relatively Since the specific heat capacity is a function of temperature and salinity, referring to previous research methods, the contribution of variation in temperature and salinity can be separately evaluated by the geostrophic currents on the isobaric surface under the condition of controlling the temperature or the salinity variation (e.g., Feng et al., 2015;Hu and Sprintall, 2017). Thus, this study calculated the geostrophic currents variation from decadal means of temperature for each decade and the mean salinity averaged over the six decades as the influence of temperature variation on geostrophic currents while the influence of salinity variation was evaluated by the geostrophic currents variation from decadal means of salinity for each decade and the mean temperature over the six decades.   Sasaki et al., 2008). The distribution of boundary currents in the tropical Pacific is prominent in Figure 2c, with the NEC bifurcating as the KC and the MC to the east of east Philippines and converging with the NGCUC, then eventually entering the NECC. The mean structure of WBCs' sections in this study was consistent with previous studies. The subsurface undercurrent system consisting of the Luzon Undercurrent (LUC), the Mindanao Undercurrent (MUC), and the North Equatorial Undercurrent (NEUC) is clearly visible. The volume transport of each current was integrated based on the vertical structure of each current shown in Figure 2. Each current was defined as follows: the NEC to the west along the section of 130°E, 8°-20°N within 1500 m deep, the MC to the south along 8°N, 126.625°-129.875°E within 1500 m, the KC to the north along 18°N, 122.375°-129.875°E within 1500 m, the NECC to the east along 130°E, 3°-10°N within 500 m, and the NGCUC to the west along 147°E, 3°-5.5°S within 1000 m. Mean volume transports are 51±4 Sv, 31±5 Sv, 17±2 Sv, 36 Sv, 7 Sv for the NEC, the MC, the KC, the NECC, the NGCUC, respectively, which are consistent with previous studies (e.g., Qu et al., 1998).

Multi-decadal mean WBCs in the tropical Pacific
The total volume transport for the MC and the KC is 48 Sv, which is 3 Sv less than that of the NEC (only 13% of the total volume transport for the NEC). The difference is mainly derived from three reasons. Firstly, the sampling stations for MC and KC are not completely close to the shore, consequently, a missing part of the two western boundary currents is omitted and the region 18°-20°N is not closed. The second is that vertical exchanges occur in the regions of NEC-MC-KC. Lastly, unexpected biases may exist in numerical calculation. Based on the aforementioned reasons, the total volume transports for the MC and the KC can be regarded as in balance with the NEC. From the comparison above, the volume transports of WBCs in the tropical Pacific obtained from the geostrophic currents based on WOA18 hydrological observation data are reasonable. Over the past 60 years, the volume transports for the NEC, the NGCUC, and the NGCUC have increased by (13±6)%, (56±23)%, (18±15)%, respectively. The KC in the origin has also increased by (10±8)%. Especially, the increase of the During the intermediate four decades from 1965 to 2004, the five currents in the tropical western Pacific showed an intensifying trend. The NEC, MC and NGCUC represented a significant enhancing trend for the three time windows, indicating that the enhancing trend is significant regardless of the selection of the decade. The differences in different algorithms of the KC and the NECC also indicate that there is certain uncertainty in estimating the trend of the two currents.

Multi-decadal trend of WBCs in tropical Pacific
Transport variation associated with the temperature variation was remarkably greater than the contribution of salinity and the multi-decadal variation of each current was almost mainly dominated by the variation of the temperature. This result is consistent with the previous study of Cheng et al. (2008) concluding that the trend of the sea surface height in the Western Pacific Warm Pool is mainly caused by thermal variations. But there is a significant contrast to the region of the ITF, where previous studies implied that the contribution of salinity variation in the trend of ITF (Hu and Sprintall, 2017)

and the sea surface height of Indian Ocean Warm Pool is quite important.
A strong increasing trend over the multiple decades was shown in the NGCUC, a total increase of 4.6 Sv over the six decades. Although the MC and the NGCUC are the sources of the NECC, however, no significant trend was detected in the NECC. Such a result may be derived from the fact that the MC and the NGCUC are also the sources of the ITF, leading to the contribution of their increase to the long-term variation of the ITF. Despite the lack of reliable estimates of the multi-decadal trend of the ITF, some previous studies addressed a strengthening trend of ITF. For example, Tillinger and Gordon (2009) found that a strengthening trend during the period from 1958 to 2007 was shown in the Indo-Pacific pressure gradient. Therefore, the increase in the volume transports for the MC and the NGCUC may balance the increase of the ITF.
The 1990s has been considered as a turning point with regard to trends in wind and circulation in the tropical Pacific (e.g., Qiu and Chen, 2012;Hu and Hu, 2014). However, Figure 3 shows that the WBCs were strengthening overall from the mid-1980s to the mid-1990s in the tropical Pacific. Employing only the decadal mean time series investigated in this paper, it is difficult to verify the consequence. As shown in Figure 3, the trends in the early 1990s showed considerable multi-decadal variations rather than entirely long-term trends.

Dynamic height field
This study examined the dynamic height (DH, cm) with temperature and salinity based on WOA18, defined as in Eq. (3), to further analyze the large-scale background of the multi-decadal trend of the WBCs in the tropical Pacific.
where g is the gravitational acceleration and P R is the reference isobaric surface (1000 dbar in this study). The horizontal gradient of the DH is one of the indicators for the field of the geostrophic current. Figure 4 displays the multi-decadal DH trend on each section and in the tropical northwestern Pacific, and the comparison to each own corresponding mean state. In general, the distribution of the DH trend is similar to the distribution of the mean state of dynamic height except for the region of the NECC, indicating that a strengthening trend in the mean state exists.
Specifically, the DH of the nearshore side decreased in the origin of the KC while increased in the offshore side, resulting in strengthened eastward DH gradient (Figure 4a). In the region of the MC, DH in the nearshore was increasing while decreasing in the offshore, resulting in strengthened westward DH gradient (Figure 4c). In the region of the NEC, the south DH was weak and even decreased while the north DH was significantly strengthening, resulting in strengthening the meridional DH gradient (Figure 4b). In the region of the NGCUC, the DH trend was almost coincident with the DH, remarkably enhancing the meridional DH gradient. In the region of the NECC, the DH gradients from the DH trend are northward and southward in the south (Figure 4b south of 5°N) and the north (5°-8°N), respectively, and are close in amplitude and offset each other accordingly. As a result, the trend of the DH gradient in the NECC was not significant. The characteristics of the spatial distribution for the multidecadal DH trend are consistent with the calculated volume transport variations in the WBCs, which corroborates the preceding discussion on the trend of the transport of the WBCs. As the variation of DH is just one of aspects of the variation of the large-scale ocean circulation, rather than the cause, it is still necessary to further investigate the reason for the multi-decadal trend of the WBCs in the tropical Pacific.

The role of wind field
Given that these currents are wind-driven circulations and the wind stress is the dominant forcing factor, it can be speculated that the strengthening trend of WBCs in the tropical Pacific over the past 60 years is caused by the variation of the wind field. Hu et al. (2020) found that the global mean tropical wind field has significantly increased since the 1950s, especially a more rapidly increasing trend over the latest 20 years. Strengthening kinetic energy input from the wind field may induce the accelerating trend of global mean ocean circulation, especially in tropical regions.
The NEC and SEC driven by trade winds converge in the western Pacific and dominate the general variation of the WBCs. Therefore, the variations of ocean circulations in the western Pacific are mainly dominated by the trade winds. For example, Duan et al. (2019a) suggested that the decadal variations of the MC is dominated by the wind fied in the region of the western tropical Pacific, with the mean zonal wind in the region (10°S-10°N, 120°-160°E) corresponding significantly with interannual and decadal variation of volume transport of the MC, indicating that variations of trade winds play an important role in the variation of the WBCs at low latitudes.
Using the ocean reanalysis dataset (ORA S3) from the European Centre for Medium-range Weather Forecasts (ECMWF), this study further analyzed the surface wind stress to investigate the role of wind field variations in the multi-decadal variations of ocean circulations in the tropical western Pacific. A hypothesis consistent with the previous studies is that the WBCs are dominated by the wind field at the same latitude Duan et al., 2019a) and the bifurcation latitude of the NEC affects the KC and the MC. Given the hypothesis, this study calculated the regional mean zonal wind stress for the NEC (8°-20°N, 120°-160°E), the KC (11°-20°N, 125°-160°E), the MC (0°-11°N, 130°-160°E), the NGCUC (0°-8°S, 140°-180°E), and the NECC (18°S-11°N, 130°E-180°), in which the north and south boundaries are the latitudes of the bifurcation point of the NEC and SEC respectively (Zhai et al., 2014).
It is evident that there is a strong correlation between the regional mean zonal wind stress and currents in the tropical western Pacific. To compare the regional mean zonal wind stress with the decadal mean transport of each current shown in Figure 5, the transport of each current needs to be multiplied by a coefficient α varying with each other for comparability (details in Figure 5). The regional mean zonal wind stress accurately describes the variation of the corresponding ocean currents, especially for the NEC, the MC and the NGCUC, indicating that the modulation of wind field lead to the decadal variations in the tropical western Pacific Ocean circulation. This study found a strengthening trend in the mean zonal wind stress in each region, where the wind field growth rates for the NEC, the MC, the KC, the NECC, the NGCUC are 4.3%, 9.4%, 3.8%, 7.7%, 7.9% per decade, respectively, indicating the trade winds in the tropical western Pacific have a strengthening trend. Figure 6 shows the spatial distribution of decadal mean zonal wind stress anomalies (minus the mean of 1959-2011). Both the decadal variation and the multi-decadal trend of the zonal wind stress are strong and significant in the region of northwest tropical Pacific, with a strengthening westward trend in the region of the NEC, which is consistent with the trend of volume transport and regional mean wind stress as shown in Figure 5. In the subtropical regions (north of 20°N), the decadal variation of zonal wind stress appears to be independent of that of tropical seas, which may explain the distinct difference between the KC and regional mean wind field in Figure 5.

Conclusion
This study employed the decadal mean hydrological observational data from the WOA18 dataset to calculate the geostrophic currents for the NEC, NECC, MC, KC in the origin and the NGCUC in the tropical Pacific. We also investigated the multi-decadal trend for each current. Furthermore, the contribution of temperature and salinity to the trend of each current was evaluated separately.
Over 63 years from 1955 to 2017, all of the currents except for the NECC showed a significant long-term strengthening trend that greater than decadal variability while the multidecadal trend for the NECC was not significant compared to the decadal variability. The sum of the multi-decadal trend of the MC and the KC in the origin was roughly consistent with that of the upstream NEC. By separating the contributions of variations of temperature and salinity to the variation of each current, this study found that the temperature variation mainly dominated the multi-decadal variation of WBCs in the tro-pical western Pacific while the contribution of the salinity variation was not significant. Analyzing the various modes of the tropical western Pacific relative to the 1000 dbar DH over 1955 to 2017, we found that the two trends of either the DH in each section of the boundary currents or the horizontal special DH were consistent with the mean state, indicating a strengthening trend of the WBCs in the tropical Pacific.
This study also investigated the important role of wind forcing in the decadal variations and trends of ocean circulations in the tropical western Pacific. The result showed a prominent correlation between the zonal wind stress and the transport of relative currents in the pivotal regions of the tropical western Pacific. It indicates that the volume transport of each current flow can be accurately captured by the areal mean zonal wind stress. During the period from the 1950s to the 2010s, a significant strengthening trend in the zonal wind stress led to a multi-decadal strengthening trend in ocean circulations in the tropical western Pacific. This study utilized areal mean zonal wind stress as an indicator that reflects the transport variation of each current, which may have important implications for understanding and predicting future variations of ocean circulation in the future.
There are still considerable uncertainties in calculated results in this study related to the WBCs in the tropical Pacific. Firstly, the deficiencies in observation data, e.g., the lack of regional data near the western boundary and the location near the core of WBCs within the western boundary region, the inhomogeneous distribution of hydrological observed data, the errors caused by the change of observation instruments, the errors from the XBT temperature observation data, and so on. Secondly, the WOA18 dataset artificially truncates each decade to obtain the decadal mean temperature and salinity field. In other words, the data at each station contains only 6 points corresponding to the six decades in this study, leading to the weak freedom of time series and an unreliable trend test accordingly. In addition, the salinity variation might be not captured precisely due to the lack of salinity observed data in the early days. Consequently, there may be Figure 5 Zonal wind stress comparing with the volume transport of currents in the tropical west Pacific. The red lines are the regional monthly mean zonal wind stress, blue box lines are the multi-decadal mean zonal wind stress anomalies, and the red dot lines are the volume transport of each current multiplied by a coefficient α for each current; 0.72×10 −3 for the NEC (a), 0.58×10 −3 for the MC (b), −3×10 −3 for the KC(c), −0.63×10 −3 for NECC (d) and 1.7×10 −3 for NGCUC (e). The first decade of the wind field is specially marked with a light color as the wind data for the first decade is available for only 5 years from 1959 to 1964. The last decade is neglected as the wind field spans from 2005 to 2011, particularly differing from the current field that spans from 2005 to 2017. some defects in the discussion on the effect of salinity variation on the transport variation of the WBCs in the tropical Pacific Ocean in this study.
Despite the considerable uncertainty, however, the WOA18 dataset covers most of the historical observational data on the hydrological environment over the tropical western Pacific and adequately represents the present quality of observation in the region. Therefore, the WOA18 data may provide the most accurate estimation of variations for the real ocean. Although we also examined other datasets of the current field from reanalysis or model-based simulations, e.g., ECMWF ORA S5 and OFES, those datasets were not adopted in this study as we found notable defects in reproducing the observed circulation structure and variations. It is a crucial step to accurately estimate the low-frequency variations and long-term trends for the large-scale ocean circulations for understanding the variabilities of ocean under the background of global warming. However, this study suggested limited evidence due to the aforementioned uncertainties in the dataset. It is expected that the academe will improve the long-term continuous observations for the ocean circulations not only to accurately evaluate the multi-decadal variations in the future and trends of large-scale ocean circulation but also further to investigate the internal physical mechanisms.
Open Access Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.