Ugo Nanni1

Florent Gimbert1

Philippe Roux2

Albanne Lecointre2

Ugo Nanni
email: ugo.nanni@univ-grenoble-alpes.fr

Observing the subglacial hydrology network and its dynamics with a dense seismic array

Keywords: Physical Sciences -Earth, Atmospheric, and Planetary Sciences Cryoseismology -Subglacial hydrology -dense seismic array -seismic noise

Subglacial water flow strongly modulates glacier basal motion, which itself strongly influences the contributions of glaciers and ice sheets to sea-level rise. However, our understanding of when and where subglacial water flow enhances or impedes glacier flow is limited due to the paucity of direct observations of subglacial drainage characteristics. Here we demonstrate that dense seismic array observations combined with an innovative systematic seismic source location technique allows for the first time the retrieval of a two-dimensional map of a subglacial drainage system, as well as its day-to-day temporal evolution. We observe with unprecedented detail when and where subglacial water flows through a cavity-like system that enhances glacier flow versus when and where water mainly flows through a channel-like system that impedes glacier flow. Most importantly, we are able to identify regions of high hydraulic connectivity within and across the cavity and channel systems, which have been identified as having a major impact on the long term glacier response to climate warming. Applying a similar seismic monitoring strategy in other glacier settings, including for ice sheets, may help to diagnose the susceptibility of their dynamics to increased meltwater input due to climate warming.

Significance Statement

Our understanding of when, where and under which conditions subglacial water flow favors or impedes glacier flow remains uncertain mainly because of sparse field observations. This strongly limits our capability to assess the susceptibility of glaciers and ice-sheets to a future increase in meltwater input due to climate warming. Here we overcome classic observational difficulties by establishing an innovative seismic-based approach. From dense-seismic array observations we retrieve with unprecedented detail the two-dimensional map of a subglacial hydrology network.

We observe its day-to-day evolution from a cavity-dominated system that favors glacier flow to a channel-dominated system that impedes glacier flow. Although our method was applied on an Alpine Glacier, it is easily adaptable and readily applicable to other settings such as Greenland or Antarctica.

Main Text

Introduction

