Julia Hidalgo

Guillaume Dumas

Valéry Masson

Gwendall Petit

Benjamin Bechtel

Erwan Bocher

Michael Foley

Robert Schoetter

Gerald Mills

Benjamin Betchtel

Comparison between local climate zones maps derived from administrative datasets and satellite observations

Keywords: Local climate zones, MApUCE, WUDAPT, urban databases Research Highlights

come

Introduction

The presence of a city disrupts local and regional weather conditions, creating a distinctive urban climate at a hierarchy of scales by altering the surface-air exchanges of heat, moisture, mass and momentum [START_REF] Oke | The urban energy balance[END_REF]. This effect is present in nearly every single meteorological variable and its magnitude depends on aspects of urban form and function, that is, the physical structure and pattern of occupation of the city. The most widely studied urban climate effect is the urban heat island (UHI), which describes the differences in surface, sub-surface and air temperatures in cities when compared to the surrounding 'natural' environment [START_REF] Oke | The energetic basis of the urban heat island[END_REF]. Other impacts include changes to flow dynamics as the overlying air adjusts to a complex, rough and heterogeneous underlying surface which influences the dispersion of air pollutants and heat emitted near the ground [START_REF] Hidalgo | Advances in Urban Climate Modelling[END_REF]. Over the last few decades, urban climate science has made significant progress in linking the properties of the urban surface cover, including its extreme spatial heterogeneity, to changes in the overlying atmosphere [START_REF] Oke | Urban Climates. In Urban Climates (p. I)[END_REF]. Significant gaps in our understanding of processes remain, but it is generally acknowledged that the outstanding issue for urban climate science is the need to transfer knowledge into urban decision-making [START_REF] Ng | Adapting Asian Cities to Climate & Urban Climatic Changes[END_REF].

One obstacle that inhibits this transfer is the absence of coherent and consistent urban databases suited to urban climate studies. These databases would contain information on both urban form (land-cover, materials and building dimensions) and function (occupation patterns) that could be used to support observational and modeling projects and inform climate mitigation and/or adaptation policies. This data gap is especially acute in the rapidly growing cities in poorer parts of the world (where the need is greatest) but even in wealthy cities, the data infrastructure that could support knowledge transfer is often incomplete, inconsistent, spatially coarse, imprecise and frequently unavailable. Ideally, the available data would be climatically meaningful and acquired using consistent methods to allow cross-city comparisons. Examples of recent efforts along these lines are the Global Human Settlement Layer at 100m of resolution and the second generation Ecoclimap database (https://opensource.umrcnrm.fr/projects/ecoclimapsg/wiki).

One approach to getting suitable data for environmental and meteorological/climatic studies is to employ the Local Climate Zone (LCZ) typology to break the urban (and natural) landscape up into 'neighbourhoods' (ideally ≥ 1 km 2 ) that are relatively homogeneous in their make-up. The LCZ scheme, which has 10 urban classes and 7 natural ones, was designed to standardize the description of observation sites used in urban heat island studies [START_REF] Stewart | Local climate zones for urban temperature studies[END_REF], hereinafter referenced as SO12), but it has two significant attributes: it is a 'universal' classification with limited cultural bias, which means that it can be applied to cities worldwide; and each LCZ type is linked to a series of variables that relate to the urban climate effect generally. These variables include the sky view factor, the aspect ratio, the mean building or tree height, the terrain roughness class, the building surface fraction, the pervious and impervious surface fractions, the thermal admittance, the albedo and the anthropogenic heat flux. The associated values are often referred to as urban canopy parameters (UCPs) as they describe the attributes of the urban surface in a manner that can be integrated into climate models. Mapping cities into LCZ types is an increasingly common practice in urban climate research, but there has been no cross-comparison of the methods in use , which include satellite images, aerial photographs, administrative databases and fieldwork.

One category of approaches can be described as 'bottom-up' as they are based on information acquired for individual cities. One approach is to sample the urban landscape using fieldwork; this requires considerable expertise as shown by [START_REF] Houet | Mapping urban climate zones and quantifying climate behaviors-an application on Toulouse urban area (France)[END_REF] and [START_REF] Leconte | Using local climate zone scheme for UHI assessment: evaluation of the method using mobile measurements[END_REF] for Toulouse and Nancy, France respectively. More commonly, administrative or topographic data (on building footprints, heights, green spaces, etc.) and Geographic Information Systems (GIS) software are used. Raster-based methods superimpose a standard grid over the urban landscape and acquire information on selected variables (e.g. building height, sky view factor) at the scale of the individual cell. Each variable is stored in a layer and the gridded layers are combined using rules to generate LCZ types. This method has been used to generate LCZ maps for Hong-Kong [START_REF] Zheng | GIS-based mapping of Local Climate Zone in the high-density city of Hong Kong[END_REF], for Nagpur in India (Kotharkan and Bagade, 2017), for three medium-sized Central European cities, Brno, Hradec Králové, and Olomouc in the Czech Republic [START_REF] Geletic | GIS-based delineation of local climate zones: The case of medium-sized Central European cities[END_REF], and for Bilbao in Spain [START_REF] Acero | Urban climate modelling. Development of climate evaluation methods for urban planning purposes[END_REF]. Vector-based methods capture the boundary of an LCZ neighbourhood and represent a more precise delineation of contiguous neighbourhood types as individual objects; [START_REF] Unger | Local Climate Zone mapping using GIS methods in Szeged[END_REF] and [START_REF] Perera | Mapping "Local Climate Zones" and Relative Warming Effects in Colombo, Sri[END_REF] have used this method for Szeged and Colombo, respectively. In all of these studies, decision-making algorithms are employed that are based on typical UCP values proposed by SO12 (referred to henceforth as SO12). For a number of reasons, implementing a workflow to generate LCZ maps from urban data is not straightforward. First, the wide ranges of UCP values associated with each LCZ type overlap with those in other types. Second, the variable needed for the LCZ types may not be in the urban dataset and other exogenous indicators must be used as substitutes; for example, Geletic and Lehnert (2015) used indicators on the Number of buildings per hectare, the Number of areas of continuous surface of crown cover or the Number of continuous fragments of all vegetation per ha. Two of the required variables (albedo and anthropogenic heat flux) are especially difficult to obtain across the urbanised landscape and consequently are rarely included in the classification workflow. Third, the sequence of steps taken to classify a neighbourhood into an LCZ type is not fixed; some begin by classifying vegetated areas (Kotharkan and Bagade, 2017) whereas others start with the built categories (i.e. LCZ types 1 to 10). The advantage of these techniques is the opportunity to make best use of available data and local expertise to create LCZ maps with known accuracy and precision. The disadvantage is the variable nature of the underlying data, which makes cross-city comparison difficult and creating an international database a formidable challenge.

Another category of LCZ approaches are 'top-down' approaches, as they are based on satellite instruments that acquire information using passive or active sensors. Passive sensors respond to natural radiation that is emitted and reflected; active sensors record radiation emitted by the sensor that is reflected back from the target surface (e.g. synthetic aperture radar). The former can provide details on the Earth's surface when not obscured and are distinguished by their spectral and spatial resolution. Perhaps the best-known project is Landsat, which was launched in 1972. In its current mode (Landsat 8), it supports two instruments, the operational land imager and thermal infrared sensor, which collect information over 11 bands at spatial resolutions ranging from 15 to 100 m.

The World Urban Database and Access Portal Tools (WUDAPT) project utilizes these freely available data, which represent the lowest level (L0) of data in its information hierarchy, to map cities into LCZ types [START_REF] Ching | World Urban Database and Access Portal Tools (WUDAPT), an urban weather, climate and environmental modeling infrastructure for the Anthropocene[END_REF].

The protocol for generating LCZ maps employs supervised classification to automatically classify Landsat scenes using training data created by an urban expert [START_REF] Bechtel | Classification of Local Climate Zones Based on Multiple Earth Observation Data[END_REF]. Various tools within WUDAPT permit these data to be downloaded in a format suited to model applications; for example, [START_REF] Brousse | WUDAPT, an efficient land use producing data tool for mesoscale models? Integration of urban LCZ in WRF over Madrid[END_REF] employed L0 data for Madrid to run the Weather Research Forecasting (WRF) model. The WUDAPT process is efficient, so cities can be mapped quickly and updated with little effort. Given the global cover provided by Landsat, this process provides the potential for creating a worldwide urban database, but the derived value is predicated on the quality of the training areas.

In WUDAPT, higher levels (Level 1, Level 2) provide more precise information on shapes (building height, street length, etc.) and functions (office, industrial, residential, etc.) for the whole city in Level 1 and with intra-city spatial variability in Level 2; but these levels are not directly achievable using the previously mentioned methods, and in this case several bottom-up methods are used. This paper contributes to the ongoing research in this area by:

1. Proposing an original method for generating LCZ maps within a GIS framework using French administrative and topographic data acquired as part of the MApUCE project (section 2). The method takes regional variations into account and defines the UCPs associated with LCZ types more precisely. Here we present the LCZ output for three case study cities: Toulouse, Paris and Nantes (section 3).

Performing a geographical analysis and quantitative comparison of

MapUCE maps with WUDAPT-L0 maps for these three cities (section 4 and section 5).

The MApUCE urban database

The Applied Modelling and Planning Laws: Climate and Energy (MApUCE) project is designed to facilitate the integration of climate information into relevant French legal documents and urban policies. This information includes quantitative data on the climate and energy consumption at the building scale. Currently these data are not available or, if they are, are scarce; for example most official weather stations are located outside cities and fine-scale energy consumption data are often considered economically sensitive information. As a result, there is little climate information for urban places and energy data are only available over long periods and large territories.

MApUCE has two objectives: to build a database and to integrate climate knowledge into policy. The database consists of indicators that are pertinent for urban climate and energy consumption studies. These indicators include information on land use, urban morphology, building type, and user behavior [START_REF] Tornay | GENIUS: A methodology to define a detailed description of buildings for urban climate and building energy consumption simulations[END_REF][START_REF] Schoetter | Parametrisation of the variety of human behaviour related to building energy consumption in TEB[END_REF]. To overcome the obstacles to knowledge integration, legal and planning documents have been analyzed by experts to identify areas where regulations or incentives could be included that would promote climate-based actions. A few 'best cases' have been studied to understand the conditions that are needed for policy innovations and to encourage interactions between actors. Finally, based on urban planning agency requirements, the Urban Climate Maps methodology (see [START_REF] Ng | The urban climatic map: A methodology for sustainable urban planning[END_REF] is being adapted to the French urban and legal context.

The MApUCE project is interdisciplinary in its structure and draws on partners in the following fields: law, urban climate, building energetics, architecture, sociology, geography and meteorology, as well as the National Federation of Urban Planning Agencies (FNAU).

The morphological database description

The MApUCE morphological database [START_REF] Bocher | A geoprocessing framework to compute urban indicators: The MApUCE tools chain[END_REF]) contains a set of indicators computed at three urban scales, corresponding to building, block and islet (Figure 1) based on input dataset provided by French Institutes (Table 1).

