Seismic risk assessment for the infrastructure in the regions adjacent to the Russian Federation Baikal–Amur Mainline based on the Unified Scaling Law for Earthquakes

The purpose of the analysis is to assess the risk of loss in performance of infrastructure facilities in the regions adjacent to the Russian Federation Baikal–Amur and Trans-Siberian Mainlines due to seismic events of maximum macroseismic intensity expected in a period of 50 years with a probability of 10%, 5% and 1%. In particular, we use earthquake data compiled at the Baikal Division of the Russian Geophysical Survey, which provides sufficiently complete earthquake determinations of M = 2.5 or larger events for the period 1994–2019 for a reliable mapping the Unified Scaling Law for Earthquakes coefficients at the seismically active cells of a regular grid. Based on these estimates, we present maps of the maximum expected magnitude in about 500, 1000 and 5000 years. Having described an anisotropic seismic effect model of seismic sites in the region, we conclude by characterizing the seismic hazard in traditional terms of macroseismic intensity and by estimating the associated seismic risk to selected infrastructures, i.e. Baikal–Amur Mainline and Trans-Siberian Railway.


Introduction
Seismic hazard and risk assessment requires an adequate understanding of the earthquake distribution in magnitude, space and time ranges. Lacking data for a period of several 1000 years makes a probabilistic approach to estimating the recurrence time of hazardous ground shaking unreliable and, as a matter of fact, misleading (Heidarzadeh and Kijko 1 3 2011; Kossobokov and Nekrasova 2012;Wyss et al. 2012). We highlight the importance of adequate rescaling of the number of earthquakes of different sizes when you zoom or expand the analysis to a smaller or larger volume, which is the critical issue for reliable seismic hazard assessment.
Neo-Deterministic Seismic Hazard Assessment (NDSHA, Fäh et al. 1993;Panza et al. 2001), which dates back to the turn of the millennium, provides an alternative to the Probabilistic Seismic Hazard Analysis (PSHA) approach. It represents the innovative scenariophysics-based multidisciplinary approach for the evaluation of seismic hazard, proven reliable by 20 years of experiments in many countries worldwide (Panza et al. 2021).
Following Bak et al. (2002) the generalization of the Gutenberg-Richter recurrence relation that accounts for naturally fractal properties of earthquake origination and distribution in space is denoted as Unified Scaling Law for Earthquakes (USLE). By applying USLE to evident patterns of distributed seismic activity one can demonstrate that the traditional estimations of seismic risk for cities and urban agglomerations are usually underestimated (Kosobokov and Mazhkenov, 1994;Kossobokov and Nekrasova, 2018a, b;Parvez et al., 2018). In fact, any estimate of seismic hazard rate (e.g. the number of earthquakes of magnitude M) depends on the size of the territory that is used for averaging and, therefore, may differ dramatically when scaled down to the proportion of the area of interest. For example, one of the very first practical conclusions drawn by Mazhkenov (1988, 1994) demonstrated that in the case of the Lake Baikal Region with an area of 1,500,000 km 2 and the fractal dimension of the set of epicentres C = 1.25, the inclusion of aseismic areas leads to underestimation of seismic activity in an area of 1000 km 2 by a factor of 15, and to its overestimation by a factor greater than 2 when a characteristic of seismic activity over an area of 1000 km 2 is computed for a grid cell of 10 km × 10 km. The up-to-date modification of the algorithm of the USLE coefficient estimations contributes to the new design of the seismic hazard and risk evaluation in the frames of the Neo-Deterministic Seismic Hazard Assessment (NDSHA).
Transportation infrastructures are complex systems of various connected components like bridges, roads, railroads and tunnels. Due to their spatial extent, they are exposed to various kinds of natural hazards such as earthquakes, tsunami, landslides, liquefaction, flooding, etc. (Adams and Heidarzadeh 2021;Zahrai and Heidarzadeh 2007;Tsuji et al. 2011). Experience from past disastrous events shows that transportation infrastructures are quite vulnerable due to the lack of redundancy, the lengthy repair time, the rerouting difficulties or the cascading failures and interdependencies. The suffered of railroad damage cites variance from temporary speed restrictions to heavy damage with repairing or failure requiring new construction or major rehabilitation (Koks et al. 2019). In case of spatially distributed network systems under natural hazards, it is needed to be concerning that vulnerability is determined by not only the hazard intensity but also the geography, topology and network flow, which includes the infrastructure itself (Hu et al. 2019). Therefore, in terms of resilience it is important to recognize and quantify the risks and global losses associated with the potential damage of transportation systems and to establish efficient risk mitigation strategies. These include, among others, enhancement of emergency preparedness, strengthening of existing structures and improvements to recovery planning.
In this paper we describe the methodology of multiscale analysis of seismic recurrence along with the algorithm for estimating scaling coefficients of USLE, and then apply it to evaluating seismic hazard and risks to regions adjacent to the Russian Federation Baikal-Amur Mainline (BAM) and the Trans-Siberian Railway (TSR). The preliminary evaluation of seismic hazard for 10-km segments of BAM and TSR and 1-km segments of the major regional pipeline in traditional terms of maximum expected magnitude in about 50, 1000 and 5000 years based on USLE, pattern recognition of earthquake-prone areas and an isotropic model of macroseismic intensity distribution is presented in (Kossobokov and Nekrasova 2021). In the following we provide a revision taking into account the anisotropic propagation model and more careful sampling of parameters and data in estimating USLE coefficients.
Finally, we evaluate the difference in USLE seismic hazard assessment maps and PSHA maps for the territory of the Lake Baikal region. Specifically, the up-to-date official General Seismic Zoning maps (version GSZ-2016) are used for comparison. The seismic hazard maps for building codes and regulations beginning with the 1997 official maps (version GSZ-97), which contributed to the Global Seismic Hazard Assessment Program (GSHAP) (Giardini 1999;Ulomov and Shumilina 1999) and its more recent generations GSZ-2010, GSZ-2012, GSZ-2014, GSZ-2015, GSZ-2016 based on the same methodology (Ulomov 2007) were prepared and modified by a large group of Russian scientists. In particular, the comparison of the two most recent versions, as well as a general description of the current methodology, mathematical and computational apparatus, is given in (Zavyalov et al. 2019).

