A revised image of the instrumental seismicity in the Lodi area (Po Plain, Italy)

. We analysed the instrumental seismicity in a sector of the Po Plain (Italy) to define the baseline for seismic monitoring of a new underground gas storage plant that will use the depleted gas reservoir of Cornegliano Laudense, near Lodi. The target area - a square approximately 80 x 80 km wide - is commonly considered aseismic. The analysed period, 10 1951-2019, includes all available instrumental data. We gathered the P- and S-phase readings collected by various agencies for more than 300 events, approximately located inside the target area. We processed the earthquakes uniformly, using absolute location algorithms and velocity models adopted by the regional and national monitoring networks. The relocated earthquake dataset depicts an image of weak and deep seismicity for this central sector of the Po Plain, which is quite different from the initial one derived from the existing earthquake catalogues. Within a distance of approximately 30 km 15 from Lodi, earthquakes are extremely rare (on average 0.5 earthquake/year, assuming a completeness magnitude Mc=2.7 from the 1980s); only 2 weak events fall at less than 15 km distance from the reservoir in the whole period 1951-2019. The strongest events instrumentally recorded are related to the seismic sequence of Caviaga in 1951 that represent the first instrumental recordings for that area. Confirming the hypocentral depths recently proposed by Caciagli et al., 2015, the events are far from the gas reservoir; we suggest common tectonic stress of the main shock of 1951 and the M4.2 earthquake 20 of Dec 17, 2020, based on the similarities in depth, location and focal mechanism. While it