Gridded population

The population data is obtained from the French National Institute for Statistics and Economic Studies (INSEE) (https://www.insee.fr/en/accueil).

Table 1. Datasets used to compute the morphological indicators in the MApUCE morphological database.

MApUCE uses vector-based geographical data on buildings, cadastral parcels and roads so that the spatial resolution is close to the scale of urban features. The database is built in three steps consisting in:pre-processing, to reduce inconsistencies and to perform islet extraction through a Voronoi tessellation algorithm. Islets, called herein RSUs, are delimited by the centre line of road surfaces.

1. the computation of 64 morphological indicators organized into 5 categories: count, area, shape, distance and others (Table 2).

2. the attribution of an urban typology class to each building. A set of Building typologies has been used as described in [START_REF] Tornay | GENIUS: A methodology to define a detailed description of buildings for urban climate and building energy consumption simulations[END_REF] to synthesize the main metrics and to illustrate the urban fabric organization. At the RSU scale, the percentage of each type is calculated and the primary and secondary types of building are identified. Figure 2 shows an example of the dominant building typology for a district in Paris.

Group

Figure 2. Dominant building type and aerial map (ESRI World Imagery) at the RSU scale for a district of Paris. Code of colours corresponds to Table 3 (Other type is represented with a transparent colour)

Case studies: urban units of Toulouse, Paris and Nantes

In MApUCE, the spatial extent of a city is defined by the notion of urban unit. In France, an urban unit (Fr: "Unité Urbaine") is a statistical area defined by the National Institute of Statistics and Economics Studies, for the measurement of contiguously built-up areas. It is defined as 'a municipality or a group of municipalities which includes a continuously built up zone (where constructions are not more than 200 meters apart) and at least 2,000 inhabitants'1 . This definition is in accordance with United Nations recommendations for the measurement of contiguous built-up areas and other comparable units used in the United States, the 'Urbanized Area' and the 'urban area' definition shared by the United Kingdom and Canada.

In this study three urban units were chosen as case studies: Paris, Toulouse and Nantes. Paris is the capital and most populous city of France; in the 2010 census it had over 10 million inhabitants across an area of 1749 km 2 .

Toulouse is the fourth-largest city in France with a population of over 800,000, occupying 1175 km 2 ; its large surface area, comparable to that of Paris but with one tenth the number of inhabitants, makes Toulouse one of the least dense cities in France. Nantes is the sixth-largest city in France with almost 600,000 inhabitants spread over 547 km 2 . The process is complicated by the nature of RSUs, which vary in size (Figure 2) and encompass one or more types of land-cover; to cope with the latter a semiautomatic method based on cluster classification has been developed. When applying a classification method, three key aspects will determine the results: the limits and specificities of the dataset (section 3.1), the hierarchy of the data treatment (section 3.2) and the definition of cut thresholds for each parameter (section 3.3). Each of these is discussed in turn.

MApUCE database specificities relevant for the LCZ classification

The various morphological indicators included in the MApUCE database exhibit different levels of precision and quality that will affect the LCZ methodology. These specificities include the following: The specificities of the MApUCE database will be taken into account in the following classification steps.



Hierarchy of the data treatment

The hierarchy proposed for the data treatment is as follows: first the RSUs with buildings are classified, followed by those that are mainly composed of water, roads and vegetation respectively.

The built-up RSUs are prioritised because of the reliability of these data and the need for an accurate discrimination among the urban LCZ types (LCZ 1 to 10).

Water cover (LCZ G) is easily obtainable as it is well defined in the MApUCE database. Impervious areas (LCZ E) are also treated before vegetation areas due to the high number of parking spaces, roads or cemeteries with trees. Vegetation areas are treated at the end, once all other covers are considered well represented. Non-classified RSUs do not clearly fit into any category, so during data processing they are called 'residue' and are assigned to LCZ D, Low Plants.

Data pre-treatment

Small and completely built-up RSUs are mostly situated in city centres; by comparison large RSUs are mostly situated in the outskirts and are often composed of a mix of vegetated and built-up areas. In order to avoid overclassification of these large RSUs into vegetated LCZs, data pre-processing is applied to all RSUs with buildings to bring out the artificial surface.

The objective is to identify the built areas. For this purpose, a buffer of 100 m around the building perimeter is applied to the building scale database. The choice of this radius is arbitrary but justified by the micro-climatic effect of single buildings [START_REF] Schmid | Spatial variability of energy fluxes in suburban terrain[END_REF]) when including the urban environment around them (car parks, roads and public places). These new polygons are intersected with the RSU limits in order to define urbanized areas inside the RSU (Figure 3.1). This makes it possible to evaluate the compactness of the built area within an RSU and better determine the LCZ density class while nevertheless retaining the MApUCE minimal surface unit, the RSU.

Some rules were applied: within the RSU, each buffer zone that doesn't touch another one from the same RSU becomes a unique built polygon, Figure 3(a). Each polygon receives a unique ID, and the Built density, the Number of buildings, the Mean building height, the mean Minimum distance, the median Minimum distance and the sum of each Building typology area are recalculated for each polygon. Any polygon that after being intersected by the RSU shape (black line in Figure 3

Characterization of built areas (LCZ 1 to 10)

To characterize these built polygons and attribute an LCZ type, eight parameters are used (Figure 4). Even if some links can be directly made between the 11 Building typologies defined in MApUCE and the ten "urbanized" LCZs (Tornay et al., 2017) --for example the "Informal building" corresponds to LCZ 7--this link is not unequivocal for the other LCZ types. Besides, the Majority Building typology cannot be the most representative within the RSU when thinking in terms of building morphology. To avoid this problem, Building typologies within an RSU are grouped together to obtain four classes based on their morphology (Figure 5).

The morphological groups are: open island buildings (Building typologies: Icio, Id and Bgh), closed island buildings (Icif and Pcif), houses (Psc, Pd and Pcio) and extended low-rise buildings. The surface fraction of these morphological groups is compared to the Majority Building typology indicator and a reclassification is applied if its value is higher. After this step, if "Extended low-rise" or "Informal building" is the main type of building, then LCZ 8 or LCZ 7 respectively are attributed. Once these two LCZs are allowed (Figure 4), the building typology is not used anymore and remaining urban polygons are separated using the Mean building height into the three height groups defined by SO12 (Low-rise: 0-10 meters, Midrise: 10-25 meters and High-rise: more than 25 meters). When exploring the relation between the LCZ classes and the Building density for the three case studies (Toulouse, Nantes and Paris), a high number of RSUs lie outside of the thresholds fixed by the SO12 LCZ type table look-up (Figure 6). These thresholds are sometimes lower and often too high. For example, this is the case for those RSUs where the Majority Building typology is "Extended low-rise building", that often have values of Building density under 30%. This is not surprising, as the generic ranges provided by SO12 are based on expert knowledge and existing inventories, and they are not meant to capture the heterogeneity in the global real world.

A semi-automatic statistical method to obtain characteristic LCZ thresholds for each city

The issue now is to find pertinent indicator thresholds to dispatch the urban polygons in three, two and two groups for the Low-, Mid-and High-rise classes defined by SO12 but taking into account the building spatial configuration (Building density, Distance between buildings) characteristic of French urban morphologies.

For that purpose a semiautomatic statistical method was developed based on cluster analysis. In this study the statistical k-means method is used. This clustering method is convenient here because the final number of clusters is already fixed by the definition of LCZ types (three groups for low-rise LCZ types and two groups for mid-rise and high-rise types). In the k-means method a distance measurement determines how the similarity of two elements is calculated when forming the clusters. Here the Euclidean distance is used.

To define the cluster structure the following indicators at the RSU scale are used: the Building density, the Mean and the Median of the Minimum Distance between buildings. To avoid differences in scale among the variables, the mean and standard deviations are used to normalize their values. The indicators characterizing the distance between buildings are needed to differentiate between LCZ 6 and LCZ 9 because two areas with the same building density can present very different spatial arrangements of buildings (Figure 7). Exploratory tests showed that, when applying the k-means method, the classification is better if we impose more than just three classes to be re-grouped afterwards. For this study 10 groups were applied for the high-rise buildings, 10 groups for the mid-rise buildings and 20 groups for the low-rise buildings. A supervised group in final 3, 2 and 2 groups for, respectively, low-, mid-and highrise classes proved a significant improvement and a finer classification, in particular for low-rise LCZ 3, 6 and 9 types.

The k-means method combined with a post supervised reclassification results in LCZ classes with thresholds that are characteristic for each site, thereby ensuring a better classification than if a single (and arbitrary) fixed threshold were applied. For example, in WUDAPT-L0, "compact" areas correspond to a Building density over 40%. The resulting mid-height groups for Nantes, Toulouse and Paris, treated with the same statistical method, are 36%, 37% and 33.8% respectively. Figure 6 presents the classification for Toulouse when using SO12 LCZ type table look-up versus the result using the proposed supervised statistical method.

Thresholds for Water, Pavement and Vegetation areas (LCZ A to G)

LCZ G is directly related to the Water fraction indicator based on the IGN BD TOPO as described in section 2. The RSU is classified here as LCZ G if the water fraction is higher than 50% of the RSU. Work is more arduous for paved surfaces, LCZ E, the road fraction is not reliable enough due to the linear definition of this indicator. It doesn't take into account the pavements or some car parks. Here the threshold value is fixed at 50% for the Road fraction indicator, considering that in reality it corresponds to over 60 to 70% of impervious surfaces in the RSU. Concerning the vegetated LCZ, with regard to the percentage of each vegetation class ("trees and bushes" and "grass"), if trees and bushes account for more than 60%, the Dense trees LCZ A class is attributed. If both vegetated classes are over 40%, Scattered trees LCZ B class is attributed. Polygons where no buildings stand and the amount of high and low vegetation isn't enough to be classified are considered residual and are attributed to the LCZ D Low Plants class.

Final LCZ MApUCE maps

