Time evolution of TGF-producing storms using ERA5, GPS and geostationary satellites observations

In this article we report the first investigation over time of the atmospheric conditions around TGFs occurrence, using GPS sensors in combination with geostationary satellite observations and ERA5 reanalyses data. The goal is to understand which characteristics are favourable to the development of these events and to investigate if any precursor signals can be expected. A total of 9 TGFs, occurred at a distance lower than 45 km from a GPS sensor, were analysed and two of them are shown here as an example analysis. Moreover, the lightning activity, collected by the World Wide Lightning Location Network (WWLLN) was used in order to identify any links and correlations with TGF occurrence and PWV trends. The combined use of GPS and the stroke rate trends identified, for all cases, a recurred pattern in which an increase of PWV is observed on a timescale of about two hours before the TGF occurrence that can be placed within the lightning peak. The temporal relation between the PWV trend and TGF occurrence is strictly related to the position of GPS sensors in relation to TGF coordinates. The life cycle of these storms observed by geostationary sensors, described TGFs producing clouds as intense with a wide range of extensions and, in all cases, the TGF is located at the edge of the convective cell. Furthermore, the satellite data give an added value in associating the GPS water vapor trend to the convective cell generating the TGF. The investigation with ERA5 reanalyses data showed that TGFs mainly occur in convective environment with not exceptional values with respect to the monthly average value of parameters measured in the same location. Moreover the analysis showed the strong potential of the use of GPS data for the troposphere characterization in areas with complex territorial morphology. This study provided indications on the dynamics of convective systems linked to TGFs and will certainly help refine our understanding on their production highlights a potential approach through the use of GPS data to explore the lightning activity trend and the TGFs occurrence. Keypoints: First use of PWV derived by GPS measurement to characterize TGF-producing storms; ERA-5 confirms TGFs are produced in convective but not extreme environments; An increase of PWV is observed on a timescale of about two hours before the TGF occurrence; TGF is located at the edge of the convective cell. Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 24 December 2020 doi:10.20944/preprints202012.0622.v1 © 2020 by the author(s). Distributed under a Creative Commons CC BY license.


