P Wave Tomography Beneath Greenland and Surrounding Regions: 2. Lower Mantle

We study the 3‐D P wave velocity (Vp) structure of the whole mantle beneath Greenland and surrounding regions using the latest P wave arrival time data. We use a new method of global seismic tomography by setting 3‐D grid nodes densely in the study volume to enhance the resolution. We invert ~5.6 million arrival times of P, pP, and PP waves from 16,257 earthquakes extracted from the ISC‐EHB catalog, which were recorded at 12,549 seismic stations in the world. Our results reveal a hot plume rising from the core‐mantle boundary beneath central Greenland, which is called the Greenland plume. This plume has the main conduit and two branches in the lower mantle, and the main conduit rises to the mantle transition zone (MTZ) beneath eastern Greenland, so it is distinguishable from the Iceland plume that appears to rise from a midmantle depth (~1,500 km) beneath Iceland. The Iceland plume itself is a powerful plume, but it may also be fed by hot mantle materials from three joints with other plumes: a branch of the Greenland plume at ~1,500 km depth, a plume beneath Western Europe at ~1,000 km depth, and the main conduit of the Greenland plume at the MTZ, leading to many active volcanoes in Iceland. We deem that the main conduit of the Greenland plume mixes with the Iceland plume incompletely and splits mainly into the Jan Mayen and Svalbard plumes in the upper mantle, which supply magma to the Jan Mayen hotspot and a geothermal area in western Svalbard, respectively.