Figure 8 shows the LCZ classification for Paris, Toulouse and Nantes based on administrative database from the MApUCE project. These maps reveal a typical European morphological structure. A denser city centre (mostly LCZ 2) is surrounded by large urban suburbs of residential housing (LCZ 6). The commercial and industrial developments (LCZ 8 and 10) are mostly concentrated along certain axes, which for the most part represent major traffic routes. However, there are also distinct differences. For instance, the dense urban core of Paris is large compared to the other two cities, while for Toulouse a large extensive development can be seen, but relevant patches of compact low rise (LCZ 3) and open mid rise (LCZ 5) can also be found adjacent to the town centre.

WUDAPT LCZ maps

The WUDAPT project gathers and organises urban data by level of detail. The lowest level of data (Level 0 or L0) breaks down urban regions into LCZ types using a supervised classification method applied to LANDSAT multi-spectral sensing data; the result of the process is a raster-based database in which each cell (100 m) has an LCZ value. Higher levels of detail in WUDAPT will provide more data on the features (e.g. buildings and trees) that make up urban areas and are needed by sophisticated climate models. As the WUDAPT-LCZ data are fundamentally based on the opinions of experts, the product is tested for reliability using a bootstrapping method [START_REF] Kaloustian | Local Climatic Zoning and Urban Heat Island in Beirut[END_REF]) that generated a set of accuracy measures [START_REF] Ching | World Urban Database and Access Portal Tools (WUDAPT), an urban weather, climate and environmental modeling infrastructure for the Anthropocene[END_REF]. While this approach can detect inconsistencies in the classification, it does not represent an independent assessment of the WUDAPT-LCZ product, which would require alternative sources of data, such as that available within MApUCE.

The structure and status of the WUDAPT project are described elsewhere [START_REF] Ching | World Urban Database and Access Portal Tools (WUDAPT), an urban weather, climate and environmental modeling infrastructure for the Anthropocene[END_REF]); here we focus on the WUDAPT-LCZ data for the three case-study cities. Figure 9 shows the results for the three case studies, Toulouse, Nantes and Paris; the class frequency distributions (not shown) are dominated by the natural classes (in particular LCZ D, low plants) but it is important to recognise that this is largely a product of the size of the area under study. However, differences between cities are evident in the urban LCZ classes as well, Table 5: 45.7% of urban classes (LCZ1 to 10) in Paris are classified as compact (LCZ 1-3), compared to 0.37% in Nantes and 0.73% in Toulouse. Nantes, due to its old industrial port close to the city centre, has a large percentage (15.9%) of warehouse and industrial areas (LCZ8 and LCZ 10), compared to 10.62% in Paris and 9.89% in Toulouse. However, these findings are somewhat distorted by the large share of sparsely built areas (LCZ 9) in the Toulouse map (27.9%), which mostly consist of natural cover and are difficult to classify [START_REF] Ching | World Urban Database and Access Portal Tools (WUDAPT), an urban weather, climate and environmental modeling infrastructure for the Anthropocene[END_REF]. Detailed cover fractions can be found in Table 5 in Appendix 1.

Comparing MApUCE-LCZ with WUDAPT-LCZ

This section compares results produced by both approaches and highlights the potential and specificities of each. To that end, WUDAPT maps are vectorized in order to obtain polygons comparable to those in MApUCE maps. As stated in section 3, boundaries of vegetated areas for the MApUCE version are distorted by the RSU cutting definition. This comparison focuses only on built areas. Two types of analysis are presented, based on the total and individual LCZ surface and on a geographical analysis based on the Building Density and Building Height class.

A. Toulouse study case

Figure 10a and Table 5 shows that total built areas are underestimated by 26.3% in the WUDAPT approach. This significant difference should be put in perspective, as the main differences appear in small suburban areas where there are few buildings and the choice of a 100-m buffer could contribute to overestimating built areas in these zones. Compact mid-rise class, LCZ 2, is under-estimated in WUDAPT map and the city centre of Toulouse is the area the most impacted by this under-representation. In WUDAPT approach LCZ 6 is also under-represented comparing to MapUCE even if it has tendency to classify the large low-rise areas (LCZ 8) into LCZ 6. This is not the case with the MApUCE method. Nevertheless, those isolated big buildings that corresponds to agricultural activities and are then surrounded by big vegetated areas are classified by MApUCE method in LCZ 8 while it would be more suitable to be classified in LCZ9. For Toulouse, the WUDAPT approach catch those subtleties, but that strongly depends on the chosen training areas.

Sparsely built, LCZ 9, is less prevalent in MApUCE representation. The use of the mean and median minimum distances between buildings and the islet referential in the classification implies an under-representation of this typology compared with the WUDAPT standards. French urban models for residential neighbourhoods are based on building distances of 5 to 15 meters. It is of course possible to find buildings that are more than 20 meters apart, but it will be sporadic. That's why the representation of LCZ 6 and 8 is higher in the MApUCE method, and the LCZ 9 class is less prevalent than in the WUDAPT method. B)

Figure 10. Geographical comparison for Toulouse based on: (A) the built areas classification obtained through MApUCE and WUDAPT-L0 methods and (B) the

Building Density and Building Height classes.

When comparing by Building height or Building density classes, some differences clearly appear (Figure 10b). The compact zones (LCZ 1, LCZ2 and LCZ3) are more prevalent on the MApUCE map and allow for the capture of small town centres, which is very important for modeling micro-climatic effects on town peripheries. The road structure to the north, south-east and south-west, appears clearly when analysing both height class maps. In this sense the MApUCE map shows better resolution and can capture the internal large avenues lined with high buildings represented by the internal red arcs in the centre of Toulouse.

B. Paris and Nantes case studies

Due to the size of the Paris Urban Unit (see figures in Appendix 2), only the city centre and Greater Paris are analysed here. As dense built areas are captured well by the WUDAPT method, and this area is highly artificial, both methods capture comparable urban extension. In terms of LCZ typologies, Paris is almost exclusively composed of Compact Mid Rise buildings, LCZ 2, overestimated of about 15.8% compared to MApUCE. Differences in LCZ 5 and LCZ 6 extension, are quite similar, under representation of WUDAPT with 36.2% and 25.2% respectively. The biggest differences being found for LCZ 3 and LCZ 4 with a high over-representation by the WUDAPT methodology, 215% and 414% respectively. In France, most LCZ 4 classes correspond to areas called in French "les grands ensembles" ("large units"), a particular type of 1960s urban planning situated near the roadside belt. Even if training areas were taken from that kind of territory, the WUDAPT method classified zones between dense and open as LCZ 4 instead of LCZ 5.

The MApUCE method shows a balanced profile with four main LCZs: Compact Mid Rise, Open Mid Rise, Large Low Rise, and Open Low rise; predominant LCZs for the WUDAPT maps are: Compact Mid Rise and Open Low Rise.

For Nantes, built areas are particularly under-represented in the WUDAPT method, with a surface difference of 37%. Nevertheless, the LCZ distribution identified by both methods is not so different. The most frequent LCZs are the same (LCZ6, LCZ8 and LCZ9) but their hierarchy slightly differs. First LCZ 6 Open Low Rise, is the most frequent in both approaches. Then, LCZ 8 Large Low Rise, and lastly, LCZ 9 Sparsely Built are the most frequent in MAPUCE while is the opposite in WUDAPT. This is explained by the fact that, for Nantes, the WUDAPT methodology identified as LCZ9 some areas classified as LCZ6 and LCZ8 by MApUCE method. Finally, density classes (Compact, Open, Sparse, Large Low Rise) or height classes (High-rise, Mid-rise, Low-rise, Large Low Rise) are correctly spatialized for this territory.

Table 5. Percentage of total and individual LCZ surface for MApUCE (M) and WUDAPT-L0 (W) maps.

Detailing LCZ indicators for French cities from the MApUCE database

In this section the analysis is focused on the LCZ ranges of urban canopy model parameters fixed on the table lookup in SO12. The objective here is to detail the typical borders of the interval of LCZ indicators for French cities.

Out of the ten indicators proposed in the table lookup in SO12, the MApUCE morphological DB calculates seven:

 the Mean building height  the Building surface fraction  the Pervious surface fraction  the Impervious surface fraction  the Aspect ratio H/W  the Roughness class  the Sky view factor However, in the morphological DB in MApUCE there are no quantitative elements with which to estimate Albedo, Surface admittance and Anthropogenic heat flux.

In order to characterize the medium-size cities of the conurbation as well, the first five indicators were computed for both the central commune and the immediate periphery of each case study. In this study, the Sky view factor was not explored.

A) Indicator calculation

Two LCZ indicators already feature on the MApUCE DB, the Mean building height and the Building surface fraction (called Building density in the MApUCE DB). The Building density corresponds in this study to the built surface with respect to the surface of the buffer polygon. This can be an entire RSU, for completely urbanized ones close to the city center, or a fraction of it for those of the rural areas that were divided into several zones during the buffering process. The Mean building height and the Building density are then recalculated for each polygon.

The Aspect ratio corresponds to the Mean building height of the polygon divided by the Mean minimum Distance between buildings. As in Stewart and Oke's study, values higher than three are truncated.

The Pervious and impervious surface fraction are calculated using the vegetation maps at 1.5m resolution combined with the buffer zones and the building footprints. In this case the "no vegetation" category is used. On the central commune and the immediate periphery where there is no unattributed land (or very little), we assume that this category accurately represents the mix of buildings and paved surfaces. To calculate the impervious surface fraction indicator, the building footprints are extracted from the buffer zones and zonal statistics are applied to the intersection of this information and the "no vegetation" category on the vegetation map. The Pervious surface fraction indicator over the buffer zone is calculated as 1-impervious surface fraction -the building footprint surface.

B) Presentation of results