Po river and its tributaries (Ogniben et al., 1975;Bigi et al., 1990). Due to the origin of its deposits, the Po Plain hosts several natural gas (i.e. methane) reservoirs that are mostly located in stratigraphic and structural traps, such as anticline structures, at depth ranging on average between 1100 m and 1800 m; the reservoir and sealing rocks are usually Pliocene sands and Lower Pleistocene clays, respectively (Fantoni and Franciosi, 2010;Boccaletti et al., 2011). The Po Plain is also one of the less seismic hazardous area of Italy (Pagani et al., 2018), but the seismic risk is not negligible (Silva et al., 2018). 35 In addition to the great exposure values, the vulnerability can also be high, due to the fact that seismic provisions in this area have become mandatory only in the last decade (e.g. DGR, 2014;Peruzza and Pessina, 2016), and earthquakes in the magnitude range 5-6 can be devastating, as demonstrated by the 2012 Emilia sequence (Liberatore et al., 2013;Centro Studi, 2014) The largest gas reservoirs of the Po Plain were discovered in the 40's of the past century by the Italian public company AGIP 40 and have extensively been used for gas production since the 50's (Dami, 1952;Bencini et al., 1999). They were the first discovery of large gas deposits in Europe; for instance, the reservoir of Caviaga, with about 12 billions Sm 3 of gas was considered a giant reservoir for the time. Once depleted after production, several reservoirs have been re-utilized as underground gas storage (UGS) to hold the natural gas imported from other sites, or countries. At present (2021), the Po the reservoir depth. Since our study area features rare and deep events, we introduce a further, larger domain, which is called 65 Outer Domain (OD); it is delimited by a spherical cap with 30 km radius, centered like ED. Fig. 2 shows the contour lines of the three domains along a vertical cross-section approximately NS-oriented. The upper crust structure, as interpreted by Martelli et al. (2016) down to a depth of approximately 10 km, is superimposed. Note that the OD has volume and surface corresponding to approximately 8 and 4 times those of the ED, respectively. This paper identifies the study area approximately in the circular region of the OD, i.e. about 30 km distance centered on the 70 UGS facility. For that area, it summarizes the main seismotectonic features and re-analyses the past instrumental seismicity, available from the national and regional catalogues from the early instrumental period (i.e., after the 2nd World War) till 2019. We deliberately excluded from our analysis the data acquired by the new local seismic network of Cornegliano Laudense, implemented in 2016 and operational since 2017 within the UGS integrated seismic and geodetic monitoring system (http://rete-cornegliano.crs.inogs.it/). The "Cornegliano Stoccaggio" network consists of 9 new permanent seismic 75 stations, each equipped by a borehole broad-band seismometer located at about 75 m depth and an accelerometer located at the ground surface. The new integrated monitoring infrastructure will be described in a forthcoming paper, together with the results of the seismic and geodetic baseline analysis and the outcomes of the first years of the UGS activities.
The result of the present study is a revised dataset of instrumental earthquake location that covers about 70 years. Although this catalogue does not feature uniform quality and completeness characteristics, it represents the most accurate and reliable 80 view of the seismicity in the Lodi area available right now. To further strengthen the relevance of our results, we have also included a brief analysis of the most recent, deep M~4 earthquake that occurred on December 17, 2020 near Milan. This study is also intended to be a critical view on the quality of earthquake catalogues in weakly seismic areas, and a stimulus for further studies, in order to improve our knowledge of the deep structures that characterize the collision between the Alps and the Apennines beneath the Po Valley. 85

Seismotectonic framework
Italy is an active tectonic province characterized by two mountain belts: the Alps and the Apennines. They were formed as the result of complex geodynamic processes linked to the progressive migration of the African plate towards the European plate. The two opposite verging fold-and-thrust belts face each other around the Po Plain that geologically coincides with the foreland of both the chains (Fig. 1); it also hosts their foredeeps, and for this reason it has undergone steady subsidence and 90 sedimentation over the past 5 Ma (Fantoni and Franciosi, 2010). The Po Plain represents the north-westernmost buried sector of the Apulian indenter too (i.e. Adria plate, Dercourt et al., 1986). New insights on the tectonic evolution of the Alps and adjacent orogens are becoming available thanks to the AlpArray initiatives (see http://www.alparray.ethz.ch/en/home/, and the special issue at https://se.copernicus.org/articles/special_issue1099.html); despite our study area is quite marginal to the main experiments (e.g. Hetényi et al., 2018), some additional constraints, on the Moho' depth for example, are given by 95 Kind et al., 2021. Global Positioning System (GPS) studies show that convergence at the latitude of northern Italy is still https://doi.org/10.5194/se-2021-35 Preprint. Discussion started: 15 April 2021 c Author(s) 2021. CC BY 4.0 License. going on at 3-8 mm/yr (Serpelloni et al., 2007), the shortening across the Po basin is estimated at 1-3 mm/yr and increases from west to east (Devoti et al., 2011).
The Po valley is formed by alluvial sediments deposited by the Po River and its tributaries coming from the southern Alpine and northern Apennine slopes; they are characterized by uneven thickness -the maximum thickness of approximately 7-8 100 km is located in the eastern Po Plain (Ghielmi et al., 2010) -and the subsidence is asymmetric. Below the alluvial and clastic Cenozoic deposits, the sedimentary sequence includes Mesozoic carbonates and a crystalline basement essentially formed by Hercynian metamorphic rocks. The main fault systems have evolved from extensional (late Triassic-early Jurassic) to compressional regimes (late Cretaceous-Cenozoic), with the inversion of tectonic mechanisms. Triassic-early Jurassic riftrelated structures were locally re-arranged by Cretacic contractional phases, with the reactivation of some of the existing 105 normal faults (Dal Piaz et al., 2004); during Miocene and Pliocene times, the basin became the foreland of the Alps and the Apennine belts and the Mesozoic rocks were deeply buried beneath the Paleogene-Neogene clastics and the associated foredeep wedges. These late sedimentary formations, consisting in alternate shale and sandstone layers, permit the accumulation of hydrocarbons derived from biogenic sources; sandstones are excellent reservoir rocks for oil and gas, thanks to their porosity, whilst shales with very low permeability act like impermeable barriers. Also the underneath Triassic 110 carbonates are great reservoir rocks, as for example in the Villafortuna-Trecate field, 30 km west of Milan. The Po Plain represents one of the major provinces of conventional hydrocarbons in continental Europe (Turrini et al., 2015).
Throughout geological time, both pre-Alpine (Mesozoic and pre-Mesozoic) and Alpine (mainly Cenozoic) orogenic phases have interacted to create the current structural and stratigraphic setting (Carminati and Doglioni, 2012). Three main outer arches have been reckoned for the Northern Apennines, buried below the Pliocene-Quaternary marine and continental 115 deposits; they are respectively from west to east ( Fig. 1): the Monferrato Arch, the Emilia Arch, and the Ferrara-Romagna Arch. On the opposite side of the Po Plain, one wide arc running from Milano to the Garda Lake is related to the Southern Alps (Vannoli et al., 2015). They host active tectonic structures capable of generating earthquakes. In most of the cases, the faults in the Po Valley are blind, i.e., they do not reach the surface and hence cannot be identified by surficial geologicalgeomorphological surveys. According to the Database of Italian Seismogenic Sources (DISS Working Group, 2018), our 120 study area does not intersect any of the seismogenic sources recognized to date. However, within a distance of approximately 10 km from it, there are the composite source ITCS002 and ITCS115 towards the NE, belonging to the outer arch of the Southern Alps (SAOA), and ITCS044 towards the S, belonging to the Emilian arch (EA). They are expected to have a maximum slip rate of 0.5 mm/yr (Maesano et al., 2015) and be capable of a maximum magnitude of 5.5-6.0.
The study area is characterized by weak to moderate seismicity. Fig. 3 shows the main historical and instrumental 125 earthquakes as reported by CPTI15 (Rovida et al., 2020;v. 3.0, Rovida et al., 2021): in the timeline of Fig. 3b during the instrumental period; 2) they occurred in an area that was thought to be aseismic, until that time, and that was the core of international interests on hydrocarbon reserves for several decades; and, 3) for decades, they have been considered anthropogenic, induced events, and they are still included in worldwide compilations of induced seismicity. The May 15 earthquake is yet considered the only M ≥ 5 gas-extraction-induced earthquake reported for Europe (Foulger et al., 2018). In their recent study, Caciagli et al. (2015) relocate the two major events (assigning an estimated Mw=5.4 to the main event of 135 May 15, 1951 and Mw=4.6 to the aftershock of May 16), by using phase readings reported in the International Seismological Summary (ISS), with the use of modern hypocentral location procedures. They provide some bounds for the uncertainty of the location parameters, and estimate the variation of the stress field due to extraction activities. As a result of that study, the epicentral location of the two events is moved far away from the most damaged area of Caviaga toward N-NW, from a minimum of 12 km to a maximum of 25 km distance (in Fig. 3 we show the macroseismic and instrumental location of the 140 the main event, namely the solution quoted as "15b" in Caciagli et al., 2015). The depth of the source is estimated at 32-34 km and 16-20 km for the main event and the aftershock, respectively, with statistical errors of 3-5 km; these are then locations compatible with faults inherited in the rocks of the crystalline basement. About the focal mechanism, the recent study does not suggest fault planes, as done by some previous analyses (e.g. Ritsema, 1961;Gasperini et al., 1985). The stress perturbation induced by the coeval gas production is estimated to be far below any value usually taken as triggering 145 threshold. Therefore, today, this earthquake is attributed to a tectonic origin, and no more to anthropogenic causes: the two 1951 events are interpreted as deep crustal events, representing active sources at the contact between the Northern Apennines and the Southern Alps beneath the Po Plain.