Methodology
The multiscale analysis involved in evaluation of the USLE coefficients at a given site combines recurrences for earthquakes obtained from the enclosed areas of different size, in order to get enough statistics on larger magnitude shocks from larger surrounding territories for a reliable confident estimation of scaling.
To obtain a seismic risk evaluation based on the concept of the USLE methodology, the following steps are performed: (a) characterization of seismic sources; (b) regional seismic effect modelling; (c) providing seismic hazard maps; (d) obtaining seismic risk estimations for infrastructure based on seismic hazard maps and the end-user hazard mitigation strategy. The following is a description of the algorithm for each of the above steps.
A preliminary data examination should be provided. It takes into account the representative nature and temporal completeness of the earthquake catalogue data, regional patterns of the distribution of active faults, parameters of the algorithm for determining the USLE coefficients and options for processing the system of active faults.

Characterization of seismic sources
The regional seismic sources are characterized by a finite set of earthquake-prone cells {c i } of linear dimension L with expected annual number of earthquakes N(M, L). Nekrasova et al. 1 3 (2015) elaborately described the Scaling Coefficients Estimation (SCE) algorithm (Nekrasova and Kossobokov 2002), which is a modified version of its prototype (Kossobokov and Mazhkenov 1988). For each c i , with reliable USLE coefficients estimation of N(M, L) can be obtained as where A-corresponds to the logarithmic estimate of seismic activity at magnitude M 0 , normalized to a unit area of 1° × 1° and a unit time of one year, B-the coefficient of magnitude balance, analogous to the b-value of the classical Gutenberg-Richter relationship (Gutenberg and Richter 1954), and C-the fractal dimension of the carrier of earthquake epicentres.
The USLE coefficients are used for estimation of the maximum magnitude Mmax expected with a p% chance of exceedance in T years following the procedures suggested in (Parvez et al. 2014). Specifically, for each c i we apply formula (1) to calculate the expected numbers of events from magnitude ranges M k in T-years and then find the magnitude The achieved values of Mmax at the entire set of grid points {c i } are used to compile a set of pairs {(c i , M i )} for a design of the expected seismic effect in the region considered in terms of macroseismic intensity.