Three levels of analysis are presented and discussed, in terms of: overall LCZ spatialization (frequency count of LCZ's), the mean indicator magnitude for the central commune and the communes situated in the immediate periphery (also called first crown), and finally the indicator variability at the LCZ scale.

In terms of overall LCZ spatialization (Table 6), it is possible to observe morphological differentiation linked to the administrative division. For example, for Paris, if the town centre is composed of 78% of LCZ 2 Compact Mid Rise, the first crown comprises 37% and the enterprise zones range from 3.8% to 9.5%. LCZ 5 and 6 account for a mere 9.1% and 1.25% respectively in the town centre, but increase to 16.3% and 19.9% in the first crown. This type of analysis can be repeated on Toulouse and Nantes. This makes it possible to observe the overwhelming presence of LCZ 6 on the peripheral territory with 78.6% of the LCZs coupled with LCZ 8 at 11.2%. For the centre of Toulouse, the LCZ 6 increases to 44.6% and the LCZ 2 appears at 25.7%. For Nantes, we find the same phenomenon with LCZ 6 and 2 for the town centre accounting for 53.9% and 15.1% and an LCZ 6 predominant with 73.1% on the periphery. It is interesting to note that the profiles of the town centre and the periphery of Toulouse and Nantes are very similar. This can be seen in the general profile presented in the previous section, but when the territory is broken up it is all the more striking. When focusing on the mean indicator magnitude for the central commune and the immediate periphery, we may observe that (Table 7):

 For the three sites and all LCZs, the Mean building height values fall on the "standard" thresholds fixed by SO12 because in our classification method the heights are fixed according to these LCZ thresholds (0-10m 10-25m and > 25m). It can also be observed that buildings are consistently higher for the city centre than for the periphery. The variability seems to be higher for Paris, in particular for LCZ 4 and 8, than for Toulouse or Nantes.

 In general terms for these French cities, values of the Building surface fraction seem to be close to the standard ones (the mean value falls on the range) for LCZ 1, 2 and 3. Values are slightly smaller (the mean value ±σ, falls on the range) for LCZ 4, 5 and 6. And values are far from the standard ones (the mean value ±σ, does not fall on the range) for LCZ 8 and 9. Evidently differences may be observed between the cities, for example, differences in values for the city centre of Toulouse and for the periphery of Nantes are higher.

 The partition of impervious and pervious surfaces in the standard look-up table seems to be quite well represented (all the values fall on or are close to the standard range) for Paris and Nantes, in particular for LCZ 2, 3, 5 and 6. Large low-rise (LCZ 8) seems to be consistently more vegetated, in particular for Toulouse.

 From the Aspect ratio values, French cities seem to be denser than the average, in particular for central towns where higher values of H/W are found. The H/W values fall on the standard threshold for just one LCZ --LCZ1, situated on the peripheral communes in Paris. While values are not so far from the standard ones in Toulouse and Nantes, values for LCZ 3, 4, 5, 6 and 8 are invariably higher for the periphery of Paris.

An example of sheets presenting ranges of urban canopy parameters according to the framework proposed by SO12 can be found in Appendix 3.

Table 7. Mean value of the indicators and its standard deviation for the central commune and the immediate periphery from the MApUCE database. Green colour indicates that mean value falls on the "Standard Threshold" (ST) proposed by SO12. Orange that the mean value ±σ, falls on the ST and red colour that the mean value ±σ, does not fall on the ST.

Indicator variability at the LCZ scale. Figure 11 clearly shows the advantage of shifting from WUDAPT level 0 to WUDAPT level 1. Effectively, in actual fact there is an intra-LCZ variability that should be described to properly evaluate the impact of the urban structure on the atmospheric boundary layer. In this example, for a neighbourhood situated in the north-western part of Toulouse, it can be clearly seen that the Mean building height, the Mean building fraction and the Impervious surface fraction for a given LCZ varies without correlation between the parameters.

Figure 11. Intra-LCZ variability in a neighbourhood of Toulouse

Conclusions

This paper aims to contribute to the scientific discussions of the pertinence and limits of the various ways ("bottom-up" vs "top-down") of obtaining Local Climate Zone maps.

From a methodological point of view the paper presents an original semiautomatic method based on cluster classification to obtain GIS-based Local Climate Zone maps. For this purpose, a rich administrative dataset developed during a national funded project, the MApUCE project, is used to obtain MApUCE-LCZ maps for three French case studies, Paris, Nantes and Toulouse. The choice was made to work with vector data that offer better comprehension of the boundaries of urban elements than rasterized data and make it easier to share results with practitioners. The work-flow in the classification is as follows: first areas with buildings are classified and then those composed mainly of water, roads and vegetation respectively. The paper also presents maps obtained through the WUDAPT level 0 methodology and comment the differences observed.

Geographical analysis based on the comparison between MApUCE and WUDAPT results for total and individual LCZ surface and Building Density and Building Height indicators shows that:

-Compared to the MApUCE approach, in general terms, WUDAPT L-0 underrepresent urbanized areas.

-Main differences appear in the small suburban areas where there are few buildings. These differences are mainly concentrated on LCZ 6, 8 and 9.

-Even if the city cores are quite well represented in the WUDAPT approach in terms of urban extension, MApUCE has better resolution for these areas and the structures of town centres are consequently better captured.

-In the WUDAPT approach, for Toulouse key differences are concentrated on LCZ 6. For Paris LCZ 4 is highly over-represented; and for Nantes, even if LCZ distribution were comparable, built areas are particularly under-represented. SO12 provided general LCZ ranges for ten urban canopy model parameters. The morphological database in MApUCE allows users to verify and detail these ranges for French cities and to analyse their intra-LCZ variability. Main findings are:

-Relatively significant morphological differentiation was observed between the central commune and the immediate periphery. The geographical domain used in numerical simulations of Urban Heat Island, and more generally urban micro-climate, used to be larger than the size of the city, in order to capture the local meteorological interactions between the city and the countryside. Therefore, it seems important to widen the urban structure analysis to peripheral cities when preparing the urban surface input data for climatic simulations.

-As already pointed out by other authors, to obtain satisfactory LCZ classification, it is important to take into account specific regional features concerning the borders of the intervals of the LCZ physical properties. When available, a high resolution administrative database is a solution; but if it is not available, other solutions exist based on manual sampling from field campaigns or crowd-sourcing and multi-source satellite images.

-Even if the LCZs are supposed to have a homogeneous climatic footprint on the near-surface atmosphere, their characteristic size is over 500x500m and, for modeling purposes, the internal variability of the building structure inside this area must be evaluated and represented in the urban surface input data. Indeed, while WUDAPT level 0 morphological parameters may be sufficient for Numerical Weather Prediction and Very-High resolution models, both typically using 1 to 4km resolutions, hectometric-scale and Large eddy Simulation models are now currently used in the urban climate community for micro-climate studies, and a finer description of the city internal structure is desirable.

The advantage of the use of the Local Climate Zone classification is the use of the same referential on the urban structure comprehension that allows comparison of morphological portraits. On this subject, Toulouse and Nantes are two cities which possess many similar morphological characteristics. The difference between the two has more to do with materials and their characteristics: white stones and light soils in Nantes, red brick and dark pavement for Toulouse, for example. It is at this level of detail that the limits of the morphological database in MApUCE or WUDAPT-L0 maps can be seen, and it is for this reason that a complementary architectural database was also developed in MApUCE. This architectural database is based on building use, building construction date and geographical location [START_REF] Tornay | GENIUS: A methodology to define a detailed description of buildings for urban climate and building energy consumption simulations[END_REF]) and makes it possible to architecturally differentiate cities and areas within a city during modeling exercises. The need for increase in architectural description arises in the WUDAPT community, and two ways are currently used to fill this gap: architectural expertise (as in Tornay et al 2017), that can be organized at global scale, and crowdsourcing, in order to gather some information for many buildings in each LCZ for many cities. Table 1. Datasets used to compute the morphological indicators in the MApUCE morphological database.

Appendix 1. WUDAPT maps meta data and evaluation results

Nantes

Figure 1 .

 1 Figure 1. spatial scales considered in the MApUCE project. Building (a), block (b)

  (b)) has no building on it, is integrated into the vegetated shape (visible in the right part of Figure 3(c)).

Figure 3 .

 3 Figure 3. Rules applied to outlined built areas inside the RSU. Black lines correspond to the RSU limits, grey stain corresponds to buffer surface and black stain to building footprint.

Figure 4 :

 4 Figure 4: Work-flow in the LCZ attribution process

Figure 5 .

 5 Figure 5. Morphological groups used to verify if the Majority Building typology corresponds to the majority morphology. Those groups are also used to identify LCZ 8 and LCZ 7.

Figure 6 .

 6 Figure 6. Classification of built areas (LCZ1 to LCZ10) for Toulouse using: left, thresholds fixed by Stewart and Oke's LCZ type table look-up; right, thresholds obtained through the supervised statistical method proposed in section 3.5.

Figure 7 .

 7 Figure 7. Example of built areas (in blue and green) with similar building density (3%) but very different spatial configuration.

Figure 8 .

 8 Figure 8. LCZ classification for (a) Nantes, (b) Toulouse and (c) Paris. Due to the larger spatial extent of the Paris Urban Unit, only the city centre and the Greater Paris area are displayed.

Figure 9 .

 9 Figure 9. WUDAPT-L0 Local Climate Zones maps at 100x100 m resolution for a) Nantes; b) Toulouse and c) Paris

Figures:Figure 1 .

 1 Figures: Figure 1. spatial scales considered in the MapUCE project. Building (a), block (b) and islet, also called Reference Spatial Unit (RSU, c) Figure 2. Dominant building type and aerial map (ESRI World Imagery) at the RSU scale for a district of Paris. Code of colours corresponds toTable 3 (Other type is represented with a transparent colour) Figure 3. Rules applied to outlined built areas inside the RSU. Black lines correspond to the RSU limits, grey stain corresponds to buffer surface and black stain to building footprint.

Figure 4 .

 4 Figure 4. Work-flow in the LCZ attribution process.

Figure 5 .

 5 Figure5. Morphological groups used to verify if the Majority Building typology corresponds to the majority morphology. These groups are also used to identify LCZ 8 and LCZ7.

Figure 6 .

 6 Figure 6. Built areas classification for Toulouse using: left, thresholds fixed by the Stewart and Oke's LCZ type table look-up; right, thresholds obtained through the proposed supervised statistical method proposed in section 3.5.

Figure 7 .

 7 Figure 7. Example of built areas with similar building density (3%) but very different spatial configuration.

Figure 8 .

 8 Figure 8. LCZ classification for Paris, Toulouse and Nantes respectively. Due to the larger spatial extent of the Paris Urban Unit, only the city centre and the Greater Paris area are displayed.

Figure 9 .

 9 Figure 9. WUDAPT-L0 Local Climate Zones maps at 100x100 m resolution for a) Nantes; b) Toulouse and c) Paris Figure 10. Geographical comparison for Toulouse based on : (A) the built areas classification obtained through MapUCE and WUDAPT-L0 methods and (B) the Building Density and Building Height classes.

Figure 11 .

 11 Figure 11. Intra-LCZ variability in a neighbourhood of Toulouse