Data revision
The initial phase of the work was the collection and comparison of instrumental earthquake catalogues that are available for 150 the area, namely instrumental earthquake catalogues released by either the main Italian institutions that manage seismic networks, or international agencies, or dedicated projects. A short description of the data providers and past initiatives is provided hereinafter.

•
The National Institute of Geophysics and Volcanology -INGV: the "Centro Nazionale Terremoti" (CNT) is the authoritative centre for earthquake monitoring in Italy. Under contract with the National Civil Protection Department, it 155 supplies the parameters of the epicentre in real time, in order to organize emergency interventions. CNT manages the National Seismic Network (RSN -Rete Sismica Nazionale), established following the catastrophic earthquake of Irpinia in 1980. The first nucleus of the network was composed by seven stations linked to an unique elaboration centre in Rome.
Today CNT manages about 500 stations and acquires the signals of several other institutions, for a unified detection, localization and magnitude determination of seismic events, and for managing the archive of recorded events, available at 160 the INGV website (http://terremoti.ingv.it/en/iside). The earthquake parameters are transmitted to international agencies too (see Margheriti et al., 2021 and references therein Accelerometrica Nazionale), which is the network dedicated to the measurement of the ground acceleration generated by moderate-to-high magnitude earthquakes. The network was acquired from the Italian electric company ENEL in 1997 and 165 over the years has been improved and upgraded to digital instruments. As of today, it is composed of almost 600 seismic stations with 3-axial accelerometers, concentrated especially in high seismicity zones. These data are complementary to the seismic monitoring done by INGV, but they are usually not taken into account for the hypocentral and magnitude determination, except in case of signal saturation of seismometers, in the near field of major earthquakes.

•
The University of Genoa -UNIGE: the Regional Seismic Network of Northwestern Italy (RSNI -Rete Sismica 170 regionale dell'Italia Nord-occidentale) belongs to the Laboratory of Seismology of the University of Genoa. It started with a single instrument at the beginning of the '60, and underwent slow developments till the '80, when the number of instruments, azimuthal coverage and transmission links became adequate to a regional network (Eva et al., 2010;Spallarossa et al., 2014). It is currently composed of 32 seismic stations useful both for monitoring and for scientific research purposes.
It releases online bulletins of automatic and revised located earthquakes for Northwestern Italy 175 (http://www.distav.unige.it/rsni/). Wood Anderson, for magnitude assessment -are still functioning, Sandron et al., 2015) were available for some decades too, and were partially integrated in the data collected by the international agencies.
The RSNI and NEI regional networks are of interest as some of their stations are close to the study area and can provide additional data to the ones detected by the national network.
The International Seismological Centre (ISC) is worth to be mentioned as well: it was set up in 1964 as a successor to the 190 International Seismological Summary (ISS), the pioneering effort in collecting, archiving and processing seismic station and network bulletins and preparing and distributing the definitive summary of world seismicity. ISC does not have its own instrumental infrastructures, and its mission is to maintain a number of important data repositories, as: the longest continuous About the Italian national earthquake catalogues, they were the outcome of dedicated projects funded within the activities of the Gruppo Nazionale per la Difesa dai Terremoti (GNDT). In the 80s, a specific project (Progetto Finalizzato Geodinamica) released the first catalogue of Italian earthquakes from 1000 to 1980 (Postpischl, 1985), mainly based on macroseismic observations and early instrumental data. Then, the very first initiative of a quality controlled and unified earthquake 200 catalogue for instrumental seismicity in Italy was promoted after the strong earthquakes of Colfiorito and Sellano in Central Italy, in 1997-98; a sub-project named "5.1.3 Catalogo Strumentale dei Terremoti" was approved and funded by GNDT. The objective of the project was the complete collection of phase readings from all the various networks and observatories existing in Italy at that time, and the uniform relocation of events, according to agreed procedures and standards. The We started the work by collecting the earthquake locations provided online and offline by the previously mentioned organizations, and falling within the OD (i.e., within 30 km distance from the reservoir). A preliminary working file has been 220 compiled by stacking the hypocentral solutions provided by INGV, UNIGE, OGS and ISC in the period they cover from 1980 to 2019. It was immediately clear that some events are present in one database only, while other earthquakes have many different solutions with large scatter in epicentral coordinates. To account for such differences, we did other searches over a larger area (Target Area, TA, white frame in Fig. 5, that envelopes the OD), and/or by individual earthquakes. The composite catalogue retains only one location for each event. Multiple locations are believed to belong to the same 225 earthquake if the origin times differ less than 20 s. In this case the following priority rules are given to the data providers: INGV>UNIGE>OGS>ISC. If an event is reported inside the TA by one agency alone we retain that location. About the OGS data (available at rts.crs.inogs.it), it is worth mentioning that a very limited number of events are located inside the target area, as it is out of the NEI authoritative region; nonetheless, phase readings for non-located events are stored in the bulletins, and have been used in the later steps. This preliminary analysis motivated us in performing a global check and relocation of events for this sector of the Po Plain.

Toward a revised instrumental catalogue
After the identification of about 340 earthquakes that potentially insist in the study area (given by the merge of catalogues previously described), we started the collection of the data needed for a uniform reprocessing of earthquake locations, i.e. Pand S-phase readings and stations data which change throughout time, and are only partially available on the web portals of 250 the data providers.
The compilation of an updated and reliable list of seismograph stations has been done in several steps, according to the needs emerged while processing the relocations; we first describe these results because they help understand the evolution in time of the seismic networks, as well as give an approximate idea of the detectability and reliability in the location of medium-tosmall magnitude events we are dealing with, in the study area. 255 The map in Fig. 6 is a simplified graphical representation of what is usually called "station book", i.e. the information about location, instrumentation installed, site response and operational conditions (ON/OFF times) usually managed by the seismic station owners. When a seismic station has a permanent installation, it can be registered with an unique code to international For the stations in grey we did not gather the details of operational conditions. The lack of stations in the Po plain is evident, especially before the mid 1990s. Around the reservoir area, this lack is confirmed at least until 2010-2015, if we exclude the Pavia station (PAV, not visible as overlapped by the symbol of the nearby EUCT station), closed presumably before 1985.
It is therefore reasonable to state that the detectability threshold of earthquakes in this area is worse than that in adjacent 270 sectors of Northern Italy, which conversely have station coverage dating back to the 80s. In addition, the completeness magnitude is strongly influenced by the level of background noise at the seismic stations, critical in this region too.
To quantify the noise level in the study area, we have analysed the continuous signal recorded at the ground surface by a set of temporary stations (13 sites) we deployed in July 2016, in the ED area, during the scouting phase carried out for designing the new monitoring seismic network for the Cornegliano Laudense UGS. In addition to our temporary stations, we have also 275 analysed the signal of the INGV permanent station of Milano (MILN). We have to keep in mind that all those stations are located in a densely populated area -MILN is in the centre of Milan -which hosts several mid-to small-size cities and industry compartments. Fig. 7 shows the median noise curves (McNamara et al., 2004) relative to 1 day of continuous recording measured at the 13 temporary sites (thin blue lines) and 1 week of signal at the MILN station (thick red curve).
The figure also shows the two curves known as New High Noise Model (NHNM) and New Low Noise Model (NLNM) 280 (lines in dark grey), respectively, obtained by Peterson (1993) from the analysis of the seismic microtremor data of 75 stations distributed in different parts of the world and generally used as a reference for estimating the quality of the signal recorded by a seismological station. For the temporary stations represented in Fig. 7, the noise curves are truncated at the period of 4 s, as the seismometer used is a Lennartz 3D-Lite with resonant period T = 1 s. We see in Fig. 7 that all the recordings carried out at the ground surface feature a seismic noise level generally near or higher than the NHNM line in the 285 short-period range.
The impact of the background noise on the minimum detectable magnitude has been evaluated by Franceschina et al. (2015) for a site located in the Po Plain within a pilot project for CO 2 geo-sequestration. In their study, the authors use both experimental data and theoretical arguments to determine the noise levels of the area, the detection capabilities of a future, possible monitoring seismic network, as well as the improvement of borehole vs. surface installations. Based on a 1-year 290 observation period, they show that the Po Plain features high seismic noise levels, although slightly smaller than our ones.
They also estimate in ML 1.9 the minimum detectable earthquake magnitude for an event-station distance of 20 km.
If we recall that the average inter-station distance existing until 2016 in our study area is greater than 30 km, we can reasonably move the detection magnitude limit about half a unit higher (ML~2.5). As the completeness magnitude (Mc) of the catalogue is higher than the detection magnitude -an event has to be detected by at least 3-4 stations for being located -295 we can therefore easily understand that the poor completeness magnitude value of the existing catalogues is the outcome of the high seismic noise existing there, combined with the low network density. Additional problems to this topic derive from https://doi.org/10.5194/se-2021-35 Preprint. Discussion started: 15 April 2021 c Author(s) 2021. CC BY 4.0 License. the use of different standards for magnitude assessment in time, problems that we will not tackle, in this paper and leave for a further analysis.

Collection of phase arrival times 300
As regards the phase arrival times, for each earthquake we gathered the following resources:

1)
the extended listing of data coming from ISC Bulletin. Each event is represented by a text file that summarizes the hypocentral solutions provided by various agencies and the PRIME one adopted by the Bulletin, the infos about magnitude, and then the list of station, phase type, arrival times, with additional metadata that concern the preferred solution (e.g. distance, azimuth, residual, station magnitude). 305 2) the listing of data coming from the INGV website, provided by tables of arrival times for each event (e.g. http://terremoti.ingv.it/event/EVENTID). These data are available jointly with the solution provided for the location, for all the events since 1985, with uniform format but different criteria adopted, for example, for assigning the uncertainties and magnitude in the various years.
3) the phase data provided by UNIGE: some phase data are available online, but an extended search on the archives 310 has been performed (Scafidi, pers. comm.), for including also non located/badly located events, that may be out of the authoritative area of the RSNI network. Readings are available in different formats, depending on the years.

4)
the phase data of OGS are available in a uniform format in the online CRS Bulletins (http://www.crs.inogs.it/bollettino/RSFVG/RSFVG.en.html) from 1977 on; in many cases, these data have never been used before in location procedures, because the events that were manually picked and reckoned outside the OGS authoritative 315 region of monitoring have not been processed till the hypocentral solution. Unlocalized phase manual readings are available in particular during the so-called "analogic" period (see Gentili et al., 2011) and after the improvement of the Veneto Region monitoring (from 2005 on; see Bragato et al., 2021). Other unpublished data derive from the seismic monitoring of the Minerbio site that was performed by OGS on behalf of A.G.I.P. S.p.A. (Rebez, 1991), or have been searched on original seismograms on paper, available in the Sgonico headquarters of OGS. Note that the readings of Trieste station (station code 320 TRS from 1904 to 1964; TRI from 1959 on) are usually available on the ISC archives, but for the earliest events they have been checked on original documents too.
Each earthquake is given an unique event identifier (OrigID), corresponding to the one assigned by ISC (if the event is present in the international Bulletin), or given by a string addressing to the year, month, day, and eventually hour minute of first arrival time, in case of multiple events during the same day. Considering the differences in data format, and reliability of 325 the final solutions, we processed the earthquakes in several time frames, as given in which no phases are listed. In practice the online dataset of arrival times begins with the aftershock of 1951: for the main 330 event, phase readings have been integrated in this study too (Caloi et al., 1956). Table 1 summarizes the number of phase arrivals by stations, as provided by the various agencies, and the total number of events processed, located and discarded. Note that some phase/station data used by INGV and UNIGE can be duplicated, and they both are usually provided to ISC Bulletin. OGS data are mainly original, as non-located phases are not provided to the international agency. ISC, collecting data worldwide, is the most represented dataset. The earthquakes not reported in 335 ISC Bulletin are usually the smallest events. Concerning magnitude, as in this study it is needed mainly for representation purposes, we avoided conversion relationships, and retained the values assigned by the original agencies, selected on the basis of priority criteria (Mw>Ms>Ml>Md). As previously said, the uniform re-assessment of magnitude is a relevant problem that falls outside the schedule of this study, and it will be the focus of future efforts.
The phase data coming from the various providers have been then transformed into a common format, and merged, for their 340 uniform processing by the selected location code, by means of a suite of original Bash scripts and Python codes, implemented for this study.

Relocation of selected events
Hypocentral location is performed by a standard code used in regional monitoring analyses, i.e. the software Hypo71 (Lee & Lahr, 1975) that uses simple 1D regional models, with fixed value of Vp/Vs. The choice of a robust method, although less 345 sophisticated than other, more recent ones, is motivated by the need of accommodating the different quality of the phase readings collected in more than 50 years and of treating all data by a common location algorithm, thus avoiding artefacts introduced by different relocation procedures. The complexity of finding robust solutions in an inverse problem, like hypocentral location, is well known in the literature: in their very recent paper Arrowsmith et al. (2020) state the superiority of global search methods which sample the full model space, on standard iterative linear least-squares approach, as originally 350 introduced by Geiger in 1912, and implemented by Hypo71 and its descendants. However, when dealing with dataset derived by the "merge" of phase readings published by several agencies throughout many decades, like the dataset of this study, there are other problematic features and uncertainties to deal with, such as: • Duplication of readings: by merging the phase readings used by different agencies we may have duplications of readings for P and S phases at a station, which can report the same, or different arrival time; 355 • Unknown or not proper phase readings; • Assessment of uncertainty given to the readings; • Differences in the phase association performed by different agencies.
We decided to face these problems as it follows: • As we have no "a priori" element to establish which are the "best" readings, in the case of multiple arrival times 360 assigned to the same station/phase, we decided to let the duplications in. In theory, the duplication of the same equation in a linear system does not affect its solution; conversely, slightly different observations may "accommodate" the uncertainty of https://doi.org/10.5194/se-2021-35 Preprint. Discussion started: 15 April 2021 c Author(s) 2021. CC BY 4.0 License. arrival times that we cannot verify on the original seismograms/time series. The evaluation of relevant discrepancies among data given for the same station is done by checking the printout of the locations, thus assigning different weights (e.g. 4, that means discard the reading, or 9, that means use the S-P time only) to the high residuals (more than 2 s for P-phases, 4 s for 365 S-phases). We thus consider the resulting statistical errors a relative measure of the goodness of the solutions within a uniformly processed dataset, and not a physical reference for the point source approximation of the hypocentre. The literature on these subjects is wide (see e.g. Husen and Hardebeck, 2010;Arrowsmith et al., 2020) and some additional comments will be given in the final sections.

•
The location algorithm uses the first arrival phase, but phase types are often not explicitly indicated by the data 370 providers (it happens mainly in the ISC dataset, and for the last decades only); we decided to exclude far stations (distances larger than 400 km, after some sensitivity tests) and improper phase readings, if explicitly labelled (e.g. Pn, Pg, P*, etc.), and to retain unknown phase types as alternative arrival times, as previously mentioned. The removal of data is driven by an "a posteriori" check on the residuals of the theoretical travel times, and giving priority to the proprietary data (i.e. readings for the stations managed by the agency). 375

•
The assessment of a weight to the picking on the basis of its reliability has no viable solution, except than that of resorting to the original seismograms. Especially in the past, it was a common practice to assign the weight without declaring the assessment criteria, and mainly referring to the residuals on theoretical travel times. We decided to adopt in the first runs the weighting codes given by the data providers, if assigned, and to set standard values (0 -full weight for P readings; 1 for S) if not given. The careful analysis of each solution led to a revision of the weights based mainly either on outliers or 380 readings that drive the algorithm to converge towards solutions with anomalous root mean squares time residual (RMS), reckoned by a generalized increase of "acceptable" residuals (less than 2 s) in most of the stations.

•
Anomalies in phase associations have been found in very few cases, for cascading events. Sometimes the revision of associations led to the analysis of other events not selected at first according to the spatial rules.
After the station list and phase arrival dataset, the last but very important component of hypocentral locations is the velocity 385 model.
We resorted to a set of 1D models already published (see Table 2), used respectively at global scale (IASPEI91; Kenneth and Engdahl, 1991), national (INGV RSN1 and RSN2, see https://emidius.mi.ingv.it/CSTI/Versione1_0/Html/LOC.HTM and Marchetti et al., 2020;CSI, Chiarabba et al., 2005, modified to accommodate the inversion in velocity profile not allowed by the location code), and regional ones (UNIGE, Spallarossa et al., 2014, adopted for a study in northwestern Italy, then 390 furtherly simplified and modified to accomplish the large interstation distance of our dataset; OGS, Bressan et al., 2003, used in northeastern Italy). We selected the solutions obtained by the model that provides the minimum time residuals (in bold in Table 2) for the final map and profiles given in Figs. 8 and 9, acknowledging the huge variability in location parameters (mainly the value of depth), that results from the different velocity models, coupled with non-optimal station distribution and large uncertainties in phase pickings. Note that the selected model (n. 4 in Tab value that is a reasonable average of lateral strongly variable depths sampled by teleseismic waves in Kind et al., 2021. Additional comments will be given in the results.

The Milan earthquake of Dec 17, 2020
On December 17th, 2020 at 15:59:22 UTC an earthquake occurred near Milan, right while we were writing this paper. The Italian agencies INGV, UNIGE, and OGS notified it in nearly real time on their institutional channels and social media, as a 400 deep event; the revised hypocentral determination given by INGV and UNIGE (http://terremoti.ingv.it/ and http://www.distav.unige.it/rsni, last accessed on 10 March 2021) assigned a local magnitude ML 3.8 and 4.2, and hypocentral depth of 55 km and 49.6 km, respectively. The event was widely felt in the Milan province, and somehow it surprised the population, convinced to live in an aseismic area. No focal mechanism has been released by any Italian or foreign agencies up to now. Since the event occurred right within our study area, we decided to pick manually the arrival 405 times and first polarities on available waveforms, to compare this earthquake with the results of our analysis. We selected the seismic stations within the cross-over distance for the area (130 km according to Spallarossa et al., 2014)  Mediterranean Very Broadband Seismographic Network -MN), from international to regional scale, equipped with seismometers or accelerometric sensors. We picked 45 P (with 31 polarities) and 37 S arrival times, weighted according to the following scheme: code 0 for uncertainty in picked time Upt ≤ 0.1 s; code 1 for 0.1< Upt ≤ 0.5 s; code 2 for 0.5< Upt ≤ 1.0 s; code 3 for 1.0< Upt ≤ 2.0 s; code 4 for Upt > 2.0 s, not used in the location process.
We used both Hypo71 to Hypoellipse (Lahr, 1999) to locate the Milan earthquake, as the latter allows us to take into account 415 the station elevation and to use velocity models featuring Vp/Vs variable with depth, or velocity inversions. We performed location tests using different velocity models: four models are those already described in section 3.1.2 (i.e. IASPEI91, INGV RSN2, UNIGE, OGS), the fifth one is the well-known AK135 model (Kennet et al., 1995) used by Caciagli et al. (2015) for relocating the 1951 Caviaga event too. Table 3 lists the hypocentral parameters obtained for each location. The most evident difference deals with the hypocentral 420 depth that ranges from 50.8 km to 59.2 km, while all the other parameters are similar: no "best location" can be defined, in our opinion, on the basis of statistical residuals. Therefore, in the followings, we adopt the hypocentral solution given as 4 in Table 3, to be coherent and uniform with the relocated earthquake catalogue.
We then calculated the fault-plane solution based on first-arrival polarities, by using the classical code by Reasenberg and Oppenheimer (1985). All the focal mechanisms obtained are represented in Table 3: note the irrelevant changes in plane 425 orientation due to the differences in hypocentral location. Note also that the solutions obtained with the same model (OGS, n. 4 and 6) are invariant with respect to the location code. A transtensive mechanism on NNW-SSE to E-W oriented structures is common to all the models: in two models (n. 3 and 4 in Table 3) an alternative compressive solution is given too, on a https://doi.org/10.5194/se-2021-35 Preprint. Discussion started: 15 April 2021 c Author(s) 2021. CC BY 4.0 License. low-angle thrust fault dipping to NW, or a nearly vertical N-S oriented fault. For the preferred solution, the Station Distribution Ratio (STDR) is always >0.63, suggesting a well-constrained result. We recognize that the data collected in this 430 study do not allow us to overcome this uncertainty, and a moment tensor analysis might be useful to obtain more robust solutions. At least a dozen of events inside the target area are still to be considered fake locations, due to improper minimum search, as it happens when the hypocentral coordinates move far away from the stations that recorded the event. The earthquakes represented in the inset, for example, are likely Swiss or French earthquakes, as testified by their first arrivals phase data.

Results
These events have been identified not exclusively considering their horizontal and vertical standard errors, but also by means of the minimum distance of the first recording station (DMIN), a parameter that usually is not taken into account in the 445 quality analyses, or it cannot be taken into account when the catalogues do not report the full solution parameters. These considerations explain the reason why in the composite catalogue (Fig. 5) we found at least two events not reckoned by national or regional agencies: it is a first warning on the uncritical utilization of earthquake catalogues in low seismicity, and badly monitored areas. Then, we removed the earthquakes with DMIN>70 km, adopting a threshold that corresponds to the 95 percentile of DMIN distances for the events in the target area: the remaining earthquakes are plotted in the cross sections 450 of Fig. 9, and in Fig. 10.
The distribution of relocated earthquakes now shows a more evident separation between the Apennines and Alpine domains than before the re-analysis. Earthquakes are in general deeper than in the composite catalogue (Fig. 5), especially towards the Apennines: if the "absolute" depths are confirmed by further studies, many earthquakes can be assigned to the lower crust (>15 km depth). For the main event of 1951 we acknowledge the consistency of our location with one of the two 455 solutions obtained by Caciagli et al. (2015), even if the ingredients and algorithms used are different. In our study, in particular, phase readings were taken from the original paper of Caloi, instead of the ISS Summary -this compilation is not included up to now in the ISC database -; our relocation uses a regional model rather than a global model; we retain only stations up to 400 km distance, discarding the worldwide stations used by Caciagli and coauthors, and uses regional  Caciagli et al. (2015) for this event; we consider the aftershock hypocentral solution poorer than the main event' one, despite the better statistical errors (see Table 3). The focal mechanism, obtained from first polarity data in this study (drawn in Table 3, too), shows multiple or poorly constrained solutions: for the main event, the compressive solution is similar to the mechanisms already published (e.g. Martelli et al., 2016), but an alternative, statistically preferable, extensional mechanism is suggested, on NW-SE to WSW-ENE oriented fault planes with 465 moderately dipping angle. Note however that the retrieved focal mechanisms do not help to fix the uncertainty in aftershock location and fault identification, or to associate the events to any of the seismogenic sources identified in the literature till now. The May 15, 1951 mechanisms, although not well constrained, have common elements with the 2020 ones, i.e a comparable unusual depth, the same ambiguity between two contrasting solutions (transtensive vs transpressive mechanism), and the maximum stress σ1 axis orientation (indicated by P in the plots in Table 3). Surprisingly, the σ1 is compatible with a 470 NW migration of Adria plate, more than to the NNE-SSW oriented convergence of the Alpine and Apennines systems.
These very preliminary considerations have to be confirmed by further analyses (e.g. full moment tensor, reflected and converted phase modelling at the Moho) that could probably help constraining the rupture location, and reducing the uncertainties.

475
Very few earthquakes are located within the 30 km distance from the reservoir (dashed black circles in Figs. 8 and 9). The subset of events that are located inside the target area and have DMIN<70 km are plotted in the double-Y graph of Fig. 10, where the blue dots represent the hypocenter-to-reservoir distance (left axis), and red crosses the magnitude assigned by the origin agency, both plotted versus time. Only two events are located within the 15 km distance (ED domain, according to the Italian Guidelines; blue half-circle in Fig. 9), i.e. the earthquakes M2.7 of Feb 2, 2002, andM3.0 of Dec 24, 1996, 480 respectively: although their location feature small standard errors, in both the cases the relevant distance from the first recording station (DMIN approximately 53 and 60 km) makes the reliability of the estimated depth still poor. Just another event is located nearby, that is the M2.7 earthquake occurred on March 25, 1983, at about 17 km distance from the reservoir (DMIN=33 km). Fig. 10 the decrease of the minimum magnitude threshold versus time, which corresponds to an increase of the 485 network sensitivity, and the unavailability of magnitude data for some events. We consider this picture a useful representation of the weak level of seismicity that characterizes the area nearby Lodi, as it: 1) helps in driving the future investigations towards individual earthquakes, or specific sequences; 2) permits a 4D representation of the volumes surroundings the investigated site, and is useful for comparison purposes too; 3) provides a visual perception of the completeness magnitude, as a statistical treatment cannot be applied to such a small and not uniform sample. 490

