Open access

Implementation of Ionospheric Asymmetry Index in TRANSMIT Prototype

Written By

M.M. Shaikh, R. Notarpietro and Bruno Nava

Published: July 17th, 2014

DOI: 10.5772/58551

Chapter metrics overview

1,947 Chapter Downloads

View Full Metrics

1. Introduction

Radio occultation (RO) [4] missions such as GPS/MET, CHAMP (pilot projects) [3, 17, 18], COSMIC and MetOP [1, 11, 19, 20] have been designed to sound the Earth’s neutral atmosphere and ionosphere via radio links between a GPS and GPS receiver on-board Low Earth Orbit (LEO) satellites. The U.S. GPS/MET experiment was the first mission, which successfully applied RO technique to the Earth atmosphere monitoring using GPS signals. Since then, RO technique has become a powerful tool to study the ionosphere [7, 10]. In neutral atmosphere, using RO technique, the bending of the signal is extracted and inverted into refractivity profiles through the Abel inversion [5, 8, 17]. In the ionosphere, where bending is negligible, carrier phase measurement and corresponding limb-TEC (LTEC in what follows) observations are used to extract electron density profiles, Ne(h), defined along the tangent points of ray paths (also known as ‘ray perigees’) between LEO and GPS (see fig. 1). The retrieval algorithm considers the time-series of LTEC below LEO orbit (between points A1 and A2 in fig. 1) observed from the same GPS satellite during an occultation event. Therefore the LTEC above LEO orbit (TEC from GPS-to-B or GPS-to-A2 in fig. 1) have to be removed before starting the inversion procedure [14] in order to take the contribution of LTEC from upper atmosphere out from the total LTEC (TEC from GPS to LEO). With the data from both pilot missions (GPS/MET, CHAMP), it was only possible to perform LTEC measurements at the highest point of LEO orbit (point B in fig. 1) so that the LTEC above LEO had to be modelled or considered constant all over the occultation [14]. With modern interpolation techniques applied in post processing of current RO missions (such as COSMIC), such above LEO orbit LTEC can be precisely estimated and removed before applying the data inversion techniques. A widely used data inversion technique used to obtain vertical electron density profile in the ionosphere is the ‘Onion-peeling’ inversion algorithm [9].

Figure 1.

Limb TEC (LTEC) measurement in the ionosphere using radio occultation technique. A1 and A2 are points defined at opposite sides of LEO orbits around ray perigees (black dots). TEC calculated between A1 and A2 is defined as ‘internal orbit LTEC’.

Onion-peeling algorithm is based on the assumption of spherical symmetry of Ne distribution in the ionosphere (Ne depends only on height). It is a very effective tool for RO data inversion in case of small horizontal gradients present in the ionosphere (particularly during undisturbed geomagnetic periods). But, for disturbed geomagnetic conditions, for example under the equatorial anomaly region, large electron density gradients may be experienced which could lead to the failure of Onion-peeling algorithm producing erroneous Ne(h) profiles as output. In the present work, a simulation study has been performed to assess the effects of the spherical symmetry assumption on the inverted electron density profiles using Onion-peeling. In order to produce synthetic background ionosphere data, we used two climatological models, NeQuick [13] and IRI [2], and a data analysis system known as ‘Multi-Instrument Data Analysis System (MIDAS)’ [12]. MIDAS use tomographic techniques to combine observation from many sources simultaneously in a single inversion, with the minimum of a priori assumptions about the form of the ionospheric electron concentration distribution. In this work, MIDAS has been used in two modes, i.e. with standard form and with RO data assimilation. In its standard form, MIDAS doesn’t use RO data as a source of ionospheric data. For the latter mode, ionospheric RO data from COSMIC mission for 26th September 2011 has been used for the assimilation. We only performed our analysis for the geomagnetic storm observed on 26th September 2011 over mid-latitudes.

In the next sections, we have presented how we assessed the problem of asymmetry in the ionosphere and its implementation in TRANSMIT prototype. In section 2, we have briefly introduced the Onion-peeling algorithm and its implementation in our work. In section 3, a discussion on the simulation results are presented. In section 4, an overview of TRANSMIT prototype with the description of processor 3C (Ionospheric asymmetry) is presented. In section 5, we have summarized the work by drawing the main conclusions and discussing prospects for the future work.


