Improving Microseismic P Wave Source Location With Multiple Seismic Arrays

Using array analysis, the direction and distance to a seismic P wave source can be determined. However, individual arrays are limited in their geographical coverage and by their resolving capability, which is determined by the array aperture, configuration, and number of stations. We demonstrate these limitations on three large seismic arrays located in Japan, Europe, and California, and find that all give a unique but imperfect insight into the P wave sources acting in the North Pacific. We then combine the data from all three arrays into one beamforming image. The combined images bring together the views offered by each array, providing a concise, comprehensive, and more robust representation of multiple source locations. Next we weight each array for distance in order to optimize the result. Being able to resolve and accurately locate source regions is an important step in being able to use seismic records to monitor ocean wave activity and track storms in real time.

Accurately locating microseism sources is essential for relating the measured microseisms to wave activity and is also important for tomography studies which may be biased by a nonisotropic distribution of sources (Ermert et al., 2016;Harmon et al., 2010;Retailleau et al., 2017;. However, individual arrays are limited by their coverage and also by their resolution in their ability to resolve closely spaced sources, as noted previously by Hillers et al. (2012) and Euler et al. (2014). In such instances, two closely spaced sources will 10.1002/2017JB015015 appear in the beamforming result as a single source centered on the average location. Using multiple seismic arrays will improve coverage and may improve the ability to resolve neighboring P wave sources. Previous studies which have combined arrays (Euler et al., 2014;Hillers et al., 2012;Landès et al., 2010;Meschede et al., 2017;Pyle et al., 2015) have mostly focused on detecting common, robust sources over monthly or seasonal time periods, but here we combine arrays with the purpose of locating multiple or closely spaced sources which are active over short (3 h) time periods. This builds on our previous work (Neale et al., 2017), where we located P wave sources in the North Pacific with the same temporal resolution using just the California array, and which found location accuracy to be up to 10 ∘ in most cases. An improvement in source location will be beneficial for possible ocean or seismic monitoring, ocean wave model validation, and data assimilation into operational wave models. We show here that a combination of arrays can improve the imaging of multiple sources not visible by individual arrays alone.

Data and Methods
2.1. Numerical Modeling of P Wave Sources P wave microseism sources were modeled following the method of Ardhuin et al. (2011), Herbers (2013), andFarra et al. (2016). The second-order pressure spectrum F P (with units of Pa 2 m 2 s) resulting from the interaction of opposing ocean waves of similar frequency f (Hasselmann, 1963) is given by where f 2 = 2f , x is location, w is the density of seawater, g is gravitational acceleration, E(x, f ) is the ocean wave frequency spectrum, and I(x, f ) is a nondimensional function that depends on the wave energy distribution M over the directions : F P was calculated on a global 0.5 ∘ longitude by 0.5 ∘ latitude grid using the numerical ocean wave model WAVEWATCH III (Tolman, 2014) with coastal reflection of R 2 = 0.1 for continents and large islands and R 2 = 0.2 for small islands (Ardhuin et al., 2011). We forced the model using 6-hourly wind and sea ice concentration from ECMWF's ERA-Interim reanalysis (Dee et al., 2011). F P was then subsampled to a 2 ∘ -by-2 ∘ grid (smoothed first by convolving with a 7-by-7 pixel Gaussian kernel) and averaged over 3 h time periods between October 2012 and March 2013.
The P wave microseism source at each grid point was calculated by multiplying the second-order pressure spectrum by the squared source site effect [2|C P | c w ] (Farra et al., 2016): We used w = 1, 020 kg m −3 and spatially varying values of crustal density c were taken from the global crustal model CRUST1.0 (http://igppweb.ucsd.edu/∼gabi/rem.html). C P is an amplification coefficient that depends on frequency, water depth, and P wave take-off-angle (Ardhuin & Herbers, 2013;Gualtieri et al., 2014). We used the formulation of Gualtieri et al. (2014) and Farra et al. (2016) to calculate C P at each frequency and take-off angle using water depth from ETOP01 bathymetry (Amante & Eakins, 2009). Bathymetry is shown in Figure 1a. Take-off angle depends on the distance between the source grid point and the array, so the site effect maps are slightly different for each array. Therefore, we use different maps of C P to calculate the modeled source relevant for each array. These maps are shown in Figure S1 of the supporting information for 0.207 Hz and are visually very similar, but do have amplitude differences of up to 30%, with the largest differences at locations closest to the array center. The maps are like those shown previously by Gualtieri et al. (2014) and Farra et al. (2016). The three separate maps of C P were averaged to calculate the modeled source for comparison with the combined array outputs, and this map for a frequency of 0.207 Hz is shown in Figure 1b.
We used the modeled source in order to evaluate the seismic observations. The modeled source may itself not be completely accurate due to potential errors in wind forcing and the resulting wave interaction terms, but should give a good indication of the location of major P wave sources as it has been found to in previous studies (Gualtieri et al., 2014;Neale et al., 2017;Obrebski et al., 2013). We considered P wave microseism sources at frequencies of 0.188 and 0.207 Hz to focus on the peak frequency of double-frequency P wave sources in the North Pacific (Meschede et al., 2017;Neale et al., 2017) and searched for time periods with multiple or closely spaced sources acting in the basin.

Seismic Data and Processing
Seismic vertical component data were downloaded for the days identified from the modeled sources. Continuous waveforms were obtained from the Southern California Earthquake Data Center (SCEDC, 2013), the Data Management Center of Japan's National Research Institute for Earth Science and Disaster Resilience (Obara et al., 2005), the European Integrated Data archive (http://www.orfeus-eu.org/eida/eida.html) and using FDSN Web Services (http://www.fdsn.org/webservices/). Figure 2 shows the networks and distribution of stations within the arrays.
Data were downsampled to 1 Hz if necessary and band-pass filtered between 0.002 and 0.400 Hz. Instrument response was removed. Earthquake events over magnitude 5 were identified using the ISC bulletin (International Seismological Centre, 2013) and removed from the data by setting a 1 h window of the waveform to zero if the RMS of that window was over 3 times the daily RMS. Any remaining bad quality data were identified by visual examination of the daily spectra and discarded. For each array separately, P wave sources were located from beamforming analysis (Farra et al., 2016;Gerstoft et al., 2006Gerstoft et al., , 2008Gualtieri et al., 2014;Obrebski et al., 2013;Rost & Thomas, 2002), using the predicted arrival times for P waves rather than a plane wave assumption. The beamformer output was calculated directly onto a 2-by-2 ∘ geographical grid based on the predicted P wave travel times, t P , from each grid point x to each station: N is the total number of stations, S j is the complex spectrum of the record at station j, f 2 is the seismic frequency, and t i refers to the start time of the Fourier transform. S j was calculated on 512 s segments overlapping by 50% and was divided by its magnitude to retain only phase. We used f 2 = 0.188 ± 0.01 and 0.207 ± 0.01 Hz. P wave travel times for a source at the surface were obtained using the global "AK135" travel time tables (Kennett et al., 1995). The angle brackets denote averaging over each frequency band and over 3 h periods to get the output at a 3 h time averaged slice t s .
Each beam projection is limited to a geographical range of 99 ∘ from the array, as beyond this the P waves will start to be effected by the Earth's outer core. In addition, P waves from locations too close to the array will be regional and the global velocity model we use is unlikely to be accurate enough. We would expect bias from upper mantle velocity heterogeneities to affect P waves originating from less than 30 ∘ . However, the region between 15 and 30 ∘ of the Japan array does contain microseism sources we are interested in imaging, so we tested the accuracy of beam projections in this area by evaluating the projection at 0.2 Hz against a known earthquake ( Figure S2). We found the maximum of the beam projection to be located within one grid point (on a 2-by-2 ∘ grid) of the earthquake location, and so we continued to project for a minimum distance of 15 ∘ . With these limits, the arrays in California and Japan both give good coverage over the North Pacific, while the European array offers coverage in the northern half of the North Pacific.
Looking at multiple images to pick out the best features from each array when the true source is unknown would be time consuming and would probably require input and judgment from a human operator. This is not ideal in the context of using the observed sources for ongoing and possibly near-real time monitoring, so in order to synthesize the results from multiple arrays we combined the data from all arrays to produce one beamforming image. Because t P of equation (4) is independent between each station and each grid point, assuming AK135 traveltimes, we simply included the stations from all arrays in the summation. While there is unlikely to be much variation in the deviation from AK135 traveltimes within an array (so that P waves propagate coherently given an expected traveltime) coherence may be broken when combining three arrays. This would degrade the beam power; however, individual peaks would remain. Additionally, because N in equation (4) is the total number of stations, not the number of stations used at grid point x, we have to adapt (4) to normalize for the number of stations within the range of each grid point, otherwise with sources being equal everywhere, the regions covered by three arrays would appear approximately 3 times stronger than those covered by one array, if the arrays contained similar numbers of stations. We therefore replace the 1∕N 2 in equation (4) with a weight w j (x) that is contained in the summation: Finally, the weight w was alternatively adjusted for each grid point (x) and each station (j) based on the distance (d) in degrees from that grid point to the station, and used in equation (5): This is a semiempirical weighting intended to boost the influence of the closest array which is expected to perform best in separating out multiple sources (e.g., as shown in Figures 8 and 9 for sources close to Japan).
Here the term ∑ 1∕d equalizes the total weighting at each grid point.

Synthetic Tests
The ability of an array to locate sources and resolve neighboring sources depends on its resolution. The angular resolution of an array is determined by the array aperture and depends on the wavelength (hence frequency) of the wave being examined (Rost & Thomas, 2002). As we are beamforming onto a geographical grid, angular resolution corresponds to a spatial resolution that also depends on distance from the source. The spatial resolution, R, of the array can be approximated by the intensity pattern width associated with a point source (Hillers et al., 2012). This is equal to R = Δ∕D where is seismic wavelength, Δ is distance from array to source, and D is the aperture of the array. For example, for the Californian array, a 0.2 Hz source at 43 ∘ N 171 ∘ E ( = 28 km, Δ = 6,100 km, D = 1,400 km) has a theoretical spatial resolution R ∼ 120 km. For the Japan array, the same source at the same location (Δ = 3,400 km, D = 610 km) has a theoretical spatial resolution of 160 km and for the Europe array (Δ = 10,000 km, D = 1,300 km) it is 210 km. In practice, the ability of an array to resolve neighboring sources will also depend on the number of stations (resolution of velocity) and station configuration (spatial aliasing) Thomas, 2002, 2009) as well as the beamforming technique used. (Hillers et al., 2012) found that neighboring sources with a separation distance much farther apart than their estimated resolution could not be resolved using conventional beamforming.
We tested the resolving capability of each array using synthetic point sources of amplitude A = 1, by replacing S j of equation (4) with Syn j , where Syn j = ∑N p k=1 Ae −i2 f 2 t P k ; that is, the sum of N p phase delays for sources located at grid points k 1 … k N p . For each array, we placed a source at 43 ∘ N 171 ∘ E and a second source at increasing distances to either the north or the east. The synthetics were tested at f 2 = 0.207 Hz. For the California array The modeled source has units of Pa 2 m 2 s and is saturated in color to 50% of the maximum value of the average modeled source. Beampower for each plot has been normalized between 0 and 1. The red circles indicate peaks in the beamformer output with the same color scale as the contours. The green lines indicate the boundaries at 15 ∘ and 99 ∘ from the array center. The number of stations in each array is given by n.
( Figure S3) the two sources separated when the second source was placed at a distance of 670-1,100 km to the north (Figures S3b and S3c) or 1,300-1,500 km to the east (Figures S3e and S3f ). The lower limit of this range, and shown in Figures S3b and S3e, is the distance for partial separation, while the upper limit, and shown in plots Figures S3c and S3f, is for full separation. For the Japan array ( Figure S4) they separated at a distance of 670-890 km to the north (Figures S4b and S4c) or 650-2,100 km to the east (Figures S4e and S4f ). For the Europe array ( Figure S5) it was 890-1,300 km to the north (Figures S5b and S5c) or 820-1,000 km to the east (Figures S5e and S5f ). It is not possible to compare the values between the different arrays in this way because of the different distances involved, with the sources placed a lot closer to the Japan array than the Europe array; however, it can be seen that these distances are all much larger than the theoretical resolutions calculated previously, and for real sources the resolution may be even lower as sources may not act as point sources. Furthermore, the synthetic tests reveal the importance of viewing direction when separating two neighboring sources. For both the California array and Japan array, north-south separated sources would be identified from their different azimuths, whereas for the Europe array north-south separated sources would be identified by their different slownesses (take-off angle). In each case, sources separated azimuthally from the array are better resolved than those separated by distance, meaning that the California and Japan arrays are better able to resolve north-south separated sources than east-west separated sources, whereas the Europe array is better able to resolve those separated east-west than north-south.
A benefit of using a combination of different arrays can be seen in the synthetic tests shown in Figure 3. Here one point source was placed at 49 ∘ N 153 ∘ E and imaged using each array individually and then combined with and without the distance weighting. Out of the three individual arrays the source was imaged best using the Japan array, as expected because of the close distance. The California array beamformer peak is much broader, and the Europe array shows more artifacts generated by beam sidelobes. Both combined outputs reduced the width of the peak to correspond much more closely with the point source that it is; however, there are some array artifacts appearing as stripes, which are stronger in the combination without distance weighting. With distance weighting, however, another feature is seen in the apparent increased amplitude in the region below 15 ∘ of the Japan array center. This is due to the California array, in particular, having a much broader peak, and in this region out of the Japan array's range this broadened peak is no longer outweighed by the low values of the Japan array as it is the region between 15 ∘ and the point source. This artifact disappears when the source is farther away from the inner boundary ( Figure S6), but because of its presence in some cases we remove the regions within 15 ∘ of an array center. These regions could still be displayed separately, but increasing amplitudes toward the outer edge would need to be treated with caution as they could indicate either a source out of range or a true extended source over the boundary.

Individual Arrays
Observed and modeled sources for each array at selected time periods are shown in Figures 4 to 11. Animations at 3-hourly intervals are also available as supporting information (Animations S1 to S8). Results are shown for either 0. Figures 6a-6d on 10 October 2012 06:00-09:00 shows a case when it is only the Europe array that observes a source in the north (again just south of the Aleutian Islands) but misses the other two sources to the east and southwest which are out of range. The California array observes both the east and southwest source, while the Japan array observes a strong southwest source and a weak blurring of the north and east sources.
A benefit of utilizing multiple arrays is therefore better coverage and the potential of observing sources that would otherwise be missed. Another, related, advantage is that different arrays are able to separate out neighboring sources in different situations. One example of this is seen on 27 October 06:00-09:00 (Figure 7), where there are two east-west aligned sources separated by 1,800 km. The Europe array distinguishes two sources, albeit with more sidelobe artifacts, but they are merged and centered on the average location between them in the outputs from California and Japan. This also indicates that the synthetic point tests overestimated the ability of the arrays to resolve neighboring sources, as both the California and Japan arrays could begin to separate the synthetic sources at distances of about 1,300 km.
In other examples, sources near the east coast of Japan were very well resolved by the Japan array but poorly or unresolved by the other arrays. This is seen in the example on 19 November at 03:00-06:00 (Figure 8). The California array merges all of the sources on a roughly north-south line east of Japan. The Europe array locates the strong source in the north more accurately than the California array, but mislocates the source farther south. The Japan array resolves the strongest source very well, and also identifies the two other source regions. On 19 January 00:00-03:00 (Figure 9) it can again be seen that the Japan array is able to resolve north-south separated sources which are merged by the California array. In other examples, however (Figures 10 and 11), it is only the California or European array that are able to separate neighboring sources.

Combined Arrays
The results presented so far have shown just how differently the same sources are imaged using different arrays. Therefore, in order to get a fuller picture of the sources acting in the ocean basin at any given time, it is beneficial to consider the results from more than one array. An ideal combined image would better represent the sources than any of the single images, but for automation of source monitoring, it would also be advantageous if it was only as good as the best single image.
Results of array combination for each of the previous examples using the real data are shown in plot (e) of Figures 4 to 11. For the first example on 13 January 2013 18:00-21:00 (Figure 4), the combined image shows better representation of the spread of the sources than any of the individual array images, with multiple peaks along the northwest to southwest line. In the second example (15 February 2013 06:00-09:00, Figure 5) the northern source seen by all arrays appears in the combined image as expected, with the east-west spread better represented, and the source southeast of Japan which was only observed by the Japan array is also observed. In the third example (10 October 2012 06:00-09:00, Figure 6) the combined image does better represent the three main source areas than any of the individual arrays but still does not have clear demarcations between the three areas. A similar situation can be seen in the fourth example on 27 October 2012 06:00-09:00 (Figure 7) where the two quite distinct source areas visible in the model are not clearly separated, with the beamformer still centered between the two sources as in the individual California and Japan images. In other examples, the combined image offers no improvement over an individual array. For example, in the fifth and sixth examples (Figures 8 and 9) the combined result is noisier and sees no extra sources than the individual Japan array. In the seventh example ( Figure 10) the combined result is noisier and sees no extra sources than the California array, while in the eighth example ( Figure 11) one of the sources observed by the California and European arrays becomes out of range in the combined image because we are not including the region within 15 ∘ of the Japan array center.
Overall, it was found that the combined image generally gave a good representation of P wave sources visible in the single arrays, balancing advantages and disadvantages of each, although in some cases offering no improvement over a single array. This meant that the combined image was sometimes less clear than the best single image, but was able in other situations to pick up sources only visible in one array, or was able to better represent the spread of sources over an area. Often the location quality of a single-array image was reduced after combining with data from stations that were much farther from the source location. This was most visible in the sources close to Japan, which were captured best by the Japan array (Figures 8  and 9). This may be because arrivals from farther away are weakened by geometric spreading and attenuation, resolution decreases with distance, and velocity deviations from AK135 could result in phases not being coherently summed. In order to retain the higher resolution of source locations nearer array centers, we adjusted the weighting w j to enhance the relative contribution of those stations that are closer to the source.

Weighting by Distance
The results of the distance weighting for the data are shown in plot (f ) of Figures 4 to 11. There beam outputs are sharper across all figures compared to the combination without distance weighting. In Figure 5 on 15 February 2013, the weaker source observed by the Japan array is now better separated from the source to the northeast and in Figure 8 on 19 November 2012 03:00-06:00 the reduction of noise means the sources are resolved in the new combined image as well as they are in the individual Japan image. In the worst case, on 13 January 2013 18:00-21:00 ( Figure 4) the combination with distance weighting degrades the output compared to the simple combination. In this case the most northern source, visible originally in the European output, is much reduced due to the European array having less influence as it is farther away.
We also tested whether the location accuracy of the strongest source was improved by combining the arrays. For each event, the location of the maximum modeled source (calculated from the averaged site effect) and the location of the corresponding beamformer peak was found. Because the strongest source may not create the strongest beamformer peak, due to array response functions and noise, differences in the size of the source 10.1002/2017JB015015 area or because the modeled source amplitudes may themselves not be accurate, this corresponding beamformer peak was not necessarily the strongest but was chosen based on proximity to the modeled source location. A very noisy output with many peaks (e.g., Figure 9d) could result in a very close by but low-valued peak being chosen which is not truly representative of the accuracy of the output, so we alleviated this problem by considering only the six beamformer peaks with the largest amplitudes. For each array and combined image, the distance between the modeled source and the corresponding beamformer peak was calculated in degrees (Table S1). Distances ranged between 0.0 and 16.4 ∘ .
The largest distances were found when an array missed a source completely, such as the California array on 10 October 2012 (distance = 12.0 ∘ ) and 19 January 2013 (distance = 8.2 ∘ ) the Japan array on 2 October 2012 (distance = 16.4 ∘ ), or when an array merged two separate sources into one averaged location, such as the Japan array on 10 October 2012 (distance = 9.9 ∘ ) or 27 October 2012 (distance = 10.3 ∘ ). In two of the four cases, the array combinations reduced these distances by better resolving the sources due to information from the other arrays. In two cases however, the distance remained high. In the array combination without distance weighting on 19 January 2013, a large number of multiple peaks meant that the peaks closest to the source did not make the top six peaks. In the array combination with distance weighting on 10 October 2012, the beamformer peaks associated with the northern source became weaker due to the reduced influence of the Europe array, again resulting in them not falling within the top six peaks. The closest beamformer peak chosen 10.1002/2017JB015015 (distance = 12.8 ∘ ) was instead associated with the sources farther south. Overall, the accuracy for each array and array combination were similar (mean distance over events = 5.0 ∘ for California, 5.7 ∘ for Japan, 3.3 ∘ for Europe, 3.6 ∘ for the array combination, and 4.0 ∘ for the array combination with weighting) especially considering our 2-by-2 ∘ grid. These errors are similar to the those expected from the use of a 1-D global velocity model (Euler et al., 2014). The main advantage of the array combination in terms of accuracy is the reduction in times when a source would be missed or merged when using only a single array.

Discussion
Accurately locating P wave microseismic sources is an important step in relating seismic records to oceanic and atmospheric conditions. The relationship has implications for using seismic records to monitor wave activity in real time and track storms across the ocean. Previous studies have shown promise for such applications; P wave microseisms have been located to specific storms such as Hurricane Katrina (Gerstoft et al., 2006), Super Typhoon Ioke (Zhang, Gerstoft, & Bromirski, 2010), and tropical cyclone Dumile (Davy et al., 2014), as well as from regions of high wave activity associated with winter storms (Euler et al., 2014;Gerstoft et al., 2008;Neale et al., 2017).
Geographical coverage of an array is determined by its location, and constraints of array aperture, configuration, and number of stations limit the ability of a single array to resolve and accurately locate microseismic sources. We have shown examples where neighboring sources become merged or sources are missed entirely by individual arrays, whereas other arrays resolve them.
With a view to using microseisms to monitor ocean wave activity and seismic sources, a single image that collates the information from each array would be more suitable than looking at many images from individual arrays. We therefore combined the stations from each array into one image and then adjusted the weighting of each station so that those closest to each grid point were weighted higher. The previous studies already mentioned have combined data from multiple arrays but with notable differences from our study, especially in the temporal resolution. Landès et al. (2010) combined three arrays located in the USA, Turkey, and Kyrgyzstan to locate P wave microseisms at 0.1-0.3 Hz, but with the aim of improving the accuracy of the main global sources over seasonal time scales. They combined their arrays by multiplying the individual output of each array. A comparison of the sources with a wave interaction model was not attempted until the later study by Hillers et al. (2012) and this again looked at sources over 13 day to seasonal time periods. Euler et al. (2014) used four arrays in Africa to locate P wave sources at 0.1-0.2 Hz, but the main difference here was that none of the arrays were deployed at the same time but were each deployed over a 1-2 year period between 1994 and 2007. They manually picked peaks in the output for each array and plotted these onto one map to find persistent sources over the entire 13 year period. They compared these source locations with a map of bathymetric excitation coefficients rather than with sources modeled from a combination of wave-wave interaction and bathymetry. Most recently, Meschede et al. (2017) identified persistent P wave source regions associated with two distinct frequency bands using five North American arrays; however, again the arrays were not all operational at the same time and only a single (the strongest) P wave source on each day was considered. Other studies have focused on higher-frequency P waves using small aperture arrays. Koper et al. (2010) used multiple arrays with apertures of 2-28 km to locate microseisms with frequencies of 0.4-4 Hz. They picked only the maxima of each array output averaged over 1 year and again plotted the P wave locations on one map to find a common source. Pyle et al. (2015) combined arrays in Asia, Australasia, and North America with apertures of 10-25 km to find robust locations of P wave sources at 0.67-1.33 Hz. They used a ranking scheme and histograms to locate sources seen by all arrays over a whole month and compared with maps of significant wave height. Our study builds on this previous work by directly comparing observed sources at three different arrays with modeled P wave sources (instead of significant wave height) over shorter (3 h) time periods, with a focus on how multiple arrays can be used to help reveal or resolve multiple sources that are not always visible from just one array.
In some cases, sources that were visible by a single array were not visible or resolved strongly in the combined image, and east-west separated sources were not well resolved. This was especially true for the Europe array, which had fewer stations and was often overpowered by the California and Japan arrays and was further reduced by the distance weighting. Because the Europe array views the North Pacific from a perpendicular viewpoint with respect to the California and Japan arrays, the influence of the Europe data could be increased by using a weighting that is dependent on azimuth as well as distance. However, because of the stronger 10.1002/2017JB015015 sidelobe artefacts of the European array, this would probably need to be combined with a check in signal to noise, whereby the array is only used in the combination if the signal to noise exceeds a certain threshold.
An algorithm such as the CLEAN algorithm used by Gal et al. (2016), which iteratively removes the array response associated with a point source, may offer further improvement on the separation of weaker sources from the stronger sources. This method is dependent on the dominant source, so by itself cannot correctly separate two closely spaced sources when the beamformer output is centered between the two sources. However, after separating or partly separating the two sources by combining multiple arrays as we have shown here, the algorithm may help separate out and resolve the sources better. Other methods that could be used may include high-resolution beamforming techniques such as sparse beamforming (Elad, 2010), MUSIC (Schmidt, 1986), and MVDR (Capon, 1969) used in ocean acoustics.
The inclusion of more stations and arrays is expected to improve the ability to separate closely spaced sources even further. For this feasibility study, we used only a subset of 202 of the 800+ stations of Japan's Hi-net network. This subset over the Chugoku district was chosen to match that used by Nishida and Takagi (2016) because of weak crustal heterogeneity there, and because the number of stations closely matches the number of stations in the California and Europe arrays. The optimal number and location of stations to use in the process would need to be the subject of future work. In our case the third array we used (Europe) was certainly not optimal evidenced by its much larger sidelobes and because it had a limited range over our area of interest (North Pacific).
Deviations from AK135 travel times will have resulted in power being mapped to incorrect grid points, so improvement of source location accuracy would further be expected if crustal and mantle structure were accounted for in the back projection. Euler et al. (2014) found that accounting for slowness biases caused by 3-D velocity structure corrected source locations by 0-4 ∘ (mostly below 2 ∘ ) which may reduce our current offsets of about 3.3-5.7 ∘ between observed and modeled sources (Table S1). Additionally, we have only been concerned with source location rather than exact amplitude. A better comparison including between amplitudes of multiple sources would be obtained if amplitude information in the seismic displacement spectra was retained. Backprojecting phases such as PP and PKP may provide further information and increased coverage (Retailleau et al., 2018).

Conclusions
We used three large seismic arrays to image P wave microseism sources in the North Pacific on occasions when multiple sources were present at any given 3 h averaged time period. The three arrays used here often gave quite different pictures of the sources acting in the North Pacific at any given time. The Europe array managed to resolve neighboring sources in the very north and those separated in an east-west direction. The California array was generally good at resolving sources in the central regions of the North Pacific but was not good at separating closely spaced sources aligned east-west and had poor resolution for sources in the west of the North Pacific near Japan. The Japan array resolved sources best in the western half of the North Pacific, the east of Japan and north into the Sea of Okhotsk. Therefore, no single array gave the best representation of sources, and a fuller understanding of the sources acting in the basin can be obtained by using data from multiple different arrays placed at different locations around the basin. Most studies, with a few exceptions such as Koper et al. (2010), Landès et al. (2010), Hillers et al. (2012), Euler et al. (2014), and Pyle et al. (2015), currently use only a single array to locate P wave sources and would benefit from including data from other arrays.
We combined the data from each array into one beamforming image. In most cases the combined images were found to give a comprehensive view of sources observed from each individual array. We also investigated a weighting term to boost the influence of arrays closest to the source in order to benefit from their higher resolution at closer distances. This mostly had the effect of sharpening the combined image, but did not perform so well when the European array contained unique information about sources in the north or east-west separated sources, as the European array was farthest from our area of interest. Nevertheless, especially after weighting by distance, the combined images were mostly as good or better than the best single image. The combined images offered improved coverage and revealed source locations that would have been missed or merged by looking at output from only one array, and often better represented the spatial spread of a source.