Meltwater produced at the surface of a glacier is mostly routed to its bed through crevasses and moulins and flows towards its terminus through the subglacial drainage system 1 . In this system, water pressure regulates the ice-bed mechanical coupling, which determines the glacier sliding speed, and therefore has major effects on the stability of glaciers and ice sheets [START_REF] Iken | The relationship between subglacial water pressure and velocity of Findelengletscher, Switzerland, during its advance and retreat[END_REF][START_REF] Davison | The influence of hydrology on the dynamics of land-terminating sectors of the Greenland ice sheet[END_REF] and their contributions to sea-level rise [START_REF] Sole | Seasonal speedup of a Greenland marine-terminating outlet glacier forced by surface melt-induced changes in subglacial hydrology[END_REF][START_REF] Ritz | Potential sea-level rise from Antarctic ice-sheet instability constrained by observations[END_REF][START_REF] King | Dynamic ice loss from the Greenland Ice Sheet driven by sustained glacier retreat[END_REF] . Subglacial water pressure shows complex dependency on the subglacial route that the water follows, and in particular on whether the water flows through a distributed and inefficient system [START_REF] Lliboutry | General Theory of Subglacial Cavitation and Sliding of Temperate Glaciers[END_REF] (e.g. cavity-dominated for hard bed glaciers) or a localized and efficient system 8 (channel-dominateddrainage system is expected to be associated with a, comparatively, lower water pressure, which promotes low glacier sliding speeds [START_REF] Röthlisberger | Water Pressure in Intra-and Subglacial Channels[END_REF][START_REF] Hubbard | Borehole water-level variations and the structure of the subglacial hydrological systelll of Haut Glacier d ' Arolla[END_REF][START_REF] Schoof | Ice-sheet acceleration driven by melt supply variability[END_REF][START_REF] Tedstone | Decadal slowdown of a land-terminating sector of the Greenland Ice Sheet despite warming[END_REF] . Where cavities are hydraulically well connected to channels they drain into the efficient drainage system, which tends to lower the overall water pressure [START_REF] Andrews | Direct observations of evolving subglacial drainage beneath the Greenland Ice Sheet[END_REF] . In addition, recent observations [START_REF] Andrews | Direct observations of evolving subglacial drainage beneath the Greenland Ice Sheet[END_REF] suggested that hydraulically isolated areas of the bed with very low permeability [START_REF] Hoffman | Greenland subglacial drainage evolution regulated by weakly connected regions of the bed[END_REF] regulate glacier basal traction during winter [START_REF] Sole | Winter motion mediates dynamic response of the Greenland Ice Sheet to warmer summers[END_REF] and over multi-annual timescales [START_REF] Williams | Dynamic response of the Greenland ice sheet to recent cooling[END_REF] . The spatial persistence of high water pressure at the glacier bed [START_REF] Rada | Channelized, distributed, and disconnected: Subglacial drainage under a valley glacier in the Yukon[END_REF] thus depends on the subglacial drainage system configuration and the hydraulic connectivity across the cavities and from the cavities to the channels. However, current observations of subglacial drainage systems are rather point-scale [START_REF] Rada | Channelized, distributed, and disconnected: Subglacial drainage under a valley glacier in the Yukon[END_REF][START_REF] Nienow | Seasonal changes in the morphology of the subglacial drainage system, Haut Glacier d'Arolla, Switzerland[END_REF] (e.g., via icedrilling) or spatially-integrated [START_REF] Tranter | Hydrochemistry as an indicator of subglacial drainage system structure: a comparison of Alpine and Sub-Polar environments[END_REF][START_REF] Chandler | Evolution of the subglacial drainage system beneath the Greenland Ice Sheet revealed by tracers[END_REF] (e.g. via dye tracing experiments or hydrochemical analysis), such that they provide only partial representations of the heterogeneous nature of the subglacial drainage system. It thus remains uncertain where and when isolated cavities, connected cavities, and/or channels operate, and therefore under what conditions (e.g., water supply rate, glacier geometry) meltwater supply to the glacial bed enhances or limits glacier flow.

Recent studies have indicated that turbulent subglacial water flow generates detectable seismic noise [START_REF] Bartholomaus | Subglacial discharge at tidewater glaciers revealed by seismic tremor[END_REF] given flow velocities of the order of a meter per second [START_REF] Gimbert | Subseasonal changes observed in subglacial channel pressure, size, and sediment transport[END_REF] , that can be used to retrieve the physical properties of subglacial channels (e.g., water pressure, conduit size) [START_REF] Gimbert | Subseasonal changes observed in subglacial channel pressure, size, and sediment transport[END_REF][START_REF] Nanni | Quantification of seasonal and diurnal dynamics of subglacial channels using seismic observations on an Alpine glacier[END_REF] as well as to estimate their location in space [START_REF] Vore | Seismic Tremor Reveals Spatial Organization and Temporal Changes of Subglacial Water System[END_REF][START_REF] Lindner | Glaciohydraulic seismic tremors on an Alpine glacier[END_REF] . Previous studies that have attempted to spatially track water flow have documented their azimuthal distribution, but not their epicentral coordinates [START_REF] Vore | Seismic Tremor Reveals Spatial Organization and Temporal Changes of Subglacial Water System[END_REF][START_REF] Lindner | Glaciohydraulic seismic tremors on an Alpine glacier[END_REF][START_REF] Burtin | Location of river-induced seismic signal from noise correlation functions[END_REF] .

Location of epicentral coordinates of simultaneously active and spatially spread noise sources indeed represents a major seismological challenge [START_REF] Rost | Array seismology: Methods and applications[END_REF] , since it requires observations of the wavefield all around the targeted sources with high enough resolution to deal with its strong incoherency [START_REF] Kuperman | Matched field acoustics[END_REF][START_REF] Almendros | Array analysis using circularwave-front geometry: An application to locate the nearby seismo-volcanic source[END_REF] . Here, we demonstrate that well-resolved maps of subglacial drainage systems can be retrieved using particularly dense seismic arrays and by adapting the matched-field processing (MFP) technique [START_REF] Kuperman | Matched field acoustics[END_REF] to the particularity of dealing with spatially spread sources.

Seismic monitoring and source location strategy

We use the 1-month seismic record from 98 three-component sensors that was acquired at the surface of glacier d'Argentière (French Alps; Fig. 1) during the onset of the 2018 melt-season (see Methods and ref. [START_REF] Gimbert | A Multi-Physics Experiment with a Temporary Dense Seismic Array on the Argentière Glacier, French Alps: The RESOLVE Project[END_REF] for detailed description of the experiment and ref. [START_REF] Nanni | Quantification of seasonal and diurnal dynamics of subglacial channels using seismic observations on an Alpine glacier[END_REF] for a broader glaciological context). At this location, subglacial water flow strongly influences glacier dynamics [START_REF] Vincent | Sliding velocity fluctuations and subglacial hydrology over the last two decades on Argentière glacier, Mont Blanc area[END_REF] and generates continuous seismic noise that is most pronounced in the [3-7] Hz frequency range [START_REF] Nanni | Quantification of seasonal and diurnal dynamics of subglacial channels using seismic observations on an Alpine glacier[END_REF] , in which seismic wavelengths are on the order of [200-500] m 32 . Our seismic array covers an area of 400 x 600 m 2 with a 40-m to 50-m sensor spacing and an aperture (i.e. largest distance between stations; Fig. 1) of ca. 720m. Such configuration allows subwavelength sampling while covering all azimuthal directions around the targeted sources, which are two necessary conditions for retrieving the epicentral coordinates of seismic sources through application of MFP [START_REF] Kuperman | Matched field acoustics[END_REF][START_REF] Almendros | Array analysis using circularwave-front geometry: An application to locate the nearby seismo-volcanic source[END_REF] (Methods). MFP consists of recursively matching the predicted and observed phase delays (Methods), and has been extensively applied to locate spatially wellseparated sources, such as those generated by hydrothermal activity [START_REF] Legaz | Self-potential and passive seismic monitoring of hydrothermal activity: A case study at Iodine Pool, Waimangu geothermal valley, New Zealand[END_REF] , oil and gas injection [START_REF] Chmiel | High-sensitivity microseismic monitoring: Automatic detection and localization of subsurface noise sources using matched-field processing and dense patch arrays[END_REF] , icequakes [START_REF] Sergeant | On the Green's function emergence from interferometry of seismic wave fields generated in high-melt glaciers: Implications for passive imaging and monitoring[END_REF] , or englacial moulins [START_REF] Sergeant | On the Green's function emergence from interferometry of seismic wave fields generated in high-melt glaciers: Implications for passive imaging and monitoring[END_REF] . However, it has been little applied to locate spatially distributed noise sources [START_REF] Chmiel | High-sensitivity microseismic monitoring: Automatic detection and localization of subsurface noise sources using matched-field processing and dense patch arrays[END_REF][START_REF] Venkatesh | Beamforming algorithm for distributed source localization and its application to jet noise[END_REF] . We adapt this technique to provide such particularity through a systematic analysis of the phase coherence over 1-s-long time windows. Within each window, we apply an efficient gradient-based minimization algorithm that allows us to locate up to 29 simultaneously active sources (Methods).We obtain each source location from maximizing the correlation between the observed and the modeled phase delays. We refer to this maximum correlation as the MFP output, which ranges from 0 to 1. An MFP output close to 1 indicates global phase coherence, as expected for a dominant and punctual source [START_REF] Gimbert | Seismological Research Letters The RESOLVE project : a multi-physics experiment with a temporary dense seismic array on the Argentière Glacier[END_REF] , while a low MFP output indicates several local phase coherences, as expected for multiple sources simultaneously acting over the area. We note that we only keep the sources that are located within 400 m of the center of our array and are associated with realistic physical properties (i.e., wave velocities, source depth; Fig. 2b,c). This yields more than a million seismic sources per day (Fig. 2a).

Retrieving the geometry of the subglacial system

In Figure 2d, e, we show the normalized spatial distribution probability of the source location obtained over the study period for high and low MFP output ranges, as indicated in Figure 2a. We observe that events associated with high (i.e., global) phase coherence (MFP output, >0.8) are mainly located at the glacier surface (Fig. 2b) where crevasses are observed (Fig. 2e), and are associated with phase velocities that are typical of surface waves [START_REF] Sergeant | On the Green ' s function emergence from interferometry of seismic wavefields generated in high-melt glaciers : implications for passive imaging and monitoring[END_REF] (ca. 1580 m.s -1 ; Fig. 2c). This is consistent with these events corresponding to crevasse-induced icequakes (Fig. 1). On the contrary, sources associated with low (i.e., local) phase coherence (MFP output, [0.05-0.3]) are preferentially located at depth near the ice-bed interface (Fig. 2b). The associated phase velocities (Fig. 2c) vary by up to ca. 3600 m.s -1 , which is consistent with body waves being generated additionally to surface waves, and thus with sources that occur at depth [START_REF] Aki | Quantitative seismology[END_REF] . The spatial distribution of these sources show two elongated regions located along the glacier centerline (one at maximum ice thickness and one 50 m further down-glacier) where hydraulic potential calculations allow materializing the likely location of subglacial channel(s) (Fig. 2d; see Methods for details on the hydraulic potential calculations). All of these concomitant features lead us to interpret these observed seismic sources as a depiction of the geometry of the subglacial drainage system. The characteristic width of the observed spots (ca. 50 m, Fig. 2d) is much greater than might be expected for a single subglacial channel (on the order of meter to few meters [START_REF] Röthlisberger | Water Pressure in Intra-and Subglacial Channels[END_REF] ), which could result from the location resolution being limited by the seismic wavelength investigated (expected to be 1/6 to 1/2 times the ca. 300m wavelength; see Methods) or from the presence of multiple channels. The absence of clear source locations in the up-glacier part of the array (Fig. 3) might be caused by reduced or less turbulent subglacial water flow there, compared to the down-glacier part, where over-deepening of the glacier bed might favor unstable and more turbulent water flow [START_REF] Lliboutry | Modifications to the theory of intraglacial waterways for the case of subglacial ones[END_REF] .

Observing the switch from a cavity-like to a channel-like drainage system

We further investigate whether with our observational technique we can retrieve the different components of a subglacial drainage system by evaluating temporal changes in the source location maps together with the subglacial flow parameters that characterize average water drainage efficiency and water pressure conditions at the glacier bed [START_REF] Gimbert | Subseasonal changes observed in subglacial channel pressure, size, and sediment transport[END_REF][START_REF] Nanni | Quantification of seasonal and diurnal dynamics of subglacial channels using seismic observations on an Alpine glacier[END_REF] (see Methods for details).

We calculate source location maps averaged over 2-day time windows (Fig. 3) and use the combined seismic signal amplitude and water discharge measurements to invert for spatiallyintegrated changes in the hydraulic pressure gradient and hydraulic radius of the subglacial drainage system (Fig. 4) [START_REF] Gimbert | Subseasonal changes observed in subglacial channel pressure, size, and sediment transport[END_REF][START_REF] Nanni | Quantification of seasonal and diurnal dynamics of subglacial channels using seismic observations on an Alpine glacier[END_REF] . All maps are obtained from sources yielding a similar phase coherence (i.e. associated with MFP output within the same narrow [0.07-0.16] range) such that they are retrieved with a similar accuracy (see Methods).

We observe that sources are spatially more distributed at the beginning of the period (until ca.

May 10; Fig. 3) than when averaging source locations over the whole period (Fig. 2d). At this time the hydraulic pressure gradient (Fig. 4a, green line) varies significantly with the discharge (from 0.1 m 3 .s -1 to 2 m 3 .s -1 ; Fig. 4b, blue line) and the hydraulic radius remains constant (Fig. 4a, purple line), which indicates low drainage efficiency [START_REF] Nanni | Quantification of seasonal and diurnal dynamics of subglacial channels using seismic observations on an Alpine glacier[END_REF] . This time period is also associated with marked glacier surface acceleration (velocity increases by up to ca. 50%) due to enhanced basal sliding [START_REF] Vincent | Sliding velocity fluctuations and subglacial hydrology over the last two decades on Argentière glacier, Mont Blanc area[END_REF] as a result of the subglacial water pressurization [START_REF] Lliboutry | General Theory of Subglacial Cavitation and Sliding of Temperate Glaciers[END_REF] . These concomitant observations provide strong support that with our source location analysis we can observe water flow in the cavity-like system. This has not been suggested in previous cryoseismic studies [START_REF] Bartholomaus | Subglacial discharge at tidewater glaciers revealed by seismic tremor[END_REF][START_REF] Gimbert | Subseasonal changes observed in subglacial channel pressure, size, and sediment transport[END_REF][START_REF] Nanni | Quantification of seasonal and diurnal dynamics of subglacial channels using seismic observations on an Alpine glacier[END_REF][START_REF] Lindner | Glaciohydraulic seismic tremors on an Alpine glacier[END_REF] on the premise that fast and turbulent water flow speeds (as m.s -1 ; i.e., seismically detectable [START_REF] Gimbert | Subseasonal changes observed in subglacial channel pressure, size, and sediment transport[END_REF] ) only occur in channels. However this is consistent with borehole field observations [START_REF] Rada | Channelized, distributed, and disconnected: Subglacial drainage under a valley glacier in the Yukon[END_REF][START_REF] Kavanaugh | Abrupt glacier motion and reorganization of basal shear stress following the establishment of a connected drainage system[END_REF] and with theoretical studies that have envisioned channel-like flow in connections between cavities (referred to as orifices) [START_REF] Kamb | Glacier surge mechanism based on linked cavity configuration of the basal water conduit system[END_REF] .

While progressing through time (from ca. May 10), the subglacial water flow localizes into a narrower zone (Fig. 3). This transition towards a channel-dominated system is well revealed by the significant increase in the two-dimensional coefficient of determination R 2 (Fig. 4a) calculated between the spatial pattern observed for May 25-26 (Fig. 3o) and each of the patterns shown in Figure 3. Concomitantly, and in contrast to the first part of the period, we observe a doubling of the hydraulic radius with a reduced and almost unvarying hydraulic pressure gradient, which indicates that the increasing drainage efficiency lowers the basal water pressure. Such an increase in drainage efficiency is consistent with the increase in water discharge not resulting in glacier acceleration, but rather in slight deceleration [START_REF] Nienow | Velocity-discharge relationships derived from dye tracer experiments in glacial meltwaters: implications for subglacial flow conditions[END_REF] . Our seismic analysis therefore provides independent observational support for the hypothesis that the development of an efficient and channel-like system reduces the basal water pressure and favors slower glacier flow.

Implications for the monitoring of hydraulic connectivity within and across subglacial drainage systems

The presence of seismically detectable turbulent water flow within the cavity-like drainage system suggests that this system is, at least in certain places, associated with high hydraulic connectivity.

The absence of seismic sources over many areas of the glacier bed indicates that these, comparatively, have lower hydraulic connectivity, and could correspond to hydraulically isolated portions of the bed with higher potential for water storage [START_REF] Andrews | Direct observations of evolving subglacial drainage beneath the Greenland Ice Sheet[END_REF][START_REF] Hoffman | Greenland subglacial drainage evolution regulated by weakly connected regions of the bed[END_REF] . This shows not only that our seismic analysis can determine when and where the transition from a cavity-like to a channel-like system occurs, but also that it can identify locations with high hydraulic connectivity within the cavity-like system and from the cavity-to the channel-like systems. This implies that we can evaluate which areas of the glacier bed are efficiently connected to channels, and thus drained when they develop, versus the areas where water is stored and is thought to regulate the basal traction over multi-year timescales [START_REF] Hoffman | Greenland subglacial drainage evolution regulated by weakly connected regions of the bed[END_REF][START_REF] Williams | Dynamic response of the Greenland ice sheet to recent cooling[END_REF] . Yielding such observations in other settings like Greenland has the potential for identifying locations where ice dynamics are expected to be particularly sensitive to the foreseen increase in melt water input rates due to climate warming [START_REF] Davison | The influence of hydrology on the dynamics of land-terminating sectors of the Greenland ice sheet[END_REF] .

Perspectives and limitations

A successful application of our methodology to other glacier systems like Ice-Sheet outlet glaciers relies on the capability to conduct subwavelength sampling over representative areas (Methods), ranging from about 1x1km 2 (e.g. Russel Glacier's tongue [START_REF] De Fleurian | A modeling study of the effect of runoff variability on the effective pressure beneath Russell Glacier, West Greenland[END_REF] ) up to about 10x10km 2 (e.g. central trunk of Pine Island Glacier [START_REF] Gillet-Chaulet | Assimilation of surface velocities acquired between 1996 and 2010 to constrain the form of the basal friction law under Pine Island Glacier[END_REF] ). These scaling constraints can be can be fulfilled through deploying c. 100 up to c. 10,000 sensors in these areas, which we foresee as achievable in the near future given the recent ease of dense seismic array installations in remote areas [START_REF] Booth | Distributed Acoustic Sensing (DAS) of Seismic Properties in a Borehole drilled on a Fast-Flowing Greenlandic Outlet Glacier[END_REF][START_REF] Soubestre | Solid Earth Network-Based Detection and Classification of Seismovolcanic Tremors : Example From the Klyuchevskoy Volcanic Group in Kamchatka[END_REF][START_REF] Moreau | Sea Ice Thickness and Elastic Properties From the Analysis of Multimodal Guided Wave Propagation Measured With a Passive Seismic Array[END_REF] .

Although the present seismic approach yields unprecedent insights on subglacial hydrology from the local (Fig. 3) up to the glacier (Fig. 4) scales, future applications will likely still strongly benefit from combination with complementary in-situ observations. Such a combination will guide the strategy for conducting local (e.g., via ice-drilling) or spatially-integrated (e.g. via dye tracing experiments or hydrochemical analysis) measurements in targeted and representative places identified from the seismic observations. Complementary measurements in those targeted places will also help enhance the level of quantitative interpretation of seismic observations through giving more quantitative constraints on the behavior of key physical variables (pressure, flow speed, connectivity) in each identified subglacial hydrology system. This integrated strategy will allow to better extrapolate physical constraints from the local up to the glacier-wide scale, which will ultimately facilitate the assimilation of seismic observations in subglacial hydrology -ice dynamics coupled models, which otherwise remain poorly constrained across scales [START_REF] Irarrazaval | Determining the evolution of an alpine glacier drainage system by solving inverse problems[END_REF] .

Conclusion

In this study, we provide well-resolved spatial observations of a localized efficient channel-like system geometry and of the extent of a distributed inefficient cavity-like system. We show that with an adapted systematic seismic investigation of low spatial phase coherences we can locate multiple seismic noise sources induced by subglacial water flow that act at the same time. This location procedure is feasible from very low water discharge (∼0.1 m 3 .s -1 ) to peak melt-season water discharge. Thus, we can simultaneously observe the distributed and localized drainage systems through time and space, and evaluate the hydraulic connectivity within and across these systems along with their changes through time. Our geophysical approach is also exportable to other glaciers in remote areas, from mountain glaciers to ice caps, especially with the current easing of dense-seismic arrays [START_REF] Moreau | Sea Ice Thickness and Elastic Properties From the Analysis of Multimodal Guided Wave Propagation Measured With a Passive Seismic Array[END_REF] and distributed acoustic sensing [START_REF] Booth | Distributed Acoustic Sensing (DAS) of Seismic Properties in a Borehole drilled on a Fast-Flowing Greenlandic Outlet Glacier[END_REF][START_REF] Walter | Distributed acoustic sensing of microseismic sources and wave propagation in glaciated terrain[END_REF] deployments. Our novel way of investigating subglacial drainage systems will allow the glaciological community to diagnose the susceptibility of ice sheets [START_REF] Hoffman | Greenland subglacial drainage evolution regulated by weakly connected regions of the bed[END_REF][START_REF] Williams | Dynamic response of the Greenland ice sheet to recent cooling[END_REF] and mountain glaciers [START_REF] Rada | Channelized, distributed, and disconnected: Subglacial drainage under a valley glacier in the Yukon[END_REF][START_REF] Irarrazaval | Determining the evolution of an alpine glacier drainage system by solving inverse problems[END_REF] to increased meltwater input due to climate warming and/or extreme melt and rainfall events. Our approach will also be particularly appropriate to study processes that generate similar spatially spread seismic noise in other environments, such as for lava flows on volcanoes [START_REF] Soubestre | Depth Migration of Seismovolcanic Tremor Sources Below the Klyuchevskoy Volcanic Group (Kamchatka) Determined From a Network-Based Analysis[END_REF] , tremors in fault zones [START_REF] Mordret | Shallow three-dimensional structure of the San Jacinto fault zone revealed from ambient noise imaging with a dense seismic array[END_REF] , and sediment transport in rivers [START_REF] Bakker | Field Application and Validation of a Seismic Bedload Transport Model[END_REF] .

Materials and Methods

Dense seismic array survey

All of the stations have a 500-Hz sampling rate and a low cut-off frequency of 4.5 Hz. The nodes were installed at a depth of 30 cm into the 4-m-thick snow cover, and had to be reinstalled on May 11 due to snow melt, which occurred at a rate of 5 cm.day -1 . We refer the reader to ref. [START_REF] Gimbert | A Multi-Physics Experiment with a Temporary Dense Seismic Array on the Argentière Glacier, French Alps: The RESOLVE Project[END_REF] for the detailed description of the experiment.

Complementary measurements

Concomitant with our seismic survey, we use continuous records of subglacial water discharge, which was measured in excavated subglacial tunnels maintained by the hydroelectric power company Emosson S.A. at ca. 600 m downstream of the center of the node-array (at 2173 m a.s.l.). We also installed four GNSS stations at the corners of the seismic array to measure glacial surface velocity, which was on the order of 0.1 m.day -1 at this time of the year and for this part of the glacier [START_REF] Vincent | Sliding velocity fluctuations and subglacial hydrology over the last two decades on Argentière glacier, Mont Blanc area[END_REF] . One week prior to the seismic deployment, we conducted a ground penetrating radar campaign over the study area with a 5-MHz signal, to improve the previous estimates of the bed topography reported by ref. [START_REF] Vincent | Glacier thickening and decay analysis from 50 years of glaciological observations performed on Glacier d'Argentière[END_REF] . Also, in September 2018, we conducted an aerial survey to derive a digital elevation model of the glacier surface using stereo-photogrammetry. Combining these two digital elevation models, we calculated the ice thickness distribution over the study area, as shown in Figure 1. The ice thickness reached up to 270 m at the center of the seismic array, with a well-marked talweg (i.e., valley-shaped bed) aligned in the glacier flow direction, and a progressive diminution of the ice thickness down-glacier both within the area covered by our array and from the location of our array towards the glacier terminus. The reader should refer to ref. [START_REF] Gimbert | A Multi-Physics Experiment with a Temporary Dense Seismic Array on the Argentière Glacier, French Alps: The RESOLVE Project[END_REF] for the detailed description of the complementary data associated with this seismic experiment.

Matched field processing

Matched field processing (MFP) consists of providing a probabilistic estimate of the epicentral coordinate of dominant source [START_REF] Kuperman | Matched field acoustics[END_REF] from evaluation of the phase coherence of the seismic signal over an array of sensors. The method consists of recursively matching the phase delays of a model-based synthetic wave-field (i.e., the 'trial source') with the phase delays observed between the sensors over the array. A condition for this method to be applicable is that spatially coherent phase information can be extracted between nearby sensors [START_REF] Kuperman | Matched field acoustics[END_REF] , which requires interstation spacing of less than one-half of the investigated wavelength. We satisfy this conditions with our 40 to 50 m interstation spacing for a 200 to 500m investigated wavelength.

We first compute the discrete Fourier transform of a given data vector 𝑑(𝑡) recorded by the 98 sensors over a frequency ω, to obtain the complex data vector 𝑑(ω) and to calculate the corresponding cross-spectral density matrix as

𝐾(ω) = 𝑑(ω)𝑑 ) (ω),
where 𝐻 is the Hermitian transpose. The cross-spectral density matrix captures the relative spatial phase difference between the sensors. We then define a set of values to be explored for trial sources. In this study, we set four degrees of freedom for the MFP processing, with a depth, range and phase velocity grid (X, Y, Z, C). An important condition for retrieving source epicentral coordinate (i.e. X, Y, Z) and not only source back azimuth is to have a source-to-array distance not greater than two to three times the array aperture (i.e. largest distance between stations) [START_REF] Almendros | Array analysis using circularwave-front geometry: An application to locate the nearby seismo-volcanic source[END_REF] .

We satisfy this condition with a c. 750 m array aperture for a c. 270m source-to-station distances [START_REF] Nanni | Quantification of seasonal and diurnal dynamics of subglacial channels using seismic observations on an Alpine glacier[END_REF] (i.e. maximum glacier thickness [START_REF] Hantz | ice permeability at depth, and water pressures at Glacier d'Argentiere, French Alps[END_REF] ). For each element 𝑎 of these four dimensions of the grid we model the Green's function replica vector 𝑑(ω, a) under the hypothesis of a homogenous medium as 𝑑(ω, a) = exp 1 𝑖ωr 5 𝑐 7 8,

where 𝑐 is the medium velocity and 𝑟 5 is the distance between each receiver and the trial source position 𝑎.

To match the observed cross-spectral density matrix with the replica vector, we calculate the Bartlett processor as

𝐵 ;5<=>?== (ω, a) = @ |𝑑(ω, a) ) 𝐾(ω)𝑑(ω, a) | B .
This operation is equivalent to cross-correlation between the observed wave-field phase and the modeled one. We refer to the 𝐵 ;5<=>?== (ω) values as the MFP output. The MFP output is calculated at specific frequencies, and ranges from 0 to 1. The closer to 1 the MFP output is, the more the modeled phase matches the observations. Sources that generate a signal resulting in similar phase coherence (i.e. MFP output value) are expected to be localized with a similar accuracy [START_REF] Kuperman | Matched field acoustics[END_REF] .

We perform source location over 1-s-long signal segments of the vertical component only. We filter the signal within the [3-7] Hz frequency range that is the most sensitive to subglacial water flow induced seismic noise [START_REF] Bartholomaus | Subglacial discharge at tidewater glaciers revealed by seismic tremor[END_REF][START_REF] Nanni | Quantification of seasonal and diurnal dynamics of subglacial channels using seismic observations on an Alpine glacier[END_REF][START_REF] Vore | Seismic Tremor Reveals Spatial Organization and Temporal Changes of Subglacial Water System[END_REF][START_REF] Lindner | Glaciohydraulic seismic tremors on an Alpine glacier[END_REF] , and coherently apply the MFP for each 0.4 Hz within this range. To maximize the efficiency of our algorithm and minimize the computational costs, we use a gradient-based minimization algorithm (i.e., Nelder-Mead optimization [START_REF] Nelder | Improved simplex method for function minimization[END_REF] ) to converge to the best match between the trial and the observed phase delays, rather than exhaustive grid-search exploration. The convergence criterion is reached when the variance of the values obtained over the last five iterations of the optimization is <1e -2 , with a maximum of 3000 iterations. Our 29 different starting points used for optimization are located 250 m below the glacier surface, and they uniformly cover an area of 800 × 800 m 2 centered on the array (Fig. 2d). We set the initial velocity to 1800 m.s -1 . The 29 point locations found per signal segment (1 s) after convergence are all located in the same place if clear global convergence exists (i.e., high MFP output) or at up to 29 different locations if up to 29 local minima exist (i.e., low MFP output). The MFP output value can also be understood as an indicator of the number of receivers over which the signal is coherent; e.g., a value of 0.1 indicates a source that generates a signal coherent over 10% of the array. We found that this approach is comparable to using sub-arrays of variable sizes to perform the MFP over different areas of the glaciers. Further details on the methods and related physics can be found in ref. [START_REF] Chmiel | Extraction of phase and group velocities from ambient surface noise in a patch-array configuration[END_REF] .

Matched field processing output statistics

We narrow the MFP output selection by keeping only: (1) the location that yields realistic phase velocities ([1200-3600] m.s-1); and ( 2) the localizations at ±400 m from the array center and ≤400 m below the glacier surface. The number of located sources after this selection (Fig. 2a) decreases as the MFP output increases, with >5000 daily sources associated with MFP output <0.05, and <1 daily source associated with MFP output >0.75. Figure 2b, c shows the distribution of the source properties (vertical position), and Figure 2e, f shows the spatial probability of the source location over the complete study period for two MFP output bands of ([0.07-0.16] and [0.75-0.99]). These locations correspond to real and detectable seismic events, since MFP output is much higher than expected for random noise (see ref. [START_REF] Gimbert | A Multi-Physics Experiment with a Temporary Dense Seismic Array on the Argentière Glacier, French Alps: The RESOLVE Project[END_REF] , Fig. 9). This two-dimensional representation is obtained by summing the total number of sources in each 8 × 8 m 2 cell of an 800 × 800 m 2 (x, y) grid centered on our seismic array.

It is important to note that at low MFP output (e.g., [0.07-0.16]) only a few sensors distinguish the source. This leads to a shorter aperture of the effective seismic array, which limits depth resolution. The trend to higher probability of the source location at depth is however observed here for tens of thousands of sources. The 2-day averaging is long enough to gather sufficient statistics, while it is also short enough to correctly investigate the temporal evolution with high enough spatial accuracy (Fig. 3).

Source location precision

When compared to the crevasse field at our study location (Fig. 1), we observe that our MFP analysis yields source locations with MFP output >0.75, a precision in range down to 10 m (Fig.

2

) for an average wavelength λ of ca. 300 m at 5 Hz 37 . While the expected Rayleigh limit gives maximum resolution of λ/4 (ca. 75 m at 5 Hz) in the far-field domains [START_REF] Rost | Array seismology: Methods and applications[END_REF] , laboratory experiments [START_REF] Pyrak-Nolte | synthetic sediments using seismic wave transmission Sample[END_REF] suggest that the seismic signal bears information of spatial structures down to λ/8 (ca. 32 m at 5

Hz) in the near-field of sources. Our unique instrumental set-up combined with our systematic analysis of phase coherence allow us to overcome the Rayleigh far-field limitations and obtain two-dimensional maps of source locations with resolution in the range of ca. 10 m for high MFP outputs. For lower MFP outputs, we expect lower resolution due to the smaller aperture of the array that is sensitive to local phase coherence. We acknowledge that the width of the area shown in Figure 2d might not represent the true width of a single subglacial channel, which is expected to be of the order of one to few meters at this location [START_REF] Nanni | Quantification of seasonal and diurnal dynamics of subglacial channels using seismic observations on an Alpine glacier[END_REF] . The observed width of c. 50 m could be due to the uncertainty for the source location or to the presence of multiple channels

Hydraulic potential calculation

We calculate the hydraulic potential 𝛷 following ref. [START_REF] Shreve | Movement of Water in Glaciers[END_REF] , as:

𝛷 = 𝜌 F 𝑔𝑧 I + 𝐾[𝜌 L 𝑔(𝑧 M -𝑧 I )],
where 𝑔 is the acceleration due to gravity (9.81 m.s -2 ), 𝜌 F is the density of water (1000 kg.m -3 ), 𝜌 L is the density of ice (917 kg.m -3 ), 𝑧 I and 𝑧 M are the elevation of the glacier bed and surface (m), respectively, and 𝐾 is the ratio of water pressure to ice overburden pressure; i.e., a uniform flotation fraction. A value of 𝐾 near 0 represents a condition where basal water pressure is negligible in comparison to ice overburden pressure, and a value near 1 represents a condition where basal water pressure is high enough to significantly counterbalance ice overburden pressure. Here we use a flotation fraction of 0.5, which is likely to represent a case where the subglacial drainage system has developed enough to efficiently drain the bed [START_REF] Shreve | Movement of Water in Glaciers[END_REF][START_REF] Willis | Subglacial drainage system structure and morphology of Brewster Glacier, New Zealand[END_REF] , and therefore reduce subglacial basal water pressure. We consider this flotation fraction to best represent the subglacial water pressure condition in our location at the end of our study period when we observe a decreasing glacier velocity while water discharge still increases (Fig. 4a). We then calculate the water flow directions based on the hydraulic potential gradients, by following the path that minimizes the gradient with a minimum upstream area of 150 m 2 for the waterway, to initiate the use of TopoToolBox, developed by ref. [START_REF] Schwanghart | Short Communication: TopoToolbox 2 -MATLAB-based software for topographic analysis and modeling in Earth surface sciences[END_REF] . We show in

Inverting hydraulic properties using water discharge and seismic power measurements

We use the theoretical framework of ref. [START_REF] Gimbert | Subseasonal changes observed in subglacial channel pressure, size, and sediment transport[END_REF] to invert for the hydraulic properties of subglacial turbulent water flow using simultaneous measurements of subglacial water discharge and seismic power. The seismic power 𝑃 is calculated at each sensor using the vertical component of ground motion within [3-7] Hz and the Welch method over 4-s time windows with 50% overlap, as reported in ref. [START_REF] Nanni | Quantification of seasonal and diurnal dynamics of subglacial channels using seismic observations on an Alpine glacier[END_REF] . The physical framework relates the changes in the measured seismic power P and discharge Q to changes in channel hydraulic radius R and hydraulic pressure gradient S. The basis of these relations is that turbulence within water flow generates frictional forces that act on the boundaries of the channels (i.e., ice, glacier bed) and create ground motion. This approach has been recently shown to be applicable over our study area [START_REF] Nanni | Quantification of seasonal and diurnal dynamics of subglacial channels using seismic observations on an Alpine glacier[END_REF] . The hydraulic radius R is defined as the ratio of the cross-sectional area of the channel flow to its wetted perimeter. The hydraulic pressure gradient S is a function of both the rate of change of water pressure and the bed slope in the flow direction. In the case of constant bed slope and channel geometry, increasing S means closed and pressurizing channel flow. Both variables are defined for subglacial conduits where turbulent water flow dominates. We calculate median values of R and S using the median seismic power P over the 98 sensors and the water discharge Q, as follows 𝑆 = 𝑆 <?V ( 𝑃 𝑃 <?V

) WX/XZ ( 𝑄 𝑄 ?V

) \]^/XZ , 𝑅 = 𝑅 <?V ( P 𝑃 <?V

) \a/bW ( 𝑄 𝑄 ?V

) \]]/bW