2. Formulation of the problem

In this work, we have applied the standard Onion-peeling algorithm to invert simulated RO data under a constraint of using ideal geometries: by considering only internal orbit ray paths (ray paths below the LEO orbit), fixed occultation planes and vertically distributed ray perigees positions. This is considered in contrast to what is shown in fig. 1 which shows an illustration of real RO event in which ray perigee positions are not vertically distributed. This happens because LEO and GPS satellites circulate in totally different orbits and independent of each other. LEO satellites, in most cases, have orbit altitudes well below 1000 km and GPS satellites have orbit altitudes approximately 20, 000 km from Earth’s surface. In a real RO event, this creates a scenario where the azimuth of the occultation plane changes with almost each ray exchanged between LEO and GPS (there are several hundred rays exchanged in a single RO event). Consequently, position of ray perigee (the tangent point of ray exchanged between LEO and GPS) is independently computed for each ray. Therefore, it is not possible to have vertically distributed ray perigee positions in a real RO event. However, in this work, rationale behind the use of this so called ‘ideal geometry’ is to focus the dependency of retrieval errors on the inversion approach only and avoiding to take into account inaccuracies due to geometry. The possibility to use external data, like vertical TEC maps [6] to improve the solution has therefore not been considered in the present work. In our analysis, we considered a LEO satellite with orbit altitude of 800 km from Earth’s surface. The background ionosphere is computed using NeQuick, IRI and MIDAS (with and without RO data assimilation).

2.1. RO data inversion using onion-peeling technique

For a given occultation event, the LTEC related to the internal orbit ray path ‘i’ can be computed considering a set of spherical shells (identified by peel ‘j’) like ‘onion shells’ [9], characterized by a constant electron density (the radius of each shell is the impact parameter of the ray). Analytically, the LTEC associated to the ith ray can be defined as:

LTECi=2j=1, Ni2LijNejE1


‘Lij’ is the length of the segment ‘i’ related to the electron density characterizing shell ‘j’

‘Nej is the electron density charactering shell ‘j’

‘LTECi is the limb-TEC value related to the ‘ith’ ray path crossing all the shells

‘N’ is total number of shells

Considering this definition, the LTEC may be easily inverted to extract the Ne characterizing each shell ‘Nej’, starting from the most external ray path. The matrix (Eq. (1)) describing the linear system of equations is triangular, and therefore it can be solved from top to bottom for LTEC inversion to extract the Ne(h) profile.

2.2. Onion-peeling derived errors

Following are the definitions of the errors and observables we have used to evaluate the impact of ionospheric asymmetry on Onion-peeling inversion. The retrieved electron density profiles are then compared with the collocated NeQuick, IRI and MIDAS (with and without RO data assimilation) ‘true’ vertical Ne(h) profiles in terms of difference on VTEC and difference on the NmF2 values. These indicators can be defined as:

ΔVTEC=h*hLEONeOnion-peeling dh- h*hLEONeNeQuickdhE2
ΔNmF2 =NmF2Onion-peeling- NmF2NeQuickE3

In [15], we defined the asymmetry level index of the ionosphere considering the degree of dissimilarity of electron density distribution along the two-halves of a given below LEO orbit ray path (as defined in section I) crossing the ionosphere from LEO satellite height down to a 100 km. In the present work, we have applied the asymmetry level index by using three different electron density distribution obtained from NeQuick, IRI and MIDAS (with and without RO data assimilation) and presented a comparative analysis of the results together with its implementation in the TRANSMIT prototype.


3. Simulation results

In this work, we have only performed our analysis over mid-latitudes where a geomagnetic storm effects were observed from UT 18:00 on 26th September 2011 to early hours of 27 September 2011. NeQuick and IRI were used with the solar flux input taken automatically by the model. For MIDAS, results are presented using both standard output and RO data assimilated output. Ionospheric RO data from COSMIC mission collected during 26 September 2011 storm (from UT 12:00 to UT 24:00) over mid-latitudes (Latitude range: 20°N-84°N/ Longitude range: 116°W – 64°E) has been used for the assimilation. Approximately, data from 80 RO events of COSMIC mission have been assimilated. We divided our analysis in two subsets, Quiet-time (from UT 13:00 to UT 14:00) before the start of geomagnetic storm and Storm-time (from UT 19:00 to UT 20:00) on 26 September, 2011.