Introduction
In 1994 a surprising observation by the National Aeronautics and Space Administration (NASA) Compton Gamma-Ray Observatory (CGRO) detected unexpected gamma-ray emissions coming from the Earth [1].
These so-called Terrestrial Gamma-ray Flashes (TGFs), produced inside storms in association with lightning with typically durations of less than 1 millisecond and energies up to few tens of megaelectronvolt [2], are the manifestation of the most energetic natural particle accelerators on Earth, strong enough to be observed by high sensitive instruments orbiting in space.
During decades this became a topic of frontier research between two disciplines: high energy physics and atmospheric physics. The first investigated the mechanisms of production, discovering that TGFs are produced by a large number of charged particles accelerated within thunderstorm and lightning intense electric fields, undergoing avalanche multiplication and subsequently emitting gamma-ray photons via bremsstrahlung [2].
On the other hand, to understand under which boundary conditions this phenomenon can trigger and propagate in air, or how rare it is, we need to involve atmospheric sciences.
Understanding the rarity, formation and evolution of these atmospheric phenomena is important in assessing the risks to which we are subjected: [7] pointed out that in the TGF production area the radiation levels are high enough to compromise health as well as electronics on board aircraft.
From the meteorological point of view, TGFs are produced by storms of all shapes and sizes, but it is still unknown why some thunderstorms produce gamma-ray bursts and others do not.
At the light of the state of the art and knowledge gaps in the understanding of correlation between TGFs and meteorological phenomena, the present study aims to investigate and to describe the atmospheric characteristics in the TGF producing both using in situ and satellite data.
The scope is to understand which conditions are favourable to the development of the TGF events and to investigate if any precursor signals can be expected.
To do that, TGFs from italian Astro-rivelatore Gamma a Immagini Leggero (AGILE) satellite, [8-10], were considered. First meteorological conditions in correspondence of TGF occurrences are studied by ERA-5 reanalysis and compared with reference values in order to evaluate the meteorological conditions leading to storms emitting TGF.
Second, the Global Positioning System precipitable water vapour (GPS-PWV) was considered as a good indicator of moisture content and matched with the strokes registered by WWLLN as well as informations obtained from geostationary satellites, allowing the time monitoring of the meteorological conditions preceding the TGF occurrence.
Several studies in the last decade had as objective to find correlations between TGFs and meteo/lightning characteristics of the associated events.
10-14 km altitude range, indicating deep convection, gave a much better match. The authors found that storm systems of all sizes could produce TGFs, showing a range in areal extent by several orders of magnitude.
[12], using radio atmospherics data from the World Wide Lightning Location Network (WWLLN), made the first analysis on the electrical evolution of a storm related to a TGF event finding a clear decline in flash rate surrounding the TGF occurrence, suggesting that TGFs occur preferentially during the declining phase of flash production.
In contrast, [13,14] showed that TGFs tend to take place during the peak of the cooling phase, when the lightning flash rate is at its maximum.
Detailed meteorological observations over 24 TGFs detected by Fermi were given by [15]. They compared the CAPE value at the TGF occurrence with the minimum, mean, and maximum CAPE value registered all days of the month at the same location and time of day The Convective Available Potential Energy (CAPE) values were calculated by the NCEP North American Regional Reanalysis (NARR) dataset [16], with spatial and temporal resolutions of ~2 km and 3h, respectively.
This study showed that the 24 TGFs originated from storms of a wide range of convective strengths, without any clear common characteristics. These results are confirmed also by the recent study of [14] that linked TGFs production to cloud instantaneous and dynamical features as extracted by visible-infrared geostationary satellite sensors.
On the other hand [17], taking advantages from satellite borne radar onboard the Global Precipitation Measurement (GPM) mission [18], performed an analysis over TGF events from a new perspective. A total of 9 TGF producing storms were analized both with active and passive instruments, finding common features: both views agree in describing TGFs producing clouds as intense thunderstorms with significant vertical development and ice content and 100% of the cases presented a cumulonimbus tower. Moreover all TGF related lightning are classified as high amplitude IC flashes.
Furthermore, as TGFs are related to lightning, it is essential to try to correlate the issue of the lightning initiation within clouds and its microphysical content processes.
Electric breakdown occurs when the electric field surpasses the so-called classical breakdown threshold Eth ~3 MVm -1 (at standard temperature and pressure). However, also in sophisticated balloon measurements, ambient field were always lower than 20% of classical breakdown [19][20][21]. The lightning initiation problem is thus the problem to start an electric discharge, when the ambient field is far below the breakdown threshold. The problem is constrained to only the start of the discharge, in the electron avalanche phase where electrons and ions grow exponentially up to the point where they produce a space charge electric field above breakdown and the field enhancement of the dielectric is not necessary anymore.
Any cloud electrification mechanism involves a small-scale process that electrifies individual hydrometeors and a process that spatially separates these charged hydrometeors by their polarity.
There are still major uncertainties about how many collisions actually occur in different regions of the cloud, what the crystal sizes and collision velocities are in these regions, what the temperatures and liquid water contents are in these regions, and what charges are actually present on the various hydrometeors throughout the cloud.
Studies aimed at obtaining such information using in situ measurements often in conjunction with radar and other observations have been conducted in various types of clouds [22][23][24]. It is possible that the primary electrification mechanism changes once a storm becomes strongly electrified. For example, collisions between ice crystals and graupel could initiate the electrification, and then the larger convective energies of the storm could continue it. It is also possible that important electrification mechanisms are still unrecognized.
So, in understanding the lightning initiation, tropospheric water vapor content and its dynamics over time could play a key role. Under this hypotheses, in this work, the search for the boundary conditions around the lightning that triggers TGF will be conducted by analizing the water vapour distribution in time.
In particular, the correlations between lightning activity and tropospheric Water Vapor (WV) content were analyzed by [25], showing that the maximum/minimum Extremely Low Frequency (ELF: 1 Hz<f<100 Hz) signal, often precede the maximum/minimum of water vapor measurements on a daily basis.
However, monitoring WV dynamics over time is very difficult. The possibility to measure precipitable WV using the global positioning system (GPS) was first explored by [26,27]. [28], combining integrated precipitable water vapour data from a GPS receiver with other meteorological data found an important trend anomaly up to about 12.5 h before the first lightning strike, able to predict lightning.
In order to identify recurrent patterns useful for improving nowcasting applications, [29] combined estimates of the WV content from GPS signal with visible /infrared measurements from the Meteosat Second Generation (MSG) and with the lightning activity collected by the ground-based lightning detection network (LINET). They found specific trends appearing before the peak of lightning activity on a timescale from 2 to 3 h.
The paper is organized as follows: in Section 2, data and instrumentations are described. In Section 3 the results from ERA5 are shown while GPS ones in Section 4. The last section reports discussion and conclusions.