Table 2 .

 2 Examples of morphological indicators for the different groups.

		Indicators	Examples
	Count	Number of	Number of buildings and blocks by RSU
		features
	Area	Areas or ratio of	Building floor; exposed building facades, road
		areas	fraction, vegetative fraction, etc. by RSU
	Shape	Volume,	Building volume, block compactness and
		compactness,	mean building height by RSU
		form factor,
		concavity, fractal
		dimension,
		contiguity
	Distance	Distance between	Distance between buildings, between
		one feature and n	buildings and roads for each RSU
		features
	Others	Passive volume,	Building or blocks main direction; Block
		Main direction,	courtyard area
		Holes area

Table 3 . Building typology used in MApUCE as described in Tornay et al. 2017. Identifiers were not translated from French to English to facilitate

 3 the identification of typologies in the MApUCE database.

	Description	Identifier and
		Legend
	Extended low-rise	Ba
	Industrial, commercial or agricultural buildings. They are characterised
	by their simple morphology and their large
	Footprint
	High-rise building	Bgh
	A building of more than 12 storeys such as an apartment tower or
	office tower
	Continuous row of mid-rise	Icif
	Perimeter islet development: Connected buildings with street front
	elevations. This typology is often present in historic centres, urban
	fabric of the industrial
	Discontinuous row of mid-rise
	A building complex in the centre of the urban islet
	7

Icio

Detached mid-rise

One or more buildings built in the centre of the islet

Id

Informal building

Ephemeral constructions, non-traced on registers (caravans, temporary prefabricated buildings, etc.)

Local

Continuous row of low-rise

Typical intermediary housing, terraced houses with patios, constructions typical of historic centres Pcif Discontinuous row of low-rise Street aligning terraced houses with gardens at the back Pcio Detached low-rise One or two-storey houses of at least four façades often located in the centre of plot of land Pd Semi-detached low-rise Town houses, terraced houses or houses detached on one side, with façades aligning the street Psc Other na

Table 4

 4 presents the characteristics of the MApUCE database for the three cities.

			Urban unit
		Toulouse	Paris	Nantes
	Number of communes	100	334	38
	Number of RSUs	15 611	61 029	15 724
	Number of
	RSUs with	12 366	47 357	10 626
	buildings
	Number of buildings	218 240	1 607 631	262 643
	Number of blocks	206 280	1 00 907	169 281
		Dominant Building Type RSU
	Ba	1355 (10.9%)	3378 (7.1%)	1039 (9.8%)
	Bgh	20 (0.2%)	321 (0.7%)	15 (0.1%)
	Icif	606 (4.9%)	7151 (15.1%)	391 (3.7%)
	Icio	1411 (11.4%)	8178 (17.3%)	846 (8.0%)

The number of RSUs and buildings in Paris is much larger than for the other cities; for this reason the results in section 5 are only shown and discussed for the city centre and Greater Paris. Toulouse has the smallest number in terms of number of buildings, but a far higher percentage categorised as detached lowrise (pd), confirming its low density. High-rise buildings (bgh) are rare in France; continuous and discontinuous rows of mid-rise buildings (icif and icio) are the most prevalent building types in city centres. Extended low-rise buildings (ba) account for around 10% and informal urbanization is almost non-existent (~1%) compared to other parts of the world.

Table 4 .

 4 Characteristics of the MApUCE database for Toulouse, Paris and Nantes3. Deriving Local Climate Zones using administrative datasetsTo classify each RSU into one of the 17 LCZ types using the MApUCE database, indicators at the RSU and building scales are used to calculate four types of land cover: built-up areas, water bodies, vegetated areas and impervious surfaces.

  The morphological indicators are calculated at the scale of the RSUs, which are roughly delimited by roads. As a result, some building types, such as the ''Informal Building'' class, are unlikely to become the dominant building type at RSU scale. and 2017. Using the methodology and software developed by[START_REF] Crombette | Traitement d'images satellitaires à Très Haute Résolution Spatiale et identification de zones à enjeux dans l'aménagement des Trames Vertes urbaines[END_REF], the fusion of panchromatic and multi-spectral band images, it provided high vegetation (trees and bushes), low vegetation (grass) and no vegetation classes. No vegetation can correspond to bare ground or impervious surfaces. The available vegetation data thus allows users to distinguish three LCZ types: dense trees (LCZ A), scattered trees (LCZ B) and low plants (LCZ D), but not bush/scrub (LCZ C) and bare ground (LCZ F). Moreover, the satellite images available are not the same for each city and this can have a strong impact on the vegetation representation. Images for Toulouse were taken in June, for Nantes in September and for Paris in May (West side) and August (East side).

	 The MApUCE database lacks information on buildings with very large
	energy throughput like power plants, steel or aluminium factories, garbage
	incineration plants, and so on. For this reason, we cannot identify heavy
	industry (LCZ 10) using the MApUCE database only.

The position of water bodies and buildings is for the most part well represented in the IGN-BD TOPO. For this reason, the indicators Building number, Building density, Building height, and Water fraction are mostly reliable and will therefore be given priority in the LCZ classification.

 The representation of impervious surfaces in the IGN-BD TOPO is not comprehensive. Main roads are included, but pavements and car parks are missing. The impervious surfaces are therefore underestimated in the MApUCE database.

 

The urban vegetation was retrieved via SPOT 6-7 satellite images provided by EQUIPEX-GEOSUD (http://ids.equipex-geosud.fr/) at

1.5mX1.5m

resolution between 2015

Table 6 .

 6 Detailed LCZ distribution based on the central commune and the immediate periphery of each case study from the MApUCE database.

Table 2 .

 2 Examples of morphological indicators for the different groups.

Table 3 .

 3 Building typology used in MApUCE as described in[START_REF] Tornay | GENIUS: A methodology to define a detailed description of buildings for urban climate and building energy consumption simulations[END_REF]. Identifiers were not translated from French to English to facilitate the identification of typologies in the MApUCE database.

Table 4 .

 4 Characteristics of the MApUCE database for Toulouse, Paris and Nantes.

Table 5 .

 5 Percentage of total and individual LCZ surface for MapUCE and WUDAPT -L0 maps.

Table 6 .

 6 Detailed LCZ distribution depending on the central commune and the immediate periphery of each case study.

Table 7 .

 7 Mean value of the indicators and its standard deviation for the central commune and the immediate periphery from the MApUCE database. Green colour indicates that mean value falls on the "Standard Threshold" (ST) proposed by SO12. Orange that the mean value ±σ, falls on the ST and red colour that the mean value ±σ, does not fall on the ST.

https://www.insee.fr/en/metadonnees/definition/c1501

Acknowledgements

This work was funded by the French National Research Agency (Agence Nationale de la Recherche) under grant #ANR-13-VBDU-004.

Appendix 2. Geographical comparison between MApUCE and WUDAPT classification for Paris and Nantes on : (A) the built area classification obtained through MApUCE and WUDAPT-L0 methods and (B) the Building Density and


Laurent Gizon
email: gizon@mps.mpg.de

Damien Fournier

Dan Yang

Aaron C Birch

Hélène Barucq

Signal and noise in helioseismic holography

Keywords: helioseismology -Sun, interior -Sun, oscillations

Context. Helioseismic holography is an imaging technique used to study heterogeneities and flows in the solar interior from observations of solar oscillations at the surface. Holograms contain noise due to the stochastic nature of solar oscillations. Aims. We provide a theoretical framework for modeling signal and noise in Porter-Bojarski helioseismic holography. Methods. The wave equation may be recast into a Helmholtz-like equation, so as to connect with the acoustics literature and define the holography Green's function in a meaningful way. Sources of wave excitation are assumed to be stationary, horizontally homogeneous, and spatially uncorrelated. Using the first Born approximation we calculate holograms in the presence of perturbations in sound-speed, density, flows, and source covariance, as well as the noise level as a function of position. This work is a direct extension of the methods used in time-distance helioseismology to model signal and noise. Results. To illustrate the theory, we compute the hologram intensity numerically for a buried sound-speed perturbation at different depths in the solar interior. The reference Green's function is obtained for a spherically-symmetric solar model using a finite-element solver in the frequency domain. Below the pupil area on the surface, we find that the spatial resolution of the hologram intensity is very close to half the local wavelength. For a sound-speed perturbation of size comparable to the local spatial resolution, the signalto-noise ratio is approximately constant with depth. Averaging the hologram intensity over a number N of frequencies above 3 mHz increases the signal-to-noise ratio by a factor nearly equal to the square root of N. This may not be the case at lower frequencies, where large variations in the holographic signal are due to the individual contributions of the long-lived modes of oscillation.

Introduction

Local helioseismology aims at studying the solar interior in three dimensions by exploiting the information contained in the waves observed at the solar surface (e.g., Gizon & Birch 2005, and references therein). Helioseismic holography is one particular approach of local helioseismology, which images subsurface scatterers by back-propagating the surface wave field to target points in the interior. Helioseismic holography is also known as Lindsey-Braun (LB) holography (Lindsey &Braun 1997, 2000a, andreferences therein). It has been used to study solar convection [START_REF] Braun | SOHO 14 Helio-and Asteroseismology: Towards a Golden Future[END_REF](Braun et al. , 2007)), active region emergence (Birch et al. 2013(Birch et al. , 2016)), sunspot subsurface structure (Braun & Birch 2008b;Birch et al. 2009), to image wave sources [START_REF] Lindsey | Solar MHD Theory and Observations: A High Spatial Resolution Perspective[END_REF], to study sunquakes caused by solar flares (Zharkov et al. 2013;Besliu-Ionescu et al. 2017), and to detect active regions on the far side of the Sun (Lindsey & Braun 2000b;Liewer et al. 2014).

In acoustics, a well-established version of holography in a medium that contains sources is Porter-Bojarski (PB) holography (Porter & Devaney 1982). PB holography uses both the wave field and its normal derivative at the surface to produce holograms (Porter 1969). PB holography was introduced in helioseismology by [START_REF] Skartlien | [END_REF]Skartlien ( , 2002)), where deterministic sources and scatterers were recovered in a solar background. Yang (2018) recently studied PB holograms in a homogeneous medium permeated by localized deterministic sources to study ghost images near the observational boundary.

In this paper we apply PB holography in a realistic helioseismological setting. First we rewrite the wave equation in Helmholtz form, in order to properly define the Green's functions that are involved in the definition of PB holograms. The background density and sound-speed are taken from a standard solar model. The model of wave excitation is described by a reasonable source covariance function, which leads to a solar-like power spectrum for acoustic oscillations.