Introduction
Greenland, the largest island in the world, spans a wide region over 20°in the central angle of the Earth. The surface-exposed rocks record~4 Gyr of Earth's history in the oldest areas (Henriksen et al., 2009), although now~80% of the land area is covered by thick ice, known as the Greenland Ice Sheet (GrIS). Seismicity in Greenland is very low, and there are no known active volcanoes; hot springs mark only shallow geothermal activity in this island (Hjartarson & Armannsson, 2010). However, the Mid-Atlantic Ridge (MAR) is located just east of Greenland, where the North American plate and the Eurasian plate separate from one another, inducing very active seismicity. The Iceland and Jan Mayen hotspots and a geothermal area in western Svalbard (Dumke et al., 2016;Portnov et al., 2016) are located along the ridge, highlighting the tectonic activity of this region (Figure 1).
The Iceland hotspot is recognized as the surface expression of the Iceland plume, which has the second largest buoyancy flux estimated among all plumes on Earth, and its evolution is believed to have played a significant role in the complex continental breakup history of the northeastern Atlantic ( Barnett-Moore et al., 2017). Attempts to reveal the shape and depth extent of the Iceland plume have been continued mainly by researchers of seismic tomography, which is a powerful tool used to obtain 3-D images of underground structures. A consensus among these studies is the existence of an upper-mantle plume that extends down to at least 400 km depth, which is clearly revealed by both local/regional tomography (e.g., Allen et al., 2002;Foulger et al., 2000Foulger et al., , 2001Hung et al., 2004;Rickers et al., 2013;Wolfe et al., 1997) and global tomography (e.g., Hosseini et al., 2020;Lei et al., 2020;Zhao, 2001Zhao, , 2004Zhao et al., 2013). However, the existence of a lower-mantle plume beneath Iceland, as well as its depth extent, have long been hotly debated. Earlier studies of global P wave tomography imaged a nearly vertical Iceland plume extending down to at least 850 km depth (Rhodes & Davies, 2001) or to the core-mantle boundary (CMB) (e.g., Bijwaard & Spakman, 1999Zhao, 2001Zhao, , 2004. On the other hand, studies of local P and S wave tomography (Foulger et al., 2000(Foulger et al., , 2001 and global S wave tomography (Ritsema et al., 1999) suggested that a low-velocity (low-V) anomaly beneath Iceland extends down to the mantle transition zone (MTZ) but no deeper. Montelli et al. (2004Montelli et al. ( , 2006 conducted global P and S wave tomography and reported a difference between the two models: the Iceland plume is absent at~1,000 km depth in both Figure 1. Surface topography and tectonic setting of the study region. The coordinate transformation is applied (see the text for details). The color scale for the topography and legends for symbols are shown at the bottom. Red triangles: active volcanoes; red circles: hot springs; yellow inverted triangles: seismic stations. The pink shaded patch denotes a geothermal area in Svalbard. GrIS, the Greenland Ice Sheet; JMFZ, the Jan Mayen fracture zone. models; however, it recovers at~1,450 km and extends to the CMB only in the S wave model. Considering these ambiguous results, Montelli et al. (2006) suggested a possibility that the Iceland plume is not continuous but pulsating. More recent studies revealed further elusive features of this plume with a complicated shape not only in the vertical direction but also in the lateral direction (e.g., Zhao et al., 2013 andHosseini et al., 2020 for global P wave tomography; Lei et al., 2020 for global S wave tomography; and Rickers et al., 2013 for regional S wave tomography). Zhao et al. (2013) revealed a prominent localized low-V zone beneath Iceland down to~500 km depth, whereas it becomes wider in the lower mantle (750-1,500 km depth) and less prominent in the lowermost mantle (2,300-2,800 km depth). Rickers et al. (2013) revealed an isolated low-V conduit beneath Iceland bending in the southeastern direction at~1,000 km depth and connecting with the lower-mantle low-V zone beneath Western Europe. Hosseini et al. (2020) revealed a continuous low-V conduit extending from the surface to the CMB, but not exactly vertical, rooted in a low-V anomaly beneath southern Greenland. Lei et al. (2020) imaged a plume that rises from the CMB right below Iceland, amplifies at 660 km depth and ascends to the surface. Large discrepancies among the latest models motivate us to investigate the whole mantle structure beneath Iceland with a much higher resolution.
The Jan Mayen hotspot is considered as the surface expression of the Jan Mayen plume (e.g., Elkins et al., 2016;Schilling et al., 1999). Geochemical studies suggested that the Jan Mayen plume is distinct from the Iceland plume. For example, the Jan Mayen and Iceland plumes are characterized by low and high 3 He/ 4 He ratios, respectively (Schilling et al., 1999). However, seismic tomography studies had long been suffered by the resolution problem to distinguish these two plumes. Recent methodological advances and high-quality seismic data acquisition enable us to resolve this issue (Rickers et al., 2013;Toyokuni et al., 2020), although there are some discrepancies among studies. For example, Rickers et al. (2013) identified an isolated low-V conduit beneath Jan Mayen bending in the southeastern direction at~1,000 km depth and connecting with the lower-mantle low-V zone beneath Western Europe. In contrast, Toyokuni et al. (2020) investigated the structure down to 700 km depth using regional P wave tomography and imaged three isolated plumes in the upper mantle beneath Iceland, Jan Mayen, and Svalbard, which seem to merge in the MTZ.
The existence of an upper-mantle plume beneath Svalbard was first proposed by Toyokuni et al. (2020), which is called the Svalbard plume. The heat source of the geothermal area in western Svalbard (Dumke et al., 2016;Portnov et al., 2016) might be this plume. Toyokuni et al. (2020) also revealed a highvelocity (high-V) anomaly beneath northeastern Greenland, considered as a remnant of ancient oceanic lithosphere and suggested that the absence of volcanoes in Svalbard might be due to the transportation of hot mantle material from the plume being hampered by the high-V body. However, the root of the Svalbard plume remains unknown due to the depth limit (~700 km) of the tomographic model of Toyokuni et al. (2020).
The lower mantle structure beneath Greenland has not been studied well. Shephard et al. (2016) reported a high-V zone at depths of 1,050-1,500 km beneath southwestern Greenland, which is found to be a common feature in many global tomography models (e.g., Ritsema et al., 2011;Simmons et al., 2010) and is interpreted as a subducted slab.
Before 2008, there were a small number of permanent seismic stations in Greenland, including only one station on the GrIS. However, an international project titled "Greenland Ice Sheet Monitoring Network (GLISN)," initiated in 2009, drastically improved the spatial coverage of seismic stations in and around Greenland (Clinton et al., 2014;Toyokuni et al., 2014). Currently, data from 34 GLISN broadband seismic stations including four on the GrIS are available (see Table S1 in the supporting information and Figure 1 of Toyokuni et al., 2020). In the past several years, seismological investigations of the structure beneath Greenland and surrounding regions have been actively conducted using the GLISN data (e.g., Darbyshire et al., 2018;Lebedev et al., 2017;Levshin et al., 2017;Mordret, 2018;Mordret et al., 2016;Pourpoint et al., 2018;Rickers et al., 2013;Toyokuni et al., 2018Toyokuni et al., , 2020. To address the issues mentioned above, the lower mantle structure beneath the MAR and Greenland should be studied in further detail. The purpose of this study is to obtain a detailed 3-D P wave velocity (Vp) model of the whole mantle, in particular, from the MTZ to the CMB beneath the MAR and Greenland, by analyzing data recorded using the latest seismic networks including GLISN, so as to improve our understanding of the deep structure and mantle dynamics of the study region ( Figure 1).

Method
In this study, we apply global tomography to reveal the whole mantle structure beneath Greenland and surrounding regions. Global tomography is a method to study the 3-D seismic velocity structure of the Earth's interior by inverting a huge number of arrival time data of earthquakes distributed all over the world (Zhao, 2015). Zhao (2001Zhao ( , 2004) developed a global tomography method using the geographical grid, that is, 3-D grid nodes are arranged along longitudes and latitudes. Zhao et al. (2013) proposed a method of flexible grid tomography by freely arranging 3-D grid nodes in the whole Earth, regardless of longitude and latitude. However, the conventional global tomography methods do not have sufficiently high resolution to reveal the detailed structure beneath a specific region because the 3-D grid is generally not dense enough, as compared with regional tomography, due to poor data coverage and limited computational resources. To resolve this problem, Zhao et al. (2017) adopted a multiscale 3-D global grid with small grid intervals in their target study region. As a result, despite treating the whole Earth, they could obtain detailed tomography of the subducting Pacific slab beneath Izu-Bonin with a resolution comparable to that of regional tomography.
In the multiscale tomography program by Zhao et al. (2017), the horizontal extent of the study region covered by a denser 3-D grid is specified by longitudinal and latitudinal ranges, which makes it difficult to arrange the grid nodes independently from the geographic coordinates. When we attempt to specify the study area located in high-latitude regions such as Greenland using longitudinal and latitudinal ranges, a fan-shaped area is cut out in the horizontal direction. Therefore, in this study, we apply a method that moves the center of the study region to the equator by a global coordinate transformation (Takenaka et al., 2017). We transform all earthquakes and seismic stations used from the equatorial to ecliptic coordinates so that the location (longitude, latitude) of a reference point in the study region is transformed to (90°, 0°). In this study, the position of the station SUMG (−38.461°, 72.574°) on the summit of the GrIS is defined as the reference point. The details of the coordinate transformation are given in Appendix A of Toyokuni et al. (2020). As a result, the study region is moved to the equatorial region. Thus, when we specify the study region by longitudinal and latitudinal ranges to make the grid interval smaller, the target region appears to be a square (or rectangle) rather than a fan on the map.
Theoretical traveltimes and raypaths are calculated using a 3-D ray tracing technique (Zhao et al., 1992) that combines the pseudobending scheme (Um & Thurber, 1987) and Snell's law. The IASP91 model (Kennett & Engdahl, 1991) is used as the 1-D initial Vp model for the tomographic inversion ( Figure S1). The tomographic inversion is conducted using the LSQR algorithm (Paige & Saunders, 1982) with damping and smoothing regularizations (Zhao, 2001(Zhao, , 2004. The optimal values of the damping parameter λ d = 15, the smoothing parameter in the vertical direction λ sv = 1.05 × 10 −2 , and the smoothing parameter in the horizontal direction λ sh = 9 × 10 −3 are adopted according to the previous study (Zhao et al., 2017). Hypocenter relocation is not conducted as a part of the tomographic inversion.

Data
To improve raypath coverage, we use not only the first P arrivals but also two kinds of later phases, that is, the depth phases (pP) and surface reflected waves (PP), from earthquakes distributed all over the world ( Figure S2). Data are collected from the ISC-EHB catalog at the International Seismological Centre (ISC)

Journal of Geophysical Research: Solid Earth
TOYOKUNI ET AL. site (http://www.isc.ac.uk/) and further selected for analysis. The P, pP, and PP arrival times from 170,435 earthquakes that occurred during 1960-2016 are obtained. The coordinate transformation is applied to all extracted longitudes and latitudes of the epicenters and seismic stations. To make the hypocentral distribution uniform, the whole crust and mantle are divided into small blocks, and only one earthquake with the largest number of data in each block is extracted. To extract as many earthquakes as possible that occurred in Greenland and surrounding regions, smaller blocks with a size of 0.1°(horizontal) × 5 km (depth) are assigned in the study region, and larger blocks with a size of 1.0°(horizontal) × 20 km (depth) are assigned in other regions. We discard data with absolute traveltime residuals exceeding ±5 s through the tomographic inversion. The number of earthquakes used finally is 16,257, and the total number of arrival time data is 5,637,957, which consist of 5,387,114 P, 155,905 pP, and 94,938 PP arrival times, recorded at 12,549 stations ( Figure 3). We note that 303 stations are located in the study region (Figure 1), including 75 stations located between 75°and 105°longitude and −16°and 16°latitude (see Table S1), which also includes 34 GLISN stations as well as other stations. Figure S3 shows the seismic ray coverage in various depth layers. The density of raypath coverage varies throughout the study region, which can be visualized by ray hit counts (HC), that is, the number of rays sampling a grid ( Figures S4-S7). To illustrate the HC distribution, areas with HC < 5, 5 ≤ HC < 50, and HC ≥ 50 are shown in black, gray, and white, respectively. White areas are where seismic rays crisscross well and the tomographic results are highly reliable. Robust results are expected in areas with a HC ≥ 5, and less reliable areas are masked in white in the tomographic images. We note that most of the study region has HC ≥ 50 at depths ≥300 km.

Results
The results are shown in map views in Figures 4 and 5. The corresponding HC distributions are shown in Figure S4. At depths of 5-80 km, only spot-like high-Vp and low-Vp zones are visible just below several stations, and in most regions, the initial Vp perturbation is not updated from zero. This is because the structure at these depths is not well resolved because the number of local earthquakes is relatively small, and so the density of seismic rays crisscrossing in the shallow regions is quite low. At depths of 120-220 km, low-Vp anomalies are remarkable beneath the Iceland and Jan Mayen hotspots, hot springs in eastern Greenland, and the geothermal area in western Svalbard. This low-Vp zone spreads further at depths of 280-490 km, and two significant low-Vp anomalies exist beneath Iceland-Jan Mayen-eastern Greenland and Svalbard. The most remarkable high-Vp zone at these depths exists beneath the northeastern coast of Greenland to the northeast offshore. In particular, at depths of 310-490 km, the high-Vp region extends inside the bend of the MAR and seems to divide the two low-Vp regions as mentioned above. At the same depths, the high-Vp region in southwestern Greenland is also remarkable. At depths of 520-670 km, the low-Vp anomaly is prominent off the southern tip of Greenland in addition to the low-Vp anomaly beneath Iceland-Jan Mayen-eastern Greenland. The low-Vp anomaly beneath Svalbard is also widened. These low-Vp zones seem connected together through low-Vp corridors. The regions just beneath inland Greenland are generally imaged as high-Vp anomalies.   At depths of 750-1,050 km, significant low-Vp regions are confined only to those beneath Iceland and eastern Greenland. As for high-Vp regions, a wide extent is observed beneath Ellesmere and Baffin Islands to southern Greenland, and other smaller areas beneath central Greenland and offshore northern Greenland to the MAR at depths of 750-850 km. However, at greater depths, these high-Vp regions become unclear except for those offshore of northern Greenland. At depths of 1,100-1,500 km, the low-Vp anomaly beneath Iceland moves toward the southeast as the depth increases. There are also spot-like low-Vp anomalies beneath northern and central Greenland and Ellesmere Island (at depths ≥1,350 km). The high-Vp anomaly is most prominent beneath western to southwestern Greenland at depths of 1,200-1,500 km.
At depths of 1,600-1,800 km, the low-Vp anomaly beneath Iceland becomes inconspicuous, and weak low-Vp zones appear in various places. The low-Vp area beneath central to southeastern Greenland and another beneath Ellesmere Island are the most prominent. High-Vp anomalies are hardly observed at these depths. At depths of 2,000-2,400 km, low-Vp anomalies located above these depths are connected, forming a wide low-Vp region running NW-SE beneath Greenland. At depths ≥2,200 km, a clear high-Vp anomaly appears beneath the northeastern coast of Greenland. At depths of 2,600-2,880 km, the low-Vp anomaly beneath Greenland tends to appear more concentrated to a central location with depth. The high-Vp anomaly beneath the northeastern coast of Greenland becomes inconspicuous, but larger amplitudes of high-Vp or low-Vp regions can be seen in various places, indicating that the structural heterogeneity is stronger at these depths compared to the shallower depths. In Figure 6, the A-A′, B-B′, and F-F′ sections show high-Vp anomalies tilting north or south from the surface to the MTZ. However, in the upper mantle above the MTZ, the regions with larger HCs are distributed with the same slope as indicated in white in Figure S5; thus, these high-Vp regions are likely to have artificial shape along the seismic rays just below the stations. In the F-F′, G-G′, and H-H′ sections through Greenland, the low-Vp anomaly that seems to be a hot plume extends from the CMB to a depth of 1,000 km. On the F-F′ section beneath western Greenland, a high-Vp zone that looks like a subducted slab exists at depths of 1,000-1,500 km. In the J-J′, K-K′, L-L′, M-M′, and N-N′ sections in the eastern part, the HC near the surface is improved (Figures S4 and S5), thereby making the tomographic patterns shallower than the MTZ reliable. In the J-J′ and K-K′ sections, a nodule-like high-Vp zone is visible from the surface to the MTZ beneath the northeastern coast of Greenland. In the K-K′, L-L′, and M-M′ sections, clear low-Vp anomalies exist beneath the Iceland and Jan Mayen hotspots. These low-Vp anomalies are separated beneath the two hotspots near the surface but are connected in the MTZ, forming a widespread low-Vp zone. In the K-K′ and L-L′ sections, the low-Vp anomaly beneath Svalbard is also remarkable; however, the anomaly seems separated from the low-Vp zones beneath Iceland and Jan Mayen due to the high-Vp body beneath northeastern Greenland.
In Figure 7, the H-H′ section passing through central Greenland and Jan Mayen shows a clear plume-like low-Vp anomaly that rises from the CMB, breaks shortly at a depth of 1,000-1,300 km, and reaches the Jan Mayen hotspot. A low-Vp anomaly with a similar shape is also visible in the J-J′, K-K′, and L-L′

10.1029/2020JB019839
Journal of Geophysical Research: Solid Earth sections beneath Iceland; however, the plume-like low-Vp zone rising from the CMB is not along the vertical cross section. It seems to be connected to the low-Vp conduit beneath Iceland (≤1,500 km depth) after meandering toward the south in the direction perpendicular to the cross section at~1,500-2,000 km depths. The K-K′ section shows that the low-Vp anomaly beneath Iceland tilts eastward from the surface to a depth of 1,000 km and continues almost vertically to a depth of 1,500 km.
In the B-B′, C-C′, and D-D′ sections in NW-SE direction (Figure 8), the plume-like low-Vp anomaly rising from the CMB beneath central Greenland is prominent. It extends to both the northwest and southeast directions at the shallow depths, but the connection with the low-Vp anomaly beneath Iceland is unclear. The C-C′ and D-D′ sections in NE-SW direction (Figure 8), passing through central Greenland and Svalbard, clearly show a plume-like low-Vp anomaly extending from the CMB to beneath Svalbard. However, it seems to be interrupted in the MTZ due to the high-Vp body beneath northeastern Greenland.
The CPU time for one tomographic inversion is approximately 34 h and 28 min using our workstation computer Xeon E5-2660 v3 (2.6 GHz, 1 core). The root-mean-square (RMS) traveltime residual is 1.60 s for the initial 1-D Vp model, but it is reduced to 1.06 s for the final 3-D Vp model.

Resolution Tests
We conducted extensive resolution tests including the checkerboard resolution tests (CRTs) (Humphreys & Clayton, 1988;Zhao et al., 2017), restoring resolution test (RRT) (Zhao et al., 1992(Zhao et al., , 2017, and synthetic resolution tests (SRTs), to evaluate the adequacy of ray coverage and spatial resolution. To conduct the CRTs, we first develop an input velocity model that contains alternate positive and negative Vp anomalies assigned to

Journal of Geophysical Research: Solid Earth
the 3-D grid nodes; the amplitude of Vp anomalies is 3% at all the nodes. The following two grid types are prepared and the inversion is conducted for each grid type: (1) CRT1: the lateral grid interval is 278 km (a great circle distance of 2.5°on the surface) inside the study region, and 556 km (a great circle distance of 5°on the surface) in other regions (Table S3), and (2) CRT2: the lateral grid interval is 167 km (a great circle distance of 1.5°on the surface) inside the study region, and 334 km (a great circle distance of 3°on the surface) in other regions (Table S4). To conduct the RRT, we highlight the patterns of the actual tomographic result when constructing the RRT input Vp model. Namely, at the grid nodes with the Vp perturbation (dVp) > +0.5% in the real tomographic model, we set dVp = +1.5%, and at the grid nodes with dVp < −0.5% in the real tomographic model, we set dVp = −1.5% in the RRT input model. The Vp anomalies at the other grid nodes are set to zero. To conduct the SRTs, we construct an input model with dVp = −1.5% representing plumes with geometrical shapes; the dVp at all other grid nodes is set to zero. Two SRT inversions with different input models are conducted: (1) SRT1: there are three continuous and tilted low-Vp conduits, which rise from the same low-Vp region on the CMB beneath central Greenland and reach the crust beneath Iceland, Jan Mayen, and Svalbard, and (2) SRT2: there are three continuous and vertical low-Vp conduits beneath Iceland, Jan Mayen, and Svalbard from the surface to the CMB. In the SRT1 input model, the horizontal sections of the three low-Vp conduits are circular, whose central position (longitude, latitude) is (98°, −6°), (99°, 1°), and (99°, 9°) at the surface of Iceland, Jan Mayen, and Svalbard, respectively. The plume radius at the surface is 1.5°, 1°, and 1°beneath Iceland, Jan Mayen, and Svalbard, respectively. These three plumes tilt and merge to the same low-Vp anomaly whose center (longitude, latitude) is (90°, 2°) at the CMB beneath Greenland. In the SRT2 input model, the locations and radiuses of three low-Vp conduits at the surface are the same as those in the SRT1, but they are vertical and completely separated in the whole mantle. Synthetic datasets for the CRT, RRT, and SRT are constructed by calculating the theoretical arrival times for each input model but with random errors added, which range between −0.2 and +0.2 s with a standard deviation of 0.1 s. In the RRT and SRT, we use the same grid distribution in the main computation (section 2).
Main features of the test results are summarized in Figure 9; the complete test results are shown in the supporting information for CRT1 ( Figures S8-S11), CRT2 (Figures S12-S15), RRT (Figures S16-S19), SRT1 ( Figures S20-S23), and SRT2 ( Figures S24-S27). For the CRT results, the recovery rate is defined as follows: On the map views of the CRT1 results (Figures 9a and S8), the output patterns are biased to either low-Vp or high-Vp, and the resolution is clearly poor at the shallow depths (15-140 km). As for the recovery rate, the black-to-gray areas with poor recovery are dominant at depths of 15-260 km, but white areas are dominant from a depth of 300 km, and almost all areas deeper than 500 km are shown in white. In the vertical cross sections ( Figures S9-S11), the depth extent of the areas with good resolution can be confirmed more clearly, displaying almost the same pattern of reliability shown by the ray HC ( Figures S5-S7). The CRT2 results show that the areas with good recovery are dominant at depths ≥300 km, as in the previous case, and can be seen in almost everywhere at depths ≥725 km (Figures 9b and S12-S15). The two CRT results show that the resolution in our study region is~200 km in the horizontal direction and the distances comparable to the vertical grid interval in the depth direction, from the MTZ to the CMB.
The RRT results show that the pattern of the input model is generally recovered, though the vertical and lateral smearing as well as spurious velocity anomalies are more prominent in and around the MTZ (Figures 9c and S16-S19). In the current dataset, the seismic ray density above the MTZ is low (Figures S3 and S4), so the near-surface effects have been imposed to the velocity anomalies in the MTZ. The SRT1 results indicate that the Vp amplitudes in the output model decrease remarkably at depths ≤100 km, though the input low-Vp amplitude is constant at all depths (Figures 9d and S20). In addition, as in the RRT case, the spurious low-Vp anomalies spread more widely in and around the MTZ (380-725 km depth) ( Figure S20). However, the pattern and amplitude of the input model are mostly recovered except around the surface, and no spurious high-Vp anomalies are generated. Furthermore, spurious bending or division of the plume does not appear on all three low-Vp conduits (Figures 9d and S21-S23). The SRT2 results also indicate that the three vertical whole-mantle plumes are well recovered, and their independency is well resolved  and S24-S27). The two SRT results show that the shapes and amplitudes of plume structures in the lower mantle can be resolved well with our dataset and method. In other words, branching and changes in the low-Vp amplitudes of the plume in the lower mantle, which appear in our tomographic results, are considered to be robust features.

Consistency With Other Models
Toyokuni et al. (2020) determined a 3-D Vp model from the surface to 700 km depth beneath Greenland and surrounding regions using regional tomography. Figure 10a shows a comparison of vertical cross sections in 10.1029/2020JB019839

Journal of Geophysical Research: Solid Earth
the N-S direction through Iceland, Jan Mayen, and Svalbard obtained by Toyokuni et al. (2020) and this study. Toyokuni et al. (2020) revealed clear separation of the Iceland, Jan Mayen, and Svalbard plumes in the upper mantle as well as connection of these three plumes in the MTZ, forming a widespread low-Vp anomaly. Toyokuni et al. (2020) also revealed that a high-Vp body existing at depths ≤500 km between Jan Mayen and Svalbard might block plume flow, resulting in insufficient magma transportation to the Svalbard region to create an active volcano. Regarding this comparison, the global tomography model (lower panel of Figure 10a) has a lower resolution at depths <~200 km but shows clear consistency with the regional model (upper panel of Figure 10a).
At depths of 310-700 km where both the regional and global models have a good resolution (Figure 4 of Toyokuni et al., 2020, and Figure 4 of this paper), the two models are generally consistent. For example, the high-Vp anomalies beneath northeastern Greenland at depths of 310-490 km and those beneath northwestern Greenland at depths of 430-700 km, and southwestern Greenland at depths of 310-700 km show good agreement (see also the left panel of Figure 10f). The difference between the two models can also be seen because of the advantages and disadvantages of both methods. For example, according to the RRT and SRT results, the global tomography tends to generate spurious structures in and around the MTZ, whereas the resolution of the regional tomography is relatively low near the boundary of the study volume. Nevertheless, the fact that the same features are revealed with the different methods and datasets indicates the robustness of the two tomography models.
Our results show a remarkable high-Vp anomaly at depths of 1,200-1,500 km beneath southwestern Greenland ( Figure 5 and the F-F′ section in Figure 6). This feature corresponds well to a possible slab remnant suggested by Shephard et al. (2016) using the existing global tomography models (section 1). This study provides a more detailed view of this high-Vp anomaly by enhancing the tomographic resolution beneath Greenland.

Volcanic and Geothermal Activities Along the MAR
This study revealed a distinct plume-like low-Vp zone rising from the CMB beneath central Greenland (Figures 10b−10e). Figure 10b shows that two branches of the low-Vp zone rise in the northwest and southeast directions, but their relationship with the Iceland plume in the upper mantle is unclear. In Figure 10c, the Iceland plume continues down to 1,500 km depth; the southeast branch of the low-Vp zone beneath Greenland is close to the Iceland plume, but their connection remains unclear. On the other hand, Figure 10d shows that the upper-mantle Jan Mayen plume is weakly inclined to the west and seems connected to the low-Vp zone beneath Greenland. Their connection at depths of 1,000-1,300 km seems rather unclear, but the similarity of the low-Vp slopes above and below these depths strongly suggests the continuity of the plume. Figure 10e also shows that the upper-mantle Svalbard plume may be continuous with the low-Vp zone beneath Greenland, although it seems to be interrupted by a remarkable high-Vp anomaly beneath northeastern Greenland in the MTZ, as observed in Toyokuni et al. (2020). Based on these results, we consider the lower-mantle low-Vp zone beneath Greenland to be a mantle plume, which may form a root of the Jan Mayen and Svalbard plumes in the upper mantle. Therefore, we call it the "Greenland plume." Hereafter we call the portion of the Greenland plume that reaches the MTZ beneath eastern Greenland the "main conduit," and the two branches shown in Figure 10b the "northwestern (NW) branch" and "southeastern (SE) branch." The branching of the Greenland plume is also clearly visible in Figure 10f.
To further examine the relationship between the Iceland plume and the SE branch of the Greenland plume, more vertical cross sections are presented in Figure 11. Figure 11a is a cross section along an irregular profile passing through North America, Greenland, Iceland, and Western Europe. This image and other cross sections (Figures 11b and 11c) show that the plume beneath Iceland dips eastward from the surface tõ 1,000 km depth, then continues vertically to 1,500 km depth but becomes indistinct below this depth. However, there seems to be weak connections both with the SE branch of the Greenland plume at 1,500 km depth ( Figure 11a) and with another plume-like low-Vp zone beneath Western Europe at 1,000 km depth (Figure 11c). Therefore, we define the distinct low-Vp conduit beneath Iceland at depths ≤1,500 km with a nearly constant diameter as the Iceland plume. The image of the Iceland plume obtained in this study supports the results of Zhao et al. (2013) showing that the Iceland plume is prominent at depths ≤1,500 km (Figure 11b). The connection of the Iceland plume with the low-V region beneath Western Europe should be further verified as a future work as the joint is located near the boundary of

10.1029/2020JB019839
Journal of Geophysical Research: Solid Earth our study region. However, this connection is also supported by the results by Rickers et al. (2013) (Figure 11c). Toyokuni et al. (2020) suggested that the Iceland, Jan Mayen, and Svalbard plumes merge in the MTZ, and therefore, the main conduit of the Greenland plume might mix with the Iceland plume. Consequently, the Iceland plume may be fed by three sources of hot mantle materials, that is, the connection to the main conduit and SE branch of the Greenland plume, and to another plume beneath Western Europe.
The SRT results show good seismic illumination of these regions, and the continuous low-Vp conduits with constant low-Vp amplitudes are well resolved from~200 km depth to the CMB (Figures 9d and 9e). Therefore, the difference in the low-Vp amplitudes along the rising Greenland plume (e.g., Figures 10d  and 10e) may reflect true differences in plume diameter and/or short interruption of the plume route, similar to the pulsating plume hypothesis suggested by Montelli et al. (2006). Figure 7 shows that the joint between the SE branch of the Greenland plume and Iceland plume is slightly bent in the lateral direction (section 4), Figure 12. Schematic diagram showing main plume features revealed in this study. The Greenland plume originates at the CMB beneath Greenland and has three branches in the lower mantle (main conduit, NW branch, and SE branch). The Iceland plume is a midmantle plume rising from~1,500 km depth. The bottom of the Iceland plume may be connected with the SE branch of the Greenland plume and a hot plume beneath Western Europe. The main conduit of the Greenland plume and the Iceland plume merge in the mantle transition zone (MTZ) and then rise as three narrow plumes toward the surface. The Jan Mayen and Svalbard plumes mainly derive from the main conduit of the Greenland plume. The high-Vp anomaly in the upper mantle beneath northeastern Greenland seems to become an obstacle to the flow of the Greenland-Svalbard plume system.

Journal of Geophysical Research: Solid Earth
which is a reliable feature confirmed by the resolution tests. We note that the shape around the bottom of the Iceland plume is complex, causing it to vary in appearance in different cross sections. Figure S28 shows vertical cross sections through the Greenland plume and the Jan Mayen hotspot. Rickers et al. (2013) suggested the continuity of the Jan Mayen plume to the low-V anomaly beneath Western Europe. However, we cannot confirm a clear connection of this plume with anything other than the Greenland plume although this should be further verified. In addition, we cannot identify any other roots of the Jan Mayen and Svalbard plumes beneath the MAR ( Figure S29).
As mentioned above, all three plumes beneath Iceland, Jan Mayen, and Svalbard may be affected by the lower-mantle Greenland plume. Combining these results, Figure 12 shows a schematic diagram of the mantle plume structures revealed by this study. The Jan Mayen and Svalbard plumes in the upper mantle are roughly considered as branches of the main conduit of the Greenland plume. We interpret that the Iceland plume is rather special because (1) the Iceland plume itself may be a powerful plume with a thick diameter and strong low-Vp amplitudes, (2) its bottom is located in the mid-lower mantle, and (3) it has three joints with other plumes. Only Iceland has many volcanoes, which may be caused by this specialty of the Iceland plume. Geochemical studies suggested differences between the Iceland and Jan Mayen plumes (e. g., Elkins et al., 2016;Schilling et al., 1999), which also support our interpretation. The source of high magma supply to Jan Mayen has long been remained unclear because the spreading rate of MAR adjacent to this island is ultraslow (e.g., Elkins et al., 2016;Mertz et al., 1991;Neumann & Schilling, 1984;Schilling et al., 1999). The discovery of the Greenland plume might be a clue to solve this puzzle.
It seems that some previous studies identified the Greenland plume as the deep portion of the Iceland plume (e.g., Bijwaard & Spakman, 1999Hosseini et al., 2020;Zhao, 2001Zhao, , 2004. Further detailed studies are needed to distinguish the Iceland and Greenland plumes. However, our present distinction seems reasonable with the following considerations: (1) the Greenland plume affects not only the Iceland hotspot but also the Jan Mayen hotspot and the geothermal area in Svalbard and (2) the Iceland plume looks continuous only at depths ≤1,500 km, whereas it becomes complex and obscure in the deeper mantle till the CMB.

Conclusions
A detailed 3-D P wave velocity (Vp) model of the whole mantle beneath Greenland and surrounding areas is obtained by inverting a large number of high-quality P wave arrival time data recorded by the updated seismograph network in the study region. The 3-D Vp structure from the MTZ to the CMB is effectively resolved. The novel tomographic model reveals the following new features.
1. The Iceland plume rises from~1,500 km depth in the lower mantle. A lower-mantle plume rising from the CMB beneath central Greenland is revealed, and we call it the Greenland plume. 2. The Greenland plume has a main conduit and two branches in the lower mantle, which we call the NW branch and SE branch. The main conduit reaches to the MTZ beneath eastern Greenland. 3. The Iceland plume itself is a powerful plume with a strong low-Vp anomaly, but it might also be fed by hot mantle materials at three joints with other plumes, that is, joints with the SE branch of the Greenland plume at~1,500 km depth, a plume beneath Western Europe at~1,000 km depth, and the main conduit of the Greenland plume at the MTZ. As a result, there are highly active volcanoes in Iceland. 4. In the upper mantle, the Greenland plume mainly splits to the Jan Mayen and Svalbard plumes, supplying hot mantle materials to feed the Jan Mayen hotspot and the geothermal area in western Svalbard. 5. A distinct high-Vp anomaly is revealed beneath southwestern Greenland at depths of 1,200-1,500 km, which may reflect a slab remnant.