AGILE MCAL
AGILE is a satellite owned and operated by the Italian Space Agency and dedicated to gamma-ray astrophysics. It was launched April 23, 2007 into a low Earth orbit (~ 550 km altitude) with an inclination of 2.5° [30]. The purpose of the mission is to provide a tool with imaging capabilities in gamma-rays with a large field of view, in order to provide better studies on galactic and extragalactic sources. Among the three instruments onboard, the Mini-Calorimeter (MCAL) is the one specific for the detection of gamma-ray transients, covering an energy range from 300 keV to 100 MeV with an absolute time accuracy of ~ 2µs [31]. Specifically, when a significant number of counts is detected in a specific time window, MCAL is triggered and data is downloaded to telemetry. Furthermore, the quasi-equatorial orbit is optimal to observe equatorial regions, where most TGFs events take place. In addition to transient sources of cosmic origin, mostly gamma-ray bursts (GRB), MCAL resulted also an optimal detector for TGF.
Additional informations on MCAL preformances are included in [10,31,32]. Additional informations on the association between TGFs data and lightning data are included in [9].
We here analyzed a total of 648 TGFs with an associated lightning sferics from March 2015 to February 2020, detected from the AGILE MCAL instrument, representing an extension of the 3rd AGILE catalog (282 TGFs in the period 2015/2018 [9]. In particular, the TGF sample is identified by the association criteria with radio sferics, detected by the WWLLN network [33], are described in [9]. Figure 1 shows the geographical distribution of the sample. Figure 2 the global distribution of longitude (a), local time (b), duration (t50) (c) and number of counts (d) for the data set used (for informations on duration and intensity procedure calculation see [9]. The global distribution of TGFs shows preferential coastal areas. The empty area covering the east side of South America represents the South Atlantic Anomaly (SAA), where the detection is not active.
It is important to underline that the peak over Africa, where the electrical activity manifests its distribution peaks on Earth, results underestimate due to the less coverage by WWLLN sensors that make difficult the recording of lightning activity.
The local time distribution shows higher rate of events occurring in the early morning as well as in the late afternoon according to [10].
The duration of TGFs, expressed in t50 parameter, shows a peak at 20-40 µs, consistent with the observations by [34] regarding TGFs with lightning sferics simultaneous association, confirming the higher WWLLN matches chance with brief TGFs durations [9,34].
Concerning the spatial accuracy on WWLLN data we here assumed an uncertainty of 15 km [33].

ERA5 Reanalyses
ERA5 reanalysis is the fifth generation of ECMWF (European Centre for Medium-Range Weather Forecasts) atmospheric reanalyses and covers the period from 1950 until five days before the real time. Reanalyses combines, optimally, model data and observations to give a complete and consistent representation of the atmosphere. This principle is called data assimilation.
Reanalysis does not have the constraint of issuing timely forecasts, so there is more time to collect observations compared to the operational analyses. In addition, when going back in time, analyses allow for the ingestion of improved versions of the original observations, giving a benefit to the quality of the reanalysis product.
The assimilation system is able to estimate biases between observations and give more weight to good-quality data compared to poor data. The laws of physics allow for estimates at locations where data coverage is low, propagating in space and time the impact of observations. ERA5 provides hourly estimates of a large number of atmospheric, land and oceanic climate variables. The data cover the Earth on a 80 km grid and resolve the atmosphere using 137 levels from the surface up to a height of 8.0 km. Data are available for surface and upper model levels, and can be interpolated on pressure, isentropic and constant potential vorticity levels. In this paper, we considered data relevant for describing atmospheric convection.

GPS
The acronym GNSS (Global Navigation Satellite Systems) defines all the constellations of artificial terrestrial satellites for positioning and navigation. In this regard, only the GPS constellation was employed in this specific study. Today, this system is used in many fields, ranging from navigation [35,36], to monitoring applications [37][38][39] and meteorology [29,38,[40][41][42].
For what concerns this latter field, it is useful to remember that in the Earth's atmosphere the signal propagation speed changes due to the physical state of the medium which is crossed. Therefore the analysis of tropospheric delay, basically caused by the presence of gas and water vapor, gets particularly interesting.
Dry air and water vapor molecules in the troposphere affect GNSS signals by lowering their propagation velocities with respect to vacuum [26,43]. This tropospheric delay can be modelled during GNSS data processing or used as source of information; in the second case, the parameter known as zenith total delay (ZTD) is estimated. The ZTD is the delay related to the zenith direction, obtained after introducing a mapping function, which depends on physical parameters, able to project into zenith direction the signal delay along each single signal path.
In order to better understand if the distribution of single signal paths, called Slant Total Delay (STD), is well representative of the area of interest, the configuration of Pierce Points (PPs) referred to the involved GPS receivers is given. The Pierce Points (PPs) are the intersections between GPS lines of sight and an ideal shell located at fix altitude; in this specific case the shell hight is set at clouds altitude.
Contributions of dry air, Zenith Hydrostatic Delay (ZHD), and water vapor, Zenith Wet Delay (ZWD), to the Zenith total delay can be separated and estimated [44], being valid the relation: where ZHD, caused by dry gases present at the troposphere, is easy to model, whereas ZWD, caused by water vapour and condensed water in form of clouds, is highly variable. Starting from meteorological data of pressure and temperature, after appropriate adjustments related to the difference in altitude between meteorological sensors and GPS receivers [45][46][47], it is possible to retrieve, by applying models based on hypotesis of standard atmosphere, e.g. [26,43], ZHD estimates and, consequently, ZWD and Precipitable Water Vapor (PWV) values, by performing the appropriate conversion.
In this paper, meteorological data of pressure and temperature are retrieved by default from the empirical GPT (Global Pressure and Temperature) model [48] and retrieved PWV data have been compared with corresponding products provided by ERA5 [49].
In this study only geodetic receivers were employed, so starting from the dual-frequency observational files (RINEX Version) collected by the devices at 30 seconds rate, the PPP technique [50], undifferenced phase observation processing, was applied using ionosphere-free combination in order to estimate ZTD values for each epoch, by daily processing sessions.
The involved ancillary products (ephemeris and clocks) were the precise products provided by International GNSS Service (IGS).
Processing were handled by the beta release of goGPS software, version 1.0, written on the basis of older releases [51] by using a new batch least-squares engine.
This section may be divided by subheadings. It should provide a concise and precise description of the experimental results, their interpretation, as well as the experimental conclusions that can be drawn.

GOES
The Geostationary Operational Environmental Satellite (GOES) program started in 1975 when the first satellite was launched under the coordination of the National Aeronautics and Space Administration (NASA) and the National Oceanic and Atmospheric Administration (NOAA). Currently, four GOES series satellites are operating: GOES-13 and GOES-14 and two GOES-R series satellites (GOES-16 and GOES-17). In this work, the GOES-16 data have been used. The GOES-16 is located at 75°W, 0° N and provides data at 10 minutes time resolution. The sub-satellite point has a spatial resolution around 1 km, which becomes coarser moving away from it. The GOES-16 was the first satellite of GOES-R series to mount the Advanced Baseline Imager (ABI) [52,53]. ABI has 16 spectral bands ranging from visible (VIS) to infrared (IR) electromagnetic spectrum in order to characterize the properties of both cloud and atmospheric gases. In this work, we limited to use only one out of the sixteen ABI channels, namely channel 13. The channel 13, being a "window" channel centered at 10.33 µm, allows to retrieve the cloud top temperature by measuring the brightness temperature (TB).
Currently, Himawari 8/9 satellites are available for operational use. In this work, the Himawari 8 data have been used. The satellite entered operational service on July 7 2015 at 140.7° E, 0° N carrying on the Advanced Himawari Imager (AHI), a visible infrared radiometer that possesses comparable capabilities to the GOES-R in terms of spatial, spectral and temporal resolutions [55]. The instrument has 16 observational bands, spanning from the visible to thermal infrared bands useful for retrieving cloud properties. In this work we used the same GOES-R channel, namely the channel centered at 10.4 µm for the retrieval of the cloud-top temperature [55].

ERA-5 results
In this section we use ERA5 reanalyses [56] to evaluate the meteorological parameters when TGFs are observed and to compare these values with their reference. The meteorological parameters considered are relevant for convection and are the Convective Available Potential Energy (CAPE), the Convection Inhibition Energy (CIN), the Total Column Water Vapour (TCWV) and the 2m dew point temperature (T2D).
To compute the values of the meteorological parameters in correspondence of the TGFs observations, we interpolated the ERA5 field to the position and time of the TGF occurrence. The spatial interpolation is bilinear starting from ERA5 fields at 0.25° horizontal resolution, while the temporal interpolation is linear, from hourly reanalyses.
Reference values for each parameter are computed considering for each TGF all the hourly values of the meteorological parameter for the month and for the position where TGF is observed. Then, all the data corresponding to TGF and reference are gathered together and represented by box plots. Figure 3 shows the results of this analysis for the CAPE (Fig. 3a), CIN (Fig. 3b), TCWV (Fig. 3c) and for the T2D (Fig. 3d).
The average value of CAPE when TGF are recorded are larger than the corresponding reference values for almost all months considered. This holds for the 25th and 75th percentile too. This result shows, as expected, that TGF occurs in convective environments. However, the values of CAPE associated with TGF are largely variable as shown on Figure 3a and are not exceptional compared to the reference values of the parameter for the locations and months where/when TGFs are observed. Indeed, the reference maximum and minimum values always include the corresponding values when TGF are observed.
CIN (Figure 3b) is larger when TGF are observed compared to reference values, even if the maximum and minimum values of the reference includes the maximum and minimum values obtained when TGF is observed. Convection is associated with small but positive value of CIN because this avoids a too fast consumption of the available potential energy that would result in shallow or no convection. TCWV (Figure 3c) is larger when TGF are recorded compared to the reference values. This occurs when convection is developing, as also shown by the GPS-ZTD analysis discussed in Section 4, but the values of TCWV associated with TGF are included in the intervals of reference values, so they are not extreme values.
Similar considerations apply for the dew point temperature at the surface. The higher values for TGF events show larger amount of water vapour at the surface compared to reference values. This is well explained by the convective environment in which TGF occur, which has higher than normal humidity at the surface.
Overall, the analysis of ERA5 data shows that, while TGF occur in convective environments, however the values of meteorological parameters describing the convective environment are not exceptional.

Case studies
In order to improve the understanding of the phenomena illustrated, a total of nine case studies were selected according to a distance-based criteria; the distances between GPS receivers location and TGF occurrence was set within 45 km. These case studies are pinpointed by orange diamonds (representing the TGFs) and black dots (representing the GPS receivers) in Figure 4 and Figure 5 (since some TGFs occurred very close to each other, even if in different dates, they overlap in the figures but are clearly reported in Table  1 and Table 2). Two out of the nine case studies (green boxes in Figure 4 and Figure 5) are more in depth discussed in sub-section 4.1 and 4.2, respectively.
In particular, the daily trend of PWV values estimated by all the available GPS receivers (within 45 km from the TGF location) and by the ERA5 reanalysis were compared.  As mentioned before, with the aim to characterize tropospheric conditions at TGF occurrence as well as the reliability limits of the PWV derivation techiques (GPS versus ERA5 based techniques) the GPS-PWV was matched with the strokes registered by WWLLN within an 1° x 1° area centered in the TGF location. In addition, the information obtained from geostationary satellites (i.e. GOES or Himawari depending on the location of the considered event) allows the time monitoring of the meteorological conditions preceding the TGF occurrence. It has to be taken into account that each case study was characterized by a distinctive geographical morphology with respect to the others.
In Table 1 and Table 2, each case study is identified through temporal (Date-Time in the first column) and spatial (φ (lat) and λ (lon) as well as altitude H in the fourth column) coordinates. Referring to each TGF occurence, nearby GPS receivers were selected and included in Tables 1 and 2 by their marker names (third column) and their coordinates (fourth column). Two additional parameters, the distance between TGF and GPS position and the correlation coefficient, computed on a daily scale, between PWV-GPS and PWV-ERA5, are placed in the last two columns. This last parameter (r) was used in order to compare the operating conditions of the two techniques.
As can be seen from the r-values, the correlation is closely linked both to the distance between TGF occurrence and GPS receiver and to the morphology of the territory where the phenomenon is located. In fact, the correlation tends to degrade as the orography becomes more complex (high altitudes) and as the distance between TGF and GPS device increases. This implies a good reliability of ERA5 in cases where the orography is simple (low altitude), although it remains a technique whose values are more smoothed, mainly because of spatial resolution (wide grid mesh), which inhibits its the ability to identify small scale variability, which is well capture by GPS. In tables the case studies described below are marked in bold.

Sumatra -March 16th, 2019
This first case study concerns an event occurred near the coast of Sumatra on March 16 th , 2019. The TGF occurred at 10:35:38 UTC at the following latitude and longitude: φ = -2.44, λ = 101.41. Two GPS receivers, LNNG and MKMK, were placed at 33 and 37 km distance from the TGF location, and at an altitude of 42m and 6m above the sea level (a.s.l.), respectively. Figure 6 shows the growth of the convective cell that will generate the TGF about two hours later. The satellite images are very useful in order to understand if the PWV variation measured by the GPS receivers is due exclusively to the precipitating system linked to the TGF. The sequence, starting 2 hours before the TGF occurrence (Figure 6a), evidences the growth of the convective cell will generate the TGF. The TBs reported in Figure 6e) (i.e. the Himawari-8 image closest to the TGF time) highlights very deep convection with values up to 205 K. The convective cell reaches its maximum vertical extension in the following 30 minutes with respect to the TGF (Figure 6e) then to move to its expiration. Furthermore, the TGF is located at the edge of the convective cell. This feature is common to almost the totality of the cases analyzed. Figure 7 and Figure 8 show a direct comparison, after the bias removal, between PWV-GPS data and PWV-ERA5 data. In particular, panel a) of each of the two figures reports the daily trend of GPS-PWV (black solid line) and ERA5-PWV at both GPS and TGF coordinates (gray and orange solid lines, respectively). The time of the TGF occurrence is marked by orange dashed vertical line. On the other hand, the panel b) shows the scatterplot between GPS-PWV and ERA5-PWV. The scatterplots show a high correlation between GPS-PWV and ERA5-PWV for both GPS receivers with rLNNG=0.87 and rMKMK=0.85, respectively (see Table 2); this means that both in this specific and similar circumstance (not particularly complex territorial morphology) it is possible to use the ERA5-PWV dataset in order to evaluate the PWV behaviour at TGF coordinates.
Comparing the GPS-PWV trend with ERA5 trend at TGF coordinates, it can be identified a slight offset between curves (Figure 7) related to the structure of the convective clouds and to the position of GPS sensors in relation to TGF coordinates.
In Figure 9 the distribution of Pierce Points (PPs), that is the intersections between GPS lines of sight and an ideal shell located at clouds altitude, referred to LNNG GPS receiver is given. The aim is to analyze the configuration of PPs with respect to GPS receiver. The analysis shows that the PPs are well distributed around the GPS receiver and close to the event; therefore the GPS-PWV trend (Figure 7) can be considered well representative of the water vapor content within the area of TGF occurrence. Finally, the GPS-PWV trend is compared with the stroke rate trend registered by the WWLLN measurements within an 1° x 1° area centered in the TGF location ( Figure 10). The TGF occurred after that a marked increase of PWV reached a local maximum at 09:30 UTC about. On the other hand, the TGF preceded the maximum stroke rate occurred at 11:00 UTC about with values reaching 16 strokes/min. This second case study consider an event occurred over the mountains of Ecuador on November 15 th , 2019.
Two TGFs, very close in time and space to each other (slight more than one minute difference and around 19 km apart), were recorded for this event. In this case, four GPS receivers are available at different distances, MZEC, BIEC, RIOP and VZCY (see Table 2). Figure 11a is the same of Figure 6 except that the data were collected by GOES-R satellite. The sequence of satellite images, starting 2 hours before the TGF occurrence (Figure 6a), shows the same features of the previous case highlighting the development of the convective cell will generate the TGF. The TBs reported in Figure 11e) (i.e. the GOES-R image closest to the TGF time) shows very deep convection with values up to 205 K even if the spatial extension of the convective cloud is quite limited, especially if compared to the convective cells covering the area. Also in this case, the TGFs are both located at the edge of the convective cell (one on the southern edge and one in northern-west edge).
For this case study, the comparison between GPS-PWV data and ERA5-PWV data showed a bad correlation with negative values of r (see Table 2). The complex topography affects the reliability of ERA5 model and the patterns, obtained by the bias removal, showed great discrepancies with respect to the GPS-PWV (e.g. Figure 12, central and right panels considering the MZEC receiver). This point highlights the strong potential of the use of GPS data, for the troposphere characterization, in areas with complex territorial morphology. Furthermore, the PWV trend for the four GPS receiver is very similar to each other (except more marked differences at the end of the day) even if the absolute values change depending on the altitude of the GPS receiver.
In Figure the distribution of Pierce Points (PPs) referred to MZEC GPS receiver is shown.
Also in this case, the PPs are well distributed nearby both the GPS receivers and TGFs. This implies that the GPS-PWV trend ( Figure 12, left panel) results particularly useful and reliable for the analysis of the event. It has to be highlighted that because of the vicinity in space and time of the two TGFs, the analysis was performed taking as refence only one out of the two TGFs.  The combined analysis of the GPS-PWV for MZEC site and the stroke rate trend ( Figure  14) shows very similar features to the previous case study. In particular, the TGFs occurred slightly before both the second stroke rate maximum and the PWV maximum. For this event, the stroke rate has higher values exceeding 20 strokes/min.