The signal is defined as the expectation value of the perturbations in hologram intensity that result from perturbations in sound speed, density, and flows with respect to the reference solar model. The corresponding sensitivity kernels are computed in the first-order Born approximation (Gizon & Birch 2002;Birch & Gizon 2007;Braun et al. 2007;Birch et al. 2011). This signal must take into account the correlations between incident and scattered wave fields, which are both connected to the sources of excitation (turbulent convection).

Random noise in holograms is due to the stochastic nature of the sources of excitation. While noise can sometimes be estimated from the data (Lindsey & Braun 1990;Braun & Birch 2008a), a theoretical understanding is useful to design holographic experiments. Here we extend to holography the noise model developed in time-distance helioseismology (Gizon & Birch 2004;Fournier et al. 2014). We do not attempt to image individual sources as in Skartlien (2002), which in our view is not a Article number, page 1 of 8 arXiv:1810.00402v1 [astro-ph.SR] 30 Sep 2018 well-posed problem (see also [START_REF] Lindsey | Solar MHD Theory and Observations: A High Spatial Resolution Perspective[END_REF], except in case of imaging the sources of sunquake waves. Instead we consider sources to be specified through a source covariance function.

Reduced wave equation

At angular frequency ω and spatial position r in the computational domain V, the propagation of acoustic waves in a 3D heterogeneous moving medium is described by the displacement vector ξ(r, ω), solution to

-(ω + iγ + iu • ∇) 2 ξ - 1 ρ ∇ ρc 2 ∇ • ξ + gravity terms = F, (1)
where ρ(r) and c(r) are the density and sound speed, and u(r) is a steady vector flow. Wave attenuation is included through the function γ(r, ω). The source term F(r, ω) is a realization from a random process; it describes the stochastic excitation of the waves by turbulent convection. Following Lamb (1909) and Deubner & Gough (1984), we consider the scalar variable

ψ = ρ 1/2 c 2 ∇ • ξ, (2)
to recast the wave equation into a Helmholtz-like equation

Lψ := -(∆ + k 2 )ψ - 2iω ρ 1/2 c ρu • ∇ ψ ρ 1/2 c = S , (3)
where S = ρ 1/2 c 2 ∇ • F is a scalar source term. The local wavenumber k(r, ω) is given by

k 2 = (ω 2 + 2iωγ) -ω 2 c c 2 , (4)
where the squared acoustic cut-off frequency is

ω 2 c = ρ 1/2 c 2 ∆(ρ -1/2 ). (5)
In obtaining Eq. (3), we ignored gravity terms and assumed slow variations of c, u, and γ compared to the wavelength (Gizon et al. 2017). The advection term is such that the corresponding operator is Hermitian symmetric for the inner product ψ 1 , ψ 2 = ψ * 1 ψ 2 dV under the conditions that the flow conserves mass and that it does not cross the boundary (u n = 0 on ∂V).

The stochastic sources of excitation are assumed to be stationary and spatially uncorrelated. They are described by a source covariance function of the form

E[S * (r, ω)S (r , ω)] = M(r, ω)δ(r -r ). (6)
To solve Eq. ( 3), one needs to specify a boundary condition at the computational boundary. Here we apply an outgoing radiation boundary condition

∂ n ψ = ik n ψ on ∂V. (7)
We apply the boundary condition (Atmo RBC 1) from Barucq et al. (2018), which assumes an exponential decay of the background density at the boundary of the domain but neglects curvature. Then, the local wavenumber k n from Eq. ( 7) is given by

k 2 n = ω 2 + 2iωγ c 2 - 1 4H 2 , ( 8
)
where H = -1/(d ln ρ/dr) is the density scale height at the boundary. The last term in Eq. ( 8) is connected to the cut-off frequency for an isothermal atmosphere (Lamb 1909), thus k n is an approximation of the wavenumber k from Eq. (4). Fournier et al. (2017) discusses several of the boundary conditions derived in Barucq et al. (2018).

Hologram and hologram intensity

The following calculations are done at constant ω, thus we drop ω from the list of function arguments when not explicitly needed. The Porter-Bojarski hologram is defined by Porter & Devaney (1982):

Φ α (x, A) := A [ψ(r)∂ n H α (r, x) -H α (r, x)∂ n ψ(r)]dr, (9)
where H α is an acoustic wave propagator for the reference medium and A is a surface on the Sun where ψ and ∂ n ψ are observed. The role of H α is to propagate the wave field backward (or forward) in time, which leads to the concept of egression (or ingression) in LB holography (Lindsey & Braun 2000a).

Several choices have been proposed in the literature for the propagators, as detailed in Table 1. These depend on the outgoing (G 0 ) and incoming (G - 0 ) Green's functions defined in a reference medium with ρ 0 , c 0 , γ 0 , and u 0 = 0:

L 0 [G 0 (r, r )] = δ(r -r ) and ∂ n G 0 = ik n G 0 on ∂V, (10) L 0 [G - 0 (r, r )] = δ(r -r ) and ∂ n G - 0 = -ik n G - 0 on ∂V, (11)
with

L 0 = -(∆ + k 2 0 ), (12)
where k 0 is k in the reference medium and k n is from Eq. ( 8). The Green's functions G * 0 or G - 0 are backward propagators (c.f. egression), while (G - 0 ) * is a forward propagator (c.f. ingression). When the surface A is closed, it is equivalent to use G * 0 and Im G 0 (Devaney & Porter 1985). Tsang et al. (1987) proposed H α = G - 0 as a backward propagator to correct for wave attenuation.

If the observations are made at the computational boundary and the wave field satisfies the same boundary condition as the Green's function, then Eq. ( 9) reads

Φ α (x, A) = A ψ(r)[∂ n H α (r, x) -ik n H α (r, x)]dr. (13)
When H α = G * 0 , we have

Φ(x, A) = -2iRe[k n ] A ψ(r)G * 0 (r, x)dr, (14)
which corresponds to the egression as defined by Lindsey & Braun (2000a). Thus LB and PB holograms are closely related, at least for the upper boundary condition employed here.

In LB holography, information is extracted from the egression-ingression correlation (wave-speed perturbations and flows) and from the egression power (source covariance). Analogously, we define the PB hologram intensity (or hologram covariance) as

I αβ (x, A, A ) = Φ * α (x, A)Φ β (x, A ). ( 15
)
For the case α = β we define

I α (x, A) = |Φ * α (x, A)| 2 . ( 16
)
Table 2. Proposed propagators and associated pupils (H α , A) and (H β , A ) for different types of perturbations. The pupil geometries are shown in Fig. 1.

Perturbations H α A H β A Source covariance Im G 0 P Im G 0 P Sound speed (G - 0 ) * P G - 0 P Flow u θ (G - 0 ) * Q south G - 0 Q north Flow u φ (G - 0 ) * Q east G - 0 Q west Fig. 1.
Pupil geometries used to compute sound-speed or source kernels (P) and flow kernels (Qs), also see Table 2.

Different choices of pupils and propagators will provide sensitivity to different quantities as shown in Table 2. Scatterers are detected by correlating forward and backward propagated holograms. Different pupil shapes will give access to different components of the flow (Table 2 and Fig. 1).

First-order perturbations with respect to a reference solar model

We wish to study how perturbations to the background medium affect holograms. Using the Born approximation, the first step is to express the perturbations to the wavefield and use this expression in Eq. ( 9) to obtain the perturbations to the hologram and hologram intensity.

Perturbations to the wavefield

We consider perturbations δc, δρ, δγ, u with respect to the reference medium defined by ρ 0 , c 0 , γ 0 , and u 0 = 0. The perturbations to the sources of excitation are described through the perturbations to the source covariance,

M(r, ω) = M 0 (r, ω) + δM(r, ω), (17)
where

δM(r, ω) = E[S * 0 (r)δS (r ) + δS * (r)S 0 (r )] = (r)M 0 (r, ω). ( 18
)
Using the Born approximation up to first order, we write the wave field as

ψ(r, ω) = ψ 0 (r, ω) + δψ(r, ω), (19)
where the zeroth-and first-order wave fields are given by

L 0 [ψ 0 ] = S 0 , (20) L 0 [δψ] = -δL[ψ 0 ] + δS . ( 21
)
The perturbed wave operator is

δL[ψ 0 ] = -δk 2 ψ 0 - 2iω ρ 1/2 0 c 0 ρ 0 u • ∇        ψ 0 ρ 1/2 0 c 0        , (22)
with

c 2 0 δk 2 = 2iω δγ -(ω 2 + 2iωγ 0 ) δc 2 c 2 0 - ∂ω 2 c ∂ρ 0 δρ. ( 23
)
In terms of the Green's function G 0 , we have

ψ(r) = ψ 0 (r) + V G 0 (r, r s )δk 2 (r s )ψ 0 (r s )dr s + V G 0 (r, s)δS (s)ds + 2iω V ρ 0 (r s ) u 0 (r s ) • ∇        G 0 (r, r s ) ρ 1/2 0 (r s )c 0 (r s )        ψ 0 (r s ) ρ 1/2 0 (r s )c 0 (r s ) dr s . ( 24
)
where the scattering location r s spans the entire volume V.

Perturbations to the hologram

For convenience, we introduce the source kernels

K α (x, s, A) = A [G(r, s)∂ n H α (r, x) -H α (r, x)∂ n G(r, s)]dr , ( 25
)
such that the PB hologram is given by

Φ α (x, A) = V K α (x, s, A)S (s)ds. ( 26
)
We denote by K α,0 the source kernel in the reference medium (when G is replaced by G 0 ).

Replacing the wavefield by its first order expansion (Eq. ( 24)) in the definition of the hologram (Eq. ( 9)), one obtains the perturbation to the hologram

Φ α (x, A) = Φ α,0 (x, A) + δΦ α (x, A), (27)
where Φ α,0 is the hologram in the background medium given by Eq. ( 26) when K α is replaced by K α,0 and δΦ α expresses the changes in the hologram due to the perturbations from the background medium

δΦ α (x, A) = V δK α (x, s, A)S (s)ds + V K α,0 (x, s, A)δS (s)ds, (28)
where

δK α (x, s, A) = V K α,0 (x, r s , A)δk 2 (r s )G 0 (r s , s)dr s + 2iω V K α,0 (x, r s , A) ρ 1/2 0 (r s )c 0 (r s ) ρ 0 u 0 • ∇        G 0 (r s , s) ρ 1/2 0 (r s )c 0 (r s )        dr s . (29)

Perturbations to hologram intensity

To first order, we write the hologram intensity in the form

I αβ (x, A, A ) = I αβ,0 (x, A, A ) + δI αβ (x, A, A ). ( 30
)
The expectation values of the zeroth-and first-order hologram intensities are

E[I αβ,0 (x, A, A )] = V K * α,0 (x, s, A)K β,0 (x, s, A )M 0 (s)ds, (31) E[δI αβ (x, A, A )] = V K * α,0 (x, s, A)δK β (x, s, A )M 0 (s)ds + V δK * α (x, s, A)K β,0 (x, s, A )M 0 (s)ds + V K * α,0 (x, s, A)K β,0 (x, s, A ) (s)M 0 (s)ds. ( 32
)
Using the definition of the source kernels, we obtain

E[I αβ,0 (x, A, A )] = C 0 (r, r ) αβ (x, A, A ), (33)
where

C 0 (r, r ) = V G * 0 (r, s)G 0 (r , s)M 0 (s)ds (34)
is the expectation value of the cross-covariance function and, for any function F(r, r ), the double brackets mean

F(r, r ) αβ (x, A, A ) = A dr A dr ∂ n H * α (r, x)F(r, r )∂ n H β (r , x) -H * α (r, x)∂ n F(r, r )∂ n H β (r , x) -∂ n H * α (r, x)∂ n F(r, r )H β (r , x) + H * α (r, x) ∂ n ∂ n F(r, r ) H β (r , x) . ( 35
)
The perturbation to the hologram intensity is given by

E[δI αβ (x, A, A )] = V (s)K αβ (x, s, A, A ) ds + V δk 2 * (r s ) K k αβ (x, r s , A, A ) + δk 2 (r s ) K k * βα (x, r s , A , A) dr s + 2iω V u(r s ) • K u αβ (x, r s , A, A ) -K u * βα (x, r s , A , A) dr s , (36)
where

K αβ (x, s, A, A ) = M 0 (s) G * 0 (r, s)G 0 (r , s) αβ (x, A, A ), (37)
K k αβ (x, r s , A, A ) = C 0 (r s , r )G * 0 (r, r s ) αβ (x, A, A ), (38)
K u αβ (x, r s , A, A ) = ∇ C 0 (r s , r ) ρ 0 (r s ) 1/2 c 0 (r s ) ρ 0 (r s ) 1/2 G * 0 (r, r s ) c 0 (r s ) αβ (x, A, A ). (39)
The kernels for δk 2 and δk 2 * can be combined using Eq. ( 23) to obtain kernels for sound-speed K c αβ , density K ρ αβ and attenuation K γ αβ . For example,

E[δI αβ (x, A, A )] = V δc(r s ) K c αβ (x, r s , A, A ) dr s (40) with K c αβ (x, r s , A, A ) = - 2(ω 2 -2iωγ) c 3 0 (r s ) K k αβ (x, r s , A, A ) - 2(ω 2 + 2iωγ) c 3 0 (r s ) K k * βα (x, r s , A , A). (41)

Choice of the source covariance

In order to be able to compute the above kernels, one still need to choose the source covariance function M 0 in order to define the reference cross-covariance C 0 using Eq. ( 34). One possibility is to place the sources at a single depth, a few hundred kilometers below the solar surface.

Another possibility is to choose a source covariance of the form

M 0 (r, ω) = Π(ω) γ(r, ω) c 2 0 (r) , (42)
where Π(ω) is the source power spectrum (see Gizon et al. 2017). This choice implies

C 0 (r, r , ω) = Π(ω) 2ω Im G 0 (r , r, ω) + surface term. ( 43
)
The surface term depends on the boundary condition. It vanishes for a Dirichlet boundary condition (free surface), while it remains for radiative boundary conditions (e.g. Sommerfeld). Below the acoustic cutoff frequency, modes are trapped well below the observational and computational boundaries and the surface term is negligible. In this paper we use Eq. ( 42) in the convection zone and switch off the sources above the photosphere. By doing so, the surface term in Eq. ( 43) vanishes.

Noise in holograms

To compute the noise level, we compute the variance of the hologram intensity in the absence of scatterers:

Var[I αβ,0 (x)] = Var K * α,0 (x, s)K β,0 (x, s )S * (s)S (s )dsds = V 4 K * α,0 (x, s 1 )K β,0 (x, s 1 )K α,0 (x, s 2 )K * β,0 (x, s 2 ) × M 4 (s 1 , s 1 , s 2 , s 2 ) ds 1 ds 1 ds 2 ds 2 - V 2 K * α,0 (x, s)K β,0 (x, s )E S * (s)S (s ) dsds 2 , (44)
where

M 4 (s 1 , s 1 , s 2 , s 2 ) = E S * (s 1 )S (s 1 )S (s 2 )S * (s 2 ) . (45)
Under the (very reasonable) assumption that S is a realization drawn from a Gaussian random process, the fourth-order moment is the sum of products of the second-order moments:

M 4 (s 1 , s 1 , s 2 , s 2 ) = E S * (s 1 )S (s 1 ) E S (s 2 )S * (s 2 ) + E S * (s 1 )S (s 2 ) E S (s 1 )S * (s 2 ) + E S * (s 1 )S * (s 2 ) E S (s 2 )S (s 1 ) .

The first term in M 4 cancels out the squared term in Eq. ( 44). The third term is zero as the frequencies are uncorrelated. Thus,

Var[I αβ,0 (x)] = V 4 K * α,0 (x, s 1 )K β,0 (x, s 1 )K α,0 (x, s 2 )K * β,0 (x, s 2 ) × E S * (s 1 )S (s 2 ) E S (s 1 )S * (s 2 ) ds 1 ds 1 ds 2 ds 2 = V |K α,0 (x, s)| 2 M(s)ds V |K β,0 (x, s)| 2 M(s)ds = E[I α,0 (x)]E[I β,0 (x)]. ( 47
)
When α = β, the standard deviation of I α,0 is equal to its expectation value. This is because the probability density function of I α,0 is a χ 2 with two degrees of freedom, i.e. an exponential distribution.

Average over frequencies and signal-to-noise ratio

In order to increase the signal-to-noise ratio, one usually averages the hologram intensity over a range of frequencies [ω 0 -∆ω/2, ω 0 + ∆ω/2]. For an observation duration T , this interval will contain N = ∆ω T/2π independent frequencies. The frequency-averaged perturbation to the hologram intensity (i.e. the signal) is denoted by

δI αβ (x) = 1 N N i=1 δI αβ (x, ω i ). ( 48
)
The variance of the noise in the average hologram intensity is then given by

Var I αβ,0 (x) = Var        1 N N i=1 I αβ,0 (x, ω i )        (49) = 1 N 2 N i=1 Var I αβ,0 (x, ω i ) (50) = 1 N Var I αβ,0 (x) , (51)
since the noise in holograms at different frequencies is uncorrelated.

The expected signal-to-noise ratio is thus

SNR(x) = E δI αβ (x) Var I αβ,0 (x) = √ N E δI αβ (x) E[I α,0 (x)] E[I β,0 (x)] . (52)
The number of available frequencies N within a fixed frequency band ∆ω is proportional to the observation duration T , hence the noise level goes like T -1/2 . Provided that the frequency interval ∆ω is small with respect to the variations of the signal, then the signal-to-noise ratio will increase like T 1/2 .

Example computations

In order to illustrate the theory, we compute holograms in the presence of sound-speed perturbations at different depths and calculate the corresponding signal-to-noise ratios.

Reference Green's function

The main input quantity required to compute holograms is the reference Green's function (Eq. 10). Here it is computed in the frequency domain using the spherically-symmetric standard solar Model S [START_REF] Christensen-Dalsgaard | [END_REF]. The wave attenuation model is taken from Gizon et al. (2017). Below 5.3 mHz, we have γ = γ 0 |ω/ω 0 | 5.77 , where γ 0 /2π = 4.29 µHz and ω 0 /2π = 3 mHz. Above 5.3 mHz, γ/2π = 125 µHz is kept constant. The radiation boundary condition defined by Eq. ( 7) is applied at the computational boundary located 500 km above the solar surface with the local wavenumber k n (where H = 105 km). The wave equation is solved using the finite-element solver Montjoie (Duruflé 2006;Gizon et al. 2017).

The reference Green's function only depends on the angular distance Θ between the two points at radii r and r . To speed up the computations, we place one of the points on the polar axis and compute the axisymmetric component of the Green's function G m=0 l (r, r , ω) at each spherical harmonic degree l, to obtain:

G 0 (r, r , ω) l max l=0 G m=0 l (r, r , ω)P l (cos Θ), (53)
where we use an approximate equality because the sum is truncated at l max = 300.

Sound-speed kernels

The sound-speed kernel is computed using Eq. ( 41) and the definition of K k αβ . One needs to evaluate two surface integrals, which can be computed analytically via a decomposition of all quantities into spherical harmonic coefficients (Fournier et al. 2018).

Figure 2 shows a sound-speed kernel K c αβ at a single frequency of 3 mHz. The pupil P is a polar cap of angular size 120 • and the wave propagators H α and H β are given in Table 2. As expected the kernel peaks around the scatterer position at z = 0.7R . A cut along the polar axis is shown in Fig. 3; the kernel width is about half the local wavelength. In addition, we find ghost values at the antipode.

Signal

At position r s = z s ẑ along the polar axis, we consider a localized increase in sound speed of 10% over a volume V s , such that the signal (Eq. 40) may be written as

E[δI αβ (x)] 0.1V s c 0 (r s ) K c αβ (x, r s ). ( 54
)
The volume V s is that of a ball of diameter λ(r s )/2 = π/[Re k(r s , ω 0 )] with ω 0 /2π = 3 mHz. This is an approxi- ] at a single frequency of 3 mHz, displayed along the z-axis (at x = zẑ). The sound-speed perturbation (see Eq. 54) is placed at two different positions along the z-axis, z s = 0.7 R (red) and 0.9 R (blue). The standard deviation of the noise Var I αβ,0 (x) is given by the black curve. Note that the jagged aspect of the curves is not due to numerical inaccuracies. Right panel: hologram intensity and noise level after averaging over 101 frequencies uniformly distributed in the interval from 2.75 to 3.25 mHz. The frequency resolution is 5 µHz, implying an observation duration of T = 55.5 h. A horizontal line segment is plotted at each depth to mark half of the local wavelength.