3.1. Quiet-time analysis

Figure 2.

Quiet-Time results for UT 19:00 and 20:00. (a) Asymmetry comparison (b) ΔVTEC comparison (c) ΔNmF2 comparison

Fig. 2(a) shows two comparative plots of asymmetry values evaluated using NeQuick, IRI and MIDAS (with and without RO data assimilation) for UT 13:00 and UT 14:00. Fig. 2(b) and 2(c) show the derived Onion-peeling inversion errors (ΔVTEC and ΔNmF2) for the same UTs, respectively. Comparing the plots, it is evident that an overall good correlation exist between evaluated asymmetry and the associated inversion errors (ΔVTEC and ΔNmF2) as function of the azimuth of occultation plane for NeQuick, IRI, and MIDAS (with and without RO data assimilation). Considering the quiet geomagnetic conditions, asymmetry evaluation is low for all three background ionosphere as expected (maximum asymmetry index values are less than 0.4). These low values can also be observed from the plots of inversion errors. Although the derived error values shown in 2(b) and 2(c) are not exactly correlated with the evaluated asymmetry for all azimuth of occultation plane, we still can find good agreement between the overall behavior of the asymmetry with the inversion errors. In fact, there is a very good agreement between them for certain range of the azimuth of occultation plane in parts. Moreover, asymmetry values together with the inversion error values are all low in numbers which is a good indication that the evaluation has been carried out during a quiet geomagnetic period. Use of MIDAS with or without RO data assimilation didn’t produce much different results in quiet-time as all values of asymmetry together with their derived inversion errors are very close to each other in both cases.

3.2. Storm-time analysis

Fig. 3 shows plots from the storm time subset. Correlation between the asymmetry values and their associated errors is evident for all electron density distributions, separately. Comparing the plots at the same UTs, it is evident that an even better correlation exists between the evaluated asymmetry and the Onion-peeling derived errors (ΔVTEC and ΔNmF2) as a function of azimuth of the occultation plane. Unlike quiet-time subset, MIDAS (with and without RO data assimilation) results show very different behavior from NeQuick and IRI. This shows a clear difference between a climatological model (NeQuick and IRI) and a data assimilated model (MIDAS). There seems to be no evidence of an active storm from NeQuick and IRI results (asymmetry range is almost similar to quiet time subset). Whereas, in MIDAS case we can easily observe large values of asymmetry and associated errors from Az=120° to 150° which provides evidence of an active storm. A higher level of correlation between the results shows that the asymmetry algorithm works very well for highly geomagnetic disturbed periods. In contrast to the quiet-time analysis, the use of RO data assimilation can be better observed in storm-time. Specifically, for UT 19:00, the difference of MIDAS as a ‘ground truth’ ionosphere with and without RO assimilation is clearly evident.

Figure 3.

Storm-Time results for UT 19:00 and 20:00. (a) Asymmetry comparison (b) ΔVTEC comparison (c) ΔNmF2 comparison


4. TRANSMIT prototype

The main deliverable of TRANSMIT project is a prototype which is being developed as a coordinated input from all level 1 partners. After discussions with the industry partners, it was concluded that that Precise Point Positioning (PPP) results would be the main output of the prototype. Both Single Frequency PPP and Double Frequency PPP are being considered as prototype outputs. To achieve this, 3 main processors were identified as follows:

  1. Processor1: S4 and TEC prediction

  2. Processor 2a: PPP mitigation

  3. Processor 2b: Improved tracking,

  4. Processor 3a: Improved MIDAS TEC

  5. Processor 3b: Ionospheric models

  6. Processor 3c: Ionospheric Asymmetry

4.1. Prototype operation

Fig. 4 shows a complete flow/block diagram of the TRANSMIT prototype. Main data flow operations through the prototype has been marked in red as 1 to 6. A brief description of each operation is given below:

Figure 4.

Prototype data flow diagram

  1. User request (a form with selection of subprocessor and case study should be available)

  2. The user parameter file is generated and sent to POLITO

  3. Automatically the processor 3C.subprocessor starts and asks for input (through INGV interface – see point 4)

  4. The 3C.subprocessor send output to INGV and SWACI (through INGV and SWACI interfaces)

  5. SWACI displays the output for the user (in the form of “surfable” maps as explained)