where the subset ref is the reference state, which is April 26 in this study. In the main text, when referring to the temporal evolution of R and S, we therefore refer to their relative changes with respect to this reference period. We evaluate P over the [3-7] Hz frequency band, as it has been shown by ref. [START_REF] Nanni | Quantification of seasonal and diurnal dynamics of subglacial channels using seismic observations on an Alpine glacier[END_REF] to be the best suited in area for studying turbulent water flow. The reader is referred to refs. [START_REF] Gimbert | Subseasonal changes observed in subglacial channel pressure, size, and sediment transport[END_REF][START_REF] Nanni | Quantification of seasonal and diurnal dynamics of subglacial channels using seismic observations on an Alpine glacier[END_REF] for more details.

We stress that the inversions of hydraulic properties depend on the average seismic amplitude, with very little dependency on the spatial variations. Therefore, our inversions are to be considered independent of the source location, even if they both arise from the seismic signals.

This also implies that our inversions represent averaged hydraulic conditions of the drainage system.

State of knowledge on the main features of channels and cavities

Cavities form in the downstream lee of bedrock bumps, due to ice sliding over the bed, and they close through ice creep. Cavities can be filled with water, which reduces the apparent rugosity of the glacier bed and weakens the ice-bed mechanical coupling [START_REF] Lliboutry | Modifications to the theory of intraglacial waterways for the case of subglacial ones[END_REF] . For soft-bed glaciers, inefficient drainage systems also include unconsolidated layers of low permeability [START_REF] De Fleurian | A double continuum hydrological model for glacier applications[END_REF] . Weakly connected cavities can have permeability that are lower by ca. 9 orders of magnitude than that of connected cavities [START_REF] Hoffman | Greenland subglacial drainage evolution regulated by weakly connected regions of the bed[END_REF] . For a similar hydraulic gradient to that of connected cavities, turbulent water flow within a weakly connected system therefore requires much higher flow velocities.