mate but much faster way to compute the effect of a perturbation of volume comparable to the highest possible holographic resolution. It has been checked that the answer does not differ significantly from the one obtained by integrating numerically the kernel over the ball of volume V s . For reference, note that λ/2 = 38 Mm at r = 0.7 R and λ/2 = 20 Mm at r = 0.9 R .

Figure 4 shows the signal E [δI αβ (x)] for sound-speed perturbations located at two different depths z s = 0.7 R (red curve) and 0.9 R (blue curve). The pupil P and the wave propagators are the same as those of Fig. 2. The left panel of Fig. 4 shows the results at a single frequency of 3 mHz. With only one frequency, the signal peaks close to the scattering location but the spatial resolution is relatively poor, with a ghost on the far side. To demonstrate the benefits of averaging, the right panel shows the signal after averaging over 101 frequencies uniformly distributed in the interval 2.75 -3.25 mHz. The frequency resolution 5 µHz corresponds to an observation duration T = 55.5 h. Averaging over frequencies improves the spatial resolution which approaches λ/2 and the ghost is suppressed. A horizontal line segment is plotted on the right panel at each depth to mark half the local wavelength.

As seen in Fig. 5 the spatial extent of the frequency-averaged kernels is approximately λ/2 in both the radial and horizontal directions, for all scattering points in the range 0.6 < z s /R < 0.98. Thus helioseismic holography is a diffractionlimited imaging technique as suggested by Lindsey & Braun (1997).