Anisotropic propagation model
To describe the seismic effect from each of the earthquake-prone cells c i , we use a model of anisotropic distribution of the macroseismic intensity outside the earthquake source zone (Nekrasova and Kossobokov 2022b). In particular, for each pair (c i , M i ) we estimate an elliptical earthquake source zone with semi-axes A(M) = 1/2 × 10 α+βM and B(M) = 1/2 × 10 γ+δM where M = M i and α, β, γ, δ are the constants characterizing typical length and width of the source zone. The constants should be regional, if available, or determined by independent studies in tectonically similar region elsewhere, e.g. (Wells and Coppersmith 1994). The direction of the elliptical shape angle φ is measured from the dominant strike ψ of active faults in the Earth's crust system, defined for the site around c i .
Following Shebalin (1968) we use here the empirical estimate of macroseismic intensity I at distance Δ from of an earthquake epicentre of magnitude M originated at depth h: where b, v and c are the empirically estimated regional constants.
We assume that seismic waves propagate uniformly from the boundary of the earthquake source zone, so that outside it macroseismic intensity follows the equation below: where Δr(M, φ) is the minimum distance from the point (Δ × cosφ, Δ × sinφ) to the boundary of the source zone, namely, to the ellipse with semi-axes A(M i ) and B(M i ) centred at c i , while within the earthquake source zone: 0, h). Note that to comply with a qualitative characterization of macroseismic intensity the value of I e is rounded to the nearest semi-integer, e.g. 7.36 to VII-VIII, 10.12 to X, etc.
We determine the semi-axis A(M i ) as the large one which holds for magnitudes M ≥ 4.9.

Providing seismic hazard map(s)
The procedure of macroseismic intensity zoning is performed as a definition of boundaries created by merging all areas with the same semi-integer value of macroseismic intensity propagated from seismic source zones with Mmax. The macroseismic intensity zones from one source are determined using Eqs. (3, 4) described in the previous paragraph.

Providing the seismic risk for infrastructure objects
Any kind of risk estimates results from a convolution of the hazard with the exposed object under consideration along with its vulnerability is natural hazard at location g, O(g) is the exposure of objects at risk at g, and V(O) is the vulnerability of objects at risk. The location g, convolution operator ⊗ and V(O) concerned depend on each certain risk estimation task determined by the end-user.