Note in
This study adds some new pieces of information to the knowledge of the seismicity of Po Plain area around the Lodi town, however we think that only an improvement of the existing observational capability will allow us to interpret which are the network that has been recently accomplished for the Cornegliano Laudense UGS will surely contribute significantly to that purpose. 495

Conclusions
The uncritical use of existing catalogues of earthquakes can lead to biased representations of the seismicity, especially in weakly seismic or poorly monitored areas. This is the case of the Lodi area, in the central sector of the Po Valley.
The collection of instrumental data allowed us to identify more than 300 earthquakes from 1951 to 2019 (Fig. 5). Through their re-processing, we have located nearly 200 events within a target area (orange rectangle in Fig. 8) defined around the 500 Cornegliano Laudense gas storage. The hypocentral solutions are obtained with a standard absolute location method that, even if non-optimal, is a suitable choice to perform a uniform processing for the whole period. Depths are influenced by the velocity model: the one we adopted for the final representation leads to the best statistical performances and depths that are compatible with the ones obtained by other recent analyses and in agreement with the results obtained for the very last earthquake occurred in 2020. The location errors are relevant for the early instrumental period and decrease in the most 505 recent years. While we recall that the statistical error cannot be considered a physical measure of hypocentral uncertainty, we underline that, in combination with a careful analysis of each event, and other solution parameters (e.g. the minimum eventstation distance, DMIN), they help detect fake locations and permit a comparative evaluation of the hypocentral solutions.
This study also shows that the detection capability in Northern Italy increased after the 1976 Friuli earthquakes, and then in the early 80s, mainly for the improvements of the regional networks; then again some betterment, with a decrease in 510 minimum magnitude threshold, is observed in 2008-2010, probably due to the growth of the national networks that followed the 2009 L'Aquila earthquake. As we did not check the magnitude assessment, however, these data are not uniform: no conclusive remarks can be drawn then, on completeness magnitude versus time.
The seismicity in this sector of the Po Plain, where the Northern Apennine and Southern Alps collide below a thick sedimentary cover, is weak and deep. This picture is quite clear, although the area is not monitored with comparable quality 515 of the adjacent regions. With regard to the Lodi area, and in particular to the Outer Domain (30 km distance, see Fig. 2) of the Cornegliano Laudense storage, the earthquakes identified in the last decades are extremely rare (on average 0.5 earthquake/yr, if we assume a completeness magnitude Mc~2.7, from the 80s). In the whole period 1951-2019, only 2 events fell within the volume with radius 15 km from the reservoir (External Domain); the reliability of their hypocentral location can be furtherly investigated by other location methods, given the limitations of the available recordings. The first 520 instrumentally sufficiently recorded events are related to the seismic sequence of Caviaga, in 1951; they represent also the strongest events occurred in that area. Our study confirms the elevated depth and large distance of these events from the Cornegliano Laudense reservoir, as already proposed by Caciagli et al. (2015). We further suggest a similarity between the focal mechanisms of the M5.4 earthquake of May 15, 1951, and  Further studies for a better knowledge and interpretation of seismogenic sources in this strategic area of the Po Plain are ongoing; the quality of the revised catalogue can be increased, for example, by approaching the hypocentral location with 530 full grid search methods, by resorting to original seismograms in some cases, and especially by facing in a uniform way the assessment of magnitude values, not tackled in this work. A better insight will be gained by improving the existing observational capability, and, in this respect, the new seismic monitoring network of Cornegliano Laudense UGS will surely contribute significantly to these purposes.     Table 3   * they refer only to the readings provided in Caloi et al. (1956) for the event of May 15, 1951, not given in the ISC archive.