Discussion and conclusions
In order to characterize the environmental and the cloud properties associated to the generation of a TGF, a total of 648 events detected by AGILE MCAL was analyzed. The investigation of the atmospheric conditions at TGF occurrence with ERA5 reanalyses shows that TGFs mainly occur in convective environment. But the values of meteorological parameters describing the convective environment are not exceptional. The CAPE measured in correspondence of a TGF occurrence does not show any particular variation with respect to the monthly average in the same location. This results is in agreement with the results obtained from the analysis over 24 TGF-producing storms conducted by [15]. On the other hand, the high values of CIN, TCWV and T2D measured at the TGF occurrence are indicative of the fact that they occurs preferentially when the updrafts reach their maximum development. The time evolution of these convective systems was made by the use of geostationary capabilities, useful to describe the development of the storms and able to establish if the PWV trend measured by the GPS sensor is due solely from the TGF associated storm.
The life cycle of the storm (from the early stages to its mature phase) through a sequence of nine snapshots, shows how, at the TGF time, cloud reaches deep convection with cloudtop temperatures dropping down to about 200 K for the totality of the observed events (see additional material). These result are in agreement with examinations over TGF-producing storms conducted by [17], where both active and passive microwave views of the considered thunderstorms agree in describing TGFs producing clouds as intense thunderstorms, with the presence of convective towers. Within the analysed structures, however, the cloud top topography shows a wide range of extensions and the TGF is located at the edge of the convective cell. This feature is common to the almost the totality of the cases analyzed.
A more in depth characterization of the atmospheric conditions at TGF occurrence was made possible by the analysis of PWV as estimated by the GPS sensors and the ERA5 dataset. To our knowledge this is the first time that this technique is applied to the characterization of TGF-producing storms. The comparison between GPS-PWV and ERA5-PWV data, for a number of case studies characterized by different topographical conditions, showed a very good correlation in case of low altitudes and bad correlation where orography is impactful.
The high correlation at low altitudes suggests how the lower resolution (in space and time) ERA5 data can be used to characterize the environmental conditions during a TGF occurrence. However, the analysis shows the strong potential of the use of GPS data for the troposphere characterization in areas with complex territorial morphology. Furthermore, the focus on the distribution of PPs with respect to GPS receivers shows that the PPs are well distributed around the GPS devices and close to the TGF location (we recall that the distance between the TGF and the GPS sensors does not exceed 45 km). Therefore the GPS-PWV trend can be considered well representative of the water vapor content within the areas of TGF occurrences.
In all the analysed case studies, the time of the TGF is within the maximum convection phase. Moreover, comparing the GPS-PWV trend with ERA5-PWV trend at TGF coordinates, it is detectable a slight offset between GPS-PWV and ERA5-PWV curves related to the structure of the convective clouds and to the position of GPS sensors in relation to TGF coordinates.
The combined use of GPS and the stroke rate trends, shows for both cases an increase of PWV about two hours before the TGF. This trend is generally confirmed also considering all the case studies shown in the additional material. Moreover, in both cases, TGFs occurred before or during the lightning peak and in time intervals where PWV reaches high values compared to the daily average.
These results suggest that the TGF occurrence and the flash rate trend is consistent with what shown in [13,14], highlighting that a TGF often occurs during the most lightning active phase of the storm. Moreover, the relation between PWV and lightning trend is in agreement with [29] where the PWV rising starts from 2 to 3 hours before the peak of lightning activity. The lightning flash rate distribution exhibiting high values (16 and 20 flashes/min), compared with [14], where the 50% events showed a flash rate less than 5 flashes/5 min with a maximum reference values of 40 flashes/5 min.
Moreover, we note that the sample presented in this work is global, including regions with a very variable lightning detection efficiency given by WWLLN, and where low, more prone to a bias towards large-scale convective systems with extended lightning activity, more easily detectable. We point out that this is an unavoidable bias affecting every global TGF sample requiring lightning association, which is in turn a mandatory requirement if a location of the source with accuracy below tens of km is needed. We note that this bias affects also other studies on the meteorological characteristics of TGF-producing thunderstorms, such as [11,12,14].
Finally, the ERA5, GPS, geostationary and stroke analysis performed in this study provided indications on the dynamics of convective systems linked to TGFs.