Data
We described in detail the USLE coefficient estimation based on the actual Lake Baikal region data in (Nekrasova and Kossobokov 2022a). Specifically, the original online version of the Baikal Division of the Geophysical Survey, Federal Research Center of the Russian Academy of Sciences catalogue data BDRGS (2020) for magnitudes equal to or more than 2.6 (energy class K ≥ 8.6, accepted in catalogue homepage) for the period 1994-2019 and within 48-58°N and 99-122°E were used for compiling the maps of the USLE coefficients. The calculations by the SCE algorithm (Nekrasova et al. 2015) performed using the hierarchy of enclosed square cells with linear sizes of 2°, 1°, 1/2°, 1/4° and 1/8° for 1813 earthquake-prone cells of a regular grid with an annual empirical probability value of seismic events of magnitude 2.6 (K = 8 or larger greater than 0.0001 and the modified Eq. (1) normalized to recurrence of strong magnitude 6.0 earthquakes, i.e. M 0 = 6.0. With the same space-time-magnitude volume as in (Nekrasova and Kossobokov 2022a) we also modified the number of intervals into which the magnitude range is subdivided. Using a magnitude step of 0.32 as the SCE algorithm parameter, we obtained reliable USLE coefficient estimates for 1715 cells, while the previously used subdivision of the magnitude scale into equal quarters of a unit allowed us to map only 1460 cells and yielded a clear "white spot" (lack) of the USLE determination values in the central part of Lake Baikal. The modified fragmentation of the magnitude MLH scale helps avoiding the bias originated by the two-step round off on conversion of the above-mentioned initial determination of the energy class K in the regional catalogue to the magnitude MLH of an earthquake. Figure 1 shows the empirical probability density distribution functions (p.d.d.f.'s) of the A, B and C coefficients obtained using the adjusted SCE parameters, as well as p.d.d.f.'s of their errors of determination σ A , σ B and σ c . The σ-values are under 0.05 for coefficients B and C and less than 0.08 for coefficient A.

3
The grey areas in Fig. 1 mark the regional empirical density distributions of the USLE coefficients estimated based on the data from 1965 to 1998 and a more robust hierarchy of square boxes with linear size from 8° down to 1/2° (Nekrasova and Kossobokov 2006). The comparison shows that distributions of A and C based on different hierarchy are within the same ranges and not affected by the magnitude steps parameter changes, while the B coefficient distribution might be strongly affected by using dM = 0.32 SCE parameter or by a change of seismic regime in the Lake Baikal region. The single sharp maximum distribution at 0.9 characterize the actual distribution of the B coefficient values in 1994-2019. Figure 2 shows the spatial distributions of the A, B and C coefficients. Figure 2A shows that the Baikal rift can be characterized as a passage from relatively low values of A in the western part (from one event with M = 6 in 1000 years), transient values of A ranges from − 2 to − 2.5 in the central part and maximum values of A (up to − 1.5, which corresponds to three earthquakes with magnitude 6.0 in 100 years) in the western part. The "white spot" zone of the SCE determinations in (Nekrasova and Kossobokov 2022a) with the modified dM = 0.32 is characterized with the values of the A coefficient about − 3 (Fig. 3A), the values of B about 0.9, which widespread from the Eastern Sayan Mountains to the Muya Ranges (Fig. 2B), and relatively high values from 1 to 1.2 of C (Fig. 2C).
We use the "Dominant directions of the local active fault system estimation" (DDLAFS) plugin (Emelyanov and Nekrasova 2022) to obtain the set of azimuths ψ 1 , ψ 2 and ψ 3 with the empirical probabilities p 1 , p 2 and p 3 for the dominant directions of the Lake Baikal Region faults reported in Active Faults of Eurasia Database (Bachmanov et al. 2017; http:// neotec. ginras. ru/ datab ase. html). Figure 3 shows approximately 2500 active faults in the Lake Baikal Region.
We used free vector and raster map data at 1:10,000,000 available from Natural Earth (https:// www. natur alear thdata. com/ downl oads/ 10m-cultu ral-vecto rs/ railr oads/) for tracing  (Nekrasova and Kossobokov 2006) from a robust hierarchy based on the data from 1965 to 1998 are marked by grey shadows Fig. 2 The spatial distributions of the USLE regional coefficients A, B and C the two mainlines BAM and TSR under study. Here we consider mainlines only; the railway network parts of the second or less relevance factor, according to the database, are not included in our analysis.

Seismic hazard maps for the Lake Baikal Region
We present the seismic hazard maps for Lake Baikal Region in terms of Mmax at earthquake-prone cells of the regular 1/8° × 1/8° grid in 50 years (Fig. 4). Specifically, we calculate the expected number of earthquakes from the magnitude range Mj in 50 years, and then find the maximum magnitude, Mmax, with the expected number 50 × N(Mmax, L) surpassing the selected chance of exceedance (in %). The objective of the analysis is to evaluate probability of damaging seismic events for the purposes of assessing loss in the performance of the major railway networks. Figure 4 shows the seismic hazard maps obtained by USLE approach (Nekrasova et al. 2015(Nekrasova et al. , 2020 for the Lake Baikal Region in terms of Mmax, at earthquake-prone cells of the regular 1/8° × 1/8° grid in 50 years with 10%, 5% and 1% chance of exceedance. Figure 5 shows the number of the regular grid cells N versus Mmax. At the level of 10% probability of exceedance 1134 of the total 1715 earthquake-prone cells featured by Mmax 5.0 or more and 269 cells by Mmax 6.0 or more, i.e. about 66.2% and 15.7% of the regional total of earthquake-prone cells, correspondingly. According to Fig. 4 the maximum values Fig. 3 Active faults of the Lake Baikal Region and strong earthquakes from (Ulomov and Medvedeva 2014) of Mmax ≥ 6.0 at 10% of exceedance extend from the northern part of the Lake Baikal rift up to the area of the Ml 7.6 Muya earthquake on June 27, 1957. The area of South Yakutia, located on the north-eastern edge of the regional seismic network and the study region, is another high-risk place further to the north-east. Similar pattern of the most seismically hazardous areas is observed in Fig. 4 at 5% and 1% probability of exceedance in 50 years. In particular, as follows from Fig. 5 for 5% of exceedance in 50 years there are about 0.9% of the cells with Mmax ≥ 6.6, 37.4% and 84.4%-with Mmax ≥ 6.0 and 5.0, respectively, and, therefore, only 15.6%-with Mmax < 5.0. Moreover, for 1% of exceedance in 50 years there are almost 7.8% of cells with Mmax ≥ 7.2, 71.0%-with 6.0 ≤ Mmax < 7.2, 21.2%with Mmax < 6.0, and just 0.6%-with Mmax < 5.0.
We prepared the set of active faults dominant directions {ψ i } at each of 1715 cells of reliable determination of the USLE coefficients. The direction ψ i is defined as the azimuth of the maximum empirical probability determined by making use of the number of active fault azimuths in non-overlapping 10-degree sectors of a 30-km radius circle centred at c i . Figure 6 demonstrates the {ψ i } directions by 60 km lines centred at {c i } and colour-coded in respect to Mmax with 10% of exceedance in 50 years shown in Fig. 4.
The set of triplets {(c i , M i , ψ i )} allows us to design the elliptical seismic source zones for the Lake Baikal Region and apply Eqs. (3, 4) with the regional constants b = 1.5, ν = 4.0, c = 4.0 to comply with the Russian building codes and recommendations (GOST 2017, Fig. 4 The seismic hazard maps for the Lake Baikal Region in terms of Mmax expected at earthquakeprone cells of the regular 1/8° × 1/8° grid in 50 years with 10%, 5% and 1% chance of exceedance Fig. 5 The number N of the regular grid earthquake-prone cells 1/8° × 1/8° in the Lake Baikal region with 10%, 5% and 1% probability of exceedance of Mmax in 50 years Appendix K, Table K). Note the results presented here are for demonstration purposes only and do not account for some very important issues of reliable seismic hazard and risk assessment for specific infrastructure facilities. We estimate macroseismic intensity fields for intensity classes VI, VII and VIII or larger and use the average coefficients α = − 2.29, β = 0.57, γ = − 1.17, δ = 0.34 from (Wells and Coppersmith 1994) to describe the elliptical seismic source zones of the expected potentially destructive earthquakes of magnitude 5.0 or larger. This resulted in considering 1136, 1448 and 1705 cells with expected 10%, 5% and 1% of exceedance of Mmax in 50 years. Figure 7 demonstrates the USLE-based seismic hazard maps for 10%, 5% and 1% probability of exceedance in 50 years in terms of macroseismic intensity VI, VII and VIII or larger.

Seismic risk for infrastructure exposures adjacent to the Baikal-Amur Mainline and Trans-Siberian Railway
We used arbitrary units of risk proportional to the macroseismic intensity to evaluate the loss in performance of infrastructure in the regions adjacent to BAM and TSR within the USLE-based seismic hazard maps (Fig. 7). Figure 8 shows the model seismic risk obtained for each 10-km segment of BAM and TSR, which statistics are summarized in Table 1.

3
The risk evaluation results show that BAM can have risk unit of VIII or larger for each of the obtained seismic hazard maps, while risk that large for TSR appears only in case of 1% chance of exceedance in 50 years (Fig. 8). Roughly the same number (84, 99 and 100) of TSR segments is at risk of VI units in each of the three seismic hazard maps, while the risk of VII units increases dramatically from 9 to 58, and then to 140 10-km segments for the 10%, 5% and 1% probability of exceedance in 50 years, correspondingly.
We do not discuss in additional detail the risk that needs to take into account the soil conditions at each site as well as the construction located along each potentially affected segment, showing only the preliminary very crude estimate of risk value proportional to the expected maximum macroseismic intensity.

Discussion
We use the rough digitized copy of the latest version of the Russian General Seismic Zoning (GSZ-2016) to compare it with the USLE-based seismic hazard maps in terms of macroseismic intensity for the Lake Baikal Region. The GSZ-2016 maps of three macroseismic intensity classes VI, VII and VIII or larger for 50 years and 10%, 5% and 1% probability of exceedance are presented in Fig. 9. The isolines of seismic hazard assessment based on the USLE approach for the same classes of macroseismic intensity and the same levels of probability are added on top (solid lines for VIII or larger, dashed lines for VII and thin solid lines for VI). The comparison reveals more accurate delineation of hazardous zones obtained by the USLE approach as an alternative to slightly modified traditional PSHA of GSZ-2016. Table 2 sums up a comparison of the area occupied by intensities VI + , VII + and VIII + ranges on the USLE-based approach and GSZ-2016 maps for the Lake Baikal region within the boundaries of Russian Federation. The percentage are given for the total area of 1,540,199 km 2 all assigned to the VI + intensity zone on the GSZ-2016/C map. Note that VIII + class includes the highest classes of the macroseismic intensity (IX-XII). One should not be surprized that in GSZ-2016 underestimation of seismic hazard for cities and urban agglomerations (Kosobokov and Mazhkenov 1994;Kossobokov and Nekrasova 2018a, b;Parvez et al. 2018) is accompanied with apparent overestimation of the hazardous area in the entire region considered.
We verified the seismic hazard maps with the data from the unified earthquake catalogue compiled by Ulomov and Medvedeva (2014) for general seismic zoning of the territory of  3) earthquake at some 10 km to the north of Baykal'sk is the last one reported. All but one epicentre (i.e. the one attributed to the MLH 6.5 earthquake in May 1827 located at 58°N and 108.5°E) fall inside the area of expected strong ground shaking on the territory of the Russian Federation (Fig. 10). Specifically, 95 earthquakes (77% of the total 124 earthquakes) are located in the area of expected intensity VIII, 22 (18%)-in area of intensity VII, and only 6 (5%) within the area of intensity VI. Finally, it should be noted that our methodology has a potential of application to other kinds of natural hazards, such as volcano eruption induced earthquakes and/or tsunamis like in the two recent cases of the 2018 Anak Krakatau and 2022 Tonga-Hunga Haʻapai eruptions (Heidarzadeh et al. 2020(Heidarzadeh et al. , 2022.

Conclusions
Using the USLE approach (Nekrasova et al. 2015;Kossobokov and Nekrasova 2018a, b, present study) can provide a good approximation of the seismic hazard map in terms of expected maximum magnitude earthquake, macroseismic intensity, peak ground acceleration and other parameters of ground shaking for the selected time interval and fixed Fig. 10 Epicentres of the magnitude MLH ≥ 5.5 earthquakes from Ulomov and Medvedeva (2014) on top the USLE-based seismic hazard map of Imax with 1% chance of exceedance in 50 years given in Fig. 7. Epicentres within the areas of macroseismic intensity VIII, VII and VI are marked with black, blue and green stars, respectively. Three epicentres fall outside the area of expected strong ground shaking (yellow stars) including the two outside the territory of the Russian Federation probability of exceedance level. The seismic risk for the different types of object with the point, linear or planar spatial shapes could be estimated for preliminary basic consideration and better understanding by decision makers.
The USLE estimations (although limited by seismic data available) addressing realistic and practical kinds of seismic risks should involve professional interpretation by experts in earthquake engineering, social sciences and economics, which is not yet done. The general levels of underestimation of seismic effect rates at exposures, as well as an overall overestimation at regional scales, are too large to be ignored in seismic risk and earthquake loss evaluations necessary for knowledgeable disaster prevention and mitigation.
The USLE suggests a possibility to improve exposure of different objects at risk by relocating them to a safer place. Vulnerability of an object at risk could be reduced by a knowledgeable choice of appropriate hazard reduction measures, either in advance of their design or in the case of an earthquake prediction alert, after simulating a comprehensive collection of earthquake scenarios of realistically distributed sizes, mechanisms and locations.
Author contributions AN performed preliminary data analysis, USLE coefficients computations, active faults system analysis and SH maps calculations. Both AN and VK collaborated on data interpretation and manuscript write up.
Funding VK and AN carried out their parts of the presented research within the framework of the Russian State Task of Scientific Research Works on "Seismic hazard assessment, development and testing of earthquake prediction methods" (No. 0143-2019-0006) IEPT RAS and AN carried out the presented research with additional support of the State Task of Scientific Research Works of IPE RAS.
Data availability statement Data sets from public sources available in January 2022 were analysed in this study. The earthquake datasets used in this study are published online by Baikal Division of the Geophysical Survey, Federal Research Center of the Russian Academy of Sciences. The Active Faults of Eurasia Database analyse in the manuscript is reported online by the Laboratory of Neotectonics and Recent Geodynamics of the Geological Institute of the Russian Academy of Sciences. Specialized catalog of earthquakes for general seismic zoning of the territory of the Russian Federation is published online by Schmidt Institute of Physics of the Earth of the Russian Academy of Sciences (IPE RAS). The vector and raster maps of railroads are available for free download from Natural Earth at www. natur alear thdata. com.

Declarations
Competing interest The authors have no relevant financial or non-financial interests to disclose.