4.2. Processor 3C: Ionospheric asymmetry

Processor 3C will deal with the ionospheric asymmetry evaluation. A complete block diagram of the processor 3C with data inputs and outputs for each block is shown in fig. 5. For each identified case study, processor 3C will exploit high correlation between asymmetry and associated errors on retrieved electron density profiles. This will be done for real RO events’ data taken from COSMIC mission. The processor will produce global asymmetry maps (not shown to the users) using different background ionosphere (provided by model data computed using NeQuick, IRI and MIDAS) for quasi-horizontal TEC observations. By querying these global maps of asymmetry, for each identified case study, processor 3C will compute the expected level of asymmetry present in the ionosphere in the geographical location of the selected RO event. Then, on the basis of the computed asymmetry, the RO event will be displayed in specified color (red, yellow or green) on a 2D map (as shown for output 1 in fig. 4). The color will be an indicator of the expected quality of RO product (considering standard Onion-peeling inversion), as shown in fig. 5. The second output (output 2 in fig. 4) will be a comparison between two RO data inversion techniques; one is the ‘standard Onion-peeling’ and the other is ‘model-aided inversion algorithm’ (we are currently working on). For model-aided data inversion, we will provide electron density profiles by taken real geometry from COSMIC mission and considering different model data as background ionosphere. Based on functionality, we have divided processor 3C in three sub-processors 3C.1, 3C.2 and 3C.3 (as shown in fig. 5). A brief description of all sub-processors are as follows:

Figure 5.

Processor 3C block diagram (see Appendix for larger view)

4.2.1. Sub-processor 3C.1: Global maps of asymmetry index

This part of processor 3C will generate global asymmetry maps for the selected case study. Ideal radio occultation (RO) geometries are taken into account. 18 global maps for each background ionosphere will be generated. Main features of sub-processor 3C.1 are as follows:

  1. Global asymmetry maps will be generated considering ideal RO geometries and three different background ionosphere computed using NeQuick, IRI and MIDAS; for each 10° azimuth of occultation plane. The asymmetry level will only be computed for one trajectory of the RO event (the one at 100 km). Users will not be able to see any of these asymmetry maps. These maps will only be used in the processing of output1 and output2 (as shown in fig. 5).

4.2.2. Sub-processor 3C.2: Electron density profile retrieval (standard techniques) effectiveness