Noise

The noise is obtained from Eq. ( 47), which requires the computation of E[I α,0 ] and E[I β,0 ] using Eq. ( 33). The reference crosscovariance C 0 is precomputed. The double surface integral is evaluated in a similar way as for the kernel computations.

For a frequency of 3 mHz the left panel of Fig. 4 (black curve) shows the noise level, together with the signal described in the previous section. The jagged aspect of the noise varia- tions with position is not due to numerical inaccuracies but to the details of the Green's function. As shown in the right panel of Fig. 4, the noise level goes down by a factor of about ten after averaging over 101 frequencies, and varies more smoothly with depth. Braun & Birch (2008a) studied the noise level in observed travel times measured from LB holography. These measurements, however, include contributions from supergranulation and so are not directly comparable to what is shown in Fig. 4.

Signal-to-noise ratio

Figure 6 shows the signal-to-noise ratio as a function of scatterer location. We recall that the sound-speed perturbation is specified by Eq. ( 54) and is the same as in Sect. frequencies in the interval 2.75 -3.25 mHz. After averaging, the signal-to-noise ratio is above 2 and is roughly independent of depth in the range 0.6 < z s /R < 0.98 for a pupil of angular size 120 • . Note that the ghost at -z s is much below the noise level. We find that both signal and noise vary rapidly with frequency for deep located sound-speed scatterers. Figure 7 shows an example of a sound-speed scatterer located at z s = 0.7 R . Strong frequency variations in signal and noise are evident for frequencies below 3.5 mHz. This can be understood as follows. Low-frequency modes have narrowly-peaked power spectra due to their long lifetimes. At these low frequencies, the amplitude of the kernel function and the noise will change rapidly when the frequency coincides with a particular mode frequency. Additionally, the kernel function may not peak exactly at the sound-speed scatterer position when only a few modes contribute to the kernel function. Figure 8 shows the signal-to-noise ratio as a function of frequency for a sound-speed scatterer located at z s = 0.9 R . For this target depth closer to the surface, the rapid variations disappear above 3 mHz, due to the larger contribution of highdegree modes which are not resolved in frequency space because of their short lifetimes.

The kernel function at frequency 2.4000 mHz for z s = 0.7R is shown in Fig. 9; this particular frequency corresponds to the peak marked in Fig. 7 with a red dot. We see that this kernel is much less localized around the scattering point than the kernel at 3 mHz (Fig. 2).

Conclusion

We derived a framework for computing the expected signal and the noise level in PB helioseismic holography. The same framework could be used to interpret LB holograms and phasesensitive holograms.

PB holography requires knowledge of the wave field, ψ = ρ 1/2 c 2 ∇ • ξ, and its normal derivative, ∂ n ψ, on the solar surface. With this definition of ψ, the Green's function used in the definition of PB holograms solves a well-defined Helmholtz-like equation, which we solve numerically (Gizon et al. 2017). The need for finite-frequency Green's functions was demonstrated in LB 8. Signal-to-noise ratio as in Fig. 7, but for a sound-speed scatterer located closer to the surface at z s = 0.9 R . holography by Pérez Hernández & González Hernández (2010). In the numerical examples shown in the previous section, we assumed that we have full knowledge of ∂ n ψ on the surface. In practice, we do not observe directly the normal derivative of the wave field; it must be approximated. According to complementary simulations (not shown here), this can be achieved by using αβ with z s = 0.7 R computed at the low frequency of 2.4000 mHz, which corresponds to the spike with a red dot in Fig. 7. This kernel displays oscillations and is not peaked as much as the 3-mHz kernel from Fig. 2. the approximate outgoing radiation condition ∂ n ψ = ik n ψ derived by Barucq et al. (2018).

We found that, for a sufficiently large pupil, scatterers can be imaged at a resolution that is very close to half the local wavelength, λ/2. This confirms the claim by Lindsey & Braun (1997, 2000a) that helioseismic holography is diffraction-limited. In that sense, helioseismic holography is superior to deep-focusing time-distance helioseismology, which gives lower spatial resolution (Munk 2001;Pourabdian et al. 2018).

For large pupils, we found that the signal-to-noise ratio in holograms does not vary much with depth in the convection zone, when a perturbation in sound-speed fills a volume corresponding to the holographic resolution.

Averaging over frequencies improves the signal-to-noise ratio. For a scatterer at the bottom of the convection zone, the signal and the noise vary smoothly with frequency above 4 mHz (see Fig. 7). At lower frequencies, however, the signal varies rapidly with frequency (due to contributions from individual long-lived p modes) and it is not obvious how the signal should be averaged. A specific analysis of low-frequency holograms is required, especially for deep scatterers.

We found that the signal-to-noise ratio in PB holography is maximum around 3.7 mHz for z s = 0.7 R (resp. at 4.3 mHz for z s = 0.9 R ). There is no indication in our calculations that there is a benefit in using the frequencies above the acoustic cutoff (unlike predictions by Ruzmaikin & Lindsey 2003, for phase sensitive holography). The signal-to-noise ratio drops to very small values above 5 mHz. One may ask if this drop is somehow compensated by the increase in spatial resolution at high frequencies. The answer is negative. Our calculations indicate that noise has a horizontal correlation length that is about half the local wavelength. Far too few independent measurements are available at high frequencies to recover a decent signal-to-noise ratio by horizontal spatial averaging.

Our synthetic data do not contain a convective background. The effect of this background on signal-to-noise ratios in holography should be studied. Future work should also investigate the performance of PB holography for target locations that are away from the axis of the pupil, especially for farside imaging applications.

Fig. 2 .

 2 Fig. 2. Meridional slice through the sound-speed kernel K c αβ (x, r s ) computed in Model S at a single frequency of 3 mHz, in units of 10 -34 kg m -3 s -3 . Both the real (left panel) and imaginary (right panel) parts of the kernel are shown. The scatterer at z s = 0.7 R is indicated by the crosses. The observation pupil P is a polar cap of full angular size 120 • . Notice the ghost values at the antipode result from the reflection of waves at the surface due to the rapid drop of the density.

Fig. 3 .

 3 Fig. 3. Cut along the z axis through the sound-speed kernel from Fig. 2. The scatterer is at z s = 0.7 R .

Fig. 4 .

 4 Fig. 4. Left panel: Hologram intensity E[δI αβ(x)] at a single frequency of 3 mHz, displayed along the z-axis (at x = zẑ). The sound-speed perturbation (see Eq. 54) is placed at two different positions along the z-axis, z s = 0.7 R (red) and 0.9 R (blue). The standard deviation of the noise Var I αβ,0 (x) is given by the black curve. Note that the jagged aspect of the curves is not due to numerical inaccuracies. Right panel: hologram intensity and noise level after averaging over 101 frequencies uniformly distributed in the interval from 2.75 to 3.25 mHz. The frequency resolution is 5 µHz, implying an observation duration of T = 55.5 h. A horizontal line segment is plotted at each depth to mark half of the local wavelength.

Fig. 5 .

 5 Fig. 5. Radial and horizontal widths of the frequency-averaged soundspeed kernel | K c | as functions of scattering position. These are close to half the local wavelength at 3 mHz (dotted line).

Fig. 6 .

 6 Figure6shows the signal-to-noise ratio as a function of scatterer location. We recall that the sound-speed perturbation is specified by Eq. (54) and is the same as in Sect. 7.3. The results are shown at a single frequency of 3 mHz and after averaging over 101

Fig. 7 .

 7 Fig. 7. Signal, noise, and signal-to-noise ratio as function of frequency for a sound-speed scatterer located at z s = 0.7 R . Here we show the result for a frequency range of 2 to 7 mHz. The rapid changes are not due to numerical inaccuracies. The red dots mark the spikes in the signal and the noise at 2.4000 mHz.

Fig.

  Fig.8. Signal-to-noise ratio as in Fig.7, but for a sound-speed scatterer located closer to the surface at z s = 0.9 R .

Fig. 9 .

 9 Fig. 9. Meridional slice of the sound-speed kernel K cαβ with z s = 0.7 R computed at the low frequency of 2.4000 mHz, which corresponds to the spike with a red dot in Fig.7. This kernel displays oscillations and is not peaked as much as the 3-mHz kernel from Fig.2.

Table 1 .

 1 Possible wave propagators.

	Wave propagator	H α	References
	Backward Backward Backward Forward	G * 0 G -0 Im G 0 Devaney & Porter (1985) Porter & Devaney (1982) Tsang et al. (1987) (G -0 ) * This work

Acknowledgements. The theoretical framework was developed by L.G. and D.F. at Mathematisches Forchungsinstitut Oberwolfach in May 2017. The numerical computations were performed by D.Y. using the Montjoie solver. D.Y. is a member of the International Max Planck Research School for Solar System Science at the University of Göttingen. We thank M. Duruflé and J. Chabassier from Inria Magique-3D for the helioseismology-related developments of Montjoie. We also thank Chris Hanson from NYUAD for the fine tuning of the model power spectrum of solar oscillations. L.G. acknowledges support from NYUAD Institute grant G1502. The computing resources were provided in part by the German Data Center for SDO, a project funded by the German Aerospace Center (DLR).