Subglacial channels can be of the R-type 8 when melted into the ice by turbulent dissipation of heat, or of the N-type [START_REF] Nye | The Flow of a Glacier in a Channel of Rectangular, Elliptic or Parabolic Cross-Section[END_REF] [START_REF] Gimbert | A Multi-Physics Experiment with a Temporary Dense Seismic Array on the Argentière Glacier, French Alps: The RESOLVE Project[END_REF] . Contact U.N. at ugo.nanni0158@gmail.com or ugo.nanni@univ-grenoble-alpes.fr for questions related to the above described datasets and codes. The closer this coefficient is to 1, the more the drainage system is observed to be channelized.

Figures and Tables

(b) Subglacial water discharge (blue line) with shaded blue area under the curve added for ease of reading. Median surface velocity over the four on-ice GNSS stations (red line) with the minimum/maximum shaded envelope.

Supplementary Figure 1.

  Figure Supp. 1 that different values of flotation fraction (𝐾 = [0.1,0.5,0.9]) result in a similar predicted waterway pattern. The expected localization of subglacial water flow thus does not depend much on the flotation fraction, likely as a result of the V-shape of the Argentière glacier being particularly pronounced.

1 .

 1 Figure 1.

Figure 2 .

 2 Figure 2.

Figure 4 .

 4 Figure 4.

  (a) Left axis shows the relative hydraulic properties, as averaged over the 98 sensors. Hydraulic pressure gradient with ±1 median absolute deviation (MAD), as the shaded envelope (green), and hydraulic radius with ±1 MAD, as the shaded envelope (purple). Both values are expressed as relative to April 26. Right axis shows temporal evolution of the spatialized seismic observations using the determination coefficient R2 of the source location pattern shown in Figure 3o (red).

  when dug into the basal sediments by the flowing water or etched into bedrock by carbonate dissolution[START_REF] Hallet | The Effect of Subglacial Chemical Processes on Glacier Sliding[END_REF] . Both types close through ice creep. When developing and reaching steady-state, a subglacial channel tends to have lower water pressure than cavities, which therefore drains the connected cavity system.The time series of physical quantities measured at Argentière glacier over the 2017 and 2018 melt-season can be found at https://doi.org/10.5281/zenodo.3701520. All of the data associated with the dense array experiment can be found at https://doi.org/10.5281/zenodo.3971815. The spatial maps of source locations, glacier geometries and hydraulic potentials as well as the code used to processed these data can be found at https://doi.org/10.5281/zenodo.4024660. Part of the dataset of the seismic signals acquired during RESOLVE-Argentière can be found at https://doi.org/10.15778/RESIF.ZO2018 in link with ref.

Acknowledgments