The main goal is to show the effectiveness of RO inversion data if the standard onion-peeling algorithm would be used. This will be done by giving a color code (green, yellow, red) to each RO event which will be based on the asymmetry evaluation for that event. The asymmetry for each RO event will be evaluated using processor 3C.1.Main features of sub-processor 3C.2 are as follows:

  1. Using real orbits of RO events available in the area defined by the case study, and the global maps of asymmetry computed by sub-processor 3C.1, processor 3C.2 will evaluate the expected level of asymmetry present in the ionosphere in the geographical location of the selected RO event. Then, on the basis of the evaluated asymmetry, the RO event will be displayed in specified color (red, yellow or green) on a 2D map (as shown for output 1 in fig. 5). The color will be an indicator of the expected quality of RO product (considering standard Onion-peeling inversion.

  2. Two electron density profile obtained considering standard retrieval algorithms will be shown. One for the best case (lowest asymmetry level among all events) and the other for the worse case (highest asymmetry level among all events).

4.2.3. Sub-processor 3C.3: Onion-peeling model aided Electron density profile retrieval

For each occultation event available in the area defined by the case study, a comparative plot will show the inverted electron density (Ne) profiles obtained using standard Onion-peeling algorithm and the model-aided inversion algorithm (in the latter case three Ne profiles will be available based on each background ionosphere evaluated using NeQuick, IRI and MIDAS). Main features of sub-processor 3C.3 are as follows:

  1. Two images of specified format containing Ne profiles of two selected events will be shown. Events with lowest and highest asymmetry will be selected (as done in sub-processor 3C.2). Comparative results related to standard and advanced retrievals will be shown as follows:

  2. One electron density profile obtained using standard Onion-peeling inversion

  3. Three electron density profiles obtained using model-aided inversion (one for each background ionosphere; NeQuick, IRI and MIDAS).


5. Conclusion & future work

In this work, we have shown the implementation of asymmetry index using three different background ionosphere computed using NeQuick, IRI and MIDAS (with and without RO data assimilation). Previously [16], while implementing the asymmetry index only with climatological models (NeQuick and IRI), it was observed that with a climatological model, it is possible to estimate asymmetry indices only for an ionosphere in ‘nominal conditions’. Indeed a climatological model cannot be used to evaluate ionospheric asymmetry for ionospheric conditions in real/near-real time as it does not support data assimilation. In this work, after completing a thorough analysis, we have observed that, in a normal solar activity condition, the climatological model may work as good as a data assimilated model in calculating ionospheric asymmetry. However, during an active solar storm, a data assimilated electron density model can outperform the climatological model and produce much improved results. We have observed a clear difference in results evaluated in two different geophysical conditions. In the storm-time case, we found asymmetry and its associated RO inversion errors much higher than in the quiet-time case. This shows a clear advantage of using data assimilated model as ‘ground truth’ in ionospheric asymmetry evaluation.

In sections 3.1 and 3.2, we presented the simulation results without highlighting how the use of MIDAS may impact the evaluation of ionospheric asymmetry. As these are our first results with only one case study using MIDAS, it would be difficult to determine exactly how much beneficial it would be to use MIDAS with RO data assimilation as it requires more computation and a higher set of input data. However, it may be concluded from the results of this study that, in case of quiet geomagnetic periods, asymmetry evaluation using MIDAS with or without RO data assimilation may not be of much difference. keeping in view of processor 3C of the TRANSMIT prototype, we are currently studying more case studies in order to analyze which mode of MIDAS would be preferable as source of background ionosphere.

In the scenario of TRANSMIT prototype implementation, we are going to evaluate global asymmetry maps (not shown to the users) with their associated inversion error (ΔVTEC and ΔNmF2) plots with background ionosphere computed using NeQuick, IRI and MIDAS. By querying these global maps of asymmetry, and for a given geometry of a real RO event, the information provided will be the prediction of the expected level of asymmetry present in the ionosphere and its potential impact on Radio Occultation inverted products using estimates of ΔNmF2 and ΔVTEC. An improved inversion technique is currently being developed which may be helpful to reduce the error in RO inversion products by tactfully removing the ionospheric symmetry hypothesis from the standard inversion techniques. This will drive the output2 of processor 3C in future.



Block Diagram of Processor 3C (Supplement of figure 5)



This research work is undertaken under the framework of the TRANSMIT ITN, funded by the Research Executive Agency within the 7th Framework Program of the European Commission, People Program, Initial Training Network, Marie Curie Actions-GA no. 264476.


  1. 1. Anthes R. A., Ector, D., Hunt D. C., Kuo Y-H., Rocken C., Schreiner W. S., Sokolovskiy S. V., Syndergaard S., Wee T-K., Zeng Z. P., Bernhardt A., Dymond K. F., Chen Y., Liu H., Manning K., Randel W. J., Trenberth K. E., Cucurull L., Healy S.B., Ho S-P., McCormick C.T., Meehan K., Thompson D.C. & Yen N.L. (2008). The COSMIC/FORMOSAT-3 Mission: Early Results. Bulletin of the American Meteorological Society, Vol. 89, pp. 313-333, doi: 10.1175/BAMS-89-3-313.
  2. 2. Bilitza D., B.W. Reinisch (2008). International Reference Ionosphere 2007: Improvements and new parameters, Advances in Space Research 42 (2008) 599–609.
  3. 3. Gorbunov, M.E. (2001). Analysis and validation of GPS/MET radio occultation data, Journal of Geophysical Research, vol. 106, No. D15, pages 17, 161-17, 169.
  4. 4. Hajj G.A., Kursinski E. R., Romans L.J., Bertiger W.I., Leroy S.S. (2001). A technical description of atmospheric sounding by GPS occultation, Journal of Atmospheric and Solar-Terrestrial Physics.
  5. 5. Healy S.B., Haase J., Lesne O. (2002). Abel transform inversion of radio occultation measurements made with a receiver inside the Earth’s atmosphere, Annales Geophysicae (2002) 20: 1253–1256 c European Geophysical Society.
  6. 6. Hernandez-Pajares M., Juan J. M., Sanz J. (2000). Improving the Abel inversion by adding ground GPS data to LEO radio occultation in ionosphere sounding, Geophysical research letters, vol. 27, no. 16, pages 2473-2476.
  7. 7. Kulikov I., A. J. Mannucci, X. Pi, C. Raymond, G. A. Hajj (2011). Electron density retrieval from occulting GNSS signals using a gradient-aided inversion technique, Advances in Space Research, 47(2), 289-295, doi:10.1016/j.asr.2010.07.002.
  8. 8. Kursinsky E.R. (1997). Observing Earth’s atmosphere with Radio Occultation measurements using the Global Positioning System, Hournal of Geophysical Research, Vol. 102, No. D19, Pages 23, 429-23, 465.
  9. 9. Leitinger R., H.P.Ladreiter, G.Kirchengast (1997). Ionosphere tomography with data from satellite reception of GNSS signals and ground reception of NNSS signals, Radio Sci., 32(4), 1657-1669.
  10. 10. Liu J. Y., C. Y. Lin, C. H. Lin, H. F. Tsai, S. C. Solomon, Y. Y. Sun, I. T. Lee, W. S. Schreiner, and Y. H. Kuo (2010). Artificial plasma cave in the low-latitude ionosphere results from the radio occultation inversion of the FORMOSAT-3/COSMIC, J. Geophys. Res., 115(A7), A07319, doi:10.1029/2009JA015079.
  11. 11. Luntama J-P., Kirchengast G., Borsche M., Foelsche U., Steiner A., Healy S., Von Engeln A, O'Clerigh E. & Marquardt C. (2008). Prospects of the EPS GRAS Mission For Operational Atmospheric Applications. Bulletin of the American Meteorological Society. Vol. 89, pp. 1863-1875, doi: 10.1175/2008BAMS2399.1.
  12. 12. Mitchell, C.N., Spencer, P. (2003). A three-dimensional time-dependent algorithm for ionospheric imaging using GPS, Annals of Geophysics, vol. 46, no. 4, pp. 678–696.
  13. 13. Nava B., Coisson P., Radicella S.M. (2008). A new version of the NeQuick ionosphere electron density model, Journal of Atmospheric and Solar-Terrestrial Physics 70 (2008) 1856-1862.
  14. 14. Schreiner W.S., Sokolovskiy S.V., Rocken C. (1999). Analysis and validation of GPS/MET radio occultation data in the ionospher, Radio Science, Volume 34, Number 4, Pages 949-966.
  15. 15. Shaikh, M.M., Notarpietro, R., Nava, B. (2013). The Impact of Spherical Symmetry Assumption on Radio Occultation Data Inversion in the Ionosphere: An Assessment Study, Advances in Space Research, doi:
  16. 16. Shaikh M. M., Notarpietro R. (2013). Evaluation of the Impact of Ionospheric Asymmetry on GNSS Radio Occultation Inversion Products using NeQuick and IRI, Proceedings of International Conference on Localization and GNSS (ICL-GNSS), ISBN: 978-1-4799-0484-6.
  17. 17. Ware R.H., Exner M., Feng D., Gorbunov M., Hardy K.R., Herman B., Kuo Y.H., Meehan T.K., Melbourne W.G., Rocken C., Schreiner W., Sokolovskiy S.V., Solheim F., Zou X., Anthes R., Businger S. & Trenberth K. (1996). GPS Sounding of the Atmosphere from Low Earth Orbit: Preliminary Results, Bull. Am. Meteorol. Soc., Vol. 77, pp. 19-40.
  18. 18. Reigber, Ch., H. Luhr, and P. Schwintzer (2000). CHAMP mission status and perspectives, Suppl. to EOS, Transactions, AGU, 81, 48, F307.
  19. 19. Dieter K., Marc Cohen, Yves Buhler, Peter Schlüssel, Rosemary Munro, Axelvon Engeln, EoinÓ Clérigh, Hans Bonekamp, Jörg Ackermann, Johannes Schmetz and Juha-Pekka Luntama (2007). An Introduction to the EUMETSAT Polar system, Bull. Amer. Meteor. Soc., 88, 1085–1096. doi:
  20. 20. Anthes, R.A., Rocken, C., and Kuo, Y.-H. (2000). Applications of COSMIC to meteorology and climate, Special issue of Terrestrial, Atmospheric and Oceanic Sciences, 11, 115-156.

Written By

M.M. Shaikh, R. Notarpietro and Bruno Nava

Published: July 17th, 2014