This work has been conducted in the framework of the RESOLVE Project (https://resolve.osug.fr/) (LabEx OSUG@2020, Investissement d'avenir -ANR10LABX56 and IDEX Université Grenoble Alpes). Most of the computations presented in this paper were performed using the GRICAD infrastructure (https://gricad.univ-grenoble-alpes.fr), which is supported by Grenoble research communities, and with the CiGri tool (https://github.com/oar-team/cigri) that was developed by Gricad, Grid5000 (https://www.grid5000.fr) and LIG (https://www.liglab.fr/). FG acknowledges support from ANR SEISMORIV (ANR-17-CE01-0008) and SAUSSURE (ANR-18-CE01-0015-01).

We thank Christian Vincent and Nathan Maier for feedback on paper writing and data interpretation. We thank Agnès Helmstetter and Benoit Urruty for fruitful discussion on data analysis. We thank Mondher Chekki for numerical support dealing with the dense seismic array big dataset. U.N. thanks Aurélien Mordret, Amandine Sergeant, Léonard Seydoux, Jean Soubestre, and Josefine Umlauft for discussions on locating spatially spread seismic noise sources, and we thank Camillo Rada for discussions on subglacial drainage. U.N. thanks Pete Akers, Jordi Bolibar, Fanny Brun, David Laliche, Fabien Moustard, Hans Segura and Jonathan Wille for fruitful discussions on this work.

Codes and data availability


Claude Bardos
email: claude.bardos@gmail.com

François Golse
email: francois.golse@polytechnique.edu

Olivier Pironneau
email: olivier.pironneau@sorbonne-universite.fr

Remarks on the Radiative Transfer Equations for Climatology

Keywords: Radiative transfer, Electronic Temperature, Integral equation, Numerical analysis AMS 3510, 35Q35, 35Q85, 80A21, 80M10

Using theoretical and numerical arguments we discuss some of the commonly accepted approximations for the radiative transfer equations in climatology.

Introduction

Satellite, atmospheric and terrestrial measurements are numerous to support the global warming. Various hypothesis are made to explain the disastrous effect of the Green-House Gases (GHG). However a rigorous proof derived from the fundamental equations of physics is not available and one of the problems is the complexity of the physical system of planet Earth in its astronomical setting around the Sun. In many references such as [8],[9], [START_REF] Chandrasekhar | Radiative Transfer[END_REF] and [START_REF] Fowler | Mathematical geoscience[END_REF] it is explained that the Sun radiates light with a heat flux Q = 1370Watt/m 2 , in the frequency range (0.5, 20) × 10 14 Hz corresponding approximately to a black body at temperature of 5800K; 70% of this light intensity reaches the ground because the atmosphere is almost transparent to this spectrum and about 30% is reflected back by the clouds or the ocean, snow, etc (albedo). The Earth behaves almost like a black body at temperature T e = 288K and as such radiates rays of frequencies ν in the infrared spectrum (0.03, 0.6) × 10 14 Hz. "The absorption coefficient variation with frequencies is such that the atmosphere is essentially transparent to solar radiation" ( [START_REF] Fowler | Mathematical geoscience[END_REF], p65). "Solar radiation that is not absorbed or reflected by the atmosphere (for example by clouds) reaches the surface of the Earth. The Earth absorbs most of the energy reaching its surface, a small fraction is reflected 4 ." Carbone dioxide renders the atmosphere opaque to infrared radiations around 20THz . Hence increasing its proportion in air increases the absorption coefficient in that range. It is believed to be one of the causes of global warming. However it must be more complex because some experimental measurements show that increasing opacity decreases the temperature at high altitude, above the clouds [START_REF] Dufresne | Greenhouse effect: The relative contributions of emission height and total absorption[END_REF].

In the framework of the radiative transfer equations, one can studies the conjectures cited above. To support theoretical claims, the numerical solver developed in [START_REF] Bardos | Radiative Transfer for the Greenhouse Effect[END_REF], [10] and [7] is used. We will address 4 questions:

• What is the effect of increasing an altitude dependent absorption coefficient?

• Is it true that sunlight crosses the Earth atmosphere unaffected?

• Is 30% of albedo equivalent to 30% reduction of the Sun radiative energy ?

• What is the effect on the temperature of increasing the absorption coefficient in a specific frequency range?

The Radiative Transfer Equations for a Stratified Atmosphere

Finding the temperature T in a fluid heated by electromagnetic radiations is a complex problem because interactions of photons with atoms in the medium involve rather intricate quantum phenomena. Assuming local thermodynamic equilibrium leads to a well-defined electronic temperature. In that case, one can write a kinetic equation for the radiative intensity I ν (x, ω, t) at time t, at position x and in the direction ω for photons of frequency ν, in terms of the temperature field T (x, t).

1 c ∂ t I ν + ω • ∇I ν + ρκ ν a ν I ν -1 4π S 2 p(ω, ω )I ν (ω )dω = ρκ ν (1 -a ν )[B ν (T ) -I ν ]. (1)
In this equation c is the speed of light, ρ is the density of the fluid, ∇ designates the gradient with respect to the position x, while

B ν (T ) = 2 ν 3 c 2 [e ν kT -1] (2)
is the Planck function at temperature T , with the Planck constant and k the Boltzmann constant. Recall the Stefan-Boltzmann identity,

∞ 0 B ν (T )dν = σT 4 , σ = 2π 4 k 4 15c 2 3 , (3)
where πσ is the Stefan-Boltzmann constant.

The intricacy of the interaction of photons with the atoms of the medium is contained in the mass-absorption κν , which is theoretically the result of vibration and rotation atomic resonance, but for practical purpose the fraction of radiative intensity at frequency ν that is absorbed by fluid per unit length. The coefficient a ν ∈ (0, 1) is the scattering albedo, and 1 4π p(ω, ω )dω is the probability that an incident ray of light with direction ω scatters in the infinitesimal element of solid angle dω centered at ω. The kinetic equation ( 1) is coupled to the temperature -or energy conservation -equations of the fluid . When thermal diffusion is small, the system decouples and energy balance becomes:

∞ 0 ρκ ν (1 -a ν ) S 2 I ν (ω)dω -4πB ν (T ) dν = 0 . (4)
When the wave source is far in the direction z the problem becomes one-dimensional in z and the radiative intensity scattered in the direction ω depends only on µ, the cosine of the angle between ω and Oz. Consequently the system becomes,

µ∂ z I ν + ρκ ν I ν = ρκ ν (1 -a ν )B ν (T ) + ρκ ν a ν 2 1 -1 p(µ, µ )I ν (z, µ )dµ , ∞ 0 ρκ ν (1 -a ν ) B ν (T ) -1 2 1 -1 I ν dµ dν = 0, z ∈ (0, H), |µ| < 1, ν ∈ R + .
(5)

For this system to be mathematically well posed, the radiation intensity I ν must be given on the domain boundary where radiation enters. For example,

I ν (H, -µ) = Q -(µ), I(0, µ) = Q + (µ), 0 < µ < 1 . (6)
2.1. Semi-analytical Solution when the scattering is isotropic Isotropic scattering is modelled by taking p ≡ 1. By introducing an optical depth (in [START_REF] Fowler | Mathematical geoscience[END_REF] and others, the signs are changed),

τ = z 0 ρ(η)dη, ρ becomes 1 and z becomes τ ∈ (0, Z), Z = H 0 ρ(η)dη. (7)
Let us denote the exponential integral by

E p (X) := 1 0 e -X/µ µ p-2 dµ . (8)
Define

J ν (τ ) = 1 2 1 -1 I ν dµ, S ν (τ ) = 1 2 1 0 e -κν τ µ Q + ν (µ) + e -κν (Z-τ ) µ Q - ν (µ) dµ. ( 9
)
The problem is equivalent to a functional integral equation:

J ν (τ ) = S ν (τ ) + κ ν 2 Z 0 E 1 (κ ν |τ -t|) ((1 -a ν )B ν (T (t)) + a ν J ν (t)) dt , ∞ 0 (1 -a ν )κ ν B ν (T (τ ))dν = ∞ 0 (1 -a ν )κ ν J ν (τ )dν,        (10)
Once J ν , T are known, I ν can be recovered by

I ν (τ, µ) =e -κν τ µ I ν (0, µ)1 µ>0 + e κν (Z-τ ) µ I ν (Z, µ)1 µ<0 + 1 µ>0 τ 0 e -κν |τ -t| |µ| κν |µ| ((1 -a ν )κ ν B ν (T (t)) + a ν κ ν J ν (t)) dt + 1 µ<0 Z τ e -κν |τ -t| |µ| κν |µ| ((1 -a ν )κ ν B ν (T (t)) + a ν κ ν J ν (t)) dt. (11)
Remark 1. Notice that two solutions (T, J ν ), (T , J ν ) with J ν ≡ J ν will have T ≡ T but not necessarily I ν ≡ I ν .

Increasing an Altitude Dependent Absorption Coeficient Decreases the Temperature in the atmosphere

The following argument goes against a common belief which claims that increasing the absorption coefficient -which GHG do indeed -implies an increase of temperature in the atmosphere. Consider the case of an altitude dependent absorption coefficient which we write as r(z)κ with κ constant.

Proposition 1. If κ is independent of ν but function of the altitude and T is monotone decreasing with altitude, then, increasing κ anywhere will decrease the temperature.

Proof The proof is a straightforward consequence of the concept of optical length (7).

Here, let τ =

µ∂ τ Ī1 + κ Ī1 = 1 2 1 -1 κ(1 -a) Ī1 + κa 1 -1 p(µ, µ ) Ī1 (µ, µ )dµ dµ, I r (H, -µ) = 0, Īr (0, µ) = Q + (µ), 0 < µ < 1. ( 12
)
where the Stefan-Boltzmann identity has been used:

σT 4 = 1 2 1 -1
Īr dµ.

Let Ī1 (τ, µ) be the solution of (12). We want to compare I r+δr with I r when δr has a small support near altitude y. By (7),

Īr+δr (z) = Ī1 ( z 0 (r(ζ) + δr(ζ))dζ) ≈ Ī1 (τ (z)) + d dτ Ī1 (τ (y))δr(y) =⇒ 1 2 1 -1 Īr+δr (z)dµ - 1 2 1 -1 Īr (z)dµ ≈ d dτ 1 2 1 -1 Ī1 (τ )dµ | τ (y) δr(y) (13)
Finally from (3),

T r+δr (z) -T r (z) ≈ 4σT 3 (y) δr(y) r(y) dT r dz (y). ( 14
)
4. Sunlight Crosses the Earth Atmosphere Unaffected?

It is true that the absorption coefficient κ ν is small in the frequency range of visible light; but does it imply that the light source can be used as a boundary condition at altitude zero instead of altitude H, the top of the troposphere? First eliminate the z dependency of the coefficients by using the optical length τ defined in (7). The following result gives an interesting answer to this question.

Proposition 2. Assume a ν = 0 and assume that for some ν * ,

:= max ν>ν * κ ν << 1 and κ ν Q 0 ν | ν>ν * = O(1).
Let I ν , T and I ν , T be solutions of (5),( 6) respectively with

Q + ν (µ) = 0, Q - ν (µ) = µQ 0 ν , 0 < µ < 1. Q + ν = µe -κν Z µ Q 0 ν , Q - ν (µ) = 0, 0 < µ < 1. ( 15
)
Then T (τ ) = T (τ ) + O( ) and

I ν -I ν = µe -κν Z µ Q 0 ν 1 µ>0 e -κν µ τ + 1 µ<0 e κν µ τ + O( ).
Proof Note that

µ∂ τ I ν + κ ν (I ν -B ν (T )) = 0, ν < ν * , I ν (0)| µ>0 = 0, I ν (Z) µ<0 = Q - ν (-µ), ν < ν * , I ν ≈ I * ν (µ) independent of τ , ν > ν * . I * ν | µ>0 = 0, I * ν | µ<0 = Q - ν (-µ), ν > ν * . ∞ 0 κ ν B ν (T )dν ≈ ν * 0 κ ν 2 1 -1 I ν dµdν + ∞ ν * κ ν 2 0 -1 Q - ν (-µ) dµdν                    (16)
Now I ν is defined by (5) but with

I ν (0)| µ>0 = Q + ν (µ), I ν (Z)| µ<0 = 0.
By the same argument, when ν > ν * , I ν ≈ I * ν (µ), independent of τ , and

µ∂ τ I ν + κ ν (I ν -B ν (T )) = 0, ν < ν * , I ν (0)| µ>0 = Q + ν (µ), I ν (Z)| µ<0 = 0, ν < ν * , ∞ 0 κ ν B ν (T )dν ≈ ν * 0 κ ν 2 1 -1 I ν dµdν + ∞ ν * κ ν 2 1 0 Q + ν (µ)dµdν          (17)
Let I ν = I ν -I ν . It holds

µ∂ τ I ν + κ ν (I ν -(B ν (T ) -B ν (T )) = 0, ν < ν * , ∞ 0 κ ν (B ν (T ) -B ν (T ))dν = ν * 0 κ ν 2 1 -1 I ν dµdν + ∞ ν * κ ν 2 1 0 (Q - ν (µ) -Q + ν (µ)) dµ dν I ν (0)| µ>0 = -Q + ν (µ), I ν (Z)| µ<0 = Q - ν (-µ), ν < ν * .                    (18)
We notice that T = T and

Q - ν (µ) = e -κν µ Z Q + ν (µ) implies µ∂ τ I ν + κ ν I ν = 0, ν < ν * , ν * 0 κ ν 2 1 -1 I ν dµ dν = ∞ ν * κ ν 2 1 0 (1 -e -κν µ Z )Q + ν (µ) dµ dν I ν (0)| µ>0 = -Q + ν (µ), I ν (Z)| µ<0 = e κν µ Z Q + ν (-µ)          (19) Let Q + ν (µ) = µQ 0 ν . Then I ν (0)| µ>0 = -µQ 0 ν , I ν (Z)| µ<0 = -µe κν µ Z Q 0 ν .
Consequently, from ( 19) and ( 11), the solution is

I ν (τ, µ) = -µQ 0 ν 1 µ>0 e -κν µ τ + 1 µ<0 e κν µ τ . =⇒ 1 -1 I ν (τ )dµ = -Q 0 ν E 3 (κ ν τ ) + Q 0 ν E 3 (κ ν τ ) = 0, (20)
We have satisfied (19) with precision O(κ 2 ν ) because

κ ν (1 -e -κν µ Z )µQ 0 ν = κ 2 ν ZQ 0 ν + O(κ 3 ν ) ≈ 0, ν > ν * Remark 2. It is easy to see that applying e -κν Z τ µQ 0 ν at τ = 0, leads to change E 3 (κ ν τ ) in (11) into E 3 (κ ν (Z + τ )).
The change in temperature is shown on Figure 1. All variables are rescaled as in [START_REF] Golse | Stratified radiative transfer in a fluid[END_REF] §7.1. In particular temperatures are in Kelvin and divided by 4780. The air density is 1.22e -z . Here Q 0 = 5.66 • 10 -5 and T sun = 1.2. Conclusion: For the specific case of sunlight on Earth It is a gross approximation to assume that sunlight crosses the atmosphere unaffected. However the theory works when T sun = 2.4 and ν * = 0.2, probably because max ν B ν (T sun ≈ 18 against 2 with T sun = 1.2.

Is the 30% Earth Albedo Equivalent to a 30% reduction of the Sun Radiative Energy?

Note that (11) can be used to find the solution of (5)-a with generalized albedo (inspired by the model proposed in [START_REF] Dufresne | Greenhouse effect: The relative contributions of emission height and total absorption[END_REF]):

I ν (0, µ) = k α k I ν (τ k , -µ) + Q 0 ν µ, I ν (Z, -µ) = 0, 0 < µ < 1.
with k = 1, 2, ..K. Indeed by (11) Altitude -km Scaled Temperature (K)

I ν (τ, µ) =e -κν τ µ k α k I ν (τ k , -µ) + Q 0 ν µ 1 µ>0 + 1 µ>0 τ 0 e -κν |τ -t| |µ| κν |µ| ((1 -a ν )κ ν B ν (T (t)) + a ν κ ν J ν (t)) dt + 1 µ<0 Z τ e -κν |τ -t| |µ| κν |µ| ((1 -a ν )κ ν B ν (T (t)) + a ν κ ν J ν (t)) dt
T 0 T 1 T 2 T 3 T 4 T 5 Figure 1: Scaled temperatures z → T (z) versus altitude, computed with 0.5(0.1 + 1 ν<3 ). The 3 curves T 1 , T 2 , T 3 corresponds to 1: Q - ν = Q 0 ν , Q + ν = 0, 2: Q - ν = 0, I ν (0, µ)| µ>0 = µe -κν Z µ Q 0 ν and to 3: Q - ν = 0, Q + ν = Q 0 .
These can be compared to the grey case T 0 with κ = 0.5. Note that T 1 is very different from T 2 . However if the Sun temperature is multiplied by 2 and κ ν = 0.5(0.1 + 1 ν<0.2 ), then the conditions in Proposition 2 are met and confirmed by the numerical simulations:

T 4 is computed like T 1 and T 5 is computed like T 2 . = 1 µ>0 e -κν τ µ k α k Z τ k e -κν (t-τ k ) µ κν µ ((1 -a ν )κ ν B ν (T (t)) + a ν κ ν J ν (t)) dt + 1 µ>0 e -κν τ µ Q 0 ν µ + τ 0 e -κν |τ -t| |µ| κν |µ| ((1 -a ν )κ ν B ν (T (t)) + a ν κ ν J ν (t)) dt + 1 µ<0 Z τ e -κν |τ -t| |µ| κν |µ| ((1 -a ν )κ ν B ν (T (t)) + a ν κ ν J ν (t)) dt (22) =⇒ J ν (τ ) = 1 2 Q 0 ν E 3 (κ ν τ ) + 1 2 Z 0 E 1 (κ ν |τ -t|) [(1 -a ν )B ν (T (t)) + a ν κ ν J ν (t)] dt , + 1 2 k α k κ ν Z τ k E 1 (κ ν (t + τ -τ k )) [(1 -a ν )κ ν B ν (T (t)) + a ν κ ν J ν (t)] dt.

Earth Albedo

Consider the case

I ν (Z, µ)| µ<0 = 0, I ν (0, µ)| µ>0 = Q 0 ν + αI ν (0, -µ) . Assume no scattering, τ 1 = 0, K = 1.
As in [7], the following iterative scheme is considered: Altitude -km Scaled Temperature (K)

J n+1 ν (τ ) = 1 2 Q 0 ν E 3 (κ ν τ ) + κ ν 2 Z 0 (E 1 (κ ν |τ -t|) + αE 1 (κ ν (t + τ ))) B ν (T n (t))dt , ∞ 0 κ ν B ν (T n+1 τ )dν = ∞ 0 κ ν J n+1 ν (τ )dν . (23
T, Q 0 T, 0.7Q 0 T T "
Figure 2: Scaled temperatures z → T (z) versus altitude. In all cases Q - ν = 0. The doted curves are T when Q + ν = λQ 0 ν with λ = 1 (top dotted curve) and λ = 0.7 (lower dotted curve); the doted line is the solution

T when I ν (0, µ)| µ>0 = 0.3I ν (0, -µ) + µe -κν Z µ Q 0 ν /0.7 . The dashed line is when I ν (0, µ)| µ>0 = 0.3I ν (0, -µ) + µQ 0 ν .

Numerical Example

Assume a ν = 0. Figure 2 compares the numerical solutions of ( 5) with κ ν = 0.5[1 ν<6 +0.1] and Q - ν = 0, Q + ν = λQ 0 ν , λ = 1 or 0.7. with the solution of (5) with β = 0 or 1 in

Q - ν = 0, I ν (0, µ) = 0.3I ν (0, -µ) + µQ 0 ν [ 1 0.7 e -κν Z µ β + 1 -β], 0 < µ < 1.
This numerical simulation is an attempt to replace the radiative source at τ = H by a radiative source at τ = 0 which gives the same result. Notice alongside that adding 0.3 albedo induces a heating of the atmosphere (i.e. comparing both computations with λ = 1, β = 1 above).

Earth Albedo, General Statement, Accommodation Coefficient

As above one considers a stratified atmosphere ( 5) with (z, µ, ν) in (0, H) × (-1, 1) × R + ; the radiation intensity I ν (z, µ) satisfies on at z = H an incoming condition

I ν (H, -µ) = Q - ν (µ), 0 < µ < 1, (24)
and on Earth, at z = 0, some "albedo" condition, denoted:

I ν, (0, µ) = A ν (I ν (0, -µ), 0 < µ < 1, (25)
with A representing the albedo effect of the earth, hence being an operator from the space of outgoing intensities into the space of incoming intensities.

Below will be given, based on energy estimates, some sufficient conditions on A for existence, stability -hence uniqueness -of the solution, beginning with the grey model and then considering a natural operator which involves the effective temperature of Earth (see [START_REF] Fowler | Mathematical geoscience[END_REF] page 66).

6.1. The energy type estimate for the grey problem For clarity we present the case without scattering (a ν = 0) and with a constant κ ν independent of the frequency ν. By rescaling z with κ, the problem can be formulated

in term of I(z, µ) = ∞ 0 I ν (z, µ, )dν as, µ∂ z I(z, µ) + I(z, µ) - 1 2 1 -1 I(z, µ )dµ = 0 , (26a)
I(H, -µ) = Q(µ) := Q - ν (µ), I(0, µ) = A(I(0, -µ)), 0 < µ < 1 . (26b)
Note that we have assumed that A commutes with the integration in ν. First observe that the operator

I → I - 1 2 1 -1 I(µ )dµ is, in the space L 2 (-1, 1
) the orthogonal projection on functions of mean value 0, because

1 -1 (I - 1 2 1 -1 I(µ )dµ )I(µ)dµ = 1 -1 (I - 1 2 1 -1 I(µ )dµ ) 2 dµ . (27)
This implies that the left hand side of ( 27) is nonnegative and equal to 0 iff I = 1 2

1 -1 I(µ )dµ is independent of µ . Multiplying the equation (26a) by I(z, µ) , integrating over (0, H) × (-1, 1) , using Green's formula and inserting in the computations the conditions (26b) one obtains:

H 0 dz 1 -1 (I - 1 2 1 -1 I(µ )dµ ) 2 dµ + 1 2 1 0 µ(I(H, -µ) 2 dµ + 1 2 1 0 µ (I(0, µ)) 2 -(A(I(0, µ)) 2 dµ = 1 2 1 0 µ(Q(µ)) 2 dµ (28)
This estimate indicates that any property which implies the positivity of

1 0 µ I(0, µ) 2 -(A(I(µ, 0)) 2 dµ
will lead to a well posed problem with a unique nonnegaative solution; for instance:

Theorem 1. Assume that the operator I → A(I) restricted to the positive cone C + (I(µ) ≥ 0) on the space (L 2 (-1, 1), µ, dµ) then for any positive µ → Q(µ) there is a unique nonnegative solution of the (26a),(26b).

Remark 3. Making use of the maximum principle one could prove similar results when the albedo operator is a contraction in the positive cone of L ∞ (0, 1) .

Remark 4. The albedo condition which is a relation, for 0 < µ < 1 , between I(0, µ) µ>0 and I(0, µ)| µ<0 can also be expressed in term of I(0, µ)| µ>0 only.

Indeed, define the operator T : I + (µ) → T (I + (µ)) = I(0, -µ)| µ>0 (defined for instance in L ∞ (0, 1)) by solving (26a) with

I(H, -µ) = Q(µ), I(0, µ) = I + (µ), 0 < µ < 1 . (29)
Then one has:

I(0, µ)| µ>0 = A(I(0, -µ)) ⇔ I(0, µ)| µ>0 = A(T (I(0, µ)| µ>0 )) . ( 30
)
With Q(µ) > 0 in L ∞ (0, 1), T is an affine operator which preserves the positivity and the monotonicity. and this leads to the simple but useful proposition, which , in turn, is the key to proof of Theorem 1.

Proposition 3. Consider two solutions {I i } i=1,2 of the (26a) with (26b) with the same Q(µ) but with two different albedo operators A i . Assume that both A i are linear contractions and one of them is a strict linear contraction, which preserve positivity. If

∀f ≥ 0 ∈ L ∞ (0, 1) ∀µ ∈ (0, 1) , A 2 (f )(µ) ≤ A 1 (f )(µ) ( 31
)
then one has the same ordering for the corresponding solutions :

∀(z, µ) ∈ (0, 1) × (0, H) I 1 (z, µ) ≤ I 2 (z, µ) . (32)
Proof Using the linearity of the operators A i one observes that R = I 1 -I 2 solves(26a) with R(H, -µ) = 0 and, for all 0 < µ < 1,

R(0, µ) - 1 2 (A 1 + A 2 )(T (R(0, µ)) = 1 2 (A 1 -A 2 )(T ((I 1 (0, µ) + I 2 (0, µ))) . (33)
Moreover, restricted to solutions which satisfy R(H, -µ) = 0, µ > 0 the operator T is a strict monotonicity preserving linear contraction; the same observation holds for

R → 1 2 (A 1 + A 2 )(T (R). As a consequence the operator R → (I -1 2 (A 1 + A 2 )(T )
) is invertible with inverse given by the Neumann series (also preserving the positivity):

R = k≥0 ( 1 2 (A 1 + A 2 )(T ))) k 1 2 (A 1 -A 2 )(T ((I 1 + I 2 )| z=0,µ>0 . (34)
I 1 and I 2 are positive intensities and for 1 2 (A 1 -A 2 )(T ((I 1 +I 2 )| z=0,µ>0 and the righthandside of (34), (31) is also true.

The frequency dependent case

Let us come back to (5) with (26b). Let us denote f + = max(f, 0). When not ambiguous, let us write I ν (µ) instead of I ν (0, µ).

Definition 1. An operator

A defined in L 1 µ (0, 1) is non-accretive if ∀I 1 , I 2 ∈ L 1 µ (0, 1), 1 0 µ (I 2 -I 1 ) + -(A(I 2 ) -A(I 1 )) + dµ ≥ 0 . ( 35
)
Remark 5. Obviously (35) follows from the same pointwise property:

∀(I 1 , I 2 ) ((A(I 2 ) -A(I 1 )) + ≤ (I 2 -I 1 ) + , µ > 0 . ( 36
)
Let us uses the symbol f

= ∞ 0 1 2 1 -1 f ν (µ)dµdν.
Then for two solutions (I 2 ν , T 2 ), (I 1 ν , T 1 ), by [START_REF] Golse | The Milne Problem for the Radiative Transfer Equations (with Frequency Dependence)[END_REF], according to Theorem 4.1 in [7] one has:

- d dz µ(I 2 ν -I 1 ν ) + = = ρκ ν (1 -a ν )((I 2 ν -I 1 ν ) -(B ν (T 2 ) -B ν (T 1 )))1 I 2 ν >I 1 ν + ρκ ν a ν (I 2 ν -I 1 ν ) -1 2 1 -1 p(µ, µ )(I 2 ν -I 1 ν )(z, µ )dµ 1 I 2 ν >I 1 ν ≥ 0 (37)
and therefore d dz µ(I 2 ν (z, µ) -I 1 ν (z, µ)) + ≤ 0 . On the other hand

I 2 ν (H, -µ) ≤ I 1 ν (H, -µ), µ > 0 =⇒ µ(I 2 ν (H, µ) -I 1 ν (H, µ)) + = 0
while the non accretivity of A gives

µ(I 2 ν (0, µ) -I 1 ν (0, µ)) + = 1 2 ∞ 0 dν 1 0 µ(A(I 2 ν (0, -µ)) -A(I 1 (0, -µ))) + -(I 2 ν -I 1 ν ) + (0, -µ))dµ ≤ 0 (38)
Therefore, unless µ(I 2 ν (0, µ) -I 1 ν (0, µ)) + = 0 there is a contradiction because it is a function which is decreasing, negative at z = 0 and zero at z = H. Then consider the difference of equations (5) with arguments (I i ν , T i ) i = 1, 2 multiplied by µ κν 1 I 2 ν >I 1 ν . Following the end of the proof of Theorem 4.1 in [7] (see also [START_REF] Golse | The Milne Problem for the Radiative Transfer Equations (with Frequency Dependence)[END_REF]), one finds by the same method that existence and uniqueness derives from

d dz µ 2 κν (I 2 ν -I 1 ν ) + = 0 , 0 < z < H .

Eventually one has:

Theorem 2. Consider (5) with incoming data 0 ≤ I ν (H, -µ) ≤ B ν (T M ), depending on a given (but possibly large) temperature T M , and with outgoing data given by an accretive operator, of the form I ν (0, µ)| µ>0 = A(I ν (0, -µ))| µ>0 . The problem has a unique well defined nonnegative solution.

Examples of Albedo operators

Below are given some examples of albedo operators which combine an accommodation parameter α and the reflection and the thermalisation effects of the Earth.

The simplest example would be given (cf. Remark 3) with 0 ≤ α ≤ by the formula, p > 0:

A(I(0, µ)) = αI(0, -µ) + (1 -α) 1 0 µ p I(0, -µ )dµ . ( 39
)
For the more realistic frequency dependent cases, one should compare the reflection with the emission from the Earth as global black body under an effective temperature T e 288 and this leads to try the albedo operator

A(I ν (0, -µ)) = αI ν (0, -µ) + (1 -α)B ν (T e ), µ > 0. ( 40
)
This obviously satisfies the hypothesis of non accretivity.

Calculus of Variations for the ν-dependent Case

The following argument sheds some light on the conditions needed to obtain a cooling (resp. heating) from a local increase of the absorption coefficient.

Conjecture 1. Let I ε ν (τ, µ), T ε (τ ) be the solution of (5) when

κ ν := κ + δκ ν , (41)
Let I 0 ν (τ, µ), T 0 (τ ) be the solution of (5) with κ ν = κ constant. When δκ ν ≈ δ ν * , the Dirac mass at ν * , and κ is small, then the sign of dT dε (τ )| ε=0 is governed by the sign of

1 2 1 -1 I 0 ν * (τ )dµ -B ν * (T 0 (τ )). (42)
7.1. Calculus of Variations Support of the Conjecture It will be convenient to define the Planck function in terms of the quantity given by the Stefan-Boltzmann law Φ := σT 4 . Henceforth we denote

B ν (T ) = b ν (Φ) .
We seek to study how the average radiation intensity or the temperature is altered if the absorption is modified on various intervals in the frequency variable. We expect that the simplest situation corresponds to absorption of the form (41) with κ > 0 independent of ν, while 0 < 1. This is of course an extremely general formulation, but in practice one could think of δκ ν := κ1 ν 1 <ν<ν 2 .

This corresponds to multiplying the absorption by (1+ ) in the frequency interval (ν 1 , ν 2 ) only, and leaving it invariant for all the other frequencies. We assume that the scattering is isotropic with constant rate λ := κa ≥ 0.

Henceforth, we consider the radiative transfer equation on (0, Z) × (-1, 1) × R,

µ∂ τ I ν + κ ν (I ν -b ν (Φ )) + λ(I ν -I ν ) = 0 , κ ν ( I ν -b ν (Φ )) = 0 , (43)
with the boundary conditions

I ν (0, µ) = µQ 0 ν , I ν (Z, -µ) = 0 , 0 < µ < 1, ( 44
)
and the notations

φ := 1 2 1 -1 φ(µ)dµ , ψ := ∞ 0 ψ(ν)dν .
First we study the case = 0: since κ is constant, one can average in frequency, and set Ī0 (τ, µ)

:= I 0 ν (τ, µ) , Ī0 (τ ) = b ν (Φ 0 ) = Φ 0 (τ ) .
One finds

µ∂ τ Ī0 (τ, µ) + κ( Ī0 (τ, µ) -Φ 0 ) + λ( Ī0 (τ, µ) -Ī0 (τ )) = 0 ,
or, equivalently

µ∂ τ Ī0 (τ, µ) + (κ + λ)( Ī0 (τ, µ) -Ī0 (τ )) = 0 , τ > 0 , Ī0 (0, µ) = µ Q 0 ν , 0 < µ < 1.
Once Ī0 is known, one recovers I 0 ν by using the semi-analytical formula (11), which in this case is

I 0 ν (τ, µ) =1 0<µ<1 µe -(κ+λ)τ µ Q 0 ν + 1 0<µ<1 τ 0 e -(κ+λ)(τ -y) µ κb ν ( Ī0 (y)) + λ I 0 ν (y) µ dy + 1 -1<µ<0 ∞ τ e -(κ+λ)(y-τ ) |µ| κb ν ( Ī0 (y)) + λ I 0 ν (y) |µ| dy . (45)
This is precisely the solution of the Milne problem in the case of radiative transfer in a grey atmosphere.

In order to study the effect of the perturbation of the absorption defined above, we seek to compute

I ν (τ, µ) = ∂I ν (τ, µ) ∂ =0 , Φ (τ ) = ∂Φ (τ ) ∂ =0 , u (τ, µ) = ∂u (τ, µ) ∂ =0 . Let ḃν (Φ) = ∂ Φ b ν (Φ).
Let us differentiate in both sides of the equations (43) at = 0: one easily finds that

µ∂ τ I ν + κ(I ν -ḃν (Φ 0 )Φ 0 ) + λ(I ν -I ν ) + δκ ν (I 0 ν -b ν (Φ 0 )) = 0 , κ I ν -κ ḃν (Φ 0 ) =1 Φ + δκ ν ( I 0 ν -b ν (Φ 0 )) = 0 ,
and we can recast the second equality as

Φ (τ ) = Ī0 (τ ) + 1 κ δκ ν ( I 0 ν (τ ) -b ν (Φ 0 (τ ))) = Ī0 (τ ) + 1 κ δκ ν ( I 0 ν (τ ) -b ν ( Ī0 (τ ))) . (46)
Averaging in frequency the RT equation, we find that

µ∂ τ Ī0 (τ, µ) + κ( Ī0 (τ, µ) -ḃν (Φ 0 (τ )) =1 Φ (τ )) + δκ ν (I 0 ν (τ, µ) -b ν (Φ 0 (τ ))) = 0 ,
and we can further eliminate Φ (τ ) between the last two equations, to find,

µ∂ τ Ī0 (τ, µ) + (κ + λ)( Ī0 (τ, µ) -Ī0 (τ )) = -δκ ν (I 0 ν (τ, µ) -I 0 ν (τ )) . ( 47
)
If δκ is small everywhere but large near ν * and if κ + λ 1 so that the righthandside dominates the seccond term on the left, then (47) becomes

µ∂ τ Ī0 (τ, µ) = -δκ(I 0 ν * (τ, µ) -I 0 ν * (τ )) .
Together with (46) it gives the sign of T (z) by,

Φ (τ ) = -δκ τ I 0 ν * (t, µ) µ dt + δκ κ ( I 0 ν * (τ ) -b ν * ( Ī0 (τ ))) . ( 48
)
when κ 1 the second terms dominate the first one and it is possible to obtain a positive sign when I 0 ν * (τ ) > b ν * ( Ī0 (τ )).

Numerical Simulations

We study the sign of (42). As it is a complex surface of ν and τ , we look at the frequencies and altitudes for which (42) changes sign when one of the two variable varies. Figure 3 shows the frequencies and altitudes where there is a change of sign, for the grey case κ = 0.5. Consider now κ ν = 0.5 + δκ1 (ν 1 ,ν 2 ) . We use an AD (automatic differentiation) extension of the computer program for (45) and display the sensitivity of T with respect to δκ for 3 ranges of support (ν 1 , ν 2 ). On Figure 4 it is seen that when ν 1 = 0.2, ν 2 = 0.3 the sensitivity of T is negative, meaning that increasing κ ν in the range (ν 1 , ν 2 ) triggers a decrease of temperature at all altitudes. Conversely with ν 1 = 0.6, ν 2 = 0.8 the sensitivity is positive above 800m meaning that an increase of κ ν in that range triggers a temperature decrease below 800m and an increase above 800m. The same is true for T 1 but at 1500m instead of 800m. 0

•10 -3

Altitude -km Temperature Sensitivities (K)

T 0 T 1 T 2
Figure 4: Scaled temperatures sensitivities z → T (z) versus altitude, computed with κ ν = 0.5 + δκ1 (ν 1 ,ν 2 ) and δκ, at δκ = 0 when case 0: ν 1 = 0.2, ν 2 = 0.3. Case 1: ν 1 = 0.3, ν 2 = 0.4. Case 2 is with ν 1 = 0.6, ν 2 = 0.8 . Notice the change of signs in the sensitivity T 2 .

[7] F. Golse and O. Pironneau. Stratified radiative transfer for multidimensional fluids.

Compte-Rendus de l'Académie des Sciences (Mécanique), to appear 2022.

[8] H. Kaper and H. Hengler. Mathematics & Climate. SIAM publications, Philadelphia., 2013.

[9] Michael F. Modest. Radiative Heat transfer. Academic Press, 2013.

[10] O. Pironneau. A fast and accurate numerical method for radiative transfer in the atmosphere. Compte Rendus de l'académie des sciences (Math), December 2021.

z 0 r

 0 (ζ)dζ. Observe that ∂ z = ∂τ ∂τ ∂ τ = r(z)∂ z ; consequently r disappears from (5), H is replaced by Z = H 0 r(ζ)dζ, and z is replaced by τ . When κ and a are constant, Īr = ∞ 0 I ν dν is solution of the "grey model":

Figure 3 :

 3 Figure 3: Switch sign points: left when ν is increased at fixed z. Right : when altitude is increase at fixed ν.

1 2 /

 2 / main . cpp : s o l v e s t h e r a d i a t i v e t r n s f e r e q u a t i o n s 3 // f o r any a b s o r p t i o n c o e f f i c i e n t kappa 4 // d o e s a u t o m a t i c d i f f e r e n t i a t i o n i f NOAUTODIF i s u n d e f i n e d 5 // Created by O l i v i e r Pironneau on 1 0 / 1 0 / 2 0 2 1 . 6 // 7 8 #i n c l u d e <i o s t r e a m > 9 #i n c l u d e <f s t r e a m > 10 #i n c l u d e <cmath> 11 #i n c l u d e <s t r i n g > 12 #i n c l u d e <time . h> 13 u s i n g namespace s t d ; 14 #d e f i n e s q r ( x ) ( x * x ) 15 16 #d e f i n e NOAUTODIF 17 18 #i f d e f i n e d ( NOAUTODIF) 19 #d e f i n e ddouble d o u b l e 20 #d e f i n e f f a b s ( x ) ( f a b s ( x ) ) 21 c o n s t d o u b l e dknu =0; 22 #e l s e 23 #i n c l u d e " ddouble . h" 24 #d e f i n e f f a b s ( x ) ( f a b s ( x . v a l [ 0 ] ) ) 25 c o n s t ddouble dknu ( 0 . , 1 . ) ; // means dknu =0.1 and AD computes d T/ d dknu 26 #e n d i f 27 c o n s t b o o l g r e y c a s e=t r u e ; // s e e b e g i n n i n g o f main ( ) f o r p a r a m e t e r s . I f f a l s e u s e below 28 c o n s t s t r i n g kappa012=" kappa0f " ; // d e f i n e s 3 r u n s with kappa vs nu i n columns 2 , 3 , 4 29 c o n s t b o o l v e r b o s e=f a l s e ; 30 c o n s t i n t Ntau =60; // nb p o i n t s i n tau 31 c o n s t i n t kmax=16; // nb f i x e d p o i n t i t e r a t i o n s 32 c o n s t d o u b l e Z=1-exp ( -12.) ; // max tau a f t e r change o f var 33 c o n s t d o u b l e SBsun =2.03 e -5 * (2 * s q r t ( 2 . ) ) / 0 . 7 ; // s c a l e d s u n l i g h t power 34 c o n s t d o u b l e Tsun = 1 . 2 0 9 ; // s c a l e d sun t e m p e r a t u r e 35 c o n s t i n t jmaxmax=600; // max o f max nb o f p o i n t s f o r i n t e g r a t i o n i n nu r a n g e 36 c o n s t i n t newton = 5 0 ; // t o compute t h e t e m p e r a t u r e from i n t k * Botlzmann= i n t k * Imean 37 c o n s t d o u b l e epsdycho =0.01 , epsnewton =1. e -12; // p r e c i s i o n f o r dychotomy b e f o r e Newton 38 c o n s t d o u b l e tmin =1. e -10; // min t i n ExpInt ( t ) 39 c o n s t d o u b l e d t t = 0 . 0 0 5 ; // min i n t e g r a t i o n s t e p s i z e i n ExpInt i n t e g r a l s // i f your c o m p i l e r has i t o r i f you can l i n k t o g s l you may adapt t h i s f u n c t i o n 64 // i t computes E1 ( t ) * B when t <18 65 d o u b l e a b s t=f f a b s ( t ) ; 66 c o n s t i n t K = 9+( abst -1) * 4 ; ; // p r e c i s i o n i n t h e e x p o n e n t i a l i n t e g r a l f u n c t i o n E1 67 c o n s t d o u b l e gaNtaua = 0 . 5 7 7 2 1 5 6 6 4 9 0 1 5 3 3 ; // s p e c i a l i n t e g r a t i o n f o r l o g ( t ) 68 i f ( abst<tmin ) r e t u r n 0 . ; // b e c a u s e i n t e g r a l 0 ˆt ( l o g x ) dx ˜0 69 i f ( abst >18) { c o u t << " v a l u e o f E 1 i s i n c o r r e c t with "<<t<<">18"<<e n d l ; r e t u r n 0 ; } 70 ddouble ak=(

  93i f s t r e a m k a p p a f i l e ( m y k a p p a f i l e ) ; 94 i n t j =-1; 95 d o u b l e kappaux =0, dummy, nuj = 0 . 0 1 ; 96 w h i l e ( ( j++<=jmaxmax ) &&(! k a p p a f i l e . e o f ( ) ) ) { 97 i f ( which==0) k a p p a f i l e >> nuj >>kappaux >> dummy >> dummy >> dummy ; 98 e l s e i f ( which==1) k a p p a f i l e >> nuj >> dummy >>kappaux >> dummy >> dummy ; 99 e l s e i f ( which==2) k a p p a f i l e >> nuj >> dummy >> dummy >>kappaux >> dummy ; 100 kappanu [ j ] = fmax ( kappaux , kappamin ) ; 101 nu [ j ]= nuj ; 102 } 103 k a p p a f i l e . c l o s e ( ) ; 104 jmax=j ; 105 r e t u r n 0 ; 106 } 107 108 d o u b l e Bsun ( c o n s t d o u b l e nu ) { r e t u r n SBsun * s q r ( nu ) * nu / ( exp ( nu/Tsun ) -1) ; } // Boltzmann 109 110 ddouble BB( c o n s t d o u b l e nu , c o n s t ddouble T) { 111 i f (T<1. e -7) r e t u r n 0 ; 112 i f ( nu<1e -10) r e t u r n T * s q r ( nu ) ; 113 r e t u r n s q r ( nu ) * nu / ( exp ( nu/T) -1) ; } // Boltzmann 114 ddouble dBB( c o n s t d o u b l e nu , c o n s t ddouble T) { 115 i f (T<1. e -7) r e t u r n 0 ; 116 i f ( nu <1. e -10) r e t u r n s q r ( nu ) ; 117 ddouble a = exp ( nu/T) ; r e t u r n a * s q r ( nu * nu / ( a -1)/T) ; 118 } 119 120 ddouble2 intBS ( c o n s t i n t jnu , c o n s t d o u b l e tau , c o n s t d o u b l e tmin , c o n s t d o u b l e tmax ) { 121 // r e t u r n s t h e c o n v o l u t i o n t-i n t e g r a l f o r t h e mean i n mu o f kappa * I n u and muˆ2 kappa * I n u 122 ddouble Imean =0, Imu2mean=0; 123 ddouble kappa=kappanu [ jnu ] ; 124 c o n s t d o u b l e dt=fmin ( dtt , nt / ( tmax-tmin ) ) ; 125 f o r ( d o u b l e t=tmin ; t<tmax ; t+=dt ) { 126 i n t i t = i n t ( ( Ntau -1) * t /Z ) ; // p a r a b o l i c l e n g t h f o r f r e q u e n c y s t e p 127 d o u b l e arnu4=a r [ it ] * s q r ( nu [ jnu ] -0 . 8 ) * s q r ( nu [ jnu ] -1 . 2 ) * ( nu [ jnu ] > 0 . 8 ) * ( nu [ jnu ] < 1 . 2 ) * 4 0 ; 128 ddouble H0 = kappa * ( BB( nu [ jnu ] , T [ i t ] ) * (1-arnu4 ) + 129 ( a i [ i t ] + 1 . 1 2 5 * arnu4 ) * G[ jnu ] [ i t ] -1.125 * arnu4 * S [ jnu ] [ i t ] ) ; 130 ddouble H2 = -0.375 * arnu4 * kappa * (G[ jnu ] [ i t ] -3 * S [ jnu ] [ i t ] ) ; 131 i f ( kappa * ( t-tau ) !=0) { 132Imean += dt * H0 * ( e x p i n t E 1 ( kappa * f a b s ( tau-t ) )+e a l b * e x p i n t E 1 ( kappa * ( tau+t ) ) ) / 2 ; 186 f o r ( i n t j =1; j <jmax ; j ++){ 187 d o u b l e dnu=nu [ j ]-nu [ j -1 ] ; // v a r i a b l e i n t e g r a l i n c r e m e n t 188 nu1=(nu [ j ]+nu [ j -1]) / 2 ; 189 r h s+=kappanu [ j ] * G[ j ] [ i ] * dnu ; 190 l e f t += kappanu [ j ] * BB( nu1 , T0) * dnu ; 191 d e r i v += kappanu [ j ] * dBB( nu1 , T0) * dnu ; 192 } 193 p r e s f u n c = rhs-l e f t ; 194 i f ( f f a b s ( d e r i v )>1e -10) T [ i ] = T0+p r e s f u n c / d e r i v ; 195 i f ( v e r b o s e ) cout<<T [ i ]<<" newton "<<p r e s f u n c <<e n d l ; 196 } 197 i f ( inewton>=newton ) c o u t << "Newton p r e c i s i o n d o u b t f u l "<<e n d l ; t getT ( ) { // r e t u r n t e m p e r a t u r e by K i r c h h o f ' s law when kappa i s c o n s t a n t 203 c o n s t d o u b l e p i = 4 * atan ( 1 . ) ; f o r ( i n t i =0; i <jmax ; i ++){ 205 ddouble r h s =0; 206 f o r ( i n t j =1; j <jmax ; j ++) 207 r h s+=G[ j ] [ i ] * ( nu [ j ]-nu [ j -1]) ; 208 T [ i ]= s q r t ( s q r t ( 1 5 * r h s ) ) / p i ; t g e t I S n u ( c o n s t i n t j ) { // r e t u r n s mean I n u & mean muˆ2 I n u 214 d o u b l e nu1 = nu [ j ] ; 215 ddouble kappa=kappanu [ j ] ; 216 i f ( kappa <0.01) kappa = 0 . 0 1 ; 217 f o r ( i n t i =0; i <Ntau ; i ++){ 218 d o u b l e x=i * Z/ ( Ntau -1) ; 219 ddouble2 aux = intBS ( j , x , 0 , Z ) ; 220 I n u t [ i ] = aux . d1 + Bsun ( nu1 ) * ( a l p h a * e x p i n t E 3 ( kappa * ( x+Z ) ) + (1a l p h a ) * e x p i n t E 3 ( kappa * ( Z-x ) ) 221 ) / 2 ; // +e a l b * a l p h a * e x p i n t E 3 ( kappa * (Z+x ) ) ) / 2 ; 222 Snut [ i ] = aux . d2 + Bsun ( nu1 ) * ( a l p h a * e x p i n t E 5 ( kappa * x ) +(1-a l p h a ) * e x p i n t E 5 ( kappa * ( Z-x ) ) 223 ) / 2 ; // +e a l b * a l p h a * ( e x p i n t E 5 ( kappa * (Z+x ) ) ) ) / 2 ; t m u l t i B l o c k ( d o u b l e i n i t T ) 229 { 230 f o r ( i n t i =0; i <Ntau ; i ++){ T [ i ]= i n i t T ; // i n i t i a l i z e 231 f o r ( i n t j =0; j<=jmax ; j ++){ G[ j ] [ i ] = 0 ; S [ j ] [ i ] = 0 ; } 232 } 233 f o r ( i n t k=0;k<kmax ; k++){ // f i x e d p o i n t l o o p : f i r s t update F 234 f o r ( i n t i =0; i <Ntau ; i ++){ 235 I n u t [ i ] = 0 ; Snut [ i ] = 0 ; 236 } 237 f o r ( i n t j =0; j <jmax ; j ++){ 238 g e t I S n u ( j ) ; 239 f o r ( i n t i =0; i <Ntau ; i ++){ 240 G[ j ] [ i ]= I n u t [ i ] ; 241 S [ j ] [ i ]= Snut [ i ] ; g r e y c a s e ) 245 #i f d e f i n e d ( NOAUTODIF) 246 genT ( ) ; // r e p l a c e by getT ( ) t o speed-up but know what you do ! 247 #e l s e 248 genT ( ) ; // Then update T . AD need genT ( ) and f a l s e with getT ( ) b e c a u s e kappa v a r i e s 249 #e n d i f ( i n t j =0; j <jmax ; j ++) 254 f o r ( i n t i =0; i <Ntau ; i ++){ normG+=f f a b s (G[ j ] [ i ] ) ; normS+=f f a b s ( S [ j ] [ i ] ) ; } 255 c o u t << "k= "<<k <<" "<<T[2]<<" "<<normG<<" "<<normS<<e n d l ; 256 // o f s t r e a m t e m p f i l e ( b a s e d i r +"tempHistScatLow"+s t d : : t o s t r i n g ( k ) +". t x t " ) ; // study c o n v e r g e n c e 257 // f o r ( i n t i =1; i <Ntau ; i ++) t e m p f i l e <<-l o g (1-i * Z/ ( Ntau -1) )<<"\t " <<T [ i ] <<e n d l ; t main ( i n t argc , c o n s t c h a r * argv [ ] ) { 263 f o r ( i n t i t =0; i t <Ntau ; i t ++){ // d e f i n e s s c a t t e r i n g v a l u e s vs tau 264 d o u b l e tm = i t * Z/Ntau ; 265 a i [ i t ] = a i s * fmax ( tm-tm1 , 0 . ) * fmax ( tm2-tm , 0 . ) * 4/ s q r ( tm2-tm1 ) ; 266 a r [ i t ] = a r s * fmax ( tm-tm2 , 0 . ) / ( Z-tm2 ) ; 267 } 268 f o r ( i n t which =0; which <1; which++){ // change t o which <1 f o r a s i n g l e computation 269 i f ( g r e y c a s e ) { 270 c o n s t d o u b l e kappa0 =0.5 , numin = 0 . 0 5 , numax = 1 5 ; 271 jmax = 4 0 0 ; 272 f o r ( i n t j =0; j <jmax ; j ++){ 273 nu [ j ]= numax/(1+( jmax-j ) * ( numax-numin ) /numin/jmax ) ; // uniform i n wavelength 274 // kappanu [ j ]= kappa0 ;//+ dknu * ( nu [ j ] > 0 . 2 ) * ( nu [ j ] < 0 . 3 ) ; 275 kappanu [ j ]= kappa0 * ( nu [ j ] <6) + 0 . 1 ; yk ap p af il e , which ) ; 279 // f o r ( i n t j =0; j <jmax ; j ++) i f ( ( nu [ j ] < 3 . / 4 ) * ( nu [ j ] > 3 . / 7 ) ) kappanu [ j ]=fmax ( 0 . 5 , 1 . 5 * kappanu [ j ] ) ; t r e a m ( b a s e d i r+m y r e s u l t t e m p e r a t u r e+t o s t r i n g ( which +10 * ( a i s >0) )+" . t x t " ) ;

Imu2mean += dt * H0 * ( e x p i n t E 3 ( kappa * f a b s ( tau-t ) )+e a l b * e x p i n t E 3 (

d o u b l e nu1 =0;

cout<<"\n i t e r a t i o n s \ t[ T ] n e a r e a r t h[ T ] f a r | | G | | and | | S | | \ n" ;

kappa * ( tau+t ) ) ) / 2 ; 134 i f (H2 !=0) { 135 Imean += dt * H2 * ( e x p i n t E 3 ( kappa * f a b s ( tau-t ) )+e a l b * e x p i n t E 3 ( kappa * ( tau+t ) ) ) / 2 ; 136 Imu2mean += dt * H2 * ( e x p i n t E 5 ( kappa * f a b s ( tau-t ) )+e a l b * e x p i n t E 5 ( kappa * ( tau+t